sindg.c 6.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256
  1. /* sindg.c
  2. *
  3. * Circular sine of angle in degrees
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * double x, y, sindg();
  10. *
  11. * y = sindg( x );
  12. *
  13. *
  14. *
  15. * DESCRIPTION:
  16. *
  17. * Range reduction is into intervals of 45 degrees.
  18. *
  19. * Two polynomial approximating functions are employed.
  20. * Between 0 and pi/4 the sine is approximated by
  21. * x + x**3 P(x**2).
  22. * Between pi/4 and pi/2 the cosine is represented as
  23. * 1 - x**2 P(x**2).
  24. *
  25. *
  26. *
  27. * ACCURACY:
  28. *
  29. * Relative error:
  30. * arithmetic domain # trials peak rms
  31. * DEC +-1000 3100 3.3e-17 9.0e-18
  32. * IEEE +-1000 30000 2.3e-16 5.6e-17
  33. *
  34. * ERROR MESSAGES:
  35. *
  36. * message condition value returned
  37. * sindg total loss x > 8.0e14 (DEC) 0.0
  38. * x > 1.0e14 (IEEE)
  39. *
  40. */
  41. /* cosdg.c
  42. *
  43. * Circular cosine of angle in degrees
  44. *
  45. *
  46. *
  47. * SYNOPSIS:
  48. *
  49. * double x, y, cosdg();
  50. *
  51. * y = cosdg( x );
  52. *
  53. *
  54. *
  55. * DESCRIPTION:
  56. *
  57. * Range reduction is into intervals of 45 degrees.
  58. *
  59. * Two polynomial approximating functions are employed.
  60. * Between 0 and pi/4 the cosine is approximated by
  61. * 1 - x**2 P(x**2).
  62. * Between pi/4 and pi/2 the sine is represented as
  63. * x + x**3 P(x**2).
  64. *
  65. *
  66. * ACCURACY:
  67. *
  68. * Relative error:
  69. * arithmetic domain # trials peak rms
  70. * DEC +-1000 3400 3.5e-17 9.1e-18
  71. * IEEE +-1000 30000 2.1e-16 5.7e-17
  72. * See also sin().
  73. *
  74. */
  75. /* Cephes Math Library Release 2.0: April, 1987
  76. * Copyright 1985, 1987 by Stephen L. Moshier
  77. * Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */
  78. #include "mconf.h"
  79. #ifdef UNK
  80. static double sincof[] = {
  81. 1.58962301572218447952E-10, -2.50507477628503540135E-8,
  82. 2.75573136213856773549E-6, -1.98412698295895384658E-4,
  83. 8.33333333332211858862E-3, -1.66666666666666307295E-1};
  84. static double coscof[] = {
  85. 1.13678171382044553091E-11, -2.08758833757683644217E-9,
  86. 2.75573155429816611547E-7, -2.48015872936186303776E-5,
  87. 1.38888888888806666760E-3, -4.16666666666666348141E-2,
  88. 4.99999999999999999798E-1};
  89. static double PI180 = 1.74532925199432957692E-2; /* pi/180 */
  90. static double lossth = 1.0e14;
  91. #endif
  92. #ifdef DEC
  93. static unsigned short sincof[] = {
  94. 0030056, 0143750, 0177170, 0073013, 0131727, 0027455, 0044510, 0132205,
  95. 0033470, 0167432, 0131752, 0042263, 0135120, 0006400, 0146776, 0174027,
  96. 0036410, 0104210, 0104207, 0137202, 0137452, 0125252, 0125252, 0125103};
  97. static unsigned short coscof[] = {
  98. 0027107, 0176030, 0153315, 0110312, 0131017, 0072476, 0007450,
  99. 0123243, 0032623, 0171174, 0070066, 0146445, 0134320, 0006400,
  100. 0147355, 0163313, 0035666, 0005540, 0133012, 0165067, 0137052,
  101. 0125252, 0125252, 0125206, 0040000, 0000000, 0000000, 0000000};
  102. static unsigned short P1[] = {0036616, 0175065, 0011224, 0164711};
  103. #define PI180 *(double *)P1
  104. static double lossth = 8.0e14;
  105. #endif
  106. #ifdef IBMPC
  107. static unsigned short sincof[] = {
  108. 0x0ec1, 0x1fcf, 0xd8fd, 0x3de5, 0x1691, 0xa929, 0xe5e5, 0xbe5a,
  109. 0x4896, 0x567d, 0x1de3, 0x3ec7, 0xdf03, 0x19bf, 0x01a0, 0xbf2a,
  110. 0xf7d0, 0x1110, 0x1111, 0x3f81, 0x5548, 0x5555, 0x5555, 0xbfc5};
  111. static unsigned short coscof[] = {
  112. 0xb219, 0x1ad9, 0xff83, 0x3da8, 0x14d4, 0xc1e5, 0xeea7,
  113. 0xbe21, 0xd9a5, 0x8e06, 0x7e4f, 0x3e92, 0xbcd9, 0x19dd,
  114. 0x01a0, 0xbefa, 0x5d47, 0x16c1, 0xc16c, 0x3f56, 0x5551,
  115. 0x5555, 0x5555, 0xbfa5, 0x0000, 0x0000, 0x0000, 0x3fe0};
  116. static unsigned short P1[] = {0x9d39, 0xa252, 0xdf46, 0x3f91};
  117. #define PI180 *(double *)P1
  118. static double lossth = 1.0e14;
  119. #endif
  120. #ifdef MIEEE
  121. static unsigned short sincof[] = {
  122. 0x3de5, 0xd8fd, 0x1fcf, 0x0ec1, 0xbe5a, 0xe5e5, 0xa929, 0x1691,
  123. 0x3ec7, 0x1de3, 0x567d, 0x4896, 0xbf2a, 0x01a0, 0x19bf, 0xdf03,
  124. 0x3f81, 0x1111, 0x1110, 0xf7d0, 0xbfc5, 0x5555, 0x5555, 0x5548};
  125. static unsigned short coscof[] = {
  126. 0x3da8, 0xff83, 0x1ad9, 0xb219, 0xbe21, 0xeea7, 0xc1e5,
  127. 0x14d4, 0x3e92, 0x7e4f, 0x8e06, 0xd9a5, 0xbefa, 0x01a0,
  128. 0x19dd, 0xbcd9, 0x3f56, 0xc16c, 0x16c1, 0x5d47, 0xbfa5,
  129. 0x5555, 0x5555, 0x5551, 0x3fe0, 0x0000, 0x0000, 0x0000};
  130. static unsigned short P1[] = {0x3f91, 0xdf46, 0xa252, 0x9d39};
  131. #define PI180 *(double *)P1
  132. static double lossth = 1.0e14;
  133. #endif
  134. #ifdef ANSIPROT
  135. extern double polevl(double, void *, int);
  136. extern double floor(double);
  137. extern double ldexp(double, int);
  138. #else
  139. double polevl(), floor(), ldexp();
  140. #endif
  141. extern double PIO4;
  142. double sindg(x) double x;
  143. {
  144. double y, z, zz;
  145. int j, sign;
  146. /* make argument positive but save the sign */
  147. sign = 1;
  148. if (x < 0) {
  149. x = -x;
  150. sign = -1;
  151. }
  152. if (x > lossth) {
  153. mtherr("sindg", TLOSS);
  154. return (0.0);
  155. }
  156. y = floor(x / 45.0); /* integer part of x/PIO4 */
  157. /* strip high bits of integer part to prevent integer overflow */
  158. z = ldexp(y, -4);
  159. z = floor(z); /* integer part of y/8 */
  160. z = y - ldexp(z, 4); /* y - 16 * (y/16) */
  161. j = z; /* convert to integer for tests on the phase angle */
  162. /* map zeros to origin */
  163. if (j & 1) {
  164. j += 1;
  165. y += 1.0;
  166. }
  167. j = j & 07; /* octant modulo 360 degrees */
  168. /* reflect in x axis */
  169. if (j > 3) {
  170. sign = -sign;
  171. j -= 4;
  172. }
  173. z = x - y * 45.0; /* x mod 45 degrees */
  174. z *= PI180; /* multiply by pi/180 to convert to radians */
  175. zz = z * z;
  176. if ((j == 1) || (j == 2)) {
  177. y = 1.0 - zz * polevl(zz, coscof, 6);
  178. } else {
  179. y = z + z * (zz * polevl(zz, sincof, 5));
  180. }
  181. if (sign < 0)
  182. y = -y;
  183. return (y);
  184. }
  185. double cosdg(x) double x;
  186. {
  187. double y, z, zz;
  188. int j, sign;
  189. /* make argument positive */
  190. sign = 1;
  191. if (x < 0)
  192. x = -x;
  193. if (x > lossth) {
  194. mtherr("cosdg", TLOSS);
  195. return (0.0);
  196. }
  197. y = floor(x / 45.0);
  198. z = ldexp(y, -4);
  199. z = floor(z); /* integer part of y/8 */
  200. z = y - ldexp(z, 4); /* y - 16 * (y/16) */
  201. /* integer and fractional part modulo one octant */
  202. j = z;
  203. if (j & 1) /* map zeros to origin */
  204. {
  205. j += 1;
  206. y += 1.0;
  207. }
  208. j = j & 07;
  209. if (j > 3) {
  210. j -= 4;
  211. sign = -sign;
  212. }
  213. if (j > 1)
  214. sign = -sign;
  215. z = x - y * 45.0; /* x mod 45 degrees */
  216. z *= PI180; /* multiply by pi/180 to convert to radians */
  217. zz = z * z;
  218. if ((j == 1) || (j == 2)) {
  219. y = z + z * (zz * polevl(zz, sincof, 5));
  220. } else {
  221. y = 1.0 - zz * polevl(zz, coscof, 6);
  222. }
  223. if (sign < 0)
  224. y = -y;
  225. return (y);
  226. }