1 // -*- compile-command: "/usr/local/wince/bin/arm-wince-pe-g++ softmath.cc -march=armv4 -mapcs-32 -malignment-traps -msoft-float -DNEWLIB -DSARM -DWIN32 -DGNUWINCE -DLINK -lm /usr/local/wince/arm-wince-pe/lib/libc.a -lgcc -lwinsock -lcoredll" -*-
2
3 // Transcendental functions implementation (not optimized at all)
4 // Use -DLINK to test or -c for object creation
5 #include "softmath.h"
6 #include <iostream>
7 #include <fstream>
8 #include <iomanip>
9 #include <cmath>
10 #include <cstdlib>
11
12 namespace std {
13 double m_ln2=0.69314718055994530942;
14 double m_twopi=2.0*M_PI;
15 double m_sqrt3=1.73205080756887729;
16
17 double m_undef=0.0/0.0;
18 double m_plusinf=1.0/0.0;
19 double m_minusinf=-1.0/0.0;
20
21 // 14!(14-k)!
22 double fact14[]={1.0,14.0,182.0,2184.0,24024.0,240240.0,2162160.0,17297280.0,121080960.0,726485760.0,3632428800.0,14529715200.0,43589145600.0,87178291200.0};
23
24 // k is an integer >=0 and < 1024
giac_gnuwince_exp2(double k)25 double giac_gnuwince_exp2(double k){
26 if (k<31.0){
27 unsigned i=0x1 << int(k);
28 return i;
29 }
30 double k1=int(k/32.0+0.001); // k1<32
31 double res=1.0;
32 k -= k1*32.0; // k<32
33 for (;k>0.0;k -=1.0)
34 res *= 2.0;
35 for (;k1>0.0; k1 -= 1.0)
36 res *= 4294967296.0;
37 return res;
38 }
39
giac_gnuwince_exp(double d)40 double giac_gnuwince_exp(double d){
41 if (d<0)
42 return 1.0/giac_gnuwince_exp(-d);
43 #ifdef DEBUG
44 ofstream of("int.txt");
45 #endif
46 double k=int(d/m_ln2+.5);
47 #ifdef DEBUG
48 of << "k:" << k << char(10) << char(13) ;
49 #endif
50 if (k>1023.0)
51 return m_plusinf;
52 d=d-k*m_ln2; // |d|<m_ln2/2
53 #ifdef DEBUG
54 of << "d:" << d << char(10) << char(13) ;
55 #endif
56 if (k<0)
57 k=1.0/giac_gnuwince_exp2(-k);
58 else
59 k=giac_gnuwince_exp2(k);
60 #ifdef DEBUG
61 of << "2^k:" << k << char(10) << char(13) ;
62 #endif
63 // Insure d is < 0, if it was > 0 compute inv(exp(-d))
64 bool inv=d>0;
65 if (inv)
66 d=-d;
67 // use Taylor expansion at order 14 since
68 // (1+(ln(2.0)/2)^14/14!)-1 -> 0.0
69 double res=d;
70 for (int i=1;i<14;++i){
71 res += fact14[i];
72 res *= d;
73 #ifdef DEBUG
74 of << "i:" << i << " ,res:" << res << char(10) << char(13) ;
75 #endif
76 }
77 res /= 87178291200.0;
78 res += 1.0;
79 #ifdef DEBUG
80 of << "res:" << res << char(10) << char(13) ;
81 #endif
82 if (inv)
83 return k/res;
84 else
85 return res*k;
86 }
87
giac_gnuwince_sinh(double d)88 double giac_gnuwince_sinh(double d){
89 double k=giac_gnuwince_exp(d);
90 return k-1.0/k;
91 }
92
giac_gnuwince_cosh(double d)93 double giac_gnuwince_cosh(double d){
94 double k=giac_gnuwince_exp(d);
95 return k+1.0/k;
96 }
97
giac_gnuwince_tanh(double d)98 double giac_gnuwince_tanh(double d){
99 if (d>700.0)
100 return 1.0;
101 if (d<-700.0)
102 return -1.0;
103 double k=giac_gnuwince_exp(d);
104 k=k*k;
105 return (k-1)/(k+1);
106 }
107
108 // 0<d<=pi/4
in_giac_gnuwince_tan(double d)109 double in_giac_gnuwince_tan(double d){
110 if (d>0.07){
111 // tan(2*x)=2*tan(x)/(1-tan(x)^2)
112 double tandsur2=in_giac_gnuwince_tan(d/2);
113 return 2.0*tandsur2/(1.0-tandsur2*tandsur2);
114 }
115 // d*(1382*d^2^5+3410*d^2^4+8415*d^2^3+20790*d^2^2+51975*d^2+155925)/155925
116 double d2=d*d;
117 double res=((((1382.0*d2+3410.0)*d2+8415.0)*d2+20790.0)*d2+51975.0)*d2+155925.0;
118 return res*d/155925.0;
119 }
120
giac_gnuwince_tan(double d)121 double giac_gnuwince_tan(double d){
122 if (d<0)
123 return -giac_gnuwince_tan(-d);
124 double k=int(d/M_PI+.5);
125 d=d-k*M_PI; // |d|<pi/2
126 bool neg=d<0,inv=false;
127 if (neg)
128 d=-d;
129 if (d>M_PI_4){
130 d=M_PI_2-d;
131 inv=true;
132 }
133 // 0<d<=pi/4
134 k= in_giac_gnuwince_tan(d);
135 if (neg)
136 k=-k;
137 if (inv)
138 k=1.0/k;
139 return k;
140 }
141
giac_gnuwince_sin(double d)142 double giac_gnuwince_sin(double d){
143 if (d<0)
144 return -giac_gnuwince_sin(-d);
145 double k=int(d/m_twopi+.5);
146 d=d-k*m_twopi;
147 bool neg=d<0;
148 if (neg)
149 d=-d;
150 if (d>M_PI_2)
151 d=M_PI-d;
152 // now |d|<pi/2
153 k=in_giac_gnuwince_tan(d/2.0);
154 k=(2.0*k)/(k*k+1.0);
155 if (neg)
156 k=-k;
157 return k;
158 }
159
giac_gnuwince_cos(double d)160 double giac_gnuwince_cos(double d){
161 if (d<0)
162 return giac_gnuwince_cos(-d);
163 double k=int(d/m_twopi+.5);
164 d=d-k*m_twopi;
165 if (d<0)
166 d=-d;
167 bool neg=false;
168 if (d>M_PI_2){
169 neg=true;
170 d=M_PI-d;
171 }
172 // now |d|<pi/2, compute sin(pi/2-d)
173 k=in_giac_gnuwince_tan((M_PI_2-d)/2.0);
174 k=(2.0*k)/(k*k+1.0);
175 if (neg)
176 k=-k;
177 return k;
178 }
179
giac_gnuwince_floor(double d)180 double giac_gnuwince_floor(double d){
181 if (d>0)
182 return int(d);
183 double k=int(d);
184 if (k==d)
185 return k;
186 else
187 return k-1;
188 }
189
giac_gnuwince_ceil(double d)190 double giac_gnuwince_ceil(double d){
191 if (d<0)
192 return int(d);
193 double k=int(d);
194 if (k==d)
195 return k;
196 else
197 return k+1;
198 }
199
200 // Taylor expansion order 23, |d|<0.27
in_giac_gnuwince_atan(double d)201 double in_giac_gnuwince_atan(double d){
202 double d2=d*d;
203 double res=((((((((((-14549535.0*d2+15935205.0)*d2-17612595.0)*d2+19684665.0)*d2-22309287.0)*d2+25741485.0)*d2-30421755.0)*d2+37182145.0)*d2-47805615.0)*d2+66927861.0)*d2-111546435.0)*d2+334639305.0;
204 return res*d/334639305.0;
205 }
206
giac_gnuwince_atan(double d)207 double giac_gnuwince_atan(double d){
208 if (d<0)
209 return -giac_gnuwince_atan(-d);
210 bool inv=false;
211 if (d>1){
212 inv=true;
213 d=1.0/d;
214 }
215 double k;
216 if (d<0.265)
217 k=in_giac_gnuwince_atan(d);
218 else {
219 if (d>0.7)
220 k=M_PI_4-in_giac_gnuwince_atan((1.0-d)/(1.0+d));
221 else
222 k=M_PI/6.0+in_giac_gnuwince_atan((m_sqrt3*d-1)/(m_sqrt3+d));
223 }
224 if (inv)
225 k=M_PI_2-k;
226 return k;
227 }
228
229 double pow2tab[]={0.13407807929942597100e155,0.11579208923731619542e78,0.34028236692093846346e39,0.18446744073709551616e20,0.42949672960000000000e10,65536.0,256.0,16.0,4.0,2.0,1.4142135623730950488,1.1892071150027210667,1.0905077326652576592};
giac_gnuwince_sqrt(double d)230 double giac_gnuwince_sqrt(double d){
231 if (d<0)
232 return m_undef;
233 if (d==0 || d==1)
234 return d;
235 if (d<1)
236 return 1.0/giac_gnuwince_sqrt(1.0/d);
237 double k=0.0,k2=512.0;
238 for (int i=0;i<8;++i,k2 /= 2.0){
239 if (d>=pow2tab[i]){
240 k += k2;
241 d /= pow2tab[i];
242 }
243 }
244 double d0=(1.0+d)/2.0,d1;
245 for (;;d0=d1){
246 d1=(d0+d/d0)/2.0;
247 if ((d0-d1)<1e-15)
248 return d1*giac_gnuwince_exp2(k/2.0);
249 }
250 }
251
giac_gnuwince_asin(double d)252 double giac_gnuwince_asin(double d){
253 if (d==1 || d==-1)
254 return d*M_PI_2;
255 double d2=d*d;
256 if (d2>1)
257 return m_undef;
258 return giac_gnuwince_atan(d/giac_gnuwince_sqrt(1-d2));
259 }
260
giac_gnuwince_acos(double d)261 double giac_gnuwince_acos(double d){
262 return M_PI_2-giac_gnuwince_asin(d);
263 }
264
giac_gnuwince_log(double d)265 double giac_gnuwince_log(double d){
266 if (d<0)
267 return m_undef;
268 if (d==0)
269 return m_minusinf;
270 if (d<1)
271 return -giac_gnuwince_log(1.0/d);
272 double k=0.0,k2=512.0;
273 // find number of powers of 2 in d
274 for (int i=0;i<12;++i,k2 /= 2.0){
275 if (d>=pow2tab[i]){
276 k += k2;
277 d /= pow2tab[i];
278 }
279 }
280 // now d>1 and d<2^(1/8) -> d-1>0 and ||<0.091
281 // Taylor expansion at order 14
282 d -= 1;
283 double res=(((((((((((((-25740.0*d+27720.0)*d-30030.0)*d+32760.0)*d-36036.0)*d+40040.0)*d-45045.0)*d+51480.0)*d-60060.0)*d+72072.0)*d-90090.0)*d+120120.0)*d-180180.0)*d+360360.0)*d;
284 return res/360360.0+k*m_ln2;
285 }
286
giac_gnuwince_log10(double d)287 double giac_gnuwince_log10(double d){
288 return giac_gnuwince_log(d)/M_LN10;
289 }
290
giac_gnuwince_acosh(double d)291 double giac_gnuwince_acosh(double d){
292 double d2=d*d;
293 if (d2<1)
294 return m_undef;
295 return giac_gnuwince_log(d+giac_gnuwince_sqrt(d2-1));
296 }
297
giac_gnuwince_asinh(double d)298 double giac_gnuwince_asinh(double d){
299 double d2=d*d;
300 return giac_gnuwince_log(d+giac_gnuwince_sqrt(d2+1));
301 }
302
giac_gnuwince_atanh(double d)303 double giac_gnuwince_atanh(double d){
304 if (d==1)
305 return m_plusinf;
306 if (d==-1)
307 return m_minusinf;
308 double d2=d*d;
309 if (d2>1)
310 return m_undef;
311 return giac_gnuwince_log(giac_gnuwince_sqrt(1-d2)/(1-d));
312 }
313
giac_gnuwince_pow(double x,double y)314 double giac_gnuwince_pow(double x,double y){
315 return giac_gnuwince_exp(y*giac_gnuwince_log(x));
316 }
317
giac_gnuwince_hypot(double x,double y)318 double giac_gnuwince_hypot(double x,double y){
319 return giac_gnuwince_sqrt(x*x+y*y);
320 }
321
giac_gnuwince_atan2(double x,double y)322 double giac_gnuwince_atan2(double x,double y){
323 if (x=0){
324 if (y>0)
325 return M_PI_2;
326 if (y<0)
327 return -M_PI_2;
328 return m_undef;
329 }
330 double res=giac_gnuwince_atan(y/x);
331 if (x>0)
332 return res;
333 if (y>0)
334 return M_PI+res;
335 else
336 return res-M_PI;
337 }
338
giac_gnuwince_exp(const complex<double> & c)339 complex<double> giac_gnuwince_exp(const complex<double> & c){
340 double t=giac_gnuwince_tan(c.imag()/2),t2=t*t;
341 return giac_gnuwince_exp(c.real())/(1.0+t2)*complex<double>(1.0-t2,(2.0*t));
342 }
343
giac_gnuwince_log(const complex<double> & c)344 complex<double> giac_gnuwince_log(const complex<double> & c){
345 double r=c.real(),i=c.imag();
346 return complex<double>(giac_gnuwince_log(r*r+i*i)/2.0,giac_gnuwince_atan2(r,i));
347 }
348
giac_gnuwince_sqrt(const complex<double> & c)349 complex<double> giac_gnuwince_sqrt(const complex<double> & c){
350 double x=c.real(),y=c.imag();
351 if (y==0)
352 return complex<double>(giac_gnuwince_sqrt(x),0);
353 double delta=giac_gnuwince_hypot(x,y);
354 double r=giac_gnuwince_sqrt((delta+x)/2.0);
355 double i=giac_gnuwince_sqrt((delta-x)/2.0);
356 return complex<double>(r,i);
357 }
358
giac_gnuwince_sin(const complex<double> & c)359 complex<double> giac_gnuwince_sin(const complex<double> & c){
360 complex<double> z=giac_gnuwince_exp(complex<double>(-c.imag(),c.real()));
361 return (z-1.0/z)/2.0;
362 }
363
giac_gnuwince_cos(const complex<double> & c)364 complex<double> giac_gnuwince_cos(const complex<double> & c){
365 complex<double> z=giac_gnuwince_exp(complex<double>(-c.imag(),c.real()));
366 return (z+1.0/z)/2.0;
367 }
368
giac_gnuwince_tan(const complex<double> & c)369 complex<double> giac_gnuwince_tan(const complex<double> & c){
370 complex<double> z=giac_gnuwince_exp(complex<double>(-c.imag(),c.real()));
371 z=z*z;
372 return (z-1.0)/(z+1.0);
373 }
374
giac_gnuwince_sinh(const complex<double> & c)375 complex<double> giac_gnuwince_sinh(const complex<double> & c){
376 complex<double> z=giac_gnuwince_exp(c);
377 return (z-1.0/z)/2.0;
378 }
379
giac_gnuwince_cosh(const complex<double> & c)380 complex<double> giac_gnuwince_cosh(const complex<double> & c){
381 complex<double> z=giac_gnuwince_exp(c);
382 return (z+1.0/z)/2.0;
383 }
384
giac_gnuwince_tanh(const complex<double> & c)385 complex<double> giac_gnuwince_tanh(const complex<double> & c){
386 complex<double> z=giac_gnuwince_exp(c);
387 z=z*z;
388 return (z-1.0)/(z+1.0);
389 }
390
391
392 } // end namespace std
393
394 #ifdef LINK
395 using namespace std;
main()396 int main(){
397 ofstream of("res.txt");
398 of << setprecision(15) ;
399 for (double x=-8.3;x<10.0;x=x+1){
400 of << "Exp" << x << char(10) << char(13) ;
401 of << std::exp(x) << " " << giac_gnuwince_exp(x) << char(10) << char(13) ;
402 }
403 for (double x=-8.3;x<10.0;x=x+1){
404 of << "Tan" << x << char(10) << char(13) ;
405 of << std::tan(x) << " " << giac_gnuwince_tan(x) << char(10) << char(13) ;
406 }
407 for (double x=-8.3;x<10.0;x=x+1){
408 of << "Sin" << x << char(10) << char(13) ;
409 of << std::sin(x) << " " << giac_gnuwince_sin(x) << char(10) << char(13) ;
410 }
411 for (double x=-8.3;x<10.0;x=x+1){
412 of << "Cos" << x << char(10) << char(13) ;
413 of << std::cos(x) << " " << giac_gnuwince_cos(x) << char(10) << char(13) ;
414 }
415 for (double x=-8.3;x<10.0;x=x+1){
416 of << "ATan" << x << char(10) << char(13) ;
417 of << std::atan(x) << " " << giac_gnuwince_atan(x) << char(10) << char(13) ;
418 }
419 for (double x=1e-10;x<1e10;x*=10){
420 of << "Natural log" << x << char(10) << char(13) ;
421 of << std::log(x) << " " << giac_gnuwince_log(x) << char(10) << char(13) ;
422 }
423 for (double x=1e-10;x<1e10;x*=10){
424 of << "Sqrt " << x << char(10) << char(13) ;
425 of << std::sqrt(x) << " " << giac_gnuwince_sqrt(x) << char(10) << char(13) ;
426 }
427 for (double x=-1;x<1;x+=0.12345){
428 of << "Asin " << x << char(10) << char(13) ;
429 of << std::asin(x) << " " << giac_gnuwince_asin(x) << char(10) << char(13) ;
430 }
431 }
432 #endif // LINK
433