12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273 |
- /*
- * (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 <errno.h>
- extern int errno;
- extern double modf();
- double
- tan(x)
- 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 (negative) x = -x;
- /* ??? avoid loss of significance, error if x is too large ??? */
- y = x * M_2_PI + 0.5;
- /* 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;
- }
|