gfs@abvax.UUCP (Greg F. Shay) (05/23/88)
Here's my program. It is basically a straightforward power series computation, albeit in integer. The only tricky part was keeping track of the decimal point and making sure the intermediate calculations did not overflow. In my first version, I left the constant expressions like (short)(SCALE*6.2831853) directly in the code, thinking the compiler would do the math and store the (short) result as a constant. To my surprise, the compiler left the computation in (@#%!* compiler!). So I resorted to sticking the 'magic numbers' in by hand. The difficulty is now that the #defines cannot be changed for any other values. By changing these constants and the #defines, 8 decimal points (#define SSHIFT 8 #define SCALE 256) works, but SSHIFT 9 overflows the intermediate computations. Greg Shay ..pyramid | ..mandrill |..!abvax!gfs ..decvax | #define SSHIFT 7 #define SCALE 128 short cosx(x) short x; { /* Return power series approximation to cos(x) */ /* using only simple math functions */ short frac,sign; int y,z; sign = 1; /* cos(-x) = cos(x) */ if (x<0) x = -x; /* cos(x+2pi) = cos(x) */ frac = x/804; /* 804=(SCALE*6.2831853) */ x = x- frac*804; /* cos(x) = cos(2pi-x) */ if (x > 402) /*402 = SCALE * 3.14159265 */ x = 804 - x; /* cos(x) = -cos(pi-x) */ if (x > 201 ) { x = 402-x; sign = -1; } z = (((int)x)*x)>>SSHIFT; /* z has x^2 */ frac = SCALE-(z>>1); y = z*z; /* y has SCALE*x^4 */ frac += y*682 >>((SSHIFT<<1)+SSHIFT); /* 682 = ((int)(.04166667*SCALE*SCALE)) */ z = y*z>>SSHIFT; /* z has SCALE*x^6 */ frac -= z*23 >>((SSHIFT<<1)+SSHIFT); /* 23 = (int)(.001388889*SCALE*SCALE) */ /* The following are the first 5 terms of the power series. Only the first */ /* three are needed at the seven decimal bit precision representation. *. /*frac = 1-x*x*.5+x*x*x*x*.041666667-x*x*x*x*x*x*.00138889+x*x*x*x*x*x*x*x* (2.4801587E-5)-x*x*x*x*x*x*x*x*x*x*(2.7557319E-7);*/ return((sign==1) ? frac : -frac); } short sinx(x) short x; { /* Return power series approximation to sin(x) */ /* Use the fact that sin(x) = cos(x-pi/2) */ return( cosx(x-((short)(SCALE*1.5707963)) )); }