asinq.c revision 1.1.1.1
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 Long double 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, see 32 <http://www.gnu.org/licenses/>. */ 33 34/* __ieee754_asin(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#include "quadmath-imp.h" 62 63static const __float128 64 one = 1, 65 huge = 1.0e+4932Q, 66 pio2_hi = 1.5707963267948966192313216916397514420986Q, 67 pio2_lo = 4.3359050650618905123985220130216759843812E-35Q, 68 pio4_hi = 7.8539816339744830961566084581987569936977E-1Q, 69 70 /* coefficient for R(x^2) */ 71 72 /* asin(x) = x + x^3 pS(x^2) / qS(x^2) 73 0 <= x <= 0.5 74 peak relative error 1.9e-35 */ 75 pS0 = -8.358099012470680544198472400254596543711E2Q, 76 pS1 = 3.674973957689619490312782828051860366493E3Q, 77 pS2 = -6.730729094812979665807581609853656623219E3Q, 78 pS3 = 6.643843795209060298375552684423454077633E3Q, 79 pS4 = -3.817341990928606692235481812252049415993E3Q, 80 pS5 = 1.284635388402653715636722822195716476156E3Q, 81 pS6 = -2.410736125231549204856567737329112037867E2Q, 82 pS7 = 2.219191969382402856557594215833622156220E1Q, 83 pS8 = -7.249056260830627156600112195061001036533E-1Q, 84 pS9 = 1.055923570937755300061509030361395604448E-3Q, 85 86 qS0 = -5.014859407482408326519083440151745519205E3Q, 87 qS1 = 2.430653047950480068881028451580393430537E4Q, 88 qS2 = -4.997904737193653607449250593976069726962E4Q, 89 qS3 = 5.675712336110456923807959930107347511086E4Q, 90 qS4 = -3.881523118339661268482937768522572588022E4Q, 91 qS5 = 1.634202194895541569749717032234510811216E4Q, 92 qS6 = -4.151452662440709301601820849901296953752E3Q, 93 qS7 = 5.956050864057192019085175976175695342168E2Q, 94 qS8 = -4.175375777334867025769346564600396877176E1Q, 95 /* 1.000000000000000000000000000000000000000E0 */ 96 97 /* asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x) 98 -0.0625 <= x <= 0.0625 99 peak relative error 3.3e-35 */ 100 rS0 = -5.619049346208901520945464704848780243887E0Q, 101 rS1 = 4.460504162777731472539175700169871920352E1Q, 102 rS2 = -1.317669505315409261479577040530751477488E2Q, 103 rS3 = 1.626532582423661989632442410808596009227E2Q, 104 rS4 = -3.144806644195158614904369445440583873264E1Q, 105 rS5 = -9.806674443470740708765165604769099559553E1Q, 106 rS6 = 5.708468492052010816555762842394927806920E1Q, 107 rS7 = 1.396540499232262112248553357962639431922E1Q, 108 rS8 = -1.126243289311910363001762058295832610344E1Q, 109 rS9 = -4.956179821329901954211277873774472383512E-1Q, 110 rS10 = 3.313227657082367169241333738391762525780E-1Q, 111 112 sS0 = -4.645814742084009935700221277307007679325E0Q, 113 sS1 = 3.879074822457694323970438316317961918430E1Q, 114 sS2 = -1.221986588013474694623973554726201001066E2Q, 115 sS3 = 1.658821150347718105012079876756201905822E2Q, 116 sS4 = -4.804379630977558197953176474426239748977E1Q, 117 sS5 = -1.004296417397316948114344573811562952793E2Q, 118 sS6 = 7.530281592861320234941101403870010111138E1Q, 119 sS7 = 1.270735595411673647119592092304357226607E1Q, 120 sS8 = -1.815144839646376500705105967064792930282E1Q, 121 sS9 = -7.821597334910963922204235247786840828217E-2Q, 122 /* 1.000000000000000000000000000000000000000E0 */ 123 124 asinr5625 = 5.9740641664535021430381036628424864397707E-1Q; 125 126 127 128__float128 129asinq (__float128 x) 130{ 131 __float128 t, w, p, q, c, r, s; 132 int32_t ix, sign, flag; 133 ieee854_float128 u; 134 135 flag = 0; 136 u.value = x; 137 sign = u.words32.w0; 138 ix = sign & 0x7fffffff; 139 u.words32.w0 = ix; /* |x| */ 140 if (ix >= 0x3fff0000) /* |x|>= 1 */ 141 { 142 if (ix == 0x3fff0000 143 && (u.words32.w1 | u.words32.w2 | u.words32.w3) == 0) 144 /* asin(1)=+-pi/2 with inexact */ 145 return x * pio2_hi + x * pio2_lo; 146 return (x - x) / (x - x); /* asin(|x|>1) is NaN */ 147 } 148 else if (ix < 0x3ffe0000) /* |x| < 0.5 */ 149 { 150 if (ix < 0x3fc60000) /* |x| < 2**-57 */ 151 { 152 math_check_force_underflow (x); 153 __float128 force_inexact = huge + x; 154 math_force_eval (force_inexact); 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