/* Parts of this software were written by Dirk Engling Those parts are considered beerware. Prost. Skol. Cheers or whatever. Original idea and the place where I heavily stole code from is http://www.quinapalus.com/efunc.html Dr Mark St John OWEN, E-mail: mail -at- quinapalus -dot- com */ #include #include #include #include #define TO_Q16(X) ((int32_t)(65536.*(X))) #define TO_Q26(X) ((int32_t)(67108864.*(X))) #define FROM_Q16(X) ((double)((X)/65536.)) #define FROM_Q26(X) ((double)((X)/67108864.)) /* It works as follows: P = pow_QX_QY_QP( x, y ) == exp( ln(x) * y ) Where if x is negative: * y must have no fractional part, * the result's sign is the lowest integral bit of y We note that in every standard Q representation ln(x) will not exceed the value 22, so that we can safely work with a Q26 representation. if ln(x) * y would overflow in that representation, so would exp( ln(x) * y ) in Q16. Finally exp() is calculated in the domain specified by script */ #if 0 The table[tm], generated from different versions of generate_table, though. +00,000000000465661287199 x = x + ( x >> 31 ) +00,000000000931322574182 x = x + ( x >> 30 ) +00,000000001862645147496 x = x + ( x >> 29 ) +00,000000003725290291523 x = x + ( x >> 28 ) +00,000000007450580569168 x = x + ( x >> 27 ) +00,000000014901161082825 x = x + ( x >> 26 ) +00,000000029802321943606 x = x + ( x >> 25 ) +00,000000059604642999034 x = x + ( x >> 24 ) +00,000000119209282445354 x = x + ( x >> 23 ) +00,000000238418550679858 x = x + ( x >> 22 ) +00,000000476837044516323 x = x + ( x >> 21 ) +00,000000953673861659188 x = x + ( x >> 20 ) +00,000001907346813825409 x = x + ( x >> 19 ) +00,000003814689989685890 x = x + ( x >> 18 ) +00,000007629365427567572 x = x + ( x >> 17 ) +00,000015258672648362398 x = x + ( x >> 16 ) +00,000030517112473186380 x = x + ( x >> 15 ) +00,000061033293680638527 x = x + ( x >> 14 ) +00,000122062862525677371 x = x + ( x >> 13 ) +00,000244110827527362707 x = x + ( x >> 12 ) +00,000488162079501351187 x = x + ( x >> 11 ) +00,000976085973055458925 x = x + ( x >> 10 ) +00,001951220131261749337 x = x + ( x >> 9 ) +00,003898640415657322889 x = x + ( x >> 8 ) +00,007782140442054948960 x = x + ( x >> 7 ) +00,015504186535965254479 x = x + ( x >> 6 ) +00,030771658666753687328 x = x + ( x >> 5 ) +00,060624621816434839938 x = x + ( x >> 4 ) +00,117783035656383455736 x = x + ( x >> 3 ) +00,223143551314209764858 x = x + ( x >> 2 ) +00,405465108108164384859 x = x + ( x >> 1 ) +00,693147180559945286227 x = ( x << 1 ) +01,386294361119890572454 x = ( x << 2 ) +02,772588722239781144907 x = ( x << 4 ) +05,545177444479562289814 x = ( x << 8 ) +11,090354888959124579628 x = ( x << 16 ) +21,487562597358305538364 x = ( x << 31 ) for negative values: -00,374693449441410697531 x-= ( x >> 2 ) + ( x >> 4 ) -00,207639364778244489562 x-= ( x >> 2 ) - ( x >> 4 ) -00,115831815525121700761 x-= ( x >> 3 ) - ( x >> 6 ) -00,060380510988907482028 x-= ( x >> 4 ) - ( x >> 8 ) -00,031748698314580298119 x-= ( x >> 5 ) -00,015996403602177928366 x-= ( x >> 6 ) + ( x >> 12 ) -00,008089270731616965762 x-= ( x >> 7 ) + ( x >> 12 ) -00,004159027401785260480 x-= ( x >> 8 ) + ( x >> 12 ) -00,003423823349553609900 x-= ( x >> 8 ) - ( x >> 11 ) -00,001832733117311761669 x-= ( x >> 9 ) - ( x >> 13 ) -00,000961766059678620302 x-= ( x >> 10 ) - ( x >> 16 ) -00,000484583944571210926 x-= ( x >> 11 ) - ( x >> 18 ) -00,000243216525424976433 x-= ( x >> 12 ) - ( x >> 20 ) #endif // input is Q16, outputQ26 static int32_t fixlog( int32_t x ) { int32_t t,y; // Assume an x in Q21 (shift by 15) and thus start with y with log(2^15) y = TO_Q26(10.39720770839918); if(x<0x00008000) x<<=16, y-= TO_Q26(11.09035488895912); /* log(65536) */ if(x<0x00800000) x<<= 8, y-= TO_Q26(5.545177444479562); /* log(256) */ if(x<0x08000000) x<<= 4, y-= TO_Q26(2.772588722239781); /* log(16) */ if(x<0x20000000) x<<= 2, y-= TO_Q26(1.386294361119891); /* log(4) */ if(x<0x40000000) x<<= 1, y-= TO_Q26(0.693147180559945); /* log(2) */ t=x+(x>>1); if((t&0x80000000)==0) x=t, y-= TO_Q26(0.405465108108164); /* log(1+1/2) */ t=x+(x>>2); if((t&0x80000000)==0) x=t, y-= TO_Q26(0.22314355131421); /* log(1+1/4) */ t=x+(x>>3); if((t&0x80000000)==0) x=t, y-= TO_Q26(0.117783035656384); /* log(1+1/8) */ t=x+(x>>4); if((t&0x80000000)==0) x=t, y-= TO_Q26(0.060624621816435); /* log(1+1/16) */ t=x+(x>>5); if((t&0x80000000)==0) x=t, y-= TO_Q26(0.030771658666754); /* log(1+1/32) */ t=x+(x>>6); if((t&0x80000000)==0) x=t, y-= TO_Q26(0.015504186535965); /* log(1+1/64) */ t=x+(x>>7); if((t&0x80000000)==0) x=t, y-= TO_Q26(0.007782140442055); /* log(1+1/128) */ t=x+(x>>8); if((t&0x80000000)==0) x=t, y-= TO_Q26(0.003898640415657); /* log(1+1/256) */ t=x+(x>>9); if((t&0x80000000)==0) x=t, y-= TO_Q26(0.001951220131262); /* log(1+1/512) */ t=x+(x>>10);if((t&0x80000000)==0) x=t, y-= TO_Q26(0.000976085973055); /* log(1+1/1024) */ t=x+(x>>11);if((t&0x80000000)==0) x=t, y-= TO_Q26(0.000488162079501); /* log(1+1/2048) */ t=x+(x>>12);if((t&0x80000000)==0) x=t, y-= TO_Q26(0.000244110827527); /* log(1+1/4096) */ t=x+(x>>13);if((t&0x80000000)==0) x=t, y-= TO_Q26(0.000122062862526); /* log(1+1/8192) */ /* Estimate residual error, log(1-x) which for small x is approx -x */ x = 2147483648-x; /* x has been promoted to Q31, substract residual error in Q26 by shifting down by 5 */ return y - ( x >> 5 ); } // input is Q26, output Q16 static uint32_t fixexp(int32_t x) { int32_t t; /* Start y with 1.0 */ uint32_t y = 0x10000; if( x >= 0 ) { t=x-TO_Q26(5.545177444479562); if(t>=0) x=t,y<<=8; /* log(256) */ t=x-TO_Q26(2.772588722239781); if(t>=0) x=t,y<<=4; /* log(16) */ t=x-TO_Q26(1.386294361119891); if(t>=0) x=t,y<<=2; /* log(4) */ t=x-TO_Q26(0.693147180559945); if(t>=0) x=t,y<<=1; /* log(2) */ t=x-TO_Q26(0.405465108108164); if(t>=0) x=t,y+=y>>1; /* log(1+1/2) */ t=x-TO_Q26(0.22314355131421); if(t>=0) x=t,y+=y>>2; /* log(1+1/4) */ t=x-TO_Q26(0.117783035656384); if(t>=0) x=t,y+=y>>3; /* log(1+1/8) */ t=x-TO_Q26(0.060624621816435); if(t>=0) x=t,y+=y>>4; /* log(1+1/16) */ t=x-TO_Q26(0.030771658666754); if(t>=0) x=t,y+=y>>5; /* log(1+1/32) */ t=x-TO_Q26(0.015504186535965); if(t>=0) x=t,y+=y>>6; /* log(1+1/64) */ t=x-TO_Q26(0.007782140442055); if(t>=0) x=t,y+=y>>7; /* log(1+1/128) */ t=x-TO_Q26(0.003898640415657); if(t>=0) x=t,y+=y>>8; /* log(1+1/256) */ t=x-TO_Q26(0.001951220131262); if(t>=0) x=t,y+=y>>9; /* log(1+1/512) */ t=x-TO_Q26(0.000976085973055); if(t>=0) x=t,y+=y>>10; /* log(1+1/1024) */ t=x-TO_Q26(0.000488162079501); if(t>=0) x=t,y+=y>>11; /* log(1+1/2048) */ t=x-TO_Q26(0.000244110827527); if(t>=0) x=t,y+=y>>12; /* log(1+1/4096) */ t=x-TO_Q26(0.000122062862526); if(t>=0) x=t,y+=y>>13; /* log(1+1/8192) */ /* from here the values in Q26 become approx 2^n, so lets check bits and fix the residual error below */ if(x&0x0001000) y+=y>>14; if(x&0x0000800) y+=y>>15; if(x&0x0000400) y+=y>>16; if(x&0x0000200) y+=y>>17; if(x&0x0000100) y+=y>>18; if(x&0x0000080) y+=y>>19; if(x&0x0000040) y+=y>>20; if(x&0x0000020) y+=y>>21; if(x&0x0000010) y+=y>>22; if(x&0x0000008) y+=y>>23; if(x&0x0000004) y+=y>>24; if(x&0x0000002) y+=y>>25; if(x&0x0000001) y+=y>>26; } else { if(x<=TO_Q26(-11.09035488895912)) return 0; t=x+TO_Q26(11.09035488895912); if(t<=0) x=t,y>>=16; t=x+TO_Q26(5.545177444479562); if(t<=0) x=t,y>>=8; t=x+TO_Q26(2.772588722239781); if(t<=0) x=t,y>>=4; t=x+TO_Q26(1.386294361119891); if(t<=0) x=t,y>>=2; t=x+TO_Q26(0.693147180559945); if(t<=0) x=t,y>>=1; t=x+TO_Q26(0.374693449441411); if(t<=0) x=t,y-=(y>>4) + (y>>2); t=x+TO_Q26(0.207639364778244); if(t<=0) x=t,y-=(y>>2) - (y>>4); t=x+TO_Q26(0.124642445207276); if(t<=0) x=t,y-=(y>>3) - (y>>7); t=x+TO_Q26(0.062457354933746); if(t<=0) x=t,y-=(y>>4) - (y>>9); t=x+TO_Q26(0.031244793038107); if(t<=0) x=t,y-=(y>>5) - (y>>11); t=x+TO_Q26(0.015748356968139); if(t<=0) x=t,y-=(y>>6); t=x+TO_Q26(0.007966216526084); if(t<=0) x=t,y-=(y>>7) + (y>>13); t=x+TO_Q26(0.004036455850488); if(t<=0) x=t,y-=(y>>8) + (y>>13); t=x+TO_Q26(0.002077351513834); if(t<=0) x=t,y-=(y>>9) + (y>>13); t=x+TO_Q26(0.001099236751907); if(t<=0) x=t,y-=(y>>10)+ (y>>13); t=x+TO_Q26(0.000610537902840); if(t<=0) x=t,y-=(y>>11)+ (y>>13); t=x+TO_Q26(0.000366278009100); if(t<=0) x=t,y-=(y>>12)+ (y>>13); t=x+TO_Q26(0.000213645867528); if(t<=0) x=t,y-=(y>>12)- (y>>15); if(x&0x0002000) y-=y>>13; if(x&0x0001000) y-=y>>14; if(x&0x0000800) y-=y>>15; if(x&0x0000400) y-=y>>16; if(x&0x0000200) y-=y>>17; if(x&0x0000100) y-=y>>18; if(x&0x0000080) y-=y>>19; if(x&0x0000040) y-=y>>20; if(x&0x0000020) y-=y>>21; if(x&0x0000010) y-=y>>22; if(x&0x0000008) y-=y>>23; if(x&0x0000004) y-=y>>24; if(x&0x0000002) y-=y>>25; if(x&0x0000001) y-=y>>26; } return y; } int fixpow(int32_t *pow, int32_t base, int32_t exponent ) { int32_t log_base, res; int64_t log_base_times_exponent; int neg = 0; if( base < 0 ) { /* negative bases only can have integer exponents */ if( exponent & 0xffff ) return -1; /* sign of power of a negative base is determined by wether exp is even */ if( exponent & 0x10000 ) neg = 1; base = abs(base); } /* To calculate pow(base,exp), we do exp( log(base) * exp ) which is mathematically the same. log_base is Q26 */ log_base = fixlog( base ); log_base_times_exponent = ( (int64_t)log_base * (int64_t)exponent ) >> 16; /* In Q0, fixexp overflows for values > 21,48756259689264 which in Q26 notation is 1442005916 */ if( log_base_times_exponent > 1442005916 ) return -2; /* In Q16, fixexp overflows for values > 10,39720770839918 */ if( log_base_times_exponent > TO_Q26(10.39720770839918) ) return -2; res = (int32_t)log_base_times_exponent; res = fixexp( res ); if( neg ) res = -res; *pow = res; return 0; } int main() { double base = -.5; double exponent = 10.; int32_t result, input, output; int error = fixpow( &result, TO_Q16(base), TO_Q16(exponent) ); // Testing fixlog for (input=1; input<=INT32_MAX/2; ++input) { int32_t expected = TO_Q26(log(FROM_Q16(input))); int32_t output = fixlog(input); if( abs(output-expected) > 8 ) printf( "%08X : %08X = %08X (%08X)\n", output, expected, abs(output-expected), input ); } // Testing fixexp for (input=-2147483648; input < TO_Q26(10.39720770839918); ++input ) { double expected = exp(FROM_Q26(input)); int32_t expected_int = TO_Q16(expected); int32_t output_int = fixexp(input); double output = FROM_Q16(output_int); if( abs(output_int-expected_int) > 8 ) { if( expected > output && output > 0 && expected / output > 1.002 ) printf( "%08X : %08X = %08X (%08d)\n", output_int, expected_int, abs(output_int-expected_int), input ); if( expected < output && expected > 0 && output / expected > 1.002 ) printf( "%08X : %08X = %08X (%08d)\n", output_int, expected_int, abs(output_int-expected_int), input ); } } printf( "pow(%lf,%lf)=%lf (%s)\n", base, exponent, FROM_Q16(result), error?"ERROR":"OK" ); return 0; }