1 /*****************************************************************************/
2 /*****************************************************************************/
3 // atanf from musl-0.9.15
4 /*****************************************************************************/
5 /*****************************************************************************/
6
7 /* origin: FreeBSD /usr/src/lib/msun/src/s_atanf.c */
8 /*
9 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
10 */
11 /*
12 * ====================================================
13 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
14 *
15 * Developed at SunPro, a Sun Microsystems, Inc. business.
16 * Permission to use, copy, modify, and distribute this
17 * software is freely granted, provided that this notice
18 * is preserved.
19 * ====================================================
20 */
21
22
23 #include "libm.h"
24
25 static const float atanhi[] = {
26 4.6364760399e-01f, /* atan(0.5)hi 0x3eed6338 */
27 7.8539812565e-01f, /* atan(1.0)hi 0x3f490fda */
28 9.8279368877e-01f, /* atan(1.5)hi 0x3f7b985e */
29 1.5707962513e+00f, /* atan(inf)hi 0x3fc90fda */
30 };
31
32 static const float atanlo[] = {
33 5.0121582440e-09f, /* atan(0.5)lo 0x31ac3769 */
34 3.7748947079e-08f, /* atan(1.0)lo 0x33222168 */
35 3.4473217170e-08f, /* atan(1.5)lo 0x33140fb4 */
36 7.5497894159e-08f, /* atan(inf)lo 0x33a22168 */
37 };
38
39 static const float aT[] = {
40 3.3333328366e-01f,
41 -1.9999158382e-01f,
42 1.4253635705e-01f,
43 -1.0648017377e-01f,
44 6.1687607318e-02f,
45 };
46
atanf(float x)47 float atanf(float x)
48 {
49 float_t w,s1,s2,z;
50 uint32_t ix,sign;
51 int id;
52
53 GET_FLOAT_WORD(ix, x);
54 sign = ix>>31;
55 ix &= 0x7fffffff;
56 if (ix >= 0x4c800000) { /* if |x| >= 2**26 */
57 if (isnan(x))
58 return x;
59 z = atanhi[3] + 0x1p-120f;
60 return sign ? -z : z;
61 }
62 if (ix < 0x3ee00000) { /* |x| < 0.4375 */
63 if (ix < 0x39800000) { /* |x| < 2**-12 */
64 if (ix < 0x00800000)
65 /* raise underflow for subnormal x */
66 FORCE_EVAL(x*x);
67 return x;
68 }
69 id = -1;
70 } else {
71 x = fabsf(x);
72 if (ix < 0x3f980000) { /* |x| < 1.1875 */
73 if (ix < 0x3f300000) { /* 7/16 <= |x| < 11/16 */
74 id = 0;
75 x = (2.0f*x - 1.0f)/(2.0f + x);
76 } else { /* 11/16 <= |x| < 19/16 */
77 id = 1;
78 x = (x - 1.0f)/(x + 1.0f);
79 }
80 } else {
81 if (ix < 0x401c0000) { /* |x| < 2.4375 */
82 id = 2;
83 x = (x - 1.5f)/(1.0f + 1.5f*x);
84 } else { /* 2.4375 <= |x| < 2**26 */
85 id = 3;
86 x = -1.0f/x;
87 }
88 }
89 }
90 /* end of argument reduction */
91 z = x*x;
92 w = z*z;
93 /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
94 s1 = z*(aT[0]+w*(aT[2]+w*aT[4]));
95 s2 = w*(aT[1]+w*aT[3]);
96 if (id < 0)
97 return x - x*(s1+s2);
98 z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
99 return sign ? -z : z;
100 }
101