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 __atanl = atanl
31
32 /*
33 * atanl(x)
34 * Table look-up algorithm
35 * By K.C. Ng, March 9, 1989
36 *
37 * Algorithm.
38 *
39 * The algorithm is based on atan(x)=atan(y)+atan((x-y)/(1+x*y)).
40 * We use poly1(x) to approximate atan(x) for x in [0,1/8] with
41 * error (relative)
42 * |(atan(x)-poly1(x))/x|<= 2^-115.94 long double
43 * |(atan(x)-poly1(x))/x|<= 2^-58.85 double
44 * |(atan(x)-poly1(x))/x|<= 2^-25.53 float
63 * (3.4) Otherwise, atan(x) = atan(inf) - poly1(1/x)
64 *
65 * (4). Now x is in (0.125, 8)
66 * Find y that match x to 4.5 bit after binary (easy).
67 * If iy is the high word of y, then
68 * single : j = (iy - 0x3e000000) >> 19
69 * double : j = (iy - 0x3fc00000) >> 16
70 * quad : j = (iy - 0x3ffc0000) >> 12
71 *
72 * Let s = (x-y)/(1+x*y). Then
73 * atan(x) = atan(y) + poly1(s)
74 * = _TBL_atanl_hi[j] + (_TBL_atanl_lo[j] + poly2(s) )
75 *
76 * Note. |s| <= 1.5384615385e-02 = 1/65. Maxium occurs at x = 1.03125
77 *
78 */
79
80 #include "libm.h"
81
82 extern const long double _TBL_atanl_hi[], _TBL_atanl_lo[];
83 static const long double
84 one = 1.0L,
85 p1 = -3.333333333333333333333333333331344526118e-0001L,
86 p2 = 1.999999999999999999999999989931277668570e-0001L,
87 p3 = -1.428571428571428571428553606221309530901e-0001L,
88 p4 = 1.111111111111111111095219842737139747418e-0001L,
89 p5 = -9.090909090909090825503603835248061123323e-0002L,
90 p6 = 7.692307692307664052130743214708925258904e-0002L,
91 p7 = -6.666666666660213835187713228363717388266e-0002L,
92 p8 = 5.882352940152439399097283359608661949504e-0002L,
93 p9 = -5.263157780447533993046614040509529668487e-0002L,
94 p10 = 4.761895816878184933175855990886788439447e-0002L,
95 p11 = -4.347345005832274022681019724553538135922e-0002L,
96 p12 = 3.983031914579635037502589204647752042736e-0002L,
97 p13 = -3.348206704469830575196657749413894897554e-0002L,
98 q1 = -3.333333333333333333333333333195273650186e-0001L,
99 q2 = 1.999999999999999999999988146114392615808e-0001L,
100 q3 = -1.428571428571428571057630319435467111253e-0001L,
101 q4 = 1.111111111111105373263048208994541544098e-0001L,
102 q5 = -9.090909090421834209167373258681021816441e-0002L,
103 q6 = 7.692305377813692706850171767150701644539e-0002L,
104 q7 = -6.660896644393861499914731734305717901330e-0002L,
105 pio2hi = 1.570796326794896619231321691639751398740e+0000L,
106 pio2lo = 4.335905065061890512398522013021675984381e-0035L;
107
108 #define i0 0
109 #define i1 3
110
111 long double
112 atanl(long double x) {
113 long double y, z, r, p, s;
114 int *px = (int *) &x, *py = (int *) &y;
115 int ix, iy, sign, j;
116
117 ix = px[i0];
118 sign = ix & 0x80000000;
119 ix ^= sign;
120
121 /* for |x| < 1/8 */
122 if (ix < 0x3ffc0000) {
123 if (ix < 0x3feb0000) { /* when |x| < 2**(-prec/6-2) */
124 if (ix < 0x3fc50000) { /* if |x| < 2**(-prec/2-2) */
125 s = one;
126 *(3 - i0 + (int *) &s) = -1; /* s = 1-ulp */
127 *(1 + (int *) &s) = -1;
128 *(2 + (int *) &s) = -1;
129 *(i0 + (int *) &s) -= 1;
130 if ((int) (s * x) < 1)
131 return (x); /* raise inexact */
132 }
133 z = x * x;
134 if (ix < 0x3fe20000) { /* if |x| < 2**(-prec/4-1) */
135 return (x + (x * z) * p1);
136 } else { /* if |x| < 2**(-prec/6-2) */
137 return (x + (x * z) * (p1 + z * p2));
138 }
139 }
140 z = x * x;
141 return (x + (x * z) * (p1 + z * (p2 + z * (p3 + z * (p4 +
142 z * (p5 + z * (p6 + z * (p7 + z * (p8 + z * (p9 +
143 z * (p10 + z * (p11 + z * (p12 + z * p13)))))))))))));
144 }
145
146 /* for |x| >= 8.0 */
147 if (ix >= 0x40020000) {
148 px[i0] = ix;
149 if (ix < 0x40050400) { /* x < 65 */
150 r = one / x;
151 z = r * r;
152 /*
153 * poly1
154 */
155 y = r * (one + z * (p1 + z * (p2 + z * (p3 +
156 z * (p4 + z * (p5 + z * (p6 + z * (p7 +
157 z * (p8 + z * (p9 + z * (p10 + z * (p11 +
158 z * (p12 + z * p13)))))))))))));
159 y -= pio2lo;
160 } else if (ix < 0x40260000) { /* x < 2**(prec/3+2) */
161 r = one / x;
162 z = r * r;
163 /*
164 * poly2
165 */
166 y = r * (one + z * (q1 + z * (q2 + z * (q3 + z * (q4 +
167 z * (q5 + z * (q6 + z * q7)))))));
168 y -= pio2lo;
169 } else if (ix < 0x40720000) { /* x < 2**(prec+2) */
170 y = one / x - pio2lo;
171 } else if (ix < 0x7fff0000) { /* x < inf */
172 y = -pio2lo;
173 } else { /* x is inf or NaN */
174 if (((ix - 0x7fff0000) | px[1] | px[2] | px[i1]) != 0)
175 return (x - x);
176 y = -pio2lo;
177 }
178
179 if (sign == 0)
180 return (pio2hi - y);
181 else
182 return (y - pio2hi);
183 }
184
185 /* now x is between 1/8 and 8 */
186 px[i0] = ix;
187 iy = (ix + 0x00000800) & 0x7ffff000;
188 py[i0] = iy;
189 py[1] = py[2] = py[i1] = 0;
190 j = (iy - 0x3ffc0000) >> 12;
191
192 if (sign == 0)
193 s = (x - y) / (one + x * y);
194 else
195 s = (y - x) / (one + x * y);
196 z = s * s;
197 if (ix == iy)
198 p = s * (one + z * (q1 + z * (q2 + z * (q3 + z * q4))));
199 else
200 p = s * (one + z * (q1 + z * (q2 + z * (q3 + z * (q4 +
201 z * (q5 + z * (q6 + z * q7)))))));
202 if (sign == 0) {
203 r = p + _TBL_atanl_lo[j];
204 return (r + _TBL_atanl_hi[j]);
205 } else {
206 r = p - _TBL_atanl_lo[j];
207 return (r - _TBL_atanl_hi[j]);
208 }
209 }
|
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 __atanl = atanl
32
33 /*
34 * atanl(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
64 * (3.4) Otherwise, atan(x) = atan(inf) - poly1(1/x)
65 *
66 * (4). Now x is in (0.125, 8)
67 * Find y that match x to 4.5 bit after binary (easy).
68 * If iy is the high word of y, then
69 * single : j = (iy - 0x3e000000) >> 19
70 * double : j = (iy - 0x3fc00000) >> 16
71 * quad : j = (iy - 0x3ffc0000) >> 12
72 *
73 * Let s = (x-y)/(1+x*y). Then
74 * atan(x) = atan(y) + poly1(s)
75 * = _TBL_atanl_hi[j] + (_TBL_atanl_lo[j] + poly2(s) )
76 *
77 * Note. |s| <= 1.5384615385e-02 = 1/65. Maxium occurs at x = 1.03125
78 *
79 */
80
81 #include "libm.h"
82
83 extern const long double _TBL_atanl_hi[], _TBL_atanl_lo[];
84 static const long double one = 1.0L,
85 p1 = -3.333333333333333333333333333331344526118e-0001L,
86 p2 = 1.999999999999999999999999989931277668570e-0001L,
87 p3 = -1.428571428571428571428553606221309530901e-0001L,
88 p4 = 1.111111111111111111095219842737139747418e-0001L,
89 p5 = -9.090909090909090825503603835248061123323e-0002L,
90 p6 = 7.692307692307664052130743214708925258904e-0002L,
91 p7 = -6.666666666660213835187713228363717388266e-0002L,
92 p8 = 5.882352940152439399097283359608661949504e-0002L,
93 p9 = -5.263157780447533993046614040509529668487e-0002L,
94 p10 = 4.761895816878184933175855990886788439447e-0002L,
95 p11 = -4.347345005832274022681019724553538135922e-0002L,
96 p12 = 3.983031914579635037502589204647752042736e-0002L,
97 p13 = -3.348206704469830575196657749413894897554e-0002L,
98 q1 = -3.333333333333333333333333333195273650186e-0001L,
99 q2 = 1.999999999999999999999988146114392615808e-0001L,
100 q3 = -1.428571428571428571057630319435467111253e-0001L,
101 q4 = 1.111111111111105373263048208994541544098e-0001L,
102 q5 = -9.090909090421834209167373258681021816441e-0002L,
103 q6 = 7.692305377813692706850171767150701644539e-0002L,
104 q7 = -6.660896644393861499914731734305717901330e-0002L,
105 pio2hi = 1.570796326794896619231321691639751398740e+0000L,
106 pio2lo = 4.335905065061890512398522013021675984381e-0035L;
107
108 #define i0 0
109 #define i1 3
110
111 long double
112 atanl(long double x)
113 {
114 long double y, z, r, p, s;
115 int *px = (int *)&x, *py = (int *)&y;
116 int ix, iy, sign, j;
117
118 ix = px[i0];
119 sign = ix & 0x80000000;
120 ix ^= sign;
121
122 /* for |x| < 1/8 */
123 if (ix < 0x3ffc0000) {
124 if (ix < 0x3feb0000) { /* when |x| < 2**(-prec/6-2) */
125 if (ix < 0x3fc50000) { /* if |x| < 2**(-prec/2-2) */
126 s = one;
127 *(3 - i0 + (int *)&s) = -1; /* s = 1-ulp */
128 *(1 + (int *)&s) = -1;
129 *(2 + (int *)&s) = -1;
130 *(i0 + (int *)&s) -= 1;
131
132 if ((int)(s * x) < 1)
133 return (x); /* raise inexact */
134 }
135
136 z = x * x;
137
138 if (ix < 0x3fe20000) /* if |x| < 2**(-prec/4-1) */
139 return (x + (x * z) * p1);
140 else /* if |x| < 2**(-prec/6-2) */
141 return (x + (x * z) * (p1 + z * p2));
142 }
143
144 z = x * x;
145 return (x + (x * z) * (p1 + z * (p2 + z * (p3 + z * (p4 + z *
146 (p5 + z * (p6 + z * (p7 + z * (p8 + z * (p9 + z *
147 (p10 + z * (p11 + z * (p12 + z * p13)))))))))))));
148 }
149
150 /* for |x| >= 8.0 */
151 if (ix >= 0x40020000) {
152 px[i0] = ix;
153
154 if (ix < 0x40050400) { /* x < 65 */
155 r = one / x;
156 z = r * r;
157
158 /*
159 * poly1
160 */
161 y = r * (one + z * (p1 + z * (p2 + z * (p3 + z * (p4 +
162 z * (p5 + z * (p6 + z * (p7 + z * (p8 + z * (p9 +
163 z * (p10 + z * (p11 + z * (p12 + z *
164 p13)))))))))))));
165 y -= pio2lo;
166 } else if (ix < 0x40260000) { /* x < 2**(prec/3+2) */
167 r = one / x;
168 z = r * r;
169
170 /*
171 * poly2
172 */
173 y = r * (one + z * (q1 + z * (q2 + z * (q3 + z * (q4 +
174 z * (q5 + z * (q6 + z * q7)))))));
175 y -= pio2lo;
176 } else if (ix < 0x40720000) { /* x < 2**(prec+2) */
177 y = one / x - pio2lo;
178 } else if (ix < 0x7fff0000) { /* x < inf */
179 y = -pio2lo;
180 } else { /* x is inf or NaN */
181 if (((ix - 0x7fff0000) | px[1] | px[2] | px[i1]) != 0)
182 return (x - x);
183
184 y = -pio2lo;
185 }
186
187 if (sign == 0)
188 return (pio2hi - y);
189 else
190 return (y - pio2hi);
191 }
192
193 /* now x is between 1/8 and 8 */
194 px[i0] = ix;
195 iy = (ix + 0x00000800) & 0x7ffff000;
196 py[i0] = iy;
197 py[1] = py[2] = py[i1] = 0;
198 j = (iy - 0x3ffc0000) >> 12;
199
200 if (sign == 0)
201 s = (x - y) / (one + x * y);
202 else
203 s = (y - x) / (one + x * y);
204
205 z = s * s;
206
207 if (ix == iy)
208 p = s * (one + z * (q1 + z * (q2 + z * (q3 + z * q4))));
209 else
210 p = s * (one + z * (q1 + z * (q2 + z * (q3 + z * (q4 + z * (q5 +
211 z * (q6 + z * q7)))))));
212
213 if (sign == 0) {
214 r = p + _TBL_atanl_lo[j];
215 return (r + _TBL_atanl_hi[j]);
216 } else {
217 r = p - _TBL_atanl_lo[j];
218 return (r - _TBL_atanl_hi[j]);
219 }
220 }
|