asinf.c 2.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186
  1. /* asinf.c
  2. *
  3. * Inverse circular sine
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * float x, y, asinf();
  10. *
  11. * y = asinf( x );
  12. *
  13. *
  14. *
  15. * DESCRIPTION:
  16. *
  17. * Returns radian angle between -pi/2 and +pi/2 whose sine is x.
  18. *
  19. * A polynomial of the form x + x**3 P(x**2)
  20. * is used for |x| in the interval [0, 0.5]. If |x| > 0.5 it is
  21. * transformed by the identity
  22. *
  23. * asin(x) = pi/2 - 2 asin( sqrt( (1-x)/2 ) ).
  24. *
  25. *
  26. * ACCURACY:
  27. *
  28. * Relative error:
  29. * arithmetic domain # trials peak rms
  30. * IEEE -1, 1 100000 2.5e-7 5.0e-8
  31. *
  32. *
  33. * ERROR MESSAGES:
  34. *
  35. * message condition value returned
  36. * asinf domain |x| > 1 0.0
  37. *
  38. */
  39. /* acosf()
  40. *
  41. * Inverse circular cosine
  42. *
  43. *
  44. *
  45. * SYNOPSIS:
  46. *
  47. * float x, y, acosf();
  48. *
  49. * y = acosf( x );
  50. *
  51. *
  52. *
  53. * DESCRIPTION:
  54. *
  55. * Returns radian angle between -pi/2 and +pi/2 whose cosine
  56. * is x.
  57. *
  58. * Analytically, acos(x) = pi/2 - asin(x). However if |x| is
  59. * near 1, there is cancellation error in subtracting asin(x)
  60. * from pi/2. Hence if x < -0.5,
  61. *
  62. * acos(x) = pi - 2.0 * asin( sqrt((1+x)/2) );
  63. *
  64. * or if x > +0.5,
  65. *
  66. * acos(x) = 2.0 * asin( sqrt((1-x)/2) ).
  67. *
  68. *
  69. * ACCURACY:
  70. *
  71. * Relative error:
  72. * arithmetic domain # trials peak rms
  73. * IEEE -1, 1 100000 1.4e-7 4.2e-8
  74. *
  75. *
  76. * ERROR MESSAGES:
  77. *
  78. * message condition value returned
  79. * acosf domain |x| > 1 0.0
  80. */
  81. /* asin.c */
  82. /*
  83. Cephes Math Library Release 2.2: June, 1992
  84. Copyright 1984, 1987, 1992 by Stephen L. Moshier
  85. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  86. */
  87. /* Single precision circular arcsine
  88. * test interval: [-0.5, +0.5]
  89. * trials: 10000
  90. * peak relative error: 6.7e-8
  91. * rms relative error: 2.5e-8
  92. */
  93. #include <math.h>
  94. extern float PIF, PIO2F;
  95. float sqrtf( float );
  96. float asinf( float xx )
  97. {
  98. float a, x, z;
  99. int sign, flag;
  100. x = xx;
  101. if( x > 0 )
  102. {
  103. sign = 1;
  104. a = x;
  105. }
  106. else
  107. {
  108. sign = -1;
  109. a = -x;
  110. }
  111. if( a > 1.0 )
  112. {
  113. mtherr( "asinf", DOMAIN );
  114. return( 0.0 );
  115. }
  116. if( a < 1.0e-4 )
  117. {
  118. z = a;
  119. goto done;
  120. }
  121. if( a > 0.5 )
  122. {
  123. z = 0.5 * (1.0 - a);
  124. x = sqrtf( z );
  125. flag = 1;
  126. }
  127. else
  128. {
  129. x = a;
  130. z = x * x;
  131. flag = 0;
  132. }
  133. z =
  134. (((( 4.2163199048E-2 * z
  135. + 2.4181311049E-2) * z
  136. + 4.5470025998E-2) * z
  137. + 7.4953002686E-2) * z
  138. + 1.6666752422E-1) * z * x
  139. + x;
  140. if( flag != 0 )
  141. {
  142. z = z + z;
  143. z = PIO2F - z;
  144. }
  145. done:
  146. if( sign < 0 )
  147. z = -z;
  148. return( z );
  149. }
  150. float acosf( float x )
  151. {
  152. if( x < -1.0 )
  153. goto domerr;
  154. if( x < -0.5)
  155. return( PIF - 2.0 * asinf( sqrtf(0.5*(1.0+x)) ) );
  156. if( x > 1.0 )
  157. {
  158. domerr: mtherr( "acosf", DOMAIN );
  159. return( 0.0 );
  160. }
  161. if( x > 0.5 )
  162. return( 2.0 * asinf( sqrtf(0.5*(1.0-x) ) ) );
  163. return( PIO2F - asinf(x) );
  164. }