@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