source: trunk/minix/lib/math/tan.c@ 11

/*
 * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
 * See the copyright notice in the ACK home directory, in the file "Copyright".
 *
 * Author: Ceriel J.H. Jacobs
 */
/* $Header: /cvsup/minix/src/lib/math/tan.c,v 1.1.1.1 2005/04/21 14:56:26 beng Exp $ */

#include        <math.h>
#include        <float.h>
#include        <errno.h>
#include        "localmath.h"

double
tan(double x)
{
        /*      Algorithm and coefficients from:
                        "Software manual for the elementary functions"
                        by W.J. Cody and W. Waite, Prentice-Hall, 1980
        */

        int negative = x < 0;
        int invert = 0;
        double  y;
        static double   p[] = {
                 1.0,
                -0.13338350006421960681e+0,
                 0.34248878235890589960e-2,
                -0.17861707342254426711e-4
        };
        static double   q[] = {
                 1.0,
                -0.46671683339755294240e+0,
                 0.25663832289440112864e-1,
                -0.31181531907010027307e-3,
                 0.49819433993786512270e-6
        };

        if (__IsNan(x)) {
                errno = EDOM;
                return x;
        }
        if (negative) x = -x;
 
        /* ??? avoid loss of significance, error if x is too large ??? */

        y = x * M_2_PI + 0.5;

        if (y >= DBL_MAX/M_PI_2) return 0.0;

        /*      Use extended precision to calculate reduced argument.
                Here we used 12 bits of the mantissa for a1.
                Also split x in integer part x1 and fraction part x2.
        */
    #define A1 1.57080078125
    #define A2 -4.454455103380768678308e-6
        {
                double x1, x2;

                modf(y, &y);
                if (modf(0.5*y, &x1)) invert = 1;
                x2 = modf(x, &x1);
                x = x1 - y * A1;
                x += x2;
                x -= y * A2;
    #undef A1
    #undef A2
        }

        /* ??? avoid underflow ??? */
        y = x * x;
        x += x * y * POLYNOM2(y, p+1);
        y = POLYNOM4(y, q);
        if (negative) x = -x;
        return invert ? -y/x : x/y;
}