5 * Common Development and Distribution License (the "License").
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
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25 /*
26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
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)
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 }
|
5 * Common Development and Distribution License (the "License").
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
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25
26 /*
27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 #pragma weak __atanf = atanf
32
33
34 /*
35 * float atanf(float x);
36 * Table look-up algorithm
37 * By K.C. Ng, March 9, 1989
38 *
39 * Algorithm.
40 *
41 * The algorithm is based on atan(x)=atan(y)+atan((x-y)/(1+x*y)).
42 * We use poly1(x) to approximate atan(x) for x in [0,1/8] with
43 * error (relative)
44 * |(atan(x)-poly1(x))/x|<= 2^-115.94 long double
45 * |(atan(x)-poly1(x))/x|<= 2^-58.85 double
46 * |(atan(x)-poly1(x))/x|<= 2^-25.53 float
47 * and use poly2(x) to approximate atan(x) for x in [0,1/65] with
48 * error (absolute)
49 * |atan(x)-poly2(x)|<= 2^-122.15 long double
50 * |atan(x)-poly2(x)|<= 2^-64.79 double
51 * |atan(x)-poly2(x)|<= 2^-35.36 float
52 * and use poly3(x) to approximate atan(x) for x in [1/8,7/16] with
53 * error (relative, on for single precision)
70 *
71 * (4). Now x is in (0.125, 8)
72 * Find y that match x to 4.5 bit after binary (easy).
73 * If iy is the high word of y, then
74 * single : j = (iy - 0x3e000000) >> 19
75 * (single is modified to (iy-0x3f000000)>>19)
76 * double : j = (iy - 0x3fc00000) >> 16
77 * quad : j = (iy - 0x3ffc0000) >> 12
78 *
79 * Let s = (x-y)/(1+x*y). Then
80 * atan(x) = atan(y) + poly1(s)
81 * = _TBL_r_atan_hi[j] + (_TBL_r_atan_lo[j] + poly2(s) )
82 *
83 * Note. |s| <= 1.5384615385e-02 = 1/65. Maxium occurs at x = 1.03125
84 *
85 */
86
87 #include "libm.h"
88
89 extern const float _TBL_r_atan_hi[], _TBL_r_atan_lo[];
90 static const float 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
103
104 float
105 atanf(float xx)
106 {
107 float x, y, z, r, p, s;
108 volatile double dummy __unused;
109 int ix, iy, sign, j;
110
111 x = xx;
112 ix = *(int *)&x;
113 sign = ix & 0x80000000;
114 ix ^= sign;
115
116 /* for |x| < 1/8 */
117 if (ix < 0x3e000000) {
118 if (ix < 0x38800000) { /* if |x| < 2**(-prec/2-2) */
119 dummy = big + x; /* get inexact flag if x != 0 */
120 #ifdef lint
121 dummy = dummy;
122 #endif
123 return (x);
124 }
125
126 z = x * x;
127
128 if (ix < 0x3c000000) { /* if |x| < 2**(-prec/4-1) */
129 x = x + (x * z) * p1;
130 return (x);
131 } else {
132 x = x + (x * z) * (p1 + z * p2);
133 return (x);
134 }
135 }
136
137 /* for |x| >= 8.0 */
138 if (ix >= 0x41000000) {
139 *(int *)&x = ix;
140
141 if (ix < 0x42820000) { /* x < 65 */
142 r = one / x;
143 z = r * r;
144 y = r * (one + z * (p1 + z * p2)); /* poly1 */
145 y -= pio2lo;
146 } else if (ix < 0x44800000) { /* x < 2**(prec/3+2) */
147 r = one / x;
148 z = r * r;
149 y = r * (one + z * q1); /* poly2 */
150 y -= pio2lo;
151 } else if (ix < 0x4c800000) { /* x < 2**(prec+2) */
152 y = one / x - pio2lo;
153 } else if (ix < 0x7f800000) { /* x < inf */
154 y = -pio2lo;
155 } else { /* x is inf or NaN */
156 if (ix > 0x7f800000)
157 return (x * x); /* - -> * for Cheetah */
158
159 y = -pio2lo;
160 }
161
162 if (sign == 0)
163 x = pio2hi - y;
164 else
165 x = y - pio2hi;
166
167 return (x);
168 }
169
170 /* now x is between 1/8 and 8 */
171 if (ix < 0x3f000000) { /* between 1/8 and 1/2 */
172 z = x * x;
173 x = x + (x * z) * (a1 + z * (a2 + z * (a3 + z * (a4 + z *
174 a5))));
175 return (x);
176 }
177
178 *(int *)&x = ix;
179 iy = (ix + 0x00040000) & 0x7ff80000;
180 *(int *)&y = iy;
181 j = (iy - 0x3f000000) >> 19;
182
183 if (ix == iy) {
184 p = x - y; /* p=0.0 */
185 } else {
186 if (sign == 0)
187 s = (x - y) / (one + x * y);
188 else
189 s = (y - x) / (one + x * y);
190
191 z = s * s;
192 p = s * (one + z * q1);
193 }
194
195 if (sign == 0) {
196 r = p + _TBL_r_atan_lo[j];
197 x = r + _TBL_r_atan_hi[j];
198 } else {
199 r = p - _TBL_r_atan_lo[j];
200 x = r - _TBL_r_atan_hi[j];
201 }
202
203 return (x);
204 }
|