sindgf.c 3.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232
  1. /* sindgf.c
  2. *
  3. * Circular sine of angle in degrees
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * float x, y, sindgf();
  10. *
  11. * y = sindgf( 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 Q(x**2).
  24. *
  25. *
  26. * ACCURACY:
  27. *
  28. * Relative error:
  29. * arithmetic domain # trials peak rms
  30. * IEEE +-3600 100,000 1.2e-7 3.0e-8
  31. *
  32. * ERROR MESSAGES:
  33. *
  34. * message condition value returned
  35. * sin total loss x > 2^24 0.0
  36. *
  37. */
  38. /* cosdgf.c
  39. *
  40. * Circular cosine of angle in degrees
  41. *
  42. *
  43. *
  44. * SYNOPSIS:
  45. *
  46. * float x, y, cosdgf();
  47. *
  48. * y = cosdgf( x );
  49. *
  50. *
  51. *
  52. * DESCRIPTION:
  53. *
  54. * Range reduction is into intervals of 45 degrees.
  55. *
  56. * Two polynomial approximating functions are employed.
  57. * Between 0 and pi/4 the cosine is approximated by
  58. * 1 - x**2 Q(x**2).
  59. * Between pi/4 and pi/2 the sine is represented as
  60. * x + x**3 P(x**2).
  61. *
  62. *
  63. * ACCURACY:
  64. *
  65. * Relative error:
  66. * arithmetic domain # trials peak rms
  67. * IEEE -8192,+8192 100,000 3.0e-7 3.0e-8
  68. */
  69. /*
  70. Cephes Math Library Release 2.2: June, 1992
  71. Copyright 1985, 1987, 1988, 1992 by Stephen L. Moshier
  72. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  73. */
  74. /* Single precision circular sine
  75. * test interval: [-pi/4, +pi/4]
  76. * trials: 10000
  77. * peak relative error: 6.8e-8
  78. * rms relative error: 2.6e-8
  79. */
  80. #include <math.h>
  81. /*static float FOPI = 1.27323954473516;*/
  82. extern float PIO4F;
  83. /* These are for a 24-bit significand: */
  84. static float T24M1 = 16777215.;
  85. static float PI180 = 0.0174532925199432957692; /* pi/180 */
  86. float sindgf( float xx )
  87. {
  88. float x, y, z;
  89. long j;
  90. int sign;
  91. sign = 1;
  92. x = xx;
  93. if( xx < 0 )
  94. {
  95. sign = -1;
  96. x = -xx;
  97. }
  98. if( x > T24M1 )
  99. {
  100. mtherr( "sindgf", TLOSS );
  101. return(0.0);
  102. }
  103. j = 0.022222222222222222222 * x; /* integer part of x/45 */
  104. y = j;
  105. /* map zeros to origin */
  106. if( j & 1 )
  107. {
  108. j += 1;
  109. y += 1.0;
  110. }
  111. j &= 7; /* octant modulo 360 degrees */
  112. /* reflect in x axis */
  113. if( j > 3)
  114. {
  115. sign = -sign;
  116. j -= 4;
  117. }
  118. x = x - y * 45.0;
  119. x *= PI180; /* multiply by pi/180 to convert to radians */
  120. z = x * x;
  121. if( (j==1) || (j==2) )
  122. {
  123. /*
  124. y = ((( 2.4462803166E-5 * z
  125. - 1.3887580023E-3) * z
  126. + 4.1666650433E-2) * z
  127. - 4.9999999968E-1) * z
  128. + 1.0;
  129. */
  130. /* measured relative error in +/- pi/4 is 7.8e-8 */
  131. y = (( 2.443315711809948E-005 * z
  132. - 1.388731625493765E-003) * z
  133. + 4.166664568298827E-002) * z * z;
  134. y -= 0.5 * z;
  135. y += 1.0;
  136. }
  137. else
  138. {
  139. /* Theoretical relative error = 3.8e-9 in [-pi/4, +pi/4] */
  140. y = ((-1.9515295891E-4 * z
  141. + 8.3321608736E-3) * z
  142. - 1.6666654611E-1) * z * x;
  143. y += x;
  144. }
  145. if(sign < 0)
  146. y = -y;
  147. return( y);
  148. }
  149. /* Single precision circular cosine
  150. * test interval: [-pi/4, +pi/4]
  151. * trials: 10000
  152. * peak relative error: 8.3e-8
  153. * rms relative error: 2.2e-8
  154. */
  155. float cosdgf( float xx )
  156. {
  157. register float x, y, z;
  158. int j, sign;
  159. /* make argument positive */
  160. sign = 1;
  161. x = xx;
  162. if( x < 0 )
  163. x = -x;
  164. if( x > T24M1 )
  165. {
  166. mtherr( "cosdgf", TLOSS );
  167. return(0.0);
  168. }
  169. j = 0.02222222222222222222222 * x; /* integer part of x/PIO4 */
  170. y = j;
  171. /* integer and fractional part modulo one octant */
  172. if( j & 1 ) /* map zeros to origin */
  173. {
  174. j += 1;
  175. y += 1.0;
  176. }
  177. j &= 7;
  178. if( j > 3)
  179. {
  180. j -=4;
  181. sign = -sign;
  182. }
  183. if( j > 1 )
  184. sign = -sign;
  185. x = x - y * 45.0; /* x mod 45 degrees */
  186. x *= PI180; /* multiply by pi/180 to convert to radians */
  187. z = x * x;
  188. if( (j==1) || (j==2) )
  189. {
  190. y = (((-1.9515295891E-4 * z
  191. + 8.3321608736E-3) * z
  192. - 1.6666654611E-1) * z * x)
  193. + x;
  194. }
  195. else
  196. {
  197. y = (( 2.443315711809948E-005 * z
  198. - 1.388731625493765E-003) * z
  199. + 4.166664568298827E-002) * z * z;
  200. y -= 0.5 * z;
  201. y += 1.0;
  202. }
  203. if(sign < 0)
  204. y = -y;
  205. return( y );
  206. }