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 |