From 2e57ec4405b2c52b48820aff75cd6d18efe2fc47 Mon Sep 17 00:00:00 2001 From: erdgeist Date: Mon, 6 Oct 2014 18:04:43 +0200 Subject: Convert hex constants to double constants and a macro to convert them to Q26 --- fixpow.c | 169 +++++++++++++++++++++++++++++++++++---------------------------- 1 file 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 #include -#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) ); -- cgit v1.2.3