/* FUNCTIONS: double sin(x) double cos(x) PURPOSE: To provide sine and cosine of angles expressed in radians. INPUT PARAMETERS: double x - The given angle in radians (maximum |x| is 3.37E+9) OUTPUT PARAMETERS: Returned - Sine(x) or Cosine(x), except where |x| > 3.37E+9. In this case, 0 is returned. DESCRIPTION: The algorithm uses a polynomial evaluation technique. See algorithm 3354 of Hart's "Computer Approximations". Note: the cosine routine simply calls the sine routine via sin(x+pi/2). WRITTEN BY: Scott D. Roth DATE: March 27, 1981 */ #define ODD(x) (((x) & 1) != 0) static double coeff[9] = { 0.15707963267948966E1d, -0.64596409750624619E0d, 0.79692626246165436E-1d, -0.46817541353026426E-2d, 0.16044118470682207E-3d, -0.35988430072086927E-5d, 0.5692134872719023E-7d, -0.66843217206396E-9d, 0.587061098171E-11d}; double sin (x) /* And "double cos (x)" follows */ double x; { double abs_x, s, c; register double *c_ptr; long int_part; register int neg_x; /* Boolean */ /* For small x, x = sin(x) */ if (1.0e-13d < (abs_x = ((neg_x = (x < 0)) ? -x : x))) if (3.37E9d < abs_x) /* Too big to figure */ x = 0; else { abs_x /= coeff[0]; /* x / (pi/2) */ if (ODD((int) (int_part = abs_x))) int_part++; c = abs_x - int_part; if (neg_x ^ ((((int) int_part) & 3) == 2)) c = -c; s = c * c; x = *(c_ptr = &coeff[8]); do x = s * x + *--c_ptr; while (c_ptr != coeff); x *= c; } return(x); } double cos (x) double x; { return(sin(x + 1.57079632679489661)); /* sin(x + pi/2) */ }