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