Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/complex/catan.c
+++ new/usr/src/lib/libm/common/complex/catan.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 __catan = catan
31 32
32 -/* INDENT OFF */
33 +
33 34 /*
34 35 * dcomplex catan(dcomplex z);
35 36 *
36 37 * If
37 38 * z = x + iy,
38 39 *
39 40 * then
40 41 * 1 ( 2x ) 1 2 2
41 42 * Re w = - arctan(-----------) = - ATAN2(2x, 1 - x - y )
42 43 * 2 ( 2 2) 2
43 44 * (1 - x - y )
44 45 *
45 46 * ( 2 2)
46 47 * 1 (x + (y+1) ) 1 4y
47 48 * Im w = - log(------------) .= --- log [ 1 + ------------- ]
48 49 * 4 ( 2 2) 4 2 2
49 50 * (x + (y-1) ) x + (y-1)
50 51 *
51 52 * 2 16 3 y
52 53 * = t - 2t + -- t - ..., where t = -----------------
53 54 * 3 x*x + (y-1)*(y-1)
54 55 *
55 56 * Note that: if catan( x, y) = ( u, v), then
56 57 * catan(-x, y) = (-u, v)
57 58 * catan( x,-y) = ( u,-v)
58 59 *
59 60 * Also, catan(x,y) = -i*catanh(-y,x), or
60 61 * catanh(x,y) = i*catan(-y,x)
61 62 * So, if catanh(y,x) = (v,u), then catan(x,y) = -i*(-v,u) = (u,v), i.e.,
62 63 * catan(x,y) = (u,v)
63 64 *
64 65 * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)):
65 66 * catan( 0 , 0 ) = (0 , 0 )
|
↓ open down ↓ |
23 lines elided |
↑ open up ↑ |
66 67 * catan( NaN, 0 ) = (NaN , 0 )
67 68 * catan( 0 , 1 ) = (0 , +inf) with divide-by-zero
68 69 * catan( inf, y ) = (pi/2 , 0 ) for finite +y
69 70 * catan( NaN, y ) = (NaN , NaN ) with invalid for finite y != 0
70 71 * catan( x , inf ) = (pi/2 , 0 ) for finite +x
71 72 * catan( inf, inf ) = (pi/2 , 0 )
72 73 * catan( NaN, inf ) = (NaN , 0 )
73 74 * catan( x , NaN ) = (NaN , NaN ) with invalid for finite x
74 75 * catan( inf, NaN ) = (pi/2 , +-0 )
75 76 */
76 -/* INDENT ON */
77 77
78 -#include "libm.h" /* atan/atan2/fabs/log/log1p */
78 +#include "libm.h" /* atan/atan2/fabs/log/log1p */
79 79 #include "complex_wrapper.h"
80 80
81 -/* INDENT OFF */
82 -static const double
83 - pi_2 = 1.570796326794896558e+00,
81 +static const double pi_2 = 1.570796326794896558e+00,
84 82 zero = 0.0,
85 83 half = 0.5,
86 84 two = 2.0,
87 85 ln2 = 6.931471805599453094172321214581765680755e-0001,
88 86 one = 1.0;
89 -/* INDENT ON */
87 +
90 88
91 89 dcomplex
92 -catan(dcomplex z) {
90 +catan(dcomplex z)
91 +{
93 92 dcomplex ans;
94 93 double x, y, ax, ay, t;
95 94 int hx, hy, ix, iy;
96 95 unsigned lx, ly;
97 96
98 97 x = D_RE(z);
99 98 y = D_IM(z);
100 99 ax = fabs(x);
101 100 ay = fabs(y);
102 101 hx = HI_WORD(x);
103 102 lx = LO_WORD(x);
104 103 hy = HI_WORD(y);
105 104 ly = LO_WORD(y);
|
↓ open down ↓ |
3 lines elided |
↑ open up ↑ |
106 105 ix = hx & 0x7fffffff;
107 106 iy = hy & 0x7fffffff;
108 107
109 108 /* x is inf or NaN */
110 109 if (ix >= 0x7ff00000) {
111 110 if (ISINF(ix, lx)) {
112 111 D_RE(ans) = pi_2;
113 112 D_IM(ans) = zero;
114 113 } else {
115 114 D_RE(ans) = x + x;
115 +
116 116 if ((iy | ly) == 0 || (ISINF(iy, ly)))
117 117 D_IM(ans) = zero;
118 118 else
119 119 D_IM(ans) = (fabs(y) - ay) / (fabs(y) - ay);
120 120 }
121 121 } else if (iy >= 0x7ff00000) {
122 122 /* y is inf or NaN */
123 123 if (ISINF(iy, ly)) {
124 124 D_RE(ans) = pi_2;
125 125 D_IM(ans) = zero;
126 126 } else {
127 127 D_RE(ans) = (fabs(x) - ax) / (fabs(x) - ax);
128 128 D_IM(ans) = y;
129 129 }
130 130 } else if ((ix | lx) == 0) {
131 - /* INDENT OFF */
131 + /* BEGIN CSTYLED */
132 132 /*
133 133 * x = 0
134 134 * 1 1
135 135 * A = --- * atan2(2x, 1-x*x-y*y) = --- atan2(0,1-|y|)
136 136 * 2 2
137 137 *
138 138 * 1 [ (y+1)*(y+1) ] 1 2 1 2y
139 139 * B = - log [ ------------ ] = - log (1+ ---) or - log(1+ ----)
140 140 * 4 [ (y-1)*(y-1) ] 2 y-1 2 1-y
141 141 */
142 - /* INDENT ON */
142 + /* END CSTYLED */
143 143 t = one - ay;
144 +
144 145 if (((iy - 0x3ff00000) | ly) == 0) {
145 146 /* y=1: catan(0,1)=(0,+inf) with 1/0 signal */
146 147 D_IM(ans) = ay / ax;
147 148 D_RE(ans) = zero;
148 149 } else if (iy >= 0x3ff00000) { /* y>1 */
149 150 D_IM(ans) = half * log1p(two / (-t));
150 151 D_RE(ans) = pi_2;
151 - } else { /* y<1 */
152 + } else { /* y<1 */
152 153 D_IM(ans) = half * log1p((ay + ay) / t);
153 154 D_RE(ans) = zero;
154 155 }
155 156 } else if (iy < 0x3e200000 || ((ix - iy) >> 20) >= 30) {
156 - /* INDENT OFF */
157 - /*
158 - * Tiny y (relative to 1+|x|)
159 - * |y| < E*(1+|x|)
160 - * where E=2**-29, -35, -60 for double, double extended, quad precision
161 - *
162 - * 1 [ x<=1: atan(x)
163 - * A = --- * atan2(2x, 1-x*x-y*y) ~ [ 1 1+x
164 - * 2 [ x>=1: - atan2(2,(1-x)*(-----))
165 - * 2 x
166 - *
167 - * y/x
168 - * B ~ t*(1-2t), where t = ----------------- is tiny
169 - * x + (y-1)*(y-1)/x
170 - */
171 - /* INDENT ON */
157 + /* BEGIN CSTYLED */
158 + /*
159 + * Tiny y (relative to 1+|x|)
160 + * |y| < E*(1+|x|)
161 + * where E=2**-29, -35, -60 for double, double extended, quad precision
162 + *
163 + * 1 [ x<=1: atan(x)
164 + * A = --- * atan2(2x, 1-x*x-y*y) ~ [ 1 1+x
165 + * 2 [ x>=1: - atan2(2,(1-x)*(-----))
166 + * 2 x
167 + *
168 + * y/x
169 + * B ~ t*(1-2t), where t = ----------------- is tiny
170 + * x + (y-1)*(y-1)/x
171 + */
172 + /* END CSTYLED */
172 173 if (ix < 0x3ff00000)
173 174 D_RE(ans) = atan(ax);
174 175 else
175 - D_RE(ans) = half * atan2(two, (one - ax) * (one +
176 - one / ax));
176 + D_RE(ans) = half * atan2(two, (one - ax) * (one + one /
177 + ax));
178 +
177 179 if ((iy | ly) == 0) {
178 180 D_IM(ans) = ay;
179 181 } else {
180 182 if (ix < 0x3e200000)
181 183 t = ay / ((ay - one) * (ay - one));
182 184 else if (ix > 0x41c00000)
183 185 t = (ay / ax) / ax;
184 186 else
185 187 t = ay / (ax * ax + (ay - one) * (ay - one));
188 +
186 189 D_IM(ans) = t * (one - (t + t));
187 190 }
188 191 } else if (iy >= 0x41c00000 && ((iy - ix) >> 20) >= 30) {
189 - /* INDENT OFF */
192 + /* BEGIN CSTYLED */
190 193 /*
191 194 * Huge y relative to 1+|x|
192 195 * |y| > Einv*(1+|x|), where Einv~2**(prec/2+3),
193 196 * 1
194 197 * A ~ --- * atan2(2x, -y*y) ~ pi/2
195 198 * 2
196 199 * y
197 200 * B ~ t*(1-2t), where t = --------------- is tiny
198 201 * (y-1)*(y-1)
199 202 */
200 - /* INDENT ON */
203 + /* END CSTYLED */
201 204 D_RE(ans) = pi_2;
202 205 t = (ay / (ay - one)) / (ay - one);
203 206 D_IM(ans) = t * (one - (t + t));
204 207 } else if (((iy - 0x3ff00000) | ly) == 0) {
205 - /* INDENT OFF */
208 + /* BEGIN CSTYLED */
206 209 /*
207 210 * y = 1
208 211 * 1 1
209 212 * A = --- * atan2(2x, -x*x) = --- atan2(2,-x)
210 213 * 2 2
211 214 *
212 215 * 1 [x*x + 4] 1 4 [ 0.5(log2-logx) if
213 216 * B = - log [-------] = - log (1+ ---) = [ |x|<E, else 0.25*
214 217 * 4 [ x*x ] 4 x*x [ log1p((2/x)*(2/x))
215 218 */
216 - /* INDENT ON */
219 + /* END CSTYLED */
217 220 D_RE(ans) = half * atan2(two, -ax);
218 - if (ix < 0x3e200000)
221 +
222 + if (ix < 0x3e200000) {
219 223 D_IM(ans) = half * (ln2 - log(ax));
220 - else {
224 + } else {
221 225 t = two / ax;
222 226 D_IM(ans) = 0.25 * log1p(t * t);
223 227 }
224 228 } else if (ix >= 0x43900000) {
225 - /* INDENT OFF */
229 + /* BEGIN CSTYLED */
226 230 /*
227 231 * Huge x:
228 232 * when |x| > 1/E^2,
229 233 * 1 pi
230 234 * A ~ --- * atan2(2x, -x*x-y*y) ~ ---
231 235 * 2 2
232 236 * y y/x
233 237 * B ~ t*(1-2t), where t = --------------- = (-------------- )/x
234 238 * x*x+(y-1)*(y-1) 1+((y-1)/x)^2
235 239 */
236 - /* INDENT ON */
240 + /* END CSTYLED */
237 241 D_RE(ans) = pi_2;
238 242 t = ((ay / ax) / (one + ((ay - one) / ax) * ((ay - one) /
239 - ax))) / ax;
243 + ax))) / ax;
240 244 D_IM(ans) = t * (one - (t + t));
241 245 } else if (ix < 0x38b00000) {
242 - /* INDENT OFF */
246 + /* BEGIN CSTYLED */
243 247 /*
244 248 * Tiny x:
245 249 * when |x| < E^4, (note that y != 1)
246 250 * 1 1
247 251 * A = --- * atan2(2x, 1-x*x-y*y) ~ --- * atan2(2x,(1-y)*(1+y))
248 252 * 2 2
249 253 *
250 254 * 1 [(y+1)*(y+1)] 1 2 1 2y
251 255 * B = - log [-----------] = - log (1+ ---) or - log(1+ ----)
252 256 * 4 [(y-1)*(y-1)] 2 y-1 2 1-y
253 257 */
254 - /* INDENT ON */
258 + /* END CSTYLED */
255 259 D_RE(ans) = half * atan2(ax + ax, (one - ay) * (one + ay));
260 +
256 261 if (iy >= 0x3ff00000)
257 262 D_IM(ans) = half * log1p(two / (ay - one));
258 263 else
259 264 D_IM(ans) = half * log1p((ay + ay) / (one - ay));
260 265 } else {
261 - /* INDENT OFF */
266 + /* BEGIN CSTYLED */
262 267 /*
263 268 * normal x,y
264 269 * 1
265 270 * A = --- * atan2(2x, 1-x*x-y*y)
266 271 * 2
267 272 *
268 273 * 1 [x*x+(y+1)*(y+1)] 1 4y
269 274 * B = - log [---------------] = - log (1+ -----------------)
270 275 * 4 [x*x+(y-1)*(y-1)] 4 x*x + (y-1)*(y-1)
271 276 */
272 - /* INDENT ON */
277 + /* END CSTYLED */
273 278 t = one - ay;
279 +
274 280 if (iy >= 0x3fe00000 && iy < 0x40000000) {
275 281 /* y close to 1 */
276 282 D_RE(ans) = half * (atan2((ax + ax), (t * (one + ay) -
277 - ax * ax)));
283 + ax * ax)));
278 284 } else if (ix >= 0x3fe00000 && ix < 0x40000000) {
279 285 /* x close to 1 */
280 - D_RE(ans) = half * atan2((ax + ax), ((one - ax) *
281 - (one + ax) - ay * ay));
282 - } else
286 + D_RE(ans) = half * atan2((ax + ax), ((one - ax) * (one +
287 + ax) - ay * ay));
288 + } else {
283 289 D_RE(ans) = half * atan2((ax + ax), ((one - ax * ax) -
284 - ay * ay));
290 + ay * ay));
291 + }
292 +
285 293 D_IM(ans) = 0.25 * log1p((4.0 * ay) / (ax * ax + t * t));
286 294 }
295 +
287 296 if (hx < 0)
288 297 D_RE(ans) = -D_RE(ans);
298 +
289 299 if (hy < 0)
290 300 D_IM(ans) = -D_IM(ans);
301 +
291 302 return (ans);
292 303 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX