diff options
-rw-r--r-- | fixpow.c | 61 |
1 files changed, 38 insertions, 23 deletions
@@ -14,7 +14,7 @@ | |||
14 | #define TO_Q16(X) ((int32_t)(65536.*(X))) | 14 | #define TO_Q16(X) ((int32_t)(65536.*(X))) |
15 | #define TO_Q26(X) ((int32_t)(67108864.*(X))) | 15 | #define TO_Q26(X) ((int32_t)(67108864.*(X))) |
16 | #define FROM_Q16(X) ((double)((X)/65536.)) | 16 | #define FROM_Q16(X) ((double)((X)/65536.)) |
17 | #define FROM_Q26(X) ((double)((x)/67108864.)) | 17 | #define FROM_Q26(X) ((double)((X)/67108864.)) |
18 | 18 | ||
19 | /* | 19 | /* |
20 | It works as follows: | 20 | It works as follows: |
@@ -118,7 +118,7 @@ static int32_t fixlog( int32_t x ) { | |||
118 | t=x+(x>>13);if((t&0x80000000)==0) x=t, y-= TO_Q26(0.000122062862526); /* log(1+1/8192) */ | 118 | t=x+(x>>13);if((t&0x80000000)==0) x=t, y-= TO_Q26(0.000122062862526); /* log(1+1/8192) */ |
119 | 119 | ||
120 | /* Estimate residual error, log(1-x) which for small x is approx -x */ | 120 | /* Estimate residual error, log(1-x) which for small x is approx -x */ |
121 | x = 0x80000000-x; | 121 | x = 2147483648-x; |
122 | 122 | ||
123 | /* x has been promoted to Q31, substract residual error in Q26 by shifting down by 5 */ | 123 | /* x has been promoted to Q31, substract residual error in Q26 by shifting down by 5 */ |
124 | return y - ( x >> 5 ); | 124 | return y - ( x >> 5 ); |
@@ -164,24 +164,25 @@ static uint32_t fixexp(int32_t x) { | |||
164 | if(x&0x0000002) y+=y>>25; | 164 | if(x&0x0000002) y+=y>>25; |
165 | if(x&0x0000001) y+=y>>26; | 165 | if(x&0x0000001) y+=y>>26; |
166 | } else { | 166 | } else { |
167 | x=-x; | 167 | if(x<=TO_Q26(-11.09035488895912)) return 0; |
168 | t=x-TO_Q26(5.545177444479562); if(t>=0) x=t,y>>=8; | 168 | t=x+TO_Q26(11.09035488895912); if(t<=0) x=t,y>>=16; |
169 | t=x-TO_Q26(2.772588722239781); if(t>=0) x=t,y>>=4; | 169 | t=x+TO_Q26(5.545177444479562); if(t<=0) x=t,y>>=8; |
170 | t=x-TO_Q26(1.386294361119891); if(t>=0) x=t,y>>=2; | 170 | t=x+TO_Q26(2.772588722239781); if(t<=0) x=t,y>>=4; |
171 | t=x-TO_Q26(0.693147180559945); if(t>=0) x=t,y>>=1; | 171 | t=x+TO_Q26(1.386294361119891); if(t<=0) x=t,y>>=2; |
172 | t=x-TO_Q26(0.374693449441411); if(t>=0) x=t,y-=(y>>4) + (y>>2); | 172 | t=x+TO_Q26(0.693147180559945); if(t<=0) x=t,y>>=1; |
173 | t=x-TO_Q26(0.207639364778244); if(t>=0) x=t,y-=(y>>2) - (y>>4); | 173 | t=x+TO_Q26(0.374693449441411); if(t<=0) x=t,y-=(y>>4) + (y>>2); |
174 | t=x-TO_Q26(0.124642445207276); if(t>=0) x=t,y-=(y>>3) - (y>>7); | 174 | t=x+TO_Q26(0.207639364778244); if(t<=0) x=t,y-=(y>>2) - (y>>4); |
175 | t=x-TO_Q26(0.062457354933746); if(t>=0) x=t,y-=(y>>4) - (y>>9); | 175 | t=x+TO_Q26(0.124642445207276); if(t<=0) x=t,y-=(y>>3) - (y>>7); |
176 | t=x-TO_Q26(0.031244793038107); if(t>=0) x=t,y-=(y>>5) - (x>>11); | 176 | t=x+TO_Q26(0.062457354933746); if(t<=0) x=t,y-=(y>>4) - (y>>9); |
177 | t=x-TO_Q26(0.015748356968139); if(t>=0) x=t,y-=(y>>6); | 177 | t=x+TO_Q26(0.031244793038107); if(t<=0) x=t,y-=(y>>5) - (y>>11); |
178 | t=x-TO_Q26(0.007966216526084); if(t>=0) x=t,y-=(y>>7) + (y>>13); | 178 | t=x+TO_Q26(0.015748356968139); if(t<=0) x=t,y-=(y>>6); |
179 | t=x-TO_Q26(0.004036455850488); if(t>=0) x=t,y-=(y>>8) + (y>>13); | 179 | t=x+TO_Q26(0.007966216526084); if(t<=0) x=t,y-=(y>>7) + (y>>13); |
180 | t=x-TO_Q26(0.002077351513834); if(t>=0) x=t,y-=(y>>9) + (y>>13); | 180 | t=x+TO_Q26(0.004036455850488); if(t<=0) x=t,y-=(y>>8) + (y>>13); |
181 | t=x-TO_Q26(0.001099236751907); if(t>=0) x=t,y-=(y>>10)+ (y>>13); | 181 | t=x+TO_Q26(0.002077351513834); if(t<=0) x=t,y-=(y>>9) + (y>>13); |
182 | t=x-TO_Q26(0.000610537902840); if(t>=0) x=t,y-=(y>>11)+ (y>>13); | 182 | t=x+TO_Q26(0.001099236751907); if(t<=0) x=t,y-=(y>>10)+ (y>>13); |
183 | t=x-TO_Q26(0.000366278009100); if(t>=0) x=t,y-=(y>>12)+ (y>>13); | 183 | t=x+TO_Q26(0.000610537902840); if(t<=0) x=t,y-=(y>>11)+ (y>>13); |
184 | t=x-TO_Q26(0.000213645867528); if(t>=0) x=t,y-=(y>>12)- (y>>15); | 184 | t=x+TO_Q26(0.000366278009100); if(t<=0) x=t,y-=(y>>12)+ (y>>13); |
185 | t=x+TO_Q26(0.000213645867528); if(t<=0) x=t,y-=(y>>12)- (y>>15); | ||
185 | if(x&0x0002000) y-=y>>13; | 186 | if(x&0x0002000) y-=y>>13; |
186 | if(x&0x0001000) y-=y>>14; | 187 | if(x&0x0001000) y-=y>>14; |
187 | if(x&0x0000800) y-=y>>15; | 188 | if(x&0x0000800) y-=y>>15; |
@@ -221,12 +222,10 @@ int fixpow(int32_t *pow, int32_t base, int32_t exponent ) | |||
221 | log_base = fixlog( base ); | 222 | log_base = fixlog( base ); |
222 | log_base_times_exponent = ( (int64_t)log_base * (int64_t)exponent ) >> 16; | 223 | log_base_times_exponent = ( (int64_t)log_base * (int64_t)exponent ) >> 16; |
223 | 224 | ||
224 | #if 0 | ||
225 | /* In Q0, fixexp overflows for values > 21,48756259689264 which | 225 | /* In Q0, fixexp overflows for values > 21,48756259689264 which |
226 | in Q26 notation is 1442005916 */ | 226 | in Q26 notation is 1442005916 */ |
227 | if( log_base_times_exponent > 1442005916 ) | 227 | if( log_base_times_exponent > 1442005916 ) |
228 | return -2; | 228 | return -2; |
229 | #endif | ||
230 | 229 | ||
231 | /* In Q16, fixexp overflows for values > 10,39720770839918 */ | 230 | /* In Q16, fixexp overflows for values > 10,39720770839918 */ |
232 | if( log_base_times_exponent > TO_Q26(10.39720770839918) ) | 231 | if( log_base_times_exponent > TO_Q26(10.39720770839918) ) |
@@ -244,9 +243,25 @@ int main() | |||
244 | { | 243 | { |
245 | double base = -.5; | 244 | double base = -.5; |
246 | double exponent = 10.; | 245 | double exponent = 10.; |
247 | int32_t result; | 246 | int32_t result, input, output; |
248 | int error = fixpow( &result, TO_Q16(base), TO_Q16(exponent) ); | 247 | int error = fixpow( &result, TO_Q16(base), TO_Q16(exponent) ); |
249 | 248 | ||
249 | // Testing fixlog | ||
250 | for (input=1; input<=INT32_MAX/2; ++input) { | ||
251 | int32_t expected = TO_Q26(log(FROM_Q16(input))); | ||
252 | int32_t output = fixlog(input); | ||
253 | if( abs(output-expected) > 8 ) | ||
254 | printf( "%08X : %08X = %08X (%08X)\n", output, expected, abs(output-expected), input ); | ||
255 | } | ||
256 | |||
257 | // Testing fixexp | ||
258 | for (input=-2147483648; input < TO_Q26(10.39720770839918); ++input ) { | ||
259 | int32_t expected = TO_Q16(exp(FROM_Q26(input))); | ||
260 | int32_t output = fixexp(input); | ||
261 | if( abs(output-expected) > 16 ) | ||
262 | printf( "%08X : %08X = %08X (%08d)\n", output, expected, abs(output-expected), input ); | ||
263 | } | ||
264 | |||
250 | printf( "pow(%lf,%lf)=%lf (%s)\n", base, exponent, FROM_Q16(result), error?"ERROR":"OK" ); | 265 | printf( "pow(%lf,%lf)=%lf (%s)\n", base, exponent, FROM_Q16(result), error?"ERROR":"OK" ); |
251 | return 0; | 266 | return 0; |
252 | } | 267 | } |