]>
Commit | Line | Data |
---|---|---|
f7eac6eb RM |
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 | |
cc3fa755 | 11 | * software is freely granted, provided that this notice |
f7eac6eb RM |
12 | * is preserved. |
13 | * ==================================================== | |
14 | */ | |
15 | ||
9b7ee67e | 16 | /* |
0ac5ae23 | 17 | Modifications for single precision expansion are |
cc7375ce | 18 | Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov> |
0ac5ae23 | 19 | and are incorporated herein by permission of the author. The author |
9cd2726c | 20 | reserves the right to distribute this material elsewhere under different |
0ac5ae23 | 21 | copying permissions. These modifications are distributed here under |
9cd2726c | 22 | the following terms: |
9b7ee67e | 23 | |
cc7375ce RM |
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 | |
59ba27a6 PE |
35 | License along with this library; if not, see |
36 | <http://www.gnu.org/licenses/>. */ | |
9b7ee67e | 37 | |
f7eac6eb RM |
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 "math.h" | |
43 | #include "math_private.h" | |
44 | ||
cc3fa755 | 45 | static const float |
f7eac6eb RM |
46 | one = 1.0000000000e+00, /* 0x3F800000 */ |
47 | huge = 1.000e+30, | |
9b7ee67e UD |
48 | |
49 | pio2_hi = 1.57079637050628662109375f, | |
50 | pio2_lo = -4.37113900018624283e-8f, | |
51 | pio4_hi = 0.785398185253143310546875f, | |
52 | ||
53 | /* asin x = x + x^3 p(x^2) | |
54 | -0.5 <= x <= 0.5; | |
55 | Peak relative error 4.8e-9 */ | |
56 | p0 = 1.666675248e-1f, | |
57 | p1 = 7.495297643e-2f, | |
58 | p2 = 4.547037598e-2f, | |
59 | p3 = 2.417951451e-2f, | |
60 | p4 = 4.216630880e-2f; | |
f7eac6eb | 61 | |
8db21882 | 62 | float __ieee754_asinf(float x) |
f7eac6eb RM |
63 | { |
64 | float t,w,p,q,c,r,s; | |
65 | int32_t hx,ix; | |
66 | GET_FLOAT_WORD(hx,x); | |
67 | ix = hx&0x7fffffff; | |
68 | if(ix==0x3f800000) { | |
69 | /* asin(1)=+-pi/2 with inexact */ | |
cc3fa755 | 70 | return x*pio2_hi+x*pio2_lo; |
f7eac6eb | 71 | } else if(ix> 0x3f800000) { /* |x|>= 1 */ |
cc3fa755 | 72 | return (x-x)/(x-x); /* asin(|x|>1) is NaN */ |
f7eac6eb RM |
73 | } else if (ix<0x3f000000) { /* |x|<0.5 */ |
74 | if(ix<0x32000000) { /* if |x| < 2**-27 */ | |
75 | if(huge+x>one) return x;/* return x with inexact if x!=0*/ | |
cc3fa755 | 76 | } else { |
f7eac6eb | 77 | t = x*x; |
9b7ee67e | 78 | w = t * (p0 + t * (p1 + t * (p2 + t * (p3 + t * p4)))); |
f7eac6eb | 79 | return x+x*w; |
cc3fa755 | 80 | } |
f7eac6eb RM |
81 | } |
82 | /* 1> |x|>= 0.5 */ | |
83 | w = one-fabsf(x); | |
9b7ee67e UD |
84 | t = w*0.5f; |
85 | p = t * (p0 + t * (p1 + t * (p2 + t * (p3 + t * p4)))); | |
f7eac6eb | 86 | s = __ieee754_sqrtf(t); |
0ac5ae23 | 87 | if(ix>=0x3F79999A) { /* if |x| > 0.975 */ |
9b7ee67e | 88 | t = pio2_hi-(2.0f*(s+s*p)-pio2_lo); |
f7eac6eb RM |
89 | } else { |
90 | int32_t iw; | |
91 | w = s; | |
92 | GET_FLOAT_WORD(iw,w); | |
93 | SET_FLOAT_WORD(w,iw&0xfffff000); | |
94 | c = (t-w*w)/(s+w); | |
9b7ee67e UD |
95 | r = p; |
96 | p = 2.0f*s*r-(pio2_lo-2.0f*c); | |
97 | q = pio4_hi-2.0f*w; | |
f7eac6eb | 98 | t = pio4_hi-(p-q); |
cc3fa755 UD |
99 | } |
100 | if(hx>0) return t; else return -t; | |
f7eac6eb | 101 | } |
0ac5ae23 | 102 | strong_alias (__ieee754_asinf, __asinf_finite) |