Branch data Line data Source code
1 : : /* e_acosf.c -- float version of e_acos.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 : : #include "cdefs-compat.h"
17 : : //__FBSDID("$FreeBSD: src/lib/msun/src/e_acosf.c,v 1.11 2008/08/03 17:39:54 das Exp $");
18 : :
19 : : #include <openlibm_math.h>
20 : :
21 : : #include "math_private.h"
22 : :
23 : : static const float
24 : : one = 1.0000000000e+00, /* 0x3F800000 */
25 : : pi = 3.1415925026e+00, /* 0x40490fda */
26 : : pio2_hi = 1.5707962513e+00; /* 0x3fc90fda */
27 : : static volatile float
28 : : pio2_lo = 7.5497894159e-08; /* 0x33a22168 */
29 : : static const float
30 : : pS0 = 1.6666586697e-01,
31 : : pS1 = -4.2743422091e-02,
32 : : pS2 = -8.6563630030e-03,
33 : : qS1 = -7.0662963390e-01;
34 : :
35 : : OLM_DLLEXPORT float
36 : 13 : __ieee754_acosf(float x)
37 : : {
38 : : float z,p,q,r,w,s,c,df;
39 : : int32_t hx,ix;
40 : 13 : GET_FLOAT_WORD(hx,x);
41 : 13 : ix = hx&0x7fffffff;
42 [ + + ]: 13 : if(ix>=0x3f800000) { /* |x| >= 1 */
43 [ + + ]: 7 : if(ix==0x3f800000) { /* |x| == 1 */
44 [ + + ]: 2 : if(hx>0) return 0.0; /* acos(1) = 0 */
45 : 1 : else return pi+(float)2.0*pio2_lo; /* acos(-1)= pi */
46 : : }
47 : 5 : return (x-x)/(x-x); /* acos(|x|>1) is NaN */
48 : : }
49 [ + + ]: 6 : if(ix<0x3f000000) { /* |x| < 0.5 */
50 [ + - ]: 3 : if(ix<=0x32800000) return pio2_hi+pio2_lo;/*if|x|<2**-26*/
51 : 0 : z = x*x;
52 : 0 : p = z*(pS0+z*(pS1+z*pS2));
53 : 0 : q = one+z*qS1;
54 : 0 : r = p/q;
55 : 0 : return pio2_hi - (x - (pio2_lo-x*r));
56 [ + + ]: 3 : } else if (hx<0) { /* x < -0.5 */
57 : 1 : z = (one+x)*(float)0.5;
58 : 1 : p = z*(pS0+z*(pS1+z*pS2));
59 : 1 : q = one+z*qS1;
60 : 1 : s = sqrtf(z);
61 : 1 : r = p/q;
62 : 1 : w = r*s-pio2_lo;
63 : 1 : return pi - (float)2.0*(s+w);
64 : : } else { /* x > 0.5 */
65 : : int32_t idf;
66 : 2 : z = (one-x)*(float)0.5;
67 : 2 : s = sqrtf(z);
68 : 2 : df = s;
69 : 2 : GET_FLOAT_WORD(idf,df);
70 : 2 : SET_FLOAT_WORD(df,idf&0xfffff000);
71 : 2 : c = (z-df*df)/(s+df);
72 : 2 : p = z*(pS0+z*(pS1+z*pS2));
73 : 2 : q = one+z*qS1;
74 : 2 : r = p/q;
75 : 2 : w = r*s+c;
76 : 2 : return (float)2.0*(df+w);
77 : : }
78 : : }
|