tan.c 1.4 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273
  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 <errno.h>
  10. extern int errno;
  11. extern double modf();
  12. double
  13. tan(x)
  14. double x;
  15. {
  16. /* Algorithm and coefficients from:
  17. "Software manual for the elementary functions"
  18. by W.J. Cody and W. Waite, Prentice-Hall, 1980
  19. */
  20. int negative = x < 0;
  21. int invert = 0;
  22. double y;
  23. static double p[] = {
  24. 1.0,
  25. -0.13338350006421960681e+0,
  26. 0.34248878235890589960e-2,
  27. -0.17861707342254426711e-4
  28. };
  29. static double q[] = {
  30. 1.0,
  31. -0.46671683339755294240e+0,
  32. 0.25663832289440112864e-1,
  33. -0.31181531907010027307e-3,
  34. 0.49819433993786512270e-6
  35. };
  36. if (negative) x = -x;
  37. /* ??? avoid loss of significance, error if x is too large ??? */
  38. y = x * M_2_PI + 0.5;
  39. /* Use extended precision to calculate reduced argument.
  40. Here we used 12 bits of the mantissa for a1.
  41. Also split x in integer part x1 and fraction part x2.
  42. */
  43. #define A1 1.57080078125
  44. #define A2 -4.454455103380768678308e-6
  45. {
  46. double x1, x2;
  47. modf(y, &y);
  48. if (modf(0.5*y, &x1)) invert = 1;
  49. x2 = modf(x, &x1);
  50. x = x1 - y * A1;
  51. x += x2;
  52. x -= y * A2;
  53. #undef A1
  54. #undef A2
  55. }
  56. /* ??? avoid underflow ??? */
  57. y = x * x;
  58. x += x * y * POLYNOM2(y, p+1);
  59. y = POLYNOM4(y, q);
  60. if (negative) x = -x;
  61. return invert ? -y/x : x/y;
  62. }