Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/R/atanf.c
+++ new/usr/src/lib/libm/common/R/atanf.c
1 1 /*
2 2 * CDDL HEADER START
3 3 *
4 4 * The contents of this file are subject to the terms of the
5 5 * Common Development and Distribution License (the "License").
6 6 * You may not use this file except in compliance with the License.
7 7 *
8 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 9 * or http://www.opensolaris.org/os/licensing.
10 10 * See the License for the specific language governing permissions
11 11 * and limitations under the License.
12 12 *
13 13 * When distributing Covered Code, include this CDDL HEADER in each
14 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
|
↓ open down ↓ |
14 lines elided |
↑ open up ↑ |
15 15 * If applicable, add the following below this CDDL HEADER, with the
16 16 * fields enclosed by brackets "[]" replaced with your own identifying
17 17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 18 *
19 19 * CDDL HEADER END
20 20 */
21 21
22 22 /*
23 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 24 */
25 +
25 26 /*
26 27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 28 * Use is subject to license terms.
28 29 */
29 30
30 31 #pragma weak __atanf = atanf
31 32
32 -/* INDENT OFF */
33 +
33 34 /*
34 35 * float atanf(float x);
35 36 * Table look-up algorithm
36 37 * By K.C. Ng, March 9, 1989
37 38 *
38 39 * Algorithm.
39 40 *
40 41 * The algorithm is based on atan(x)=atan(y)+atan((x-y)/(1+x*y)).
41 42 * We use poly1(x) to approximate atan(x) for x in [0,1/8] with
42 43 * 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
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
46 47 * and use poly2(x) to approximate atan(x) for x in [0,1/65] with
47 48 * error (absolute)
48 49 * |atan(x)-poly2(x)|<= 2^-122.15 long double
49 50 * |atan(x)-poly2(x)|<= 2^-64.79 double
50 51 * |atan(x)-poly2(x)|<= 2^-35.36 float
51 52 * and use poly3(x) to approximate atan(x) for x in [1/8,7/16] with
52 53 * error (relative, on for single precision)
53 - * |(atan(x)-poly1(x))/x|<= 2^-25.53 float
54 + * |(atan(x)-poly1(x))/x|<= 2^-25.53 float
54 55 *
55 56 * Here poly1-3 are odd polynomial with the following form:
56 57 * x + x^3*(a1+x^2*(a2+...))
57 58 *
58 59 * (0). Purge off Inf and NaN and 0
59 60 * (1). Reduce x to positive by atan(x) = -atan(-x).
60 61 * (2). For x <= 1/8, use
61 62 * (2.1) if x < 2^(-prec/2-2), atan(x) = x with inexact
62 63 * (2.2) Otherwise
63 64 * atan(x) = poly1(x)
64 65 * (3). For x >= 8 then
65 66 * (3.1) if x >= 2^(prec+2), atan(x) = atan(inf) - pio2lo
66 67 * (3.2) if x >= 2^(prec/3+2), atan(x) = atan(inf) - 1/x
67 68 * (3.3) if x > 65, atan(x) = atan(inf) - poly2(1/x)
68 69 * (3.4) Otherwise, atan(x) = atan(inf) - poly1(1/x)
69 70 *
70 71 * (4). Now x is in (0.125, 8)
71 72 * Find y that match x to 4.5 bit after binary (easy).
72 73 * If iy is the high word of y, then
73 74 * single : j = (iy - 0x3e000000) >> 19
74 75 * (single is modified to (iy-0x3f000000)>>19)
75 76 * double : j = (iy - 0x3fc00000) >> 16
76 77 * quad : j = (iy - 0x3ffc0000) >> 12
77 78 *
78 79 * Let s = (x-y)/(1+x*y). Then
|
↓ open down ↓ |
15 lines elided |
↑ open up ↑ |
79 80 * atan(x) = atan(y) + poly1(s)
80 81 * = _TBL_r_atan_hi[j] + (_TBL_r_atan_lo[j] + poly2(s) )
81 82 *
82 83 * Note. |s| <= 1.5384615385e-02 = 1/65. Maxium occurs at x = 1.03125
83 84 *
84 85 */
85 86
86 87 #include "libm.h"
87 88
88 89 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 */
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 103
104 104 float
105 -atanf(float xx) {
105 +atanf(float xx)
106 +{
106 107 float x, y, z, r, p, s;
107 108 volatile double dummy __unused;
108 109 int ix, iy, sign, j;
109 110
110 111 x = xx;
111 - ix = *(int *) &x;
112 + ix = *(int *)&x;
112 113 sign = ix & 0x80000000;
113 114 ix ^= sign;
114 115
115 116 /* for |x| < 1/8 */
116 117 if (ix < 0x3e000000) {
117 - if (ix < 0x38800000) { /* if |x| < 2**(-prec/2-2) */
118 + if (ix < 0x38800000) { /* if |x| < 2**(-prec/2-2) */
118 119 dummy = big + x; /* get inexact flag if x != 0 */
119 120 #ifdef lint
120 121 dummy = dummy;
121 122 #endif
122 123 return (x);
123 124 }
125 +
124 126 z = x * x;
127 +
125 128 if (ix < 0x3c000000) { /* if |x| < 2**(-prec/4-1) */
126 129 x = x + (x * z) * p1;
127 130 return (x);
128 131 } else {
129 132 x = x + (x * z) * (p1 + z * p2);
130 133 return (x);
131 134 }
132 135 }
133 136
134 137 /* for |x| >= 8.0 */
135 138 if (ix >= 0x41000000) {
136 - *(int *) &x = ix;
137 - if (ix < 0x42820000) { /* x < 65 */
139 + *(int *)&x = ix;
140 +
141 + if (ix < 0x42820000) { /* x < 65 */
138 142 r = one / x;
139 143 z = r * r;
140 144 y = r * (one + z * (p1 + z * p2)); /* poly1 */
141 145 y -= pio2lo;
142 - } else if (ix < 0x44800000) { /* x < 2**(prec/3+2) */
146 + } else if (ix < 0x44800000) { /* x < 2**(prec/3+2) */
143 147 r = one / x;
144 148 z = r * r;
145 149 y = r * (one + z * q1); /* poly2 */
146 150 y -= pio2lo;
147 - } else if (ix < 0x4c800000) { /* x < 2**(prec+2) */
151 + } else if (ix < 0x4c800000) { /* x < 2**(prec+2) */
148 152 y = one / x - pio2lo;
149 - } else if (ix < 0x7f800000) { /* x < inf */
153 + } else if (ix < 0x7f800000) { /* x < inf */
150 154 y = -pio2lo;
151 - } else { /* x is inf or NaN */
152 - if (ix > 0x7f800000) {
155 + } else { /* x is inf or NaN */
156 + if (ix > 0x7f800000)
153 157 return (x * x); /* - -> * for Cheetah */
154 - }
158 +
155 159 y = -pio2lo;
156 160 }
157 161
158 162 if (sign == 0)
159 163 x = pio2hi - y;
160 164 else
161 165 x = y - pio2hi;
166 +
162 167 return (x);
163 168 }
164 169
165 -
166 170 /* now x is between 1/8 and 8 */
167 - if (ix < 0x3f000000) { /* between 1/8 and 1/2 */
171 + if (ix < 0x3f000000) { /* between 1/8 and 1/2 */
168 172 z = x * x;
169 - x = x + (x * z) * (a1 + z * (a2 + z * (a3 + z * (a4 +
170 - z * a5))));
173 + x = x + (x * z) * (a1 + z * (a2 + z * (a3 + z * (a4 + z *
174 + a5))));
171 175 return (x);
172 176 }
173 - *(int *) &x = ix;
177 +
178 + *(int *)&x = ix;
174 179 iy = (ix + 0x00040000) & 0x7ff80000;
175 - *(int *) &y = iy;
180 + *(int *)&y = iy;
176 181 j = (iy - 0x3f000000) >> 19;
177 182
178 - if (ix == iy)
179 - p = x - y; /* p=0.0 */
180 - else {
183 + if (ix == iy) {
184 + p = x - y; /* p=0.0 */
185 + } else {
181 186 if (sign == 0)
182 187 s = (x - y) / (one + x * y);
183 188 else
184 189 s = (y - x) / (one + x * y);
190 +
185 191 z = s * s;
186 192 p = s * (one + z * q1);
187 193 }
194 +
188 195 if (sign == 0) {
189 196 r = p + _TBL_r_atan_lo[j];
190 197 x = r + _TBL_r_atan_hi[j];
191 198 } else {
192 199 r = p - _TBL_r_atan_lo[j];
193 200 x = r - _TBL_r_atan_hi[j];
194 201 }
202 +
195 203 return (x);
196 204 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX