Initial version. Works. Untested. Unreviewed

2 files changed, 269 insertions, 0 deletions

diff --git a/fixpow.c b/fixpow.c
new file mode 100644
index 0000000..be787b2
--- /dev/null
+++ b/fixpow.c
@@ -0,0 +1,235 @@
+/* Parts of this software was written by Dirk Engling <erdgeist@erdgeist.org>
+   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 <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#define TO_Q16(X)    ((int32_t)(65536.f*(X)))
+#define FROM_Q16(X)  ((double)((X)/65536.f))
+
+// 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 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 )
+
+#endif
+
+// input is Q16, outputQ26
+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;
+  x = 0x80000000-x;
+  y-= x>>5;
+  return y;
+}
+
+// input is Q26, output Q16
+static uint32_t fixexp(int32_t x) {
+  int32_t t;
+  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&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 {
+      x=-x;
+      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-0x017FAFA3; if(t>=0) x=t,y-=(y>>2) + (y>>4);
+      t=x-0x00D49F69; if(t>=0) x=t,y-=(y>>2) - (y>>4);
+      t=x-0x00769C9D; if(t>=0) x=t,y-=(y>>3) - (y>>6);
+      t=x-0x003DD463; if(t>=0) x=t,y-=(y>>4) - (y>>8);
+      t=x-0x002082BB; if(t>=0) x=t,y-=y>>5;
+      t=x-0x0010615C; if(t>=0) x=t,y-=(y>>6) + (y>>12);
+      t=x-0x0008488D; if(t>=0) x=t,y-=(y>>7) + (y>>12);
+      t=x-0x00044243; if(t>=0) x=t,y-=(y>>8) + (y>>12);
+      t=x-0x00038188; if(t>=0) x=t,y-=(y>>8) - (y>>11);
+      t=x-0x0001E070; if(t>=0) x=t,y-=(y>>9) - (y>>13);
+      t=x-0x0000FC1F; if(t>=0) x=t,y-=(y>>10)- (y>>16);
+      t=x-0x00007F07; if(t>=0) x=t,y-=(y>>11)- (y>>18);
+      t=x-0x00003FC1; if(t>=0) x=t,y-=(y>>12)- (y>>20);
+      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 )
+{
+  int neg = 0;
+  int32_t log_base, res;
+  int64_t log_base_times_exponent;
+
+  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;
+
+  // fixexp overflows for values > 21,48756259689264 which
+  // in Q26 notation is 1442005916
+  if( log_base_times_exponent > 1442005916 )
+    return -2;
+
+  res = (int32_t)log_base_times_exponent;
+  res = fixexp( res );
+  if( neg ) res = -res;
+  *pow = res;
+
+  return 0;
+}
+
+int main()
+{
+  double base     = -.5f;
+  double exponent = 10.f;
+  int32_t result;
+  int error = fixpow( &result, TO_Q16(base), TO_Q16(exponent) );
+
+  printf( "pow(%lf,%lf)=%lf (%s)\n", base, exponent, FROM_Q16(result), error?"ERROR":"OK" );
+  return 0;
+}
diff --git a/generate_table.c b/generate_table.c
new file mode 100644
index 0000000..c7a9d0c
--- /dev/null
+++ b/generate_table.c
@@ -0,0 +1,34 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+// s(-1) n(14) off(+1)   x =   x - ( x >> 14 )
+// s( 1) n(-4) off( 0)   x =       ( x <<  4 )
+// s( 1) n( 8) off(-1)   x = - x + ( x <<  8 )
+
+int main() {
+  int n, m, nsign;
+
+  // loop over all constructs in the form +-1*x^+-2n, n=-31..+31
+  for (n=-31; n<= 31; ++n )
+    for ( nsign=-1; nsign<=1; nsign+=2 )
+    {
+      // The one term only case
+      double v = (double)nsign * pow( 2.f, (double)n );
+      if( v > 0.f )
+          printf( "%0+25.21lf %08X x =             %s ( x %s %2d )%s\n", log(v), abs((int)(67108864.f*log(v))), nsign==-1?"-":" ", n>=0?"<<":">>", (int)abs(n), " !RECOMMENDED!" );
+
+      // Loop over second term
+      for (m=-31; m<=31; ++m )
+      {
+        double v = pow( 2.f, (double)m ) + (double)nsign * pow( 2.f, (double)n );
+        if( v > 0.f ) {
+          printf( "%0+25.21lf %08X x = ( x %s %2d ) %s ( x %s %2d )%s\n", log(v), abs((int)(67108864.f*log(v))), m>=0?"<<":">>", (int)abs(m), nsign==-1?"-":"+", n>=0?"<<":">>", (int)abs(n), n*m?"":" !RECOMMENDED!" );
+          if( v < 1.f )
+              printf( "%0+25.21lf %08X x-= ( x %s %2d ) %s ( x %s %2d )\n", log(1.0f - v), abs((int)(67108864.f*log(1.0f - v))), m>=0?"<<":">>", (int)abs(m), nsign==-1?"-":"+", n>=0?"<<":">>", (int)abs(n) );
+        }
+      }
+    }
+
+  return 0;
+}