Need faster Exponentiate

General FreeBASIC programming questions.
Provoni
Posts: 513
Joined: Jan 05, 2014 12:33
Location: Belgium

Need faster Exponentiate

Post by Provoni »

Hey all,

I need a faster exponentiate for my program, is that possible, perhaps an approximation?

Thanks

https://www.freebasic.net/wiki/wikka.ph ... ponentiate

Test environment:

Code: Select all

'test environment

randomize timer

screenres 640,480,32

dim as integer i,maxrn=10000000
dim as double o,rn(maxrn)

for i=1 to maxrn
	rn(i)=rnd
next i

sleep 100
dim as double t1=timer

for i=1 to maxrn
	o+=rn(i)^1.1 'this simulates what I need though the 1.1 could be 0.9 or 1.123 etc
next i

print o,timer-t1

sleep
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Need faster Exponentiate

Post by dodicat »

slight speed improvement

Code: Select all

'test environment


#define pwr(x,n) exp((n)*log((x)))


randomize timer

screenres 640,480,32

dim as integer i,maxrn=10000000
dim as double o,rn(maxrn)

for i=1 to maxrn
   rn(i)=rnd
next i

sleep 100
dim as double t1=timer

for i=1 to maxrn
   'o+=rn(i)^1.1 'this simulates what I need though the 1.1 could be 0.9 or 1.123 etc
   o+=pwr(rn(i),1.1)
next i

print o,timer-t1

sleep
print pwr(5,9.3)

'y=x^n
'log y=log x^n=n*log x
'y= exp(n*log x) 
Provoni
Posts: 513
Joined: Jan 05, 2014 12:33
Location: Belgium

Re: Need faster Exponentiate

Post by Provoni »

Thanks dodicat, it is a noteworthy speed improvement for me.

I also found that using numbers like ^0.25, ^0.50, ^0.75, ^1.50 have almost no calculation cost with the standard Exponentiate function.
Provoni
Posts: 513
Joined: Jan 05, 2014 12:33
Location: Belgium

Re: Need faster Exponentiate

Post by Provoni »

Here is a power approximation function that I found and converted:

Code: Select all

function powa1(byval a as double,byval b as double)as double
	
	dim as double rv,ln,am1
	ln=log(a)
	am1=b-1.0
	rv=a*ln*am1
	ln*=ln
	am1*=am1
	rv+=.5*a*ln*am1
	rv+=a
	return rv
	
end function
I would like to convert this one to FreeBASIC https://martin.ankerl.com/2007/10/04/op ... a-and-c-c/ but I do not understand some parts of the code.
Provoni
Posts: 513
Joined: Jan 05, 2014 12:33
Location: Belgium

Re: Need faster Exponentiate

Post by Provoni »

Code: Select all

double Fast_Pow(double a, double b) //fastpower originally developed by Martin Ankerl
 {
  int tmp = (*(1 + (int *)&a));
  int tmp2 = (int)(b * (tmp - 1072632447) + 1072632447);
  double p = 0.0;
  *(1 + (int * )&p) = tmp2;
  //p = p * a / 2.71828F ; failed attempt to auto correct the accuracy
  return tmp;
 }
paul doe
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: Need faster Exponentiate

Post by paul doe »

Provoni wrote:...
I would like to convert this one to FreeBASIC https://martin.ankerl.com/2007/10/04/op ... a-and-c-c/ but I do not understand some parts of the code.
Here:

Code: Select all

function approximatePow( _
  byval a as double, _
  byval b as double ) _
  as double
 
  dim as longint _
    tmp => clng( *cptr( longint ptr, @a ) shr 32 ), _
    tmp2 => clng( b * ( tmp - 1072632447 ) + 1072632447 ), _
    tmp3 => tmp2 shl 32
  
  return( *cptr( _
    double ptr, _
    @tmp3 ) )
end function

? approximatePow( 3.14, 2.32 )
? 3.14 ^ 2.23

sleep()
Didn't test the equivalence with the C++ version yet...
Last edited by paul doe on May 13, 2019 2:35, edited 1 time in total.
paul doe
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: Need faster Exponentiate

Post by paul doe »

This is the result of the C++ version:

Code: Select all

