12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576 |
- /*
- * (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
- */
- /* $Id$ */
- #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;
- }
|