11210 libm should be cstyle(1ONBLD) clean

1 /*
3  *
4  * The contents of this file are subject to the terms of the
6  * You may not use this file except in compliance with the License.
7  *
8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9  * or http://www.opensolaris.org/os/licensing.
10  * See the License for the specific language governing permissions
11  * and limitations under the License.
12  *
13  * When distributing Covered Code, include this CDDL HEADER in each
15  * If applicable, add the following below this CDDL HEADER, with the
16  * fields enclosed by brackets "[]" replaced with your own identifying
18  *
20  */
21
22 /*
24  */

25 /*
27  * Use is subject to license terms.
28  */
29
30 #pragma weak __atanf = atanf
31
32 /* INDENT OFF */
33 /*
34  * float atanf(float x);
35  * Table look-up algorithm
36  * By K.C. Ng, March 9, 1989
37  *
38  * Algorithm.
39  *
40  * The algorithm is based on atan(x)=atan(y)+atan((x-y)/(1+x*y)).
41  * We use poly1(x) to approximate atan(x) for x in [0,1/8] with
42  * error (relative)
43  *      |(atan(x)-poly1(x))/x|<= 2^-115.94   long double
44  *      |(atan(x)-poly1(x))/x|<= 2^-58.85    double
45  *      |(atan(x)-poly1(x))/x|<= 2^-25.53    float
46  * and use poly2(x) to approximate atan(x) for x in [0,1/65] with
47  * error (absolute)
48  *      |atan(x)-poly2(x)|<= 2^-122.15               long double
49  *      |atan(x)-poly2(x)|<= 2^-64.79                double
50  *      |atan(x)-poly2(x)|<= 2^-35.36                float
51  * and use poly3(x) to approximate atan(x) for x in [1/8,7/16] with
52  * error (relative, on for single precision)
53  *      |(atan(x)-poly1(x))/x|<= 2^-25.53    float
54  *
55  * Here poly1-3 are odd polynomial with the following form:
56  *              x + x^3*(a1+x^2*(a2+...))
57  *
58  * (0). Purge off Inf and NaN and 0
59  * (1). Reduce x to positive by atan(x) = -atan(-x).
60  * (2). For x <= 1/8, use
61  *      (2.1) if x < 2^(-prec/2-2), atan(x) =  x  with inexact
62  *      (2.2) Otherwise
63  *              atan(x) = poly1(x)
64  * (3). For x >= 8 then
65  *      (3.1) if x >= 2^(prec+2),   atan(x) = atan(inf) - pio2lo
66  *      (3.2) if x >= 2^(prec/3+2), atan(x) = atan(inf) - 1/x
67  *      (3.3) if x >  65,           atan(x) = atan(inf) - poly2(1/x)
68  *      (3.4) Otherwise,            atan(x) = atan(inf) - poly1(1/x)
69  *
70  * (4). Now x is in (0.125, 8)
71  *      Find y that match x to 4.5 bit after binary (easy).
72  *      If iy is the high word of y, then
73  *              single : j = (iy - 0x3e000000) >> 19
74  *              (single is modified to (iy-0x3f000000)>>19)
75  *              double : j = (iy - 0x3fc00000) >> 16
76  *              quad   : j = (iy - 0x3ffc0000) >> 12
77  *
78  *      Let s = (x-y)/(1+x*y). Then
79  *      atan(x) = atan(y) + poly1(s)
80  *              = _TBL_r_atan_hi[j] + (_TBL_r_atan_lo[j] + poly2(s) )
81  *
82  *      Note. |s| <= 1.5384615385e-02 = 1/65. Maxium occurs at x = 1.03125
83  *
84  */
85
86 #include "libm.h"
87
88 extern const float _TBL_r_atan_hi[], _TBL_r_atan_lo[];
89 static const float
90         big     =   1.0e37F,
91         one     =   1.0F,
92         p1      =  -3.333185951111688247225368498733544672172e-0001F,
93         p2      =   1.969352894213455405211341983203180636021e-0001F,
94         q1      =  -3.332921964095646819563419704110132937456e-0001F,
95         a1      =  -3.333323465223893614063523351509338934592e-0001F,
96         a2      =   1.999425625935277805494082274808174062403e-0001F,
97         a3      =  -1.417547090509737780085769846290301788559e-0001F,
98         a4      =   1.016250813871991983097273733227432685084e-0001F,
99         a5      =  -5.137023693688358515753093811791755221805e-0002F,
100         pio2hi  =   1.570796371e+0000F,
101         pio2lo  =  -4.371139000e-0008F;
102 /* INDENT ON */
103
104 float
105 atanf(float xx) {

106         float x, y, z, r, p, s;
107         volatile double dummy __unused;
108         int ix, iy, sign, j;
109
110         x = xx;
111         ix = *(int *) &x;
112         sign = ix & 0x80000000;
113         ix ^= sign;
114
115         /* for |x| < 1/8 */
116         if (ix < 0x3e000000) {
117                 if (ix < 0x38800000) {       /* if |x| < 2**(-prec/2-2) */
118                         dummy = big + x;        /* get inexact flag if x != 0 */
119 #ifdef lint
120                         dummy = dummy;
121 #endif
122                         return (x);
123                 }

124                 z = x * x;

125                 if (ix < 0x3c000000) {       /* if |x| < 2**(-prec/4-1) */
126                         x = x + (x * z) * p1;
127                         return (x);
128                 } else {
129                         x = x + (x * z) * (p1 + z * p2);
130                         return (x);
131                 }
132         }
133
134         /* for |x| >= 8.0 */
135         if (ix >= 0x41000000) {
136                 *(int *) &x = ix;

137                 if (ix < 0x42820000) {       /* x <  65 */
138                         r = one / x;
139                         z = r * r;
140                         y = r * (one + z * (p1 + z * p2));      /* poly1 */
141                         y -= pio2lo;
142                 } else if (ix < 0x44800000) {        /* x <  2**(prec/3+2) */
143                         r = one / x;
144                         z = r * r;
145                         y = r * (one + z * q1); /* poly2 */
146                         y -= pio2lo;
147                 } else if (ix < 0x4c800000) {        /* x <  2**(prec+2) */
148                         y = one / x - pio2lo;
149                 } else if (ix < 0x7f800000) {        /* x <  inf */
150                         y = -pio2lo;
151                 } else {                /* x is inf or NaN */
152                         if (ix > 0x7f800000) {
153                                 return (x * x); /* - -> * for Cheetah */
154                         }
155                         y = -pio2lo;
156                 }
157
158                 if (sign == 0)
159                         x = pio2hi - y;
160                 else
161                         x = y - pio2hi;

162                 return (x);
163         }
164
165
166         /* now x is between 1/8 and 8 */
167         if (ix < 0x3f000000) {       /* between 1/8 and 1/2 */
168                 z = x * x;
169                 x = x + (x * z) * (a1 + z * (a2 + z * (a3 + z * (a4 +
170                         z * a5))));
171                 return (x);
172         }
173         *(int *) &x = ix;

174         iy = (ix + 0x00040000) & 0x7ff80000;
175         *(int *) &y = iy;
176         j = (iy - 0x3f000000) >> 19;
177
178         if (ix == iy)
179                 p = x - y;      /* p=0.0 */
180         else {
181                 if (sign == 0)
182                         s = (x - y) / (one + x * y);
183                 else
184                         s = (y - x) / (one + x * y);

185                 z = s * s;
186                 p = s * (one + z * q1);
187         }

188         if (sign == 0) {
189                 r = p + _TBL_r_atan_lo[j];
190                 x = r + _TBL_r_atan_hi[j];
191         } else {
192                 r = p - _TBL_r_atan_lo[j];
193                 x = r - _TBL_r_atan_hi[j];
194         }

195         return (x);
196 }
--- EOF ---