1 files changed, 93 insertions, 76 deletions

diff --git a/fixpow.c b/fixpow.c
index d81881f..cb687b4 100644
--- a/fixpow.c
+++ b/fixpow.c
@@ -11,23 +11,27 @@
 #include <stdlib.h>
 #include <math.h>
 
-#define TO_Q16(X)    ((int32_t)(65536.f*(X)))
-#define FROM_Q16(X)  ((double)((X)/65536.f))
+#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 )
+/*
+  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
+  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
+  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
@@ -72,19 +76,19 @@
 
 for negative values:
 
--00,374693449441410697531 017FAFA3 x-= ( x >>  2 ) + ( x >>  4 )
--00,207639364778244489562 00D49F69 x-= ( x >>  2 ) - ( x >>  4 )
--00,115831815525121700761 00769C9D x-= ( x >>  3 ) - ( x >>  6 )
--00,060380510988907482028 003DD463 x-= ( x >>  4 ) - ( x >>  8 )
--00,031748698314580298119 002082BB x-=               ( x >>  5 )
--00,015996403602177928366 0010615C x-= ( x >>  6 ) + ( x >> 12 )
--00,008089270731616965762 0008488D x-= ( x >>  7 ) + ( x >> 12 )
--00,004159027401785260480 00044243 x-= ( x >>  8 ) + ( x >> 12 )
--00,003423823349553609900 00038188 x-= ( x >>  8 ) - ( x >> 11 )
--00,001832733117311761669 0001E070 x-= ( x >>  9 ) - ( x >> 13 )
--00,000961766059678620302 0000FC1F x-= ( x >> 10 ) - ( x >> 16 )
--00,000484583944571210926 00007F07 x-= ( x >> 11 ) - ( x >> 18 )
--00,000243216525424976433 00003FC1 x-= ( x >> 12 ) - ( x >> 20 )
+-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
 
@@ -92,53 +96,60 @@ for negative values:
 static int32_t fixlog( int32_t x ) {
   int32_t t,y;
 
-  y = 0x2996BD9E; // ln(2^15) << 26
-  if(x<0x00008000) x<<=16,               y-= 0x2C5C85FD;
-  if(x<0x00800000) x<<= 8,               y-= 0x162E42FE;
-  if(x<0x08000000) x<<= 4,               y-= 0x0B17217F;
-  if(x<0x20000000) x<<= 2,               y-= 0x058B90BF;
-  if(x<0x40000000) x<<= 1,               y-= 0x02C5C85F;
-  t=x+(x>>1); if((t&0x80000000)==0) x=t, y-= 0x019F323E;
-  t=x+(x>>2); if((t&0x80000000)==0) x=t, y-= 0x00E47FBE;
-  t=x+(x>>3); if((t&0x80000000)==0) x=t, y-= 0x00789C1D;
-  t=x+(x>>4); if((t&0x80000000)==0) x=t, y-= 0x003E1461;
-  t=x+(x>>5); if((t&0x80000000)==0) x=t, y-= 0x001F829B;
-  t=x+(x>>6); if((t&0x80000000)==0) x=t, y-= 0x000FE054;
-  t=x+(x>>7); if((t&0x80000000)==0) x=t, y-= 0x0007F80A;
-  t=x+(x>>8); if((t&0x80000000)==0) x=t, y-= 0x0003FE01;
-  t=x+(x>>9); if((t&0x80000000)==0) x=t, y-= 0x0001FF80;
-  t=x+(x>>10);if((t&0x80000000)==0) x=t, y-= 0x0000FFE0;
-  t=x+(x>>11);if((t&0x80000000)==0) x=t, y-= 0x00007FF8;
-  t=x+(x>>12);if((t&0x80000000)==0) x=t, y-= 0x00003FFE;
-  t=x+(x>>13);if((t&0x80000000)==0) x=t, y-= 0x00001FFF;
+  // 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 = 0x80000000-x;
-  y-= x>>5;
-  return y;
+
+  /* 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 & 0x80000000) == 0 ) {
-      t=x-0x162E42FE; if(t>=0) x=t,y<<=8;
-      t=x-0x0B17217F; if(t>=0) x=t,y<<=4;
-      t=x-0x058B90BF; if(t>=0) x=t,y<<=2;
-      t=x-0x02C5C85F; if(t>=0) x=t,y<<=1;
-      t=x-0x019F323E; if(t>=0) x=t,y+=y>>1;
-      t=x-0x00E47FBE; if(t>=0) x=t,y+=y>>2;
-      t=x-0x00789C1D; if(t>=0) x=t,y+=y>>3;
-      t=x-0x003E1461; if(t>=0) x=t,y+=y>>4;
-      t=x-0x001F829B; if(t>=0) x=t,y+=y>>5;
-      t=x-0x000FE054; if(t>=0) x=t,y+=y>>6;
-      t=x-0x0007F80A; if(t>=0) x=t,y+=y>>7;
-      t=x-0x0003FE01; if(t>=0) x=t,y+=y>>8;
-      t=x-0x0001FF80; if(t>=0) x=t,y+=y>>9;
-      t=x-0x0000FFE0; if(t>=0) x=t,y+=y>>10;
-      t=x-0x00007FF8; if(t>=0) x=t,y+=y>>11;
-      t=x-0x00003FFE; if(t>=0) x=t,y+=y>>12;
-      t=x-0x00001FFF; if(t>=0) x=t,y+=y>>13;
+  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;
@@ -192,28 +203,34 @@ static uint32_t fixexp(int32_t x) {
 
 int fixpow(int32_t *pow, int32_t base, int32_t exponent )
 {
-  int neg = 0;
   int32_t log_base, res;
   int64_t log_base_times_exponent;
+  int neg = 0;
 
   if( base < 0 ) {
-     // negative bases only can have integer exponents
+     /* 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
+    /* 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
+  /* 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;
 
-  // fixexp overflows for values > 21,48756259689264 which
-  // in Q26 notation is 1442005916
+#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) )
+    return -2;
 
   res = (int32_t)log_base_times_exponent;
   res = fixexp( res );
@@ -225,8 +242,8 @@ int fixpow(int32_t *pow, int32_t base, int32_t exponent )
 
 int main()
 {
-  double base     = -.5f;
-  double exponent = 10.f;
+  double base     = -.5;
+  double exponent = 10.;
   int32_t result;
   int error = fixpow( &result, TO_Q16(base), TO_Q16(exponent) );