asinq.c 7.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255
  1. /*
  2. * ====================================================
  3. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  4. *
  5. * Developed at SunPro, a Sun Microsystems, Inc. business.
  6. * Permission to use, copy, modify, and distribute this
  7. * software is freely granted, provided that this notice
  8. * is preserved.
  9. * ====================================================
  10. */
  11. /*
  12. __float128 expansions are
  13. Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
  14. and are incorporated herein by permission of the author. The author
  15. reserves the right to distribute this material elsewhere under different
  16. copying permissions. These modifications are distributed here under the
  17. following terms:
  18. This library is free software; you can redistribute it and/or
  19. modify it under the terms of the GNU Lesser General Public
  20. License as published by the Free Software Foundation; either
  21. version 2.1 of the License, or (at your option) any later version.
  22. This library is distributed in the hope that it will be useful,
  23. but WITHOUT ANY WARRANTY; without even the implied warranty of
  24. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  25. Lesser General Public License for more details.
  26. You should have received a copy of the GNU Lesser General Public
  27. License along with this library; if not, write to the Free Software
  28. Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
  29. /* asinq(x)
  30. * Method :
  31. * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
  32. * we approximate asin(x) on [0,0.5] by
  33. * asin(x) = x + x*x^2*R(x^2)
  34. * Between .5 and .625 the approximation is
  35. * asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
  36. * For x in [0.625,1]
  37. * asin(x) = pi/2-2*asin(sqrt((1-x)/2))
  38. * Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
  39. * then for x>0.98
  40. * asin(x) = pi/2 - 2*(s+s*z*R(z))
  41. * = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
  42. * For x<=0.98, let pio4_hi = pio2_hi/2, then
  43. * f = hi part of s;
  44. * c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z)
  45. * and
  46. * asin(x) = pi/2 - 2*(s+s*z*R(z))
  47. * = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
  48. * = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
  49. *
  50. * Special cases:
  51. * if x is NaN, return x itself;
  52. * if |x|>1, return NaN with invalid signal.
  53. *
  54. */
  55. #include "quadmath-imp.h"
  56. static const __float128
  57. one = 1.0Q,
  58. huge = 1.0e+4932Q,
  59. pio2_hi = 1.5707963267948966192313216916397514420986Q,
  60. pio2_lo = 4.3359050650618905123985220130216759843812E-35Q,
  61. pio4_hi = 7.8539816339744830961566084581987569936977E-1Q,
  62. /* coefficient for R(x^2) */
  63. /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
  64. 0 <= x <= 0.5
  65. peak relative error 1.9e-35 */
  66. pS0 = -8.358099012470680544198472400254596543711E2Q,
  67. pS1 = 3.674973957689619490312782828051860366493E3Q,
  68. pS2 = -6.730729094812979665807581609853656623219E3Q,
  69. pS3 = 6.643843795209060298375552684423454077633E3Q,
  70. pS4 = -3.817341990928606692235481812252049415993E3Q,
  71. pS5 = 1.284635388402653715636722822195716476156E3Q,
  72. pS6 = -2.410736125231549204856567737329112037867E2Q,
  73. pS7 = 2.219191969382402856557594215833622156220E1Q,
  74. pS8 = -7.249056260830627156600112195061001036533E-1Q,
  75. pS9 = 1.055923570937755300061509030361395604448E-3Q,
  76. qS0 = -5.014859407482408326519083440151745519205E3Q,
  77. qS1 = 2.430653047950480068881028451580393430537E4Q,
  78. qS2 = -4.997904737193653607449250593976069726962E4Q,
  79. qS3 = 5.675712336110456923807959930107347511086E4Q,
  80. qS4 = -3.881523118339661268482937768522572588022E4Q,
  81. qS5 = 1.634202194895541569749717032234510811216E4Q,
  82. qS6 = -4.151452662440709301601820849901296953752E3Q,
  83. qS7 = 5.956050864057192019085175976175695342168E2Q,
  84. qS8 = -4.175375777334867025769346564600396877176E1Q,
  85. /* 1.000000000000000000000000000000000000000E0 */
  86. /* asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
  87. -0.0625 <= x <= 0.0625
  88. peak relative error 3.3e-35 */
  89. rS0 = -5.619049346208901520945464704848780243887E0Q,
  90. rS1 = 4.460504162777731472539175700169871920352E1Q,
  91. rS2 = -1.317669505315409261479577040530751477488E2Q,
  92. rS3 = 1.626532582423661989632442410808596009227E2Q,
  93. rS4 = -3.144806644195158614904369445440583873264E1Q,
  94. rS5 = -9.806674443470740708765165604769099559553E1Q,
  95. rS6 = 5.708468492052010816555762842394927806920E1Q,
  96. rS7 = 1.396540499232262112248553357962639431922E1Q,
  97. rS8 = -1.126243289311910363001762058295832610344E1Q,
  98. rS9 = -4.956179821329901954211277873774472383512E-1Q,
  99. rS10 = 3.313227657082367169241333738391762525780E-1Q,
  100. sS0 = -4.645814742084009935700221277307007679325E0Q,
  101. sS1 = 3.879074822457694323970438316317961918430E1Q,
  102. sS2 = -1.221986588013474694623973554726201001066E2Q,
  103. sS3 = 1.658821150347718105012079876756201905822E2Q,
  104. sS4 = -4.804379630977558197953176474426239748977E1Q,
  105. sS5 = -1.004296417397316948114344573811562952793E2Q,
  106. sS6 = 7.530281592861320234941101403870010111138E1Q,
  107. sS7 = 1.270735595411673647119592092304357226607E1Q,
  108. sS8 = -1.815144839646376500705105967064792930282E1Q,
  109. sS9 = -7.821597334910963922204235247786840828217E-2Q,
  110. /* 1.000000000000000000000000000000000000000E0 */
  111. asinr5625 = 5.9740641664535021430381036628424864397707E-1Q;
  112. __float128
  113. asinq (__float128 x)
  114. {
  115. __float128 t = 0;
  116. __float128 w, p, q, c, r, s;
  117. int32_t ix, sign, flag;
  118. ieee854_float128 u;
  119. flag = 0;
  120. u.value = x;
  121. sign = u.words32.w0;
  122. ix = sign & 0x7fffffff;
  123. u.words32.w0 = ix; /* |x| */
  124. if (ix >= 0x3fff0000) /* |x|>= 1 */
  125. {
  126. if (ix == 0x3fff0000
  127. && (u.words32.w1 | u.words32.w2 | u.words32.w3) == 0)
  128. /* asin(1)=+-pi/2 with inexact */
  129. return x * pio2_hi + x * pio2_lo;
  130. return (x - x) / (x - x); /* asin(|x|>1) is NaN */
  131. }
  132. else if (ix < 0x3ffe0000) /* |x| < 0.5 */
  133. {
  134. if (ix < 0x3fc60000) /* |x| < 2**-57 */
  135. {
  136. if (huge + x > one)
  137. return x; /* return x with inexact if x!=0 */
  138. }
  139. else
  140. {
  141. t = x * x;
  142. /* Mark to use pS, qS later on. */
  143. flag = 1;
  144. }
  145. }
  146. else if (ix < 0x3ffe4000) /* 0.625 */
  147. {
  148. t = u.value - 0.5625;
  149. p = ((((((((((rS10 * t
  150. + rS9) * t
  151. + rS8) * t
  152. + rS7) * t
  153. + rS6) * t
  154. + rS5) * t
  155. + rS4) * t
  156. + rS3) * t
  157. + rS2) * t
  158. + rS1) * t
  159. + rS0) * t;
  160. q = ((((((((( t
  161. + sS9) * t
  162. + sS8) * t
  163. + sS7) * t
  164. + sS6) * t
  165. + sS5) * t
  166. + sS4) * t
  167. + sS3) * t
  168. + sS2) * t
  169. + sS1) * t
  170. + sS0;
  171. t = asinr5625 + p / q;
  172. if ((sign & 0x80000000) == 0)
  173. return t;
  174. else
  175. return -t;
  176. }
  177. else
  178. {
  179. /* 1 > |x| >= 0.625 */
  180. w = one - u.value;
  181. t = w * 0.5;
  182. }
  183. p = (((((((((pS9 * t
  184. + pS8) * t
  185. + pS7) * t
  186. + pS6) * t
  187. + pS5) * t
  188. + pS4) * t
  189. + pS3) * t
  190. + pS2) * t
  191. + pS1) * t
  192. + pS0) * t;
  193. q = (((((((( t
  194. + qS8) * t
  195. + qS7) * t
  196. + qS6) * t
  197. + qS5) * t
  198. + qS4) * t
  199. + qS3) * t
  200. + qS2) * t
  201. + qS1) * t
  202. + qS0;
  203. if (flag) /* 2^-57 < |x| < 0.5 */
  204. {
  205. w = p / q;
  206. return x + x * w;
  207. }
  208. s = sqrtq (t);
  209. if (ix >= 0x3ffef333) /* |x| > 0.975 */
  210. {
  211. w = p / q;
  212. t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
  213. }
  214. else
  215. {
  216. u.value = s;
  217. u.words32.w3 = 0;
  218. u.words32.w2 = 0;
  219. w = u.value;
  220. c = (t - w * w) / (s + w);
  221. r = p / q;
  222. p = 2.0 * s * r - (pio2_lo - 2.0 * c);
  223. q = pio4_hi - 2.0 * w;
  224. t = pio4_hi - (p - q);
  225. }
  226. if ((sign & 0x80000000) == 0)
  227. return t;
  228. else
  229. return -t;
  230. }