Code: Select all
'=======================================================================
'an approach which applies the (usually quite fast) runtime-Sqr-function (multiple times)
'the optional i_0_11 allows to adjust the amount of iterations - and thus the precision...
'with its default: i_0_11=6 the error is +-1% ... with 7 it is half that, a.s.o. up to 11
'-----------------------------------------------------------------------
' Written by OSchmidt, modified by Richard
' WARNING: errors rise rapidly as x approaches zero
'-----------------------------------------------------------------------
Function FastPow( Byval x As Single, Byval n As Single) As Single
Static As Single L( 0 To 11 ) = { 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625,_
7.8125e-3, 3.90625e-3, 1.953125e-3, 9.765625e-4, 4.8828125e-4, 2.44140625e-4 }
If x = 0 Then Return 0 ' special case for 0-values of x
Dim As Long i, i_0_11 = 6 ' this parameter has been internalised for testing
Dim As Single r = 1
If n > 0 Then ' handle positive exponents
While n >= 1
r *= x
n -= 1
Wend
For i = 0 To i_0_11
x = Sqr(x)
If n >= L(i) Then r *= x: n -= L(i)
If n <= L(i_0_11) Then Exit For ' early exit (n can't get any smaller)
Next
Else
If n = 0 Then Return 1 ' special case, avoids n=0 crash
r = 1 / FastPow( x, -n ) ' handle negative exponents
End If
Return r
End Function
'-----------------------------------------------------------------------
' Compare the single precision x^b with algorithm under test
' reference ref time test test time speed ratio Error bounds
' single_pwr 69.8ns FastPow11 154.7ns 221.67% 0.03125%
' single_pwr 69.8ns FastPow10 141.5ns 202.68% 0.0625%
' single_pwr 69.7ns FastPow9 127.0ns 182.23% 0.125%
' single_pwr 69.7ns FastPow8 112.3ns 161.07% 0.25%
' single_pwr 69.7ns FastPow7 99.8ns 143.11% 0.5%
' single_pwr 69.7ns FastPow6 87.6ns 125.68% 1.%
' single_pwr 69.7ns FastPow5 75.5ns 108.31% 2.%
' single_pwr 69.8ns single_pwr 69.8ns 100.00% 0.00001% XXXXXXXX
' single_pwr 69.7ns FastPow4 63.3ns 90.81% 4.%
' single_pwr 69.8ns FastPow3 52.0ns 74.42% 8.%
' single_pwr 69.8ns Raise 41.7ns 59.72% 0.50% XXXXXXXXXXXX
' single_pwr 69.8ns FastPow2 41.4ns 59.36% 16.%
' single_pwr 69.9ns FastPow1 32.5ns 46.52% 32.%
'=======================================================================
' Approximate the power function; written by Richard
'-----------------------------------------------------------------------
' Approximate x ^ b, for x > 0; |b| < 60. Avoid infinite or denormalised x values
Function Raise( Byval x As Single, Byval b As Single ) As Single
' to raise x to power b, y = x^b, y = Exp2( b * Log2( x ) )
' IEEE 754 Single Seeeeeeeefffffffffffffffffffffff ' Sign_Log_Linear format
#define expo_bias &b00111111000000000000000000000000uL ' bias = 126 * 2^23
#define frac_mask &b00000000011111111111111111111111uL ' positive fraction
#define frac_bits 23 ' bits used to store fraction, ignoring implicit msb
#define bias 126 ' exponent that puts single into range from 0.500 to 0.999
Dim As Long Ptr fp = Cptr( Long Ptr, @x ) ' pointer to bit pattern of x
'-------- approximate Log2(x) --------------------------------------
Dim As Long expo = ( *fp Shr frac_bits ) - bias ' extract exponent
*fp = ( *fp And frac_mask ) Or expo_bias ' bias x into range 0.5 to 0.999
#Define a0 -3.5353058e0 ' log2 polynomial approximation
#Define a1 8.2616852e0 ' over the range from 0.5 to 1.0
#Define a2 -8.7235965e0 ' maximum error is 1.54e-4 = 0.000154
#Define a3 5.3685806e0 ' Zeros are located at
#Define a4 -1.3713634e0 ' 0.500, 0.575, 0.720, 0.885, 1.000
x = b * ( Csng( expo ) + a0 + ( a1 + ( a2 + ( a3 + a4 * x ) * x ) * x ) * x )
'-------- approximate 2^x, is antiLog2( x ) ------------------------
expo = Int( x ) ' this integer part adjusts the exponent of Alog2(x)
x -= expo ' x reduced to fraction in range 0.000 thru 0.999
x = 1e0 + ( 0.6607687e0 + 0.3392313e0 * x ) * x ' approx 2^x over range 0 to 1
*fp += ( expo Shl frac_bits ) ' restore early integer to biased log exponent
Return x
End Function
'-----------------------------------------------------------------------
' Maximum error is < ± 0.30% for b < ± 2.5
' Maximum error is < ± 0.35% for b < ± 7.0
' Maximum error is < ± 0.50% for b < ± 22.
' Maximum error is < ± 0.92% for b < ± 60.
' Beyond which the exponent overflows, so it all falls apart
'
'=======================================================================
' set up graphics screen
Dim As Integer sx, sy, sz
Screeninfo sx, sy, sz
sx -= 32 ' shrink
sy -= 24
sz = 4 ' 16 colours
Screenres sx, sy, sz
'--------------------------------------
' graphics scale
Dim As Single x, y, z, t, q
Dim As Single lhs = 0.0, rhs = 4.2, lo = -0.5, hi = 0.5
Window ( lhs, lo ) - ( rhs, hi )
Line( lhs, 0 )-( rhs, 0 ), 1
For x = 1 To 4
Line( x, lo )-( x, hi ), 1
Next x
'=======================================================================
' plot the error curves
Dim As Single max = -1e6, min = 1e6, scale = 40
Dim As Integer i, k = 7
Print
Print " Powers = ";
'--------------------------------------
For q = -4.5 To 4.51 Step 0.1
k += 1 ' 16 colours
If k > 15 Then k = 8
Color k
Print Using "####.##"; q;
'----------------------------------
For x = 0 To rhs Step ( rhs - lhs ) / ( sx * 2 )
y = x^q ' correct
'--------------------------------------
z = FastPow( x, q ) ' Raise( x, q ) ' ' approximation function to plot
'--------------------------------------
t = ( z - y ) / y ' t is the error
If min > t Then min = t
If max < t Then max = t
Pset( x, scale * t ), k
Next x
Next q
Line( lhs, 0.005 * scale )-( rhs, 0.005 * scale ), 1
Line( lhs, -0.005 * scale )-( rhs, -0.005 * scale ), 1
Draw String ( 0.01, 0.0049 * scale ), " +0.5 %", 15
Draw String ( 0.01, -.0052 * scale ), " -0.5 %", 15
Line( lhs, max * scale )-( rhs, max * scale ), 2
Line( lhs, min * scale )-( rhs, min * scale ), 2
'=======================================================================
For x = 0 To 4
Draw String ( x+.01, -0.005), Str(x), 15
Next x
Color 15
Print
Print
Print Using " Error bounds ####.### % ####.### %"; max * 100; min * 100
'=======================================================================
Sleep
'=======================================================================