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