1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
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;
}
|
