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 __atan = atan
31
32 /* INDENT OFF */
33 /*
34 * atan(x)
35 * Accurate Table look-up algorithm with polynomial approximation in
36 * partially product form.
37 *
38 * -- K.C. Ng, October 17, 2004
39 *
40 * Algorithm
41 *
42 * (1). Purge off Inf and NaN and 0
43 * (2). Reduce x to positive by atan(x) = -atan(-x).
44 * (3). For x <= 1/8 and let z = x*x, return
45 * (2.1) if x < 2^(-prec/2), atan(x) = x with inexact flag raised
46 * (2.2) if x < 2^(-prec/4-1), atan(x) = x+(x/3)(x*x)
47 * (2.3) if x < 2^(-prec/6-2), atan(x) = x+(z-5/3)(z*x/5)
48 * (2.4) Otherwise
49 * atan(x) = poly1(x) = x + A * B,
50 * where
51 * A = (p1*x*z) * (p2+z(p3+z))
52 * B = (p4+z)+z*z) * (p5+z(p6+z))
62 * poly2(r) = (q1*r) * (q2+z(q3+z)) * (q4+z),
63 * its domain is [0, 0.0154]; and its remez absolute
64 * approximation error is bounded by
65 * |atan(x)-poly2(x)|<= 2^-59.45
66 *
67 * (5). Now x is in (0.125, 8).
68 * Recall identity
69 * atan(x) = atan(y) + atan((x-y)/(1+x*y)).
70 * Let j = (ix - 0x3fc00000) >> 16, 0 <= j < 96, where ix is the high
71 * part of x in IEEE double format. Then
72 * atan(x) = atan(y[j]) + poly2((x-y[j])/(1+x*y[j]))
73 * where y[j] are carefully chosen so that it matches x to around 4.5
74 * bits and at the same time atan(y[j]) is very close to an IEEE double
75 * floating point number. Calculation indicates that
76 * max|(x-y[j])/(1+x*y[j])| < 0.0154
77 * j,x
78 *
79 * Accuracy: Maximum error observed is bounded by 0.6 ulp after testing
80 * more than 10 million random arguments
81 */
82 /* INDENT ON */
83
84 #include "libm.h"
85 #include "libm_protos.h"
86
87 extern const double _TBL_atan[];
88 static const double g[] = {
89 /* one = */ 1.0,
90 /* p1 = */ 8.02176624254765935351230154992663301527500152588e-0002,
91 /* p2 = */ 1.27223421700559402580665846471674740314483642578e+0000,
92 /* p3 = */ -1.20606901800503640842521235754247754812240600586e+0000,
93 /* p4 = */ -2.36088967922325565496066701598465442657470703125e+0000,
94 /* p5 = */ 1.38345799501389166152875986881554126739501953125e+0000,
95 /* p6 = */ 1.06742368078953453469637224770849570631980895996e+0000,
96 /* q1 = */ -1.42796626333911796935538518482644576579332351685e-0001,
97 /* q2 = */ 3.51427110447873227059810477159863497078605962912e+0000,
98 /* q3 = */ 5.92129112708164262457444237952586263418197631836e-0001,
99 /* q4 = */ -1.99272234785683144409063061175402253866195678711e+0000,
100 /* pio2hi */ 1.570796326794896558e+00,
101 /* pio2lo */ 6.123233995736765886e-17,
102 /* t1 = */ -0.333333333333333333333333333333333,
103 /* t2 = */ 0.2,
104 /* t3 = */ -1.666666666666666666666666666666666,
105 };
106
107 #define one g[0]
108 #define p1 g[1]
109 #define p2 g[2]
110 #define p3 g[3]
111 #define p4 g[4]
112 #define p5 g[5]
113 #define p6 g[6]
114 #define q1 g[7]
115 #define q2 g[8]
116 #define q3 g[9]
117 #define q4 g[10]
118 #define pio2hi g[11]
119 #define pio2lo g[12]
120 #define t1 g[13]
121 #define t2 g[14]
122 #define t3 g[15]
123
124
125 double
126 atan(double x) {
127 double y, z, r, p, s;
128 int ix, lx, hx, j;
129
130 hx = ((int *) &x)[HIWORD];
131 lx = ((int *) &x)[LOWORD];
132 ix = hx & ~0x80000000;
133 j = ix >> 20;
134
135 /* for |x| < 1/8 */
136 if (j < 0x3fc) {
137 if (j < 0x3f5) { /* when |x| < 2**(-prec/6-2) */
138 if (j < 0x3e3) { /* if |x| < 2**(-prec/2-2) */
139 return ((int) x == 0 ? x : one);
140 }
141 if (j < 0x3f1) { /* if |x| < 2**(-prec/4-1) */
142 return (x + (x * t1) * (x * x));
143 } else { /* if |x| < 2**(-prec/6-2) */
144 z = x * x;
145 s = t2 * x;
146 return (x + (t3 + z) * (s * z));
147 }
148 }
149 z = x * x; s = p1 * x;
150 return (x + ((s * z) * (p2 + z * (p3 + z))) *
151 (((p4 + z) + z * z) * (p5 + z * (p6 + z))));
152 }
153
154 /* for |x| >= 8.0 */
155 if (j >= 0x402) {
156 if (j < 0x436) {
157 r = one / x;
158 if (hx >= 0) {
159 y = pio2hi; p = pio2lo;
160 } else {
161 y = -pio2hi; p = -pio2lo;
162 }
163 if (ix < 0x40504000) { /* x < 65 */
164 z = r * r;
165 s = p1 * r;
166 return (y + ((p - r) - ((s * z) *
167 (p2 + z * (p3 + z))) *
168 (((p4 + z) + z * z) *
169 (p5 + z * (p6 + z)))));
170 } else if (j < 0x412) {
171 z = r * r;
172 return (y + (p - ((q1 * r) * (q4 + z)) *
173 (q2 + z * (q3 + z))));
174 } else
175 return (y + (p - r));
176 } else {
177 if (j >= 0x7ff) /* x is inf or NaN */
178 if (((ix - 0x7ff00000) | lx) != 0)
179 #if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN)
180 return (ix >= 0x7ff80000 ? x : x - x);
181 /* assumes sparc-like QNaN */
182 #else
183 return (x - x);
184 #endif
185 y = -pio2lo;
186 return (hx >= 0 ? pio2hi - y : y - pio2hi);
187 }
188 } else { /* now x is between 1/8 and 8 */
189 double *w, w0, w1, s, z;
190 w = (double *) _TBL_atan + (((ix - 0x3fc00000) >> 16) << 1);
191 w0 = (hx >= 0)? w[0] : -w[0];
192 s = (x - w0) / (one + x * w0);
193 w1 = (hx >= 0)? w[1] : -w[1];
194 z = s * s;
195 return (((q1 * s) * (q4 + z)) * (q2 + z * (q3 + z)) + w1);
196 }
197 }
|
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 __atan = atan
32
33
34 /*
35 * atan(x)
36 * Accurate Table look-up algorithm with polynomial approximation in
37 * partially product form.
38 *
39 * -- K.C. Ng, October 17, 2004
40 *
41 * Algorithm
42 *
43 * (1). Purge off Inf and NaN and 0
44 * (2). Reduce x to positive by atan(x) = -atan(-x).
45 * (3). For x <= 1/8 and let z = x*x, return
46 * (2.1) if x < 2^(-prec/2), atan(x) = x with inexact flag raised
47 * (2.2) if x < 2^(-prec/4-1), atan(x) = x+(x/3)(x*x)
48 * (2.3) if x < 2^(-prec/6-2), atan(x) = x+(z-5/3)(z*x/5)
49 * (2.4) Otherwise
50 * atan(x) = poly1(x) = x + A * B,
51 * where
52 * A = (p1*x*z) * (p2+z(p3+z))
53 * B = (p4+z)+z*z) * (p5+z(p6+z))
63 * poly2(r) = (q1*r) * (q2+z(q3+z)) * (q4+z),
64 * its domain is [0, 0.0154]; and its remez absolute
65 * approximation error is bounded by
66 * |atan(x)-poly2(x)|<= 2^-59.45
67 *
68 * (5). Now x is in (0.125, 8).
69 * Recall identity
70 * atan(x) = atan(y) + atan((x-y)/(1+x*y)).
71 * Let j = (ix - 0x3fc00000) >> 16, 0 <= j < 96, where ix is the high
72 * part of x in IEEE double format. Then
73 * atan(x) = atan(y[j]) + poly2((x-y[j])/(1+x*y[j]))
74 * where y[j] are carefully chosen so that it matches x to around 4.5
75 * bits and at the same time atan(y[j]) is very close to an IEEE double
76 * floating point number. Calculation indicates that
77 * max|(x-y[j])/(1+x*y[j])| < 0.0154
78 * j,x
79 *
80 * Accuracy: Maximum error observed is bounded by 0.6 ulp after testing
81 * more than 10 million random arguments
82 */
83
84 #include "libm.h"
85 #include "libm_protos.h"
86
87 extern const double _TBL_atan[];
88
89 static const double g[] = {
90 /* one = */
91 1.0,
92 /* p1 = */8.02176624254765935351230154992663301527500152588e-0002,
93 /* p2 = */1.27223421700559402580665846471674740314483642578e+0000,
94 /* p3 = */-1.20606901800503640842521235754247754812240600586e+0000,
95 /* p4 = */-2.36088967922325565496066701598465442657470703125e+0000,
96 /* p5 = */1.38345799501389166152875986881554126739501953125e+0000,
97 /* p6 = */1.06742368078953453469637224770849570631980895996e+0000,
98 /* q1 = */ -1.42796626333911796935538518482644576579332351685e-0001,
99 /* q2 = */ 3.51427110447873227059810477159863497078605962912e+0000,
100 /* q3 = */ 5.92129112708164262457444237952586263418197631836e-0001,
101 /* q4 = */ -1.99272234785683144409063061175402253866195678711e+0000,
102 /* pio2hi */ 1.570796326794896558e+00,
103 /* pio2lo */ 6.123233995736765886e-17,
104 /* t1 = */ -0.333333333333333333333333333333333,
105 /* t2 = */ 0.2,
106 /* t3 = */ -1.666666666666666666666666666666666,
107 };
108
109 #define one g[0]
110 #define p1 g[1]
111 #define p2 g[2]
112 #define p3 g[3]
113 #define p4 g[4]
114 #define p5 g[5]
115 #define p6 g[6]
116 #define q1 g[7]
117 #define q2 g[8]
118 #define q3 g[9]
119 #define q4 g[10]
120 #define pio2hi g[11]
121 #define pio2lo g[12]
122 #define t1 g[13]
123 #define t2 g[14]
124 #define t3 g[15]
125
126 double
127 atan(double x)
128 {
129 double y, z, r, p, s;
130 int ix, lx, hx, j;
131
132 hx = ((int *)&x)[HIWORD];
133 lx = ((int *)&x)[LOWORD];
134 ix = hx & ~0x80000000;
135 j = ix >> 20;
136
137 /* for |x| < 1/8 */
138 if (j < 0x3fc) {
139 if (j < 0x3f5) { /* when |x| < 2**(-prec/6-2) */
140 if (j < 0x3e3) /* if |x| < 2**(-prec/2-2) */
141 return ((int)x == 0 ? x : one);
142
143 if (j < 0x3f1) { /* if |x| < 2**(-prec/4-1) */
144 return (x + (x * t1) * (x * x));
145 } else { /* if |x| < 2**(-prec/6-2) */
146 z = x * x;
147 s = t2 * x;
148 return (x + (t3 + z) * (s * z));
149 }
150 }
151
152 z = x * x;
153 s = p1 * x;
154 return (x + ((s * z) * (p2 + z * (p3 + z))) *
155 (((p4 + z) + z * z) * (p5 + z * (p6 + z))));
156 }
157
158 /* for |x| >= 8.0 */
159 if (j >= 0x402) {
160 if (j < 0x436) {
161 r = one / x;
162
163 if (hx >= 0) {
164 y = pio2hi;
165 p = pio2lo;
166 } else {
167 y = -pio2hi;
168 p = -pio2lo;
169 }
170
171 if (ix < 0x40504000) { /* x < 65 */
172 z = r * r;
173 s = p1 * r;
174 return (y + ((p - r) - ((s * z) * (p2 + z *
175 (p3 + z))) * (((p4 + z) + z * z) *
176 (p5 + z * (p6 + z)))));
177 } else if (j < 0x412) {
178 z = r * r;
179 return (y + (p - ((q1 * r) * (q4 + z)) * (q2 +
180 z * (q3 + z))));
181 } else {
182 return (y + (p - r));
183 }
184 } else {
185 if (j >= 0x7ff) /* x is inf or NaN */
186 if (((ix - 0x7ff00000) | lx) != 0)
187 #if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN)
188 return (ix >= 0x7ff80000 ? x : x - x);
189
190 /* assumes sparc-like QNaN */
191 #else
192 return (x - x);
193 #endif
194 y = -pio2lo;
195 return (hx >= 0 ? pio2hi - y : y - pio2hi);
196 }
197 } else { /* now x is between 1/8 and 8 */
198 double *w, w0, w1, s, z;
199
200 w = (double *)_TBL_atan + (((ix - 0x3fc00000) >> 16) << 1);
201 w0 = (hx >= 0) ? w[0] : -w[0];
202 s = (x - w0) / (one + x * w0);
203 w1 = (hx >= 0) ? w[1] : -w[1];
204 z = s * s;
205 return (((q1 * s) * (q4 + z)) * (q2 + z * (q3 + z)) + w1);
206 }
207 }
|