| OLD | NEW |
| 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_asinl.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_asinl.c */ |
| 2 /* | 2 /* |
| 3 * ==================================================== | 3 * ==================================================== |
| 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
| 5 * | 5 * |
| 6 * Developed at SunSoft, a Sun Microsystems, Inc. business. | 6 * Developed at SunSoft, a Sun Microsystems, Inc. business. |
| 7 * Permission to use, copy, modify, and distribute this | 7 * Permission to use, copy, modify, and distribute this |
| 8 * software is freely granted, provided that this notice | 8 * software is freely granted, provided that this notice |
| 9 * is preserved. | 9 * is preserved. |
| 10 * ==================================================== | 10 * ==================================================== |
| 11 */ | 11 */ |
| 12 /* | 12 /* |
| 13 * See comments in asin.c. | 13 * See comments in asin.c. |
| 14 * Converted to long double by David Schultz <das@FreeBSD.ORG>. | 14 * Converted to long double by David Schultz <das@FreeBSD.ORG>. |
| 15 */ | 15 */ |
| 16 | 16 |
| 17 #include "libm.h" | 17 #include "libm.h" |
| 18 | 18 |
| 19 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 | 19 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 |
| 20 long double asinl(long double x) | 20 long double asinl(long double x) { |
| 21 { | 21 return asin(x); |
| 22 » return asin(x); | |
| 23 } | 22 } |
| 24 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 | 23 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 |
| 25 #include "__invtrigl.h" | 24 #include "__invtrigl.h" |
| 26 #if LDBL_MANT_DIG == 64 | 25 #if LDBL_MANT_DIG == 64 |
| 27 #define CLOSETO1(u) (u.i.m>>56 >= 0xf7) | 26 #define CLOSETO1(u) (u.i.m >> 56 >= 0xf7) |
| 28 #define CLEARBOTTOM(u) (u.i.m &= -1ULL << 32) | 27 #define CLEARBOTTOM(u) (u.i.m &= -1ULL << 32) |
| 29 #elif LDBL_MANT_DIG == 113 | 28 #elif LDBL_MANT_DIG == 113 |
| 30 #define CLOSETO1(u) (u.i.top >= 0xee00) | 29 #define CLOSETO1(u) (u.i.top >= 0xee00) |
| 31 #define CLEARBOTTOM(u) (u.i.lo = 0) | 30 #define CLEARBOTTOM(u) (u.i.lo = 0) |
| 32 #endif | 31 #endif |
| 33 | 32 |
| 34 long double asinl(long double x) | 33 long double asinl(long double x) { |
| 35 { | 34 union ldshape u = {x}; |
| 36 » union ldshape u = {x}; | 35 long double z, r, s; |
| 37 » long double z, r, s; | 36 uint16_t e = u.i.se & 0x7fff; |
| 38 » uint16_t e = u.i.se & 0x7fff; | 37 int sign = u.i.se >> 15; |
| 39 » int sign = u.i.se >> 15; | |
| 40 | 38 |
| 41 » if (e >= 0x3fff) { /* |x| >= 1 or nan */ | 39 if (e >= 0x3fff) { /* |x| >= 1 or nan */ |
| 42 » » /* asin(+-1)=+-pi/2 with inexact */ | 40 /* asin(+-1)=+-pi/2 with inexact */ |
| 43 » » if (x == 1 || x == -1) | 41 if (x == 1 || x == -1) |
| 44 » » » return x*pio2_hi + 0x1p-120f; | 42 return x * pio2_hi + 0x1p-120f; |
| 45 » » return 0/(x-x); | 43 return 0 / (x - x); |
| 46 » } | 44 } |
| 47 » if (e < 0x3fff - 1) { /* |x| < 0.5 */ | 45 if (e < 0x3fff - 1) { /* |x| < 0.5 */ |
| 48 » » if (e < 0x3fff - (LDBL_MANT_DIG+1)/2) { | 46 if (e < 0x3fff - (LDBL_MANT_DIG + 1) / 2) { |
| 49 » » » /* return x with inexact if x!=0 */ | 47 /* return x with inexact if x!=0 */ |
| 50 » » » FORCE_EVAL(x + 0x1p120f); | 48 FORCE_EVAL(x + 0x1p120f); |
| 51 » » » return x; | 49 return x; |
| 52 » » } | 50 } |
| 53 » » return x + x*__invtrigl_R(x*x); | 51 return x + x * __invtrigl_R(x * x); |
| 54 » } | 52 } |
| 55 » /* 1 > |x| >= 0.5 */ | 53 /* 1 > |x| >= 0.5 */ |
| 56 » z = (1.0 - fabsl(x))*0.5; | 54 z = (1.0 - fabsl(x)) * 0.5; |
| 57 » s = sqrtl(z); | 55 s = sqrtl(z); |
| 58 » r = __invtrigl_R(z); | 56 r = __invtrigl_R(z); |
| 59 » if (CLOSETO1(u)) { | 57 if (CLOSETO1(u)) { |
| 60 » » x = pio2_hi - (2*(s+s*r)-pio2_lo); | 58 x = pio2_hi - (2 * (s + s * r) - pio2_lo); |
| 61 » } else { | 59 } else { |
| 62 » » long double f, c; | 60 long double f, c; |
| 63 » » u.f = s; | 61 u.f = s; |
| 64 » » CLEARBOTTOM(u); | 62 CLEARBOTTOM(u); |
| 65 » » f = u.f; | 63 f = u.f; |
| 66 » » c = (z - f*f)/(s + f); | 64 c = (z - f * f) / (s + f); |
| 67 » » x = 0.5*pio2_hi-(2*s*r - (pio2_lo-2*c) - (0.5*pio2_hi-2*f)); | 65 x = 0.5 * pio2_hi - |
| 68 » } | 66 (2 * s * r - (pio2_lo - 2 * c) - (0.5 * pio2_hi - 2 * f)); |
| 69 » return sign ? -x : x; | 67 } |
| 68 return sign ? -x : x; |
| 70 } | 69 } |
| 71 #endif | 70 #endif |
| OLD | NEW |