-1375731712
--------------------------------
Process exited after 0.07882 seconds with return value 0
Presione una tecla para continuar . . .
Seems fine to me (if even more off than the final approximated version).
Provoni
Posts: 513
Joined: Jan 05, 2014 12:33
Location: Belgium

Re: Need faster Exponentiate

Post by Provoni »

Thanks paul doe, it is very fast but the approximation is not very good though it may be good enough for my needs.
paul doe
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: Need faster Exponentiate

Post by paul doe »

Provoni wrote:Thanks paul doe, it is very fast but the approximation is not very good though it may be good enough for my needs.
Indeed, the error is huge (and it'll only get worse the more you use them). Ahh, these kind of codes brings back so many fond memories... ;)
Provoni
Posts: 513
Joined: Jan 05, 2014 12:33
Location: Belgium

Re: Need faster Exponentiate

Post by Provoni »

paul doe wrote:
Provoni wrote:Thanks paul doe, it is very fast but the approximation is not very good though it may be good enough for my needs.
Indeed, the error is huge (and it'll only get worse the more you use them). Ahh, these kind of codes brings back so many fond memories... ;)
I saw this on YouTube a while ago. Of course I also played Wolfenstein 3D, Doom 1 & 2, Duke3D, Quake 1, 2 & 3 back in the days. :-)
paul doe
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: Need faster Exponentiate

Post by paul doe »

Provoni wrote:I saw this on YouTube a while ago. Of course I also played Wolfenstein 3D, Doom 1 & 2, Duke3D, Quake 1, 2 & 3 back in the days. :-)
Yep, those were the days (I remember how we had to make our own wiring to play deathmatch in DooM haha). Doom II deathmatch just doesn't get old ;)
Richard
Posts: 3096
Joined: Jan 15, 2007 20:44
Location: Australia

Re: Need faster Exponentiate

Post by Richard »

Code: Select all

'=======================================================================
' approximate a ^ b
Function approx_power( Byval a As Double, Byval b As Double ) As Double
    Dim As Long Ptr fp = 1 + Cptr( Long Ptr, @a )
    Dim As Long tmp = ( *fp - 1072632447 ) * b
    *fp = tmp + 1072632447
    Return a
End Function

'=======================================================================
' How and Why It Works.
' To accurately compute  y = a^b;  for non-integer b;  with a > 0;
'   take the log(a), multiply it by b, then look up the antilog, which is
'   the Exp() function, to give; y = Exp( b * Log( a ) )
' For a computer using base 2 it is convenient to use; y = 2^( b * Log2(a) )
' The IEEE double precision floating point format has an 11 bit, base two,
'   integer exponent. That will become the integer part of the log2(x) function.
'   It is simple to extract the exponent, multiply it by b, then return it to
'   the exponent, but that would be a clunky approximation because the exponent
'   is an integer count of powers of two with poor resolution.
' So that the multiply will change the result for small changes of a or power b,
'   the leading 20 mantissa bits are included as a false fractional extension to
'   the logarithmic exponent. That is the source of error in this approximation.
'   Note that the implicit 1 msb is missing from the mantissa.
' The exponent offset bias of 1022 must be removed by subtraction before the
'   multiplication. In the same operation an amount can be deducted from the
'   mantissa to significantly reduce the approximation error. The approximation
'   is computed from only the top 32 bits of the 64 bit double precision format.
' &bSeeeeeeeeeeeffffffffffffffffffff  Sign bit, exponent and fraction format
' &b00111111111000000000000000000000    1071644672 = 1022 * 2^20, shifted bias
' &b00000000000011110001001001111111        987775 = fudge amount
' &b00111111111011110001001001111111    1072632447 = sum is the magic number
'=======================================================================
Richard
Posts: 3096
Joined: Jan 15, 2007 20:44
Location: Australia

Re: Need faster Exponentiate

Post by Richard »

@Provoni
The 12% errors in the above crude approximation can be reduced by taking the log2 of the mantissa and the antilog2 of the result. That can be done with simple polynomials faster than the intrinsic function calls. It is important to use only primitive +, –, * and bitwise operators. Avoid division, or decisions that break the instruction flow.

