Adding test code for fixexp. We're accurate to 22 significant bits

1 files 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;
 }