| e_asin.c (117912) | e_asin.c (141296) |
|---|---|
| 1/* @(#)e_asin.c 5.1 93/09/24 */ | 1 2/* @(#)e_asin.c 1.3 95/01/18 */ |
| 2/* 3 * ==================================================== 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5 * | 3/* 4 * ==================================================== 5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6 * |
| 6 * Developed at SunPro, a Sun Microsystems, Inc. business. | 7 * Developed at SunSoft, a Sun Microsystems, Inc. business. |
| 7 * Permission to use, copy, modify, and distribute this | 8 * Permission to use, copy, modify, and distribute this |
| 8 * software is freely granted, provided that this notice | 9 * software is freely granted, provided that this notice |
| 9 * is preserved. 10 * ==================================================== 11 */ 12 13#ifndef lint | 10 * is preserved. 11 * ==================================================== 12 */ 13 14#ifndef lint |
| 14static char rcsid[] = "$FreeBSD: head/lib/msun/src/e_asin.c 117912 2003-07-23 04:53:47Z peter $"; | 15static char rcsid[] = "$FreeBSD: head/lib/msun/src/e_asin.c 141296 2005-02-04 18:26:06Z das $"; |
| 15#endif 16 17/* __ieee754_asin(x) | 16#endif 17 18/* __ieee754_asin(x) |
| 18 * Method : | 19 * Method : |
| 19 * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ... 20 * we approximate asin(x) on [0,0.5] by 21 * asin(x) = x + x*x^2*R(x^2) 22 * where | 20 * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ... 21 * we approximate asin(x) on [0,0.5] by 22 * asin(x) = x + x*x^2*R(x^2) 23 * where |
| 23 * R(x^2) is a rational approximation of (asin(x)-x)/x^3 | 24 * R(x^2) is a rational approximation of (asin(x)-x)/x^3 |
| 24 * and its remez error is bounded by 25 * |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75) 26 * 27 * For x in [0.5,1] 28 * asin(x) = pi/2-2*asin(sqrt((1-x)/2)) 29 * Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2; 30 * then for x>0.98 31 * asin(x) = pi/2 - 2*(s+s*z*R(z)) --- 41 unchanged lines hidden (view full) --- 73 int32_t hx,ix; 74 GET_HIGH_WORD(hx,x); 75 ix = hx&0x7fffffff; 76 if(ix>= 0x3ff00000) { /* |x|>= 1 */ 77 u_int32_t lx; 78 GET_LOW_WORD(lx,x); 79 if(((ix-0x3ff00000)|lx)==0) 80 /* asin(1)=+-pi/2 with inexact */ | 25 * and its remez error is bounded by 26 * |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75) 27 * 28 * For x in [0.5,1] 29 * asin(x) = pi/2-2*asin(sqrt((1-x)/2)) 30 * Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2; 31 * then for x>0.98 32 * asin(x) = pi/2 - 2*(s+s*z*R(z)) --- 41 unchanged lines hidden (view full) --- 74 int32_t hx,ix; 75 GET_HIGH_WORD(hx,x); 76 ix = hx&0x7fffffff; 77 if(ix>= 0x3ff00000) { /* |x|>= 1 */ 78 u_int32_t lx; 79 GET_LOW_WORD(lx,x); 80 if(((ix-0x3ff00000)|lx)==0) 81 /* asin(1)=+-pi/2 with inexact */ |
| 81 return x*pio2_hi+x*pio2_lo; 82 return (x-x)/(x-x); /* asin(|x|>1) is NaN */ | 82 return x*pio2_hi+x*pio2_lo; 83 return (x-x)/(x-x); /* asin(|x|>1) is NaN */ |
| 83 } else if (ix<0x3fe00000) { /* |x|<0.5 */ 84 if(ix<0x3e400000) { /* if |x| < 2**-27 */ 85 if(huge+x>one) return x;/* return x with inexact if x!=0*/ | 84 } else if (ix<0x3fe00000) { /* |x|<0.5 */ 85 if(ix<0x3e400000) { /* if |x| < 2**-27 */ 86 if(huge+x>one) return x;/* return x with inexact if x!=0*/ |
| 86 } else | 87 } else |
| 87 t = x*x; 88 p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); 89 q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); 90 w = p/q; 91 return x+x*w; 92 } 93 /* 1> |x|>= 0.5 */ 94 w = one-fabs(x); 95 t = w*0.5; 96 p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); 97 q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); | 88 t = x*x; 89 p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); 90 q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); 91 w = p/q; 92 return x+x*w; 93 } 94 /* 1> |x|>= 0.5 */ 95 w = one-fabs(x); 96 t = w*0.5; 97 p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); 98 q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); |
| 98 s = __ieee754_sqrt(t); | 99 s = sqrt(t); |
| 99 if(ix>=0x3FEF3333) { /* if |x| > 0.975 */ 100 w = p/q; 101 t = pio2_hi-(2.0*(s+s*w)-pio2_lo); 102 } else { 103 w = s; 104 SET_LOW_WORD(w,0); 105 c = (t-w*w)/(s+w); 106 r = p/q; 107 p = 2.0*s*r-(pio2_lo-2.0*c); 108 q = pio4_hi-2.0*w; 109 t = pio4_hi-(p-q); | 100 if(ix>=0x3FEF3333) { /* if |x| > 0.975 */ 101 w = p/q; 102 t = pio2_hi-(2.0*(s+s*w)-pio2_lo); 103 } else { 104 w = s; 105 SET_LOW_WORD(w,0); 106 c = (t-w*w)/(s+w); 107 r = p/q; 108 p = 2.0*s*r-(pio2_lo-2.0*c); 109 q = pio4_hi-2.0*w; 110 t = pio4_hi-(p-q); |
| 110 } 111 if(hx>0) return t; else return -t; | 111 } 112 if(hx>0) return t; else return -t; |
| 112} | 113} |