From da540be956937b49beba8e02569f16b80561c5a8 Mon Sep 17 00:00:00 2001 From: erdgeist Date: Tue, 7 Oct 2014 04:51:49 +0200 Subject: Adding test code for fixexp. We're accurate to 22 significant bits --- fixpow.c | 61 ++++++++++++++++++++++++++++++++++++++----------------------- 1 file changed, 38 insertions(+), 23 deletions(-) diff --git a/fixpow.c b/fixpow.c index 3dc4815..3db6df3 100644 --- a/fixpow.c +++ b/fixpow.c @@ -14,7 +14,7 @@ #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.)) +#define FROM_Q26(X) ((double)((X)/67108864.)) /* It works as follows: @@ -118,7 +118,7 @@ static int32_t fixlog( int32_t x ) { 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 = 0x80000000-x; + x = 2147483648-x; /* x has been promoted to Q31, substract residual error in Q26 by shifting down by 5 */ return y - ( x >> 5 ); @@ -164,24 +164,25 @@ static uint32_t fixexp(int32_t x) { if(x&0x0000002) y+=y>>25; if(x&0x0000001) y+=y>>26; } else { - x=-x; - 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) - (x>>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<=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; @@ -221,12 +222,10 @@ int fixpow(int32_t *pow, int32_t base, int32_t exponent ) log_base = fixlog( base ); log_base_times_exponent = ( (int64_t)log_base * (int64_t)exponent ) >> 16; -#if 0 /* In Q0, fixexp overflows for values > 21,48756259689264 which in Q26 notation is 1442005916 */ if( log_base_times_exponent > 1442005916 ) return -2; -#endif /* In Q16, fixexp overflows for values > 10,39720770839918 */ if( log_base_times_exponent > TO_Q26(10.39720770839918) ) @@ -244,9 +243,25 @@ int main() { double base = -.5; double exponent = 10.; - int32_t result; + 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 ) { + int32_t expected = TO_Q16(exp(FROM_Q26(input))); + int32_t output = fixexp(input); + if( abs(output-expected) > 16 ) + printf( "%08X : %08X = %08X (%08d)\n", output, expected, abs(output-expected), input ); + } + printf( "pow(%lf,%lf)=%lf (%s)\n", base, exponent, FROM_Q16(result), error?"ERROR":"OK" ); return 0; } -- cgit v1.2.3