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