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 18 the 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 /* acosq(x) 35 * Method : 36 * acos(x) = pi/2 - asin(x) acoshq(__float128 x)37 * acos(-x) = pi/2 + asin(x) 38 * For |x| <= 0.375 39 * acos(x) = pi/2 - asin(x) 40 * Between .375 and .5 the approximation is 41 * acos(0.4375 + x) = acos(0.4375) + x P(x) / Q(x) 42 * Between .5 and .625 the approximation is 43 * acos(0.5625 + x) = acos(0.5625) + x rS(x) / sS(x) 44 * For x > 0.625, 45 * acos(x) = 2 asin(sqrt((1-x)/2)) 46 * computed with an extended precision square root in the leading term. 47 * For x < -0.625 48 * acos(x) = pi - 2 asin(sqrt((1-|x|)/2)) 49 * 50 * Special cases: 51 * if x is NaN, return x itself; 52 * if |x|>1, return NaN with invalid signal. 53 * 54 * Functions needed: sqrtq. 55 */ 56 57 #include "quadmath-imp.h" 58 59 static const __float128 60 one = 1, 61 pio2_hi = 1.5707963267948966192313216916397514420986Q, 62 pio2_lo = 4.3359050650618905123985220130216759843812E-35Q, 63 64 /* acos(0.5625 + x) = acos(0.5625) + x rS(x) / sS(x) 65 -0.0625 <= x <= 0.0625 66 peak relative error 3.3e-35 */ 67 68 rS0 = 5.619049346208901520945464704848780243887E0Q, 69 rS1 = -4.460504162777731472539175700169871920352E1Q, 70 rS2 = 1.317669505315409261479577040530751477488E2Q, 71 rS3 = -1.626532582423661989632442410808596009227E2Q, 72 rS4 = 3.144806644195158614904369445440583873264E1Q, 73 rS5 = 9.806674443470740708765165604769099559553E1Q, 74 rS6 = -5.708468492052010816555762842394927806920E1Q, 75 rS7 = -1.396540499232262112248553357962639431922E1Q, 76 rS8 = 1.126243289311910363001762058295832610344E1Q, 77 rS9 = 4.956179821329901954211277873774472383512E-1Q, 78 rS10 = -3.313227657082367169241333738391762525780E-1Q, 79 80 sS0 = -4.645814742084009935700221277307007679325E0Q, 81 sS1 = 3.879074822457694323970438316317961918430E1Q, 82 sS2 = -1.221986588013474694623973554726201001066E2Q, 83 sS3 = 1.658821150347718105012079876756201905822E2Q, 84 sS4 = -4.804379630977558197953176474426239748977E1Q, 85 sS5 = -1.004296417397316948114344573811562952793E2Q, 86 sS6 = 7.530281592861320234941101403870010111138E1Q, 87 sS7 = 1.270735595411673647119592092304357226607E1Q, 88 sS8 = -1.815144839646376500705105967064792930282E1Q, 89 sS9 = -7.821597334910963922204235247786840828217E-2Q, 90 /* 1.000000000000000000000000000000000000000E0 */ 91 92 acosr5625 = 9.7338991014954640492751132535550279812151E-1Q, 93 pimacosr5625 = 2.1682027434402468335351320579240000860757E0Q, 94 95 /* acos(0.4375 + x) = acos(0.4375) + x rS(x) / sS(x) 96 -0.0625 <= x <= 0.0625 97 peak relative error 2.1e-35 */ 98 99 P0 = 2.177690192235413635229046633751390484892E0Q, 100 P1 = -2.848698225706605746657192566166142909573E1Q, 101 P2 = 1.040076477655245590871244795403659880304E2Q, 102 P3 = -1.400087608918906358323551402881238180553E2Q, 103 P4 = 2.221047917671449176051896400503615543757E1Q, 104 P5 = 9.643714856395587663736110523917499638702E1Q, 105 P6 = -5.158406639829833829027457284942389079196E1Q, 106 P7 = -1.578651828337585944715290382181219741813E1Q, 107 P8 = 1.093632715903802870546857764647931045906E1Q, 108 P9 = 5.448925479898460003048760932274085300103E-1Q, 109 P10 = -3.315886001095605268470690485170092986337E-1Q, 110 Q0 = -1.958219113487162405143608843774587557016E0Q, 111 Q1 = 2.614577866876185080678907676023269360520E1Q, 112 Q2 = -9.990858606464150981009763389881793660938E1Q, 113 Q3 = 1.443958741356995763628660823395334281596E2Q, 114 Q4 = -3.206441012484232867657763518369723873129E1Q, 115 Q5 = -1.048560885341833443564920145642588991492E2Q, 116 Q6 = 6.745883931909770880159915641984874746358E1Q, 117 Q7 = 1.806809656342804436118449982647641392951E1Q, 118 Q8 = -1.770150690652438294290020775359580915464E1Q, 119 Q9 = -5.659156469628629327045433069052560211164E-1Q, 120 /* 1.000000000000000000000000000000000000000E0 */ 121 122 acosr4375 = 1.1179797320499710475919903296900511518755E0Q, 123 pimacosr4375 = 2.0236129215398221908706530535894517323217E0Q, 124 125 /* asin(x) = x + x^3 pS(x^2) / qS(x^2) 126 0 <= x <= 0.5 127 peak relative error 1.9e-35 */ 128 pS0 = -8.358099012470680544198472400254596543711E2Q, 129 pS1 = 3.674973957689619490312782828051860366493E3Q, 130 pS2 = -6.730729094812979665807581609853656623219E3Q, 131 pS3 = 6.643843795209060298375552684423454077633E3Q, 132 pS4 = -3.817341990928606692235481812252049415993E3Q, 133 pS5 = 1.284635388402653715636722822195716476156E3Q, 134 pS6 = -2.410736125231549204856567737329112037867E2Q, 135 pS7 = 2.219191969382402856557594215833622156220E1Q, 136 pS8 = -7.249056260830627156600112195061001036533E-1Q, 137 pS9 = 1.055923570937755300061509030361395604448E-3Q, 138 139 qS0 = -5.014859407482408326519083440151745519205E3Q, 140 qS1 = 2.430653047950480068881028451580393430537E4Q, 141 qS2 = -4.997904737193653607449250593976069726962E4Q, 142 qS3 = 5.675712336110456923807959930107347511086E4Q, 143 qS4 = -3.881523118339661268482937768522572588022E4Q, 144 qS5 = 1.634202194895541569749717032234510811216E4Q, 145 qS6 = -4.151452662440709301601820849901296953752E3Q, 146 qS7 = 5.956050864057192019085175976175695342168E2Q, 147 qS8 = -4.175375777334867025769346564600396877176E1Q; 148 /* 1.000000000000000000000000000000000000000E0 */ 149 150 __float128 151 acosq (__float128 x) 152 { 153 __float128 z, r, w, p, q, s, t, f2; 154 int32_t ix, sign; 155 ieee854_float128 u; 156 157 u.value = x; 158 sign = u.words32.w0; 159 ix = sign & 0x7fffffff; 160 u.words32.w0 = ix; /* |x| */ 161 if (ix >= 0x3fff0000) /* |x| >= 1 */ 162 { 163 if (ix == 0x3fff0000 164 && (u.words32.w1 | u.words32.w2 | u.words32.w3) == 0) 165 { /* |x| == 1 */ 166 if ((sign & 0x80000000) == 0) 167 return 0.0; /* acos(1) = 0 */ 168 else 169 return (2.0 * pio2_hi) + (2.0 * pio2_lo); /* acos(-1)= pi */ 170 } 171 return (x - x) / (x - x); /* acos(|x| > 1) is NaN */ 172 } 173 else if (ix < 0x3ffe0000) /* |x| < 0.5 */ 174 { 175 if (ix < 0x3f8e0000) /* |x| < 2**-113 */ 176 return pio2_hi + pio2_lo; 177 if (ix < 0x3ffde000) /* |x| < .4375 */ 178 { 179 /* Arcsine of x. */ 180 z = x * x; 181 p = (((((((((pS9 * z 182 + pS8) * z 183 + pS7) * z 184 + pS6) * z 185 + pS5) * z 186 + pS4) * z 187 + pS3) * z 188 + pS2) * z 189 + pS1) * z 190 + pS0) * z; 191 q = (((((((( z 192 + qS8) * z 193 + qS7) * z 194 + qS6) * z 195 + qS5) * z 196 + qS4) * z 197 + qS3) * z 198 + qS2) * z 199 + qS1) * z 200 + qS0; 201 r = x + x * p / q; 202 z = pio2_hi - (r - pio2_lo); 203 return z; 204 } 205 /* .4375 <= |x| < .5 */ 206 t = u.value - 0.4375Q; 207 p = ((((((((((P10 * t 208 + P9) * t 209 + P8) * t 210 + P7) * t 211 + P6) * t 212 + P5) * t 213 + P4) * t 214 + P3) * t 215 + P2) * t 216 + P1) * t 217 + P0) * t; 218 219 q = (((((((((t 220 + Q9) * t 221 + Q8) * t 222 + Q7) * t 223 + Q6) * t 224 + Q5) * t 225 + Q4) * t 226 + Q3) * t 227 + Q2) * t 228 + Q1) * t 229 + Q0; 230 r = p / q; 231 if (sign & 0x80000000) 232 r = pimacosr4375 - r; 233 else 234 r = acosr4375 + r; 235 return r; 236 } 237 else if (ix < 0x3ffe4000) /* |x| < 0.625 */ 238 { 239 t = u.value - 0.5625Q; 240 p = ((((((((((rS10 * t 241 + rS9) * t 242 + rS8) * t 243 + rS7) * t 244 + rS6) * t 245 + rS5) * t 246 + rS4) * t 247 + rS3) * t 248 + rS2) * t 249 + rS1) * t 250 + rS0) * t; 251 252 q = (((((((((t 253 + sS9) * t 254 + sS8) * t 255 + sS7) * t 256 + sS6) * t 257 + sS5) * t 258 + sS4) * t 259 + sS3) * t 260 + sS2) * t 261 + sS1) * t 262 + sS0; 263 if (sign & 0x80000000) 264 r = pimacosr5625 - p / q; 265 else 266 r = acosr5625 + p / q; 267 return r; 268 } 269 else 270 { /* |x| >= .625 */ 271 z = (one - u.value) * 0.5; 272 s = sqrtq (z); 273 /* Compute an extended precision square root from 274 the Newton iteration s -> 0.5 * (s + z / s). 275 The change w from s to the improved value is 276 w = 0.5 * (s + z / s) - s = (s^2 + z)/2s - s = (z - s^2)/2s. 277 Express s = f1 + f2 where f1 * f1 is exactly representable. 278 w = (z - s^2)/2s = (z - f1^2 - 2 f1 f2 - f2^2)/2s . 279 s + w has extended precision. */ 280 u.value = s; 281 u.words32.w2 = 0; 282 u.words32.w3 = 0; 283 f2 = s - u.value; 284 w = z - u.value * u.value; 285 w = w - 2.0 * u.value * f2; 286 w = w - f2 * f2; 287 w = w / (2.0 * s); 288 /* Arcsine of s. */ 289 p = (((((((((pS9 * z 290 + pS8) * z 291 + pS7) * z 292 + pS6) * z 293 + pS5) * z 294 + pS4) * z 295 + pS3) * z 296 + pS2) * z 297 + pS1) * z 298 + pS0) * z; 299 q = (((((((( z 300 + qS8) * z 301 + qS7) * z 302 + qS6) * z 303 + qS5) * z 304 + qS4) * z 305 + qS3) * z 306 + qS2) * z 307 + qS1) * z 308 + qS0; 309 r = s + (w + s * p / q); 310 311 if (sign & 0x80000000) 312 w = pio2_hi + (pio2_lo - r); 313 else 314 w = r; 315 return 2.0 * w; 316 } 317 } 318