sin.c 1.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899
  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. static double
  13. sinus(x, cos_flag)
  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. static double r[] = {
  21. -0.16666666666666665052e+0,
  22. 0.83333333333331650314e-2,
  23. -0.19841269841201840457e-3,
  24. 0.27557319210152756119e-5,
  25. -0.25052106798274584544e-7,
  26. 0.16058936490371589114e-9,
  27. -0.76429178068910467734e-12,
  28. 0.27204790957888846175e-14
  29. };
  30. double xsqr;
  31. double y;
  32. int neg = 0;
  33. if (x < 0) {
  34. x = -x;
  35. neg = 1;
  36. }
  37. if (cos_flag) {
  38. neg = 0;
  39. y = M_PI_2 + x;
  40. }
  41. else y = x;
  42. /* ??? avoid loss of significance, if y is too large, error ??? */
  43. y = y * M_1_PI + 0.5;
  44. /* Use extended precision to calculate reduced argument.
  45. Here we used 12 bits of the mantissa for a1.
  46. Also split x in integer part x1 and fraction part x2.
  47. */
  48. #define A1 3.1416015625
  49. #define A2 -8.908910206761537356617e-6
  50. {
  51. double x1, x2;
  52. modf(y, &y);
  53. if (modf(0.5*y, &x1)) neg = !neg;
  54. if (cos_flag) y -= 0.5;
  55. x2 = modf(x, &x1);
  56. x = x1 - y * A1;
  57. x += x2;
  58. x -= y * A2;
  59. #undef A1
  60. #undef A2
  61. }
  62. if (x < 0) {
  63. neg = !neg;
  64. x = -x;
  65. }
  66. /* ??? avoid underflow ??? */
  67. y = x * x;
  68. x += x * y * POLYNOM7(y, r);
  69. return neg ? -x : x;
  70. }
  71. double
  72. sin(x)
  73. double x;
  74. {
  75. return sinus(x, 0);
  76. }
  77. double
  78. cos(x)
  79. double x;
  80. {
  81. if (x < 0) x = -x;
  82. return sinus(x, 1);
  83. }