tan.c 1.5 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576
  1. /*
  2. * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
  3. * See the copyright notice in the ACK home directory, in the file "Copyright".
  4. *
  5. * Author: Ceriel J.H. Jacobs
  6. */
  7. /* $Id$ */
  8. #include <math.h>
  9. #include <float.h>
  10. #include <errno.h>
  11. #include "localmath.h"
  12. double
  13. tan(double x)
  14. {
  15. /* Algorithm and coefficients from:
  16. "Software manual for the elementary functions"
  17. by W.J. Cody and W. Waite, Prentice-Hall, 1980
  18. */
  19. int negative = x < 0;
  20. int invert = 0;
  21. double y;
  22. static double p[] = {
  23. 1.0,
  24. -0.13338350006421960681e+0,
  25. 0.34248878235890589960e-2,
  26. -0.17861707342254426711e-4
  27. };
  28. static double q[] = {
  29. 1.0,
  30. -0.46671683339755294240e+0,
  31. 0.25663832289440112864e-1,
  32. -0.31181531907010027307e-3,
  33. 0.49819433993786512270e-6
  34. };
  35. if (__IsNan(x)) {
  36. errno = EDOM;
  37. return x;
  38. }
  39. if (negative) x = -x;
  40. /* ??? avoid loss of significance, error if x is too large ??? */
  41. y = x * M_2_PI + 0.5;
  42. if (y >= DBL_MAX/M_PI_2) return 0.0;
  43. /* Use extended precision to calculate reduced argument.
  44. Here we used 12 bits of the mantissa for a1.
  45. Also split x in integer part x1 and fraction part x2.
  46. */
  47. #define A1 1.57080078125
  48. #define A2 -4.454455103380768678308e-6
  49. {
  50. double x1, x2;
  51. modf(y, &y);
  52. if (modf(0.5*y, &x1)) invert = 1;
  53. x2 = modf(x, &x1);
  54. x = x1 - y * A1;
  55. x += x2;
  56. x -= y * A2;
  57. #undef A1
  58. #undef A2
  59. }
  60. /* ??? avoid underflow ??? */
  61. y = x * x;
  62. x += x * y * POLYNOM2(y, p+1);
  63. y = POLYNOM4(y, q);
  64. if (negative) x = -x;
  65. return invert ? -y/x : x/y;
  66. }