sindg.c

Go to the documentation of this file.

/*                                                     sindg.c
 *
 *     Circular sine of angle in degrees
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, sindg();
 *
 * y = sindg( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Range reduction is into intervals of 45 degrees.
 *
 * Two polynomial approximating functions are employed.
 * Between 0 and pi/4 the sine is approximated by
 *      x  +  x**3 P(x**2).
 * Between pi/4 and pi/2 the cosine is represented as
 *      1  -  x**2 P(x**2).
 *
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain      # trials      peak         rms
 *    IEEE      +-1000       30000      2.3e-16      5.6e-17
 *
 * ERROR MESSAGES:
 *
 *   message           condition        value returned
 * sindg total loss   x > 1.0e14 (IEEE)     0.0
 *
 */
/*                                                     cosdg.c
 *
 *      Circular cosine of angle in degrees
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, cosdg();
 *
 * y = cosdg( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Range reduction is into intervals of 45 degrees.
 *
 * Two polynomial approximating functions are employed.
 * Between 0 and pi/4 the cosine is approximated by
 *      1  -  x**2 P(x**2).
 * Between pi/4 and pi/2 the sine is represented as
 *      x  +  x**3 P(x**2).
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain      # trials      peak         rms
 *    IEEE     +-1000        30000       2.1e-16     5.7e-17
 *  See also sin().
 *
 */

/* Cephes Math Library Release 2.0:  April, 1987
 * Copyright 1985, 1987 by Stephen L. Moshier
 * Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */
 
#include "mconf.h"
 
static double sincof[] = {
    1.58962301572218447952E-10,
    -2.50507477628503540135E-8,
    2.75573136213856773549E-6,
    -1.98412698295895384658E-4,
    8.33333333332211858862E-3,
    -1.66666666666666307295E-1
};
 
static double coscof[] = {
    1.13678171382044553091E-11,
    -2.08758833757683644217E-9,
    2.75573155429816611547E-7,
    -2.48015872936186303776E-5,
    1.38888888888806666760E-3,
    -4.16666666666666348141E-2,
    4.99999999999999999798E-1
};
 
static double PI180 = 1.74532925199432957692E-2;        /* pi/180 */
static double lossth = 1.0e14;
 
double sindg(double x)
{
    double y, z, zz;
    int j, sign;
 
    /* make argument positive but save the sign */
    sign = 1;
    if (x < 0) {
        x = -x;
        sign = -1;
    }
 
    if (x > lossth) {
        sf_error("sindg", SF_ERROR_NO_RESULT, NULL);
        return (0.0);
    }
 
    y = floor(x / 45.0);        /* integer part of x/M_PI_4 */
 
    /* strip high bits of integer part to prevent integer overflow */
    z = ldexp(y, -4);
    z = floor(z);               /* integer part of y/8 */
    z = y - ldexp(z, 4);        /* y - 16 * (y/16) */
 
    j = z;                      /* convert to integer for tests on the phase angle */
    /* map zeros to origin */
    if (j & 1) {
        j += 1;
        y += 1.0;
    }
    j = j & 07;                 /* octant modulo 360 degrees */
    /* reflect in x axis */
    if (j > 3) {
        sign = -sign;
        j -= 4;
    }
 
    z = x - y * 45.0;           /* x mod 45 degrees */
    z *= PI180;                 /* multiply by pi/180 to convert to radians */
    zz = z * z;
 
    if ((j == 1) || (j == 2)) {
        y = 1.0 - zz * polevl(zz, coscof, 6);
    }
    else {
        y = z + z * (zz * polevl(zz, sincof, 5));
    }
 
    if (sign < 0)
        y = -y;
 
    return (y);
}
 
 
double cosdg(double x)
{
    double y, z, zz;
    int j, sign;
 
    /* make argument positive */
    sign = 1;
    if (x < 0)
        x = -x;
 
    if (x > lossth) {
        sf_error("cosdg", SF_ERROR_NO_RESULT, NULL);
        return (0.0);
    }
 
    y = floor(x / 45.0);
    z = ldexp(y, -4);
    z = floor(z);               /* integer part of y/8 */
    z = y - ldexp(z, 4);        /* y - 16 * (y/16) */
 
    /* integer and fractional part modulo one octant */
    j = z;
    if (j & 1) {                /* map zeros to origin */
        j += 1;
        y += 1.0;
    }
    j = j & 07;
    if (j > 3) {
        j -= 4;
        sign = -sign;
    }
 
    if (j > 1)
        sign = -sign;
 
    z = x - y * 45.0;           /* x mod 45 degrees */
    z *= PI180;                 /* multiply by pi/180 to convert to radians */
 
    zz = z * z;
 
    if ((j == 1) || (j == 2)) {
        y = z + z * (zz * polevl(zz, sincof, 5));
    }
    else {
        y = 1.0 - zz * polevl(zz, coscof, 6);
    }
 
    if (sign < 0)
        y = -y;
 
    return (y);
}
 
 
/* Degrees, minutes, seconds to radians: */
 
/* 1 arc second, in radians = 4.848136811095359935899141023579479759563533023727e-6 */
static double P64800 =
    4.848136811095359935899141023579479759563533023727e-6;
 
double radian(double d, double m, double s)
{
    return (((d * 60.0 + m) * 60.0 + s) * P64800);
}