Here are times and errors for routines that work for Single and for Double inputs.
Double intrinsic a^b. Time = 100.0 Accuracy is about 16 digits.
Single intrinsic a^b. Time = 87.5 Accuracy is about 7 digits.
Double approximation. Time = 65.3 Maximum error is about ±0.5% of output.
Single approximation. Time = 55.97 Maximum error is about ±0.5% of output.
Crude approximation. Time = 7.10 Error is about ±12.% of output.

Here is the Double routine.

Code: Select all

Function Double_power( Byval x As Double, Byval b As Double ) As Double
    #define expo_bias &b00111111111000000000000000000000uL ' bias = 1022 * 2^20
    #define expo_mask &b01111111111100000000000000000000uL ' exponent field
    #define mant_mask &b00000000000011111111111111111111uL ' mantissa field, positive
    Dim As Long Ptr fp = 1 + Cptr( Long Ptr, @x ) ' make pointer to bit pattern of x
    Dim As Long expo = ( ( *fp And expo_mask ) Shr 20 ) - 1022  ' extract exponent
    *fp And= mant_mask  ' kill the exponent bits
    *fp Or=  expo_bias  ' bias exponent so input is between 0.5000 and 0.9999
    x = expo - 3.153635 + x * ( 6.095832 + x * ( -4.207588 + x * 1.266028 ) ) ' Log2(x)
    x *= b          ' raise Log to power b,  x = b * Log(x)
    expo = Int( x ) ' integer part will become exponent of antilog2( x )
    x -= expo       ' remove to make x range 0.0000 to 0.9999
    x = 1.00 + x * ( 0.6557133 + x * 0.3442867 )   ' my quick 2^x over 0 to 1
    *fp += ( expo Shl 20 )  ' add the integer to biased log exponent
    Return x        ' approximate x ^ b
End Function        ' maximum error is about 0.5 percent of output
and here is the faster single routine using the same approximations.

Code: Select all

Function single_power( Byval x As Single, Byval b As Single ) As Single
    ' IEEE 754 Single   Seeeeeeeefffffffffffffffffffffff  ' single format
    #define expo_bias &b00111111000000000000000000000000uL ' bias = 126 * 2^23
    #define expo_mask &b01111111100000000000000000000000uL ' exponent field
    #define mant_mask &b00000000011111111111111111111111uL ' positive mantissa
    #define mant_bits 23
    #define bias 126
    Dim As Long Ptr fp = Cptr( Long Ptr, @x )   ' make ptr to bit pattern x
    Dim As Long expo = ( *fp Shr mant_bits ) - bias  ' extract exponent
    *fp And= mant_mask  ' kill the exponent bits
    *fp Or=  expo_bias  ' bias exponent so input is between 0.5000 and 0.9999
    x = expo - 3.153635 + x * ( 6.095832 + x * ( -4.207588 + x * 1.266028 ) ) ' Log2(x)
    x *= b          ' raise Log to power b,  x = b * Log(x)
    expo = Int( x ) ' integer part will become exponent of antilog2( x )
    x -= expo       ' remove to make x range 0.0000 to 0.9999
    x = 1.00 + x * ( 0.6557133 + x * 0.3442867 )   ' quick 2^x over 0 to 1
    *fp += ( expo Shl mant_bits )  ' add the integer to biased log exponent
    Return x        ' approximate x ^ b
End Function        ' maximum error is about 0.5 percent of output
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Need faster Exponentiate

Post by dodicat »

Did anybody here actually finish Return to castle wolfenstein?
I got to the village, caused some carnage and went into the pit.
There I got stuck.
My other machine overheats badly now, so I don't have enough time to traverse the pit.
The horror of the ever speeding fan on top of everything else!
This machine won't run wolfenstein.(some opengl hiccup)
I have the option now of re-greasing my two pentium xeon's on my other box or going back to Ebay for another machine.
I fear I shall never complete wolfenstein, I have had the game now for 15 years.
Tourist Trap
Posts: 2958
Joined: Jun 02, 2015 16:24

Re: Need faster Exponentiate

Post by Tourist Trap »

dodicat wrote:Did anybody here actually finish Return to castle wolfenstein?
You are repaid from your investment for such a longevity :D
I can't help, apart for the scene in the mountain at start, I never played tge campaign, only the multiplayer mode (as a medic).
Not a bad game, was famous for its flames effects - maybe what causes trouble to your fanes?
Post Reply