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, write to the Free Software 32 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ 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 * 40 * For x in [0.5,1] 41 * asin(x) = pi/2-2*asin(sqrt((1-x)/2)) 42 * Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2; 43 * then for x>0.98 44 * asin(x) = pi/2 - 2*(s+s*z*R(z)) 45 * = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo) 46 * For x<=0.98, let pio4_hi = pio2_hi/2, then 47 * f = hi part of s; 48 * c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z) 49 * and 50 * asin(x) = pi/2 - 2*(s+s*z*R(z)) 51 * = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo) 52 * = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c)) 53 * 54 * Special cases: 55 * if x is NaN, return x itself; 56 * if |x|>1, return NaN with invalid signal. 57 * 58 */ 59 60 61#include "math.h" 62#include "math_private.h" 63 64#ifdef __STDC__ 65static const long double 66#else 67static long double 68#endif 69 one = 1.0L, 70 huge = 1.0e+4932L, 71 pio2_hi = 1.5707963267948966192021943710788178805159986950457096099853515625L, 72 pio2_lo = 2.9127320560933561582586004641843300502121E-20L, 73 pio4_hi = 7.8539816339744830960109718553940894025800E-1L, 74 75 /* coefficient for R(x^2) */ 76 77 /* asin(x) = x + x^3 pS(x^2) / qS(x^2) 78 0 <= x <= 0.5 79 peak relative error 1.9e-21 */ 80 pS0 = -1.008714657938491626019651170502036851607E1L, 81 pS1 = 2.331460313214179572063441834101394865259E1L, 82 pS2 = -1.863169762159016144159202387315381830227E1L, 83 pS3 = 5.930399351579141771077475766877674661747E0L, 84 pS4 = -6.121291917696920296944056882932695185001E-1L, 85 pS5 = 3.776934006243367487161248678019350338383E-3L, 86 87 qS0 = -6.052287947630949712886794360635592886517E1L, 88 qS1 = 1.671229145571899593737596543114258558503E2L, 89 qS2 = -1.707840117062586426144397688315411324388E2L, 90 qS3 = 7.870295154902110425886636075950077640623E1L, 91 qS4 = -1.568433562487314651121702982333303458814E1L; 92 /* 1.000000000000000000000000000000000000000E0 */ 93 94#ifdef __STDC__ 95long double 96__ieee754_asinl (long double x) 97#else 98double 99__ieee754_asinl (x) 100 long double x; 101#endif 102{ 103 long double t, w, p, q, c, r, s; 104 int32_t ix; 105 u_int32_t se, i0, i1, k; 106 107 GET_LDOUBLE_WORDS (se, i0, i1, x); 108 ix = se & 0x7fff; 109 ix = (ix << 16) | (i0 >> 16); 110 if (ix >= 0x3fff8000) 111 { /* |x|>= 1 */ 112 if (ix == 0x3fff8000 && ((i0 - 0x80000000) | i1) == 0) 113 /* asin(1)=+-pi/2 with inexact */ 114 return x * pio2_hi + x * pio2_lo; 115 return (x - x) / (x - x); /* asin(|x|>1) is NaN */ 116 } 117 else if (ix < 0x3ffe8000) 118 { /* |x|<0.5 */ 119 if (ix < 0x3fde8000) 120 { /* if |x| < 2**-33 */ 121 if (huge + x > one) 122 return x; /* return x with inexact if x!=0 */ 123 } 124 else 125 { 126 t = x * x; 127 p = 128 t * (pS0 + 129 t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5))))); 130 q = qS0 + t * (qS1 + t * (qS2 + t * (qS3 + t * (qS4 + t)))); 131 w = p / q; 132 return x + x * w; 133 } 134 } 135 /* 1> |x|>= 0.5 */ 136 w = one - fabsl (x); 137 t = w * 0.5; 138 p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5))))); 139 q = qS0 + t * (qS1 + t * (qS2 + t * (qS3 + t * (qS4 + t)))); 140 s = __ieee754_sqrtl (t); 141 if (ix >= 0x3ffef999) 142 { /* if |x| > 0.975 */ 143 w = p / q; 144 t = pio2_hi - (2.0 * (s + s * w) - pio2_lo); 145 } 146 else 147 { 148 GET_LDOUBLE_WORDS (k, i0, i1, s); 149 i1 = 0; 150 SET_LDOUBLE_WORDS (w,k,i0,i1); 151 c = (t - w * w) / (s + w); 152 r = p / q; 153 p = 2.0 * s * r - (pio2_lo - 2.0 * c); 154 q = pio4_hi - 2.0 * w; 155 t = pio4_hi - (p - q); 156 } 157 if ((se & 0x8000) == 0) 158 return t; 159 else 160 return -t; 161} 162