summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--fixpow.c235
-rw-r--r--generate_table.c34
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 @@
1/* Parts of this software was written by Dirk Engling <erdgeist@erdgeist.org>
2 Those parts are considered beerware. Prost. Skol. Cheers or whatever.
3
4 Original idea and the place where I heavily stole code from is
5 http://www.quinapalus.com/efunc.html
6 Dr Mark St John OWEN, E-mail: mail -at- quinapalus -dot- com
7*/
8
9#include <stdint.h>
10#include <stdio.h>
11#include <stdlib.h>
12#include <math.h>
13
14#define TO_Q16(X) ((int32_t)(65536.f*(X)))
15#define FROM_Q16(X) ((double)((X)/65536.f))
16
17// It works as follows:
18// P = pow_QX_QY_QP( x, y ) == exp( ln(x) * y )
19
20// Where if x is negative:
21// * y must have no fractional part,
22// * the result's sign is the lowest integral bit of y
23
24// We note that in every standard Q representation ln(x) will not exceed
25// the value 22, so that we can safely work with a Q26 representation.
26//
27// if ln(x) * y would overflow in that representation, so would
28// exp( ln(x) * y ) in Q16.
29//
30// Finally exp() is calculated in the domain specified by script
31
32#if 0
33 The table[tm], generated from different versions of
34 generate_table, though.
35+00,000000000465661287199 x = x + ( x >> 31 )
36+00,000000000931322574182 x = x + ( x >> 30 )
37+00,000000001862645147496 x = x + ( x >> 29 )
38+00,000000003725290291523 x = x + ( x >> 28 )
39+00,000000007450580569168 x = x + ( x >> 27 )
40+00,000000014901161082825 x = x + ( x >> 26 )
41+00,000000029802321943606 x = x + ( x >> 25 )
42+00,000000059604642999034 x = x + ( x >> 24 )
43+00,000000119209282445354 x = x + ( x >> 23 )
44+00,000000238418550679858 x = x + ( x >> 22 )
45+00,000000476837044516323 x = x + ( x >> 21 )
46+00,000000953673861659188 x = x + ( x >> 20 )
47+00,000001907346813825409 x = x + ( x >> 19 )
48+00,000003814689989685890 x = x + ( x >> 18 )
49+00,000007629365427567572 x = x + ( x >> 17 )
50+00,000015258672648362398 x = x + ( x >> 16 )
51+00,000030517112473186380 x = x + ( x >> 15 )
52+00,000061033293680638527 x = x + ( x >> 14 )
53+00,000122062862525677371 x = x + ( x >> 13 )
54+00,000244110827527362707 x = x + ( x >> 12 )
55+00,000488162079501351187 x = x + ( x >> 11 )
56+00,000976085973055458925 x = x + ( x >> 10 )
57+00,001951220131261749337 x = x + ( x >> 9 )
58+00,003898640415657322889 x = x + ( x >> 8 )
59+00,007782140442054948960 x = x + ( x >> 7 )
60+00,015504186535965254479 x = x + ( x >> 6 )
61+00,030771658666753687328 x = x + ( x >> 5 )
62+00,060624621816434839938 x = x + ( x >> 4 )
63+00,117783035656383455736 x = x + ( x >> 3 )
64+00,223143551314209764858 x = x + ( x >> 2 )
65+00,405465108108164384859 x = x + ( x >> 1 )
66+00,693147180559945286227 x = ( x << 1 )
67+01,386294361119890572454 x = ( x << 2 )
68+02,772588722239781144907 x = ( x << 4 )
69+05,545177444479562289814 x = ( x << 8 )
70+11,090354888959124579628 x = ( x << 16 )
71+21,487562597358305538364 x = ( x << 31 )
72
73for negative values:
74
75-00,374693449441410697531 017FAFA3 x-= ( x >> 2 ) + ( x >> 4 )
76-00,207639364778244489562 00D49F69 x-= ( x >> 2 ) - ( x >> 4 )
77-00,115831815525121700761 00769C9D x-= ( x >> 3 ) - ( x >> 6 )
78-00,060380510988907482028 003DD463 x-= ( x >> 4 ) - ( x >> 8 )
79-00,031748698314580298119 002082BB x-= ( x >> 5 )
80-00,015996403602177928366 0010615C x-= ( x >> 6 ) + ( x >> 12 )
81-00,008089270731616965762 0008488D x-= ( x >> 7 ) + ( x >> 12 )
82-00,004159027401785260480 00044243 x-= ( x >> 8 ) + ( x >> 12 )
83-00,003423823349553609900 00038188 x-= ( x >> 8 ) - ( x >> 11 )
84-00,001832733117311761669 0001E070 x-= ( x >> 9 ) - ( x >> 13 )
85-00,000961766059678620302 0000FC1F x-= ( x >> 10 ) - ( x >> 16 )
86-00,000484583944571210926 00007F07 x-= ( x >> 11 ) - ( x >> 18 )
87-00,000243216525424976433 00003FC1 x-= ( x >> 12 ) - ( x >> 20 )
88
89#endif
90
91// input is Q16, outputQ26
92static int32_t fixlog( int32_t x ) {
93 int32_t t,y;
94
95 y = 0x2996BD9E; // ln(2^15) << 26
96 if(x<0x00008000) x<<=16, y-= 0x2C5C85FD;
97 if(x<0x00800000) x<<= 8, y-= 0x162E42FE;
98 if(x<0x08000000) x<<= 4, y-= 0x0B17217F;
99 if(x<0x20000000) x<<= 2, y-= 0x058B90BF;
100 if(x<0x40000000) x<<= 1, y-= 0x02C5C85F;
101 t=x+(x>>1); if((t&0x80000000)==0) x=t, y-= 0x019F323E;
102 t=x+(x>>2); if((t&0x80000000)==0) x=t, y-= 0x00E47FBE;
103 t=x+(x>>3); if((t&0x80000000)==0) x=t, y-= 0x00789C1D;
104 t=x+(x>>4); if((t&0x80000000)==0) x=t, y-= 0x003E1461;
105 t=x+(x>>5); if((t&0x80000000)==0) x=t, y-= 0x001F829B;
106 t=x+(x>>6); if((t&0x80000000)==0) x=t, y-= 0x000FE054;
107 t=x+(x>>7); if((t&0x80000000)==0) x=t, y-= 0x0007F80A;
108 t=x+(x>>8); if((t&0x80000000)==0) x=t, y-= 0x0003FE01;
109 t=x+(x>>9); if((t&0x80000000)==0) x=t, y-= 0x0001FF80;
110 t=x+(x>>10);if((t&0x80000000)==0) x=t, y-= 0x0000FFE0;
111 t=x+(x>>11);if((t&0x80000000)==0) x=t, y-= 0x00007FF8;
112 t=x+(x>>12);if((t&0x80000000)==0) x=t, y-= 0x00003FFE;
113 t=x+(x>>13);if((t&0x80000000)==0) x=t, y-= 0x00001FFF;
114 x = 0x80000000-x;
115 y-= x>>5;
116 return y;
117}
118
119// input is Q26, output Q16
120static uint32_t fixexp(int32_t x) {
121 int32_t t;
122 uint32_t y = 0x10000;
123
124 if( (x & 0x80000000) == 0 ) {
125 t=x-0x162E42FE; if(t>=0) x=t,y<<=8;
126 t=x-0x0B17217F; if(t>=0) x=t,y<<=4;
127 t=x-0x058B90BF; if(t>=0) x=t,y<<=2;
128 t=x-0x02C5C85F; if(t>=0) x=t,y<<=1;
129 t=x-0x019F323E; if(t>=0) x=t,y+=y>>1;
130 t=x-0x00E47FBE; if(t>=0) x=t,y+=y>>2;
131 t=x-0x00789C1D; if(t>=0) x=t,y+=y>>3;
132 t=x-0x003E1461; if(t>=0) x=t,y+=y>>4;
133 t=x-0x001F829B; if(t>=0) x=t,y+=y>>5;
134 t=x-0x000FE054; if(t>=0) x=t,y+=y>>6;
135 t=x-0x0007F80A; if(t>=0) x=t,y+=y>>7;
136 t=x-0x0003FE01; if(t>=0) x=t,y+=y>>8;
137 t=x-0x0001FF80; if(t>=0) x=t,y+=y>>9;
138 t=x-0x0000FFE0; if(t>=0) x=t,y+=y>>10;
139 t=x-0x00007FF8; if(t>=0) x=t,y+=y>>11;
140 t=x-0x00003FFE; if(t>=0) x=t,y+=y>>12;
141 t=x-0x00001FFF; if(t>=0) x=t,y+=y>>13;
142 if(x&0x0001000) y+=y>>14;
143 if(x&0x0000800) y+=y>>15;
144 if(x&0x0000400) y+=y>>16;
145 if(x&0x0000200) y+=y>>17;
146 if(x&0x0000100) y+=y>>18;
147 if(x&0x0000080) y+=y>>19;
148 if(x&0x0000040) y+=y>>20;
149 if(x&0x0000020) y+=y>>21;
150 if(x&0x0000010) y+=y>>22;
151 if(x&0x0000008) y+=y>>23;
152 if(x&0x0000004) y+=y>>24;
153 if(x&0x0000002) y+=y>>25;
154 if(x&0x0000001) y+=y>>26;
155 } else {
156 x=-x;
157 t=x-0x162E42FE; if(t>=0) x=t,y>>=8;
158 t=x-0x0B17217F; if(t>=0) x=t,y>>=4;
159 t=x-0x058B90BF; if(t>=0) x=t,y>>=2;
160 t=x-0x02C5C85F; if(t>=0) x=t,y>>=1;
161 t=x-0x017FAFA3; if(t>=0) x=t,y-=(y>>2) + (y>>4);
162 t=x-0x00D49F69; if(t>=0) x=t,y-=(y>>2) - (y>>4);
163 t=x-0x00769C9D; if(t>=0) x=t,y-=(y>>3) - (y>>6);
164 t=x-0x003DD463; if(t>=0) x=t,y-=(y>>4) - (y>>8);
165 t=x-0x002082BB; if(t>=0) x=t,y-=y>>5;
166 t=x-0x0010615C; if(t>=0) x=t,y-=(y>>6) + (y>>12);
167 t=x-0x0008488D; if(t>=0) x=t,y-=(y>>7) + (y>>12);
168 t=x-0x00044243; if(t>=0) x=t,y-=(y>>8) + (y>>12);
169 t=x-0x00038188; if(t>=0) x=t,y-=(y>>8) - (y>>11);
170 t=x-0x0001E070; if(t>=0) x=t,y-=(y>>9) - (y>>13);
171 t=x-0x0000FC1F; if(t>=0) x=t,y-=(y>>10)- (y>>16);
172 t=x-0x00007F07; if(t>=0) x=t,y-=(y>>11)- (y>>18);
173 t=x-0x00003FC1; if(t>=0) x=t,y-=(y>>12)- (y>>20);
174 if(x&0x0002000) y-=y>>13;
175 if(x&0x0001000) y-=y>>14;
176 if(x&0x0000800) y-=y>>15;
177 if(x&0x0000400) y-=y>>16;
178 if(x&0x0000200) y-=y>>17;
179 if(x&0x0000100) y-=y>>18;
180 if(x&0x0000080) y-=y>>19;
181 if(x&0x0000040) y-=y>>20;
182 if(x&0x0000020) y-=y>>21;
183 if(x&0x0000010) y-=y>>22;
184 if(x&0x0000008) y-=y>>23;
185 if(x&0x0000004) y-=y>>24;
186 if(x&0x0000002) y-=y>>25;
187 if(x&0x0000001) y-=y>>26;
188 }
189
190 return y;
191}
192
193int fixpow(int32_t *pow, int32_t base, int32_t exponent )
194{
195 int neg = 0;
196 int32_t log_base, res;
197 int64_t log_base_times_exponent;
198
199 if( base < 0 ) {
200 // negative bases only can have integer exponents
201 if( exponent & 0xffff ) return -1;
202 // sign of power of a negative base is determined by
203 // wether exp is even
204 if( exponent & 0x10000 ) neg = 1;
205 base = abs(base);
206 }
207
208 // To calculate pow(base,exp), we do exp( log(base) * exp )
209 // which is mathematically the same. log_base is Q26
210 log_base = fixlog( base );
211 log_base_times_exponent = ( (int64_t)log_base * (int64_t)exponent ) >> 16;
212
213 // fixexp overflows for values > 21,48756259689264 which
214 // in Q26 notation is 1442005916
215 if( log_base_times_exponent > 1442005916 )
216 return -2;
217
218 res = (int32_t)log_base_times_exponent;
219 res = fixexp( res );
220 if( neg ) res = -res;
221 *pow = res;
222
223 return 0;
224}
225
226int main()
227{
228 double base = -.5f;
229 double exponent = 10.f;
230 int32_t result;
231 int error = fixpow( &result, TO_Q16(base), TO_Q16(exponent) );
232
233 printf( "pow(%lf,%lf)=%lf (%s)\n", base, exponent, FROM_Q16(result), error?"ERROR":"OK" );
234 return 0;
235}
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 @@
1#include <math.h>
2#include <stdio.h>
3#include <stdlib.h>
4
5// s(-1) n(14) off(+1) x = x - ( x >> 14 )
6// s( 1) n(-4) off( 0) x = ( x << 4 )
7// s( 1) n( 8) off(-1) x = - x + ( x << 8 )
8
9int main() {
10 int n, m, nsign;
11
12 // loop over all constructs in the form +-1*x^+-2n, n=-31..+31
13 for (n=-31; n<= 31; ++n )
14 for ( nsign=-1; nsign<=1; nsign+=2 )
15 {
16 // The one term only case
17 double v = (double)nsign * pow( 2.f, (double)n );
18 if( v > 0.f )
19 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!" );
20
21 // Loop over second term
22 for (m=-31; m<=31; ++m )
23 {
24 double v = pow( 2.f, (double)m ) + (double)nsign * pow( 2.f, (double)n );
25 if( v > 0.f ) {
26 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!" );
27 if( v < 1.f )
28 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) );
29 }
30 }
31 }
32
33 return 0;
34}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-21 14:46:27 +0000