1 | /* e_asinf.c -- float version of e_asin.c. |
2 | * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
3 | */ |
4 | |
5 | /* |
6 | * ==================================================== |
7 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
8 | * |
9 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
10 | * Permission to use, copy, modify, and distribute this |
11 | * software is freely granted, provided that this notice |
12 | * is preserved. |
13 | * ==================================================== |
14 | */ |
15 | |
16 | /* |
17 | Modifications for single precision expansion are |
18 | Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov> |
19 | and are incorporated herein by permission of the author. The author |
20 | reserves the right to distribute this material elsewhere under different |
21 | copying permissions. These modifications are distributed here under |
22 | the following terms: |
23 | |
24 | This library is free software; you can redistribute it and/or |
25 | modify it under the terms of the GNU Lesser General Public |
26 | License as published by the Free Software Foundation; either |
27 | version 2.1 of the License, or (at your option) any later version. |
28 | |
29 | This library is distributed in the hope that it will be useful, |
30 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
31 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
32 | Lesser General Public License for more details. |
33 | |
34 | You should have received a copy of the GNU Lesser General Public |
35 | License along with this library; if not, see |
36 | <http://www.gnu.org/licenses/>. */ |
37 | |
38 | #if defined(LIBM_SCCS) && !defined(lint) |
39 | static char rcsid[] = "$NetBSD: e_asinf.c,v 1.5 1995/05/12 04:57:25 jtc Exp $" ; |
40 | #endif |
41 | |
42 | #include <float.h> |
43 | #include <math.h> |
44 | #include <math_private.h> |
45 | |
46 | static const float |
47 | one = 1.0000000000e+00, /* 0x3F800000 */ |
48 | huge = 1.000e+30, |
49 | |
50 | pio2_hi = 1.57079637050628662109375f, |
51 | pio2_lo = -4.37113900018624283e-8f, |
52 | pio4_hi = 0.785398185253143310546875f, |
53 | |
54 | /* asin x = x + x^3 p(x^2) |
55 | -0.5 <= x <= 0.5; |
56 | Peak relative error 4.8e-9 */ |
57 | p0 = 1.666675248e-1f, |
58 | p1 = 7.495297643e-2f, |
59 | p2 = 4.547037598e-2f, |
60 | p3 = 2.417951451e-2f, |
61 | p4 = 4.216630880e-2f; |
62 | |
63 | float __ieee754_asinf(float x) |
64 | { |
65 | float t,w,p,q,c,r,s; |
66 | int32_t hx,ix; |
67 | GET_FLOAT_WORD(hx,x); |
68 | ix = hx&0x7fffffff; |
69 | if(ix==0x3f800000) { |
70 | /* asin(1)=+-pi/2 with inexact */ |
71 | return x*pio2_hi+x*pio2_lo; |
72 | } else if(ix> 0x3f800000) { /* |x|>= 1 */ |
73 | return (x-x)/(x-x); /* asin(|x|>1) is NaN */ |
74 | } else if (ix<0x3f000000) { /* |x|<0.5 */ |
75 | if(ix<0x32000000) { /* if |x| < 2**-27 */ |
76 | math_check_force_underflow (x); |
77 | if(huge+x>one) return x;/* return x with inexact if x!=0*/ |
78 | } else { |
79 | t = x*x; |
80 | w = t * (p0 + t * (p1 + t * (p2 + t * (p3 + t * p4)))); |
81 | return x+x*w; |
82 | } |
83 | } |
84 | /* 1> |x|>= 0.5 */ |
85 | w = one-fabsf(x); |
86 | t = w*0.5f; |
87 | p = t * (p0 + t * (p1 + t * (p2 + t * (p3 + t * p4)))); |
88 | s = __ieee754_sqrtf(t); |
89 | if(ix>=0x3F79999A) { /* if |x| > 0.975 */ |
90 | t = pio2_hi-(2.0f*(s+s*p)-pio2_lo); |
91 | } else { |
92 | int32_t iw; |
93 | w = s; |
94 | GET_FLOAT_WORD(iw,w); |
95 | SET_FLOAT_WORD(w,iw&0xfffff000); |
96 | c = (t-w*w)/(s+w); |
97 | r = p; |
98 | p = 2.0f*s*r-(pio2_lo-2.0f*c); |
99 | q = pio4_hi-2.0f*w; |
100 | t = pio4_hi-(p-q); |
101 | } |
102 | if(hx>0) return t; else return -t; |
103 | } |
104 | strong_alias (__ieee754_asinf, __asinf_finite) |
105 | |