Code: Select all
/*
* cbm_sine.c
*
* written 2008-08-21 by Michael Kircher
*
* Analysis of Chebyshev approximation of SIN(x) in CBM BASIC
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.141592653589793
double horner(double x, double *coef, unsigned int n)
{
double result=coef[n];
while(n>0) result=result*x+coef[--n];
return(result);
}
double cbm_sine(double x)
{
double u;
double coef[6]=
{/* Coefficients for interpolation polynomial, as given by cbm2dbl. */
6.283185307,
-41.3417021,
81.60522369,
-76.70417026,
42.00779712,
-14.38139067,
};
/* scale argument, and reduce range to [0,1[ */
u=x/(2*PI);
u=u-floor(u);
/*
* Employ symmetries to bring the argument into the range [-0.25,0.25].
* A similar, but more complicated routine also is used in the BASIC ROM.
*/
if(u>0.75) u=u-1.0;
if(u>0.25) u=0.5-u;
/*
* Finally, do "series" evaluation in u², and
* multiply with u to get odd powers of u.
*/
return(u*horner(u*u,coef,5));
}
/*
* Aside: COS(x)=SIN(x+pi/2), TAN(x)=SIN(x)/COS(x)
*/
double cbm2dbl(unsigned char cbm_exponent, unsigned long cbm_mantissa)
{
double result;
if(0 == cbm_exponent && 0 == cbm_mantissa) return(0.0);
result=ldexp(0.5+(cbm_mantissa & 0x7FFFFFFF)/4294967296.0,cbm_exponent-128);
return((cbm_mantissa & 0x80000000)?-result:result);
}
int main(int argc, char *argv[])
{
int i;
double x;
double cbm_result,ieee_result;
double abs_diff;
/*
* Convert the coefficients used for the Chebyshev approximation
* of the sine function, from CBM float to IEEE double. Only the
* first 10 significant digits are shown, as CBM float doesn't
* provide more precision anyway.
*/
printf("%.10G\n",cbm2dbl(0x83,0x490FDAA2));
printf("%.10G\n",cbm2dbl(0x86,0xA55DE728));
printf("%.10G\n",cbm2dbl(0x87,0x2335DFE1));
printf("%.10G\n",cbm2dbl(0x87,0x99688901));
printf("%.10G\n",cbm2dbl(0x86,0x2807FBF8));
printf("%.10G\n",cbm2dbl(0x84,0xE61A2D1B));
/* first test */
printf("%.9G\n",cbm_sine(0.25*PI)); /* = 0.707106781, O.K. */
/*
* Check accuracy against IEEE sin(x) function,
* in full double precision, for some arguments
* in the range 0 to pi/2:
*/
abs_diff=0.0;
for(i=0; i<=1000; i++)
{
x=PI/2.0*(i/1000.0);
cbm_result=cbm_sine(x);
ieee_result= sin (x);
printf("x=%f,\tSIN(x)=%.15E,\tsin(x)=%.15E\n",x,cbm_result,ieee_result);
if(abs_diff<fabs(cbm_result-ieee_result)) abs_diff=fabs(cbm_result-ieee_result);
}
printf("Maximum absolute difference: %E\n",abs_diff); /* =8.150047E-011,
which is quite good for a 5th degree polynomial. CBM BASIC however must do
the same calculation with CBM floats, which results in additional rounding
errors. Thus the last 1 or 2 bits of the mantissa may not be correct, but
this is still accurate to 9 decimal digits. */
exit(EXIT_SUCCESS);
}
