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)) ));
}