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
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_synonyms.h"
86 #include "libm_protos.h"
87
88 extern const double _TBL_atan[];
89 static const double g[] = {
90 /* one = */ 1.0,
91 /* p1 = */ 8.02176624254765935351230154992663301527500152588e-0002,
92 /* p2 = */ 1.27223421700559402580665846471674740314483642578e+0000,
93 /* p3 = */ -1.20606901800503640842521235754247754812240600586e+0000,
94 /* p4 = */ -2.36088967922325565496066701598465442657470703125e+0000,
95 /* p5 = */ 1.38345799501389166152875986881554126739501953125e+0000,
96 /* p6 = */ 1.06742368078953453469637224770849570631980895996e+0000,
97 /* q1 = */ -1.42796626333911796935538518482644576579332351685e-0001,
98 /* q2 = */ 3.51427110447873227059810477159863497078605962912e+0000,
99 /* q3 = */ 5.92129112708164262457444237952586263418197631836e-0001,
100 /* q4 = */ -1.99272234785683144409063061175402253866195678711e+0000,
101 /* pio2hi */ 1.570796326794896558e+00,
102 /* pio2lo */ 6.123233995736765886e-17,
103 /* t1 = */ -0.333333333333333333333333333333333,
104 /* t2 = */ 0.2,
105 /* t3 = */ -1.666666666666666666666666666666666,
|
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
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,
|