Print this page
5261 libm should stop using synonyms.h
5298 fabs is 0-sized, confuses dis(1) and others
Reviewed by: Josef 'Jeff' Sipek <jeffpc@josefsipek.net>
Approved by: Gordon Ross <gwr@nexenta.com>
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/complex/catanl.c
+++ new/usr/src/lib/libm/common/complex/catanl.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.
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
|
↓ open down ↓ |
19 lines elided |
↑ open up ↑ |
20 20 */
21 21
22 22 /*
23 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 24 */
25 25 /*
26 26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 27 * Use is subject to license terms.
28 28 */
29 29
30 -#pragma weak catanl = __catanl
30 +#pragma weak __catanl = catanl
31 31
32 32 /* INDENT OFF */
33 33 /*
34 34 * ldcomplex catanl(ldcomplex z);
35 35 *
36 36 * Atan(z) return A + Bi where,
37 37 * 1
38 38 * A = --- * atan2(2x, 1-x*x-y*y)
39 39 * 2
40 40 *
41 41 * 1 [ x*x + (y+1)*(y+1) ] 1 4y
42 42 * B = --- log [ ----------------- ] = - log (1+ -----------------)
43 43 * 4 [ x*x + (y-1)*(y-1) ] 4 x*x + (y-1)*(y-1)
44 44 *
45 45 * 2 16 3 y
46 46 * = t - 2t + -- t - ..., where t = -----------------
47 47 * 3 x*x + (y-1)*(y-1)
48 48 * Proof:
49 49 * Let w = atan(z=x+yi) = A + B i. Then tan(w) = z.
50 50 * Since sin(w) = (exp(iw)-exp(-iw))/(2i), cos(w)=(exp(iw)+exp(-iw))/(2),
51 51 * Let p = exp(iw), then z = tan(w) = ((p-1/p)/(p+1/p))/i, or
52 52 * iz = (p*p-1)/(p*p+1), or, after simplification,
53 53 * p*p = (1+iz)/(1-iz) ... (1)
54 54 * LHS of (1) = exp(2iw) = exp(2i(A+Bi)) = exp(-2B)*exp(2iA)
55 55 * = exp(-2B)*(cos(2A)+i*sin(2A)) ... (2)
56 56 * 1-y+ix (1-y+ix)*(1+y+ix) 1-x*x-y*y + 2xi
57 57 * RHS of (1) = ------ = ----------------- = --------------- ... (3)
58 58 * 1+y-ix (1+y)**2 + x**2 (1+y)**2 + x**2
59 59 *
60 60 * Comparing the real and imaginary parts of (2) and (3), we have:
61 61 * cos(2A) : 1-x*x-y*y = sin(2A) : 2x
62 62 * and hence
63 63 * tan(2A) = 2x/(1-x*x-y*y), or
64 64 * A = 0.5 * atan2(2x, 1-x*x-y*y) ... (4)
65 65 *
66 66 * For the imaginary part B, Note that |p*p| = exp(-2B), and
67 67 * |1+iz| |i-z| hypot(x,(y-1))
68 68 * |----| = |---| = --------------
69 69 * |1-iz| |i+z| hypot(x,(y+1))
70 70 * Thus
71 71 * x*x + (y+1)*(y+1)
72 72 * exp(4B) = -----------------, or
73 73 * x*x + (y-1)*(y-1)
74 74 *
75 75 * 1 [x^2+(y+1)^2] 1 4y
76 76 * B = - log [-----------] = - log(1+ -------------) ... (5)
77 77 * 4 [x^2+(y-1)^2] 4 x^2+(y-1)^2
78 78 *
79 79 * QED.
80 80 *
81 81 * Note that: if catan( x, y) = ( u, v), then
82 82 * catan(-x, y) = (-u, v)
83 83 * catan( x,-y) = ( u,-v)
84 84 *
85 85 * Also, catan(x,y) = -i*catanh(-y,x), or
86 86 * catanh(x,y) = i*catan(-y,x)
87 87 * So, if catanh(y,x) = (v,u), then catan(x,y) = -i*(-v,u) = (u,v), i.e.,
88 88 * catan(x,y) = (u,v)
89 89 *
90 90 * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)):
91 91 * catan( 0 , 0 ) = (0 , 0 )
92 92 * catan( NaN, 0 ) = (NaN , 0 )
93 93 * catan( 0 , 1 ) = (0 , +inf) with divide-by-zero
94 94 * catan( inf, y ) = (pi/2 , 0 ) for finite +y
95 95 * catan( NaN, y ) = (NaN , NaN ) with invalid for finite y != 0
96 96 * catan( x , inf ) = (pi/2 , 0 ) for finite +x
97 97 * catan( inf, inf ) = (pi/2 , 0 )
98 98 * catan( NaN, inf ) = (NaN , 0 )
99 99 * catan( x , NaN ) = (NaN , NaN ) with invalid for finite x
100 100 * catan( inf, NaN ) = (pi/2 , +-0 )
101 101 */
102 102 /* INDENT ON */
103 103
104 104 #include "libm.h" /* atan2l/atanl/fabsl/isinfl/iszerol/log1pl/logl */
105 105 #include "complex_wrapper.h"
106 106 #include "longdouble.h"
107 107
108 108 /* INDENT OFF */
109 109 static const long double
110 110 zero = 0.0L,
111 111 one = 1.0L,
112 112 two = 2.0L,
113 113 half = 0.5L,
114 114 ln2 = 6.931471805599453094172321214581765680755e-0001L,
115 115 pi_2 = 1.570796326794896619231321691639751442098584699687552910487472L,
116 116 #if defined(__x86)
117 117 E = 2.910383045673370361328125000000000000000e-11L, /* 2**-35 */
118 118 Einv = 3.435973836800000000000000000000000000000e+10L; /* 2**+35 */
119 119 #else
120 120 E = 8.673617379884035472059622406959533691406e-19L, /* 2**-60 */
121 121 Einv = 1.152921504606846976000000000000000000000e18L; /* 2**+60 */
122 122 #endif
123 123 /* INDENT ON */
124 124
125 125 ldcomplex
126 126 catanl(ldcomplex z) {
127 127 ldcomplex ans;
128 128 long double x, y, t1, ax, ay, t;
129 129 int hx, hy, ix, iy;
130 130
131 131 x = LD_RE(z);
132 132 y = LD_IM(z);
133 133 ax = fabsl(x);
134 134 ay = fabsl(y);
135 135 hx = HI_XWORD(x);
136 136 hy = HI_XWORD(y);
137 137 ix = hx & 0x7fffffff;
138 138 iy = hy & 0x7fffffff;
139 139
140 140 /* x is inf or NaN */
141 141 if (ix >= 0x7fff0000) {
142 142 if (isinfl(x)) {
143 143 LD_RE(ans) = pi_2;
144 144 LD_IM(ans) = zero;
145 145 } else {
146 146 LD_RE(ans) = x + x;
147 147 if (iszerol(y) || (isinfl(y)))
148 148 LD_IM(ans) = zero;
149 149 else
150 150 LD_IM(ans) = (fabsl(y) - ay) / (fabsl(y) - ay);
151 151 }
152 152 } else if (iy >= 0x7fff0000) {
153 153 /* y is inf or NaN */
154 154 if (isinfl(y)) {
155 155 LD_RE(ans) = pi_2;
156 156 LD_IM(ans) = zero;
157 157 } else {
158 158 LD_RE(ans) = (fabsl(x) - ax) / (fabsl(x) - ax);
159 159 LD_IM(ans) = y;
160 160 }
161 161 } else if (iszerol(x)) {
162 162 /* INDENT OFF */
163 163 /*
164 164 * x = 0
165 165 * 1 1
166 166 * A = --- * atan2(2x, 1-x*x-y*y) = --- atan2(0,1-|y|)
167 167 * 2 2
168 168 *
169 169 * 1 [ (y+1)*(y+1) ] 1 2 1 2y
170 170 * B = - log [ ----------- ] = - log (1+ ---) or - log(1+ ----)
171 171 * 4 [ (y-1)*(y-1) ] 2 y-1 2 1-y
172 172 */
173 173 /* INDENT ON */
174 174 t = one - ay;
175 175 if (ay == one) {
176 176 /* y=1: catan(0,1)=(0,+inf) with 1/0 signal */
177 177 LD_IM(ans) = ay / ax;
178 178 LD_RE(ans) = zero;
179 179 } else if (ay > one) { /* y>1 */
180 180 LD_IM(ans) = half * log1pl(two / (-t));
181 181 LD_RE(ans) = pi_2;
182 182 } else { /* y<1 */
183 183 LD_IM(ans) = half * log1pl((ay + ay) / t);
184 184 LD_RE(ans) = zero;
185 185 }
186 186 } else if (ay < E * (one + ax)) {
187 187 /* INDENT OFF */
188 188 /*
189 189 * Tiny y (relative to 1+|x|)
190 190 * |y| < E*(1+|x|)
191 191 * where E=2**-29, -35, -60 for double, extended, quad precision
192 192 *
193 193 * 1 [x<=1: atan(x)
194 194 * A = - * atan2(2x,1-x*x-y*y) ~ [ 1 1+x
195 195 * 2 [x>=1: - atan2(2,(1-x)*(-----))
196 196 * 2 x
197 197 *
198 198 * y/x
199 199 * B ~ t*(1-2t), where t = ----------------- is tiny
200 200 * x + (y-1)*(y-1)/x
201 201 *
202 202 * y
203 203 * (when x < 2**-60, t = ----------- )
204 204 * (y-1)*(y-1)
205 205 */
206 206 /* INDENT ON */
207 207 if (ay == zero)
208 208 LD_IM(ans) = ay;
209 209 else {
210 210 t1 = ay - one;
211 211 if (ix < 0x3fc30000)
212 212 t = ay / (t1 * t1);
213 213 else if (ix > 0x403b0000)
214 214 t = (ay / ax) / ax;
215 215 else
216 216 t = ay / (ax * ax + t1 * t1);
217 217 LD_IM(ans) = t * (one - two * t);
218 218 }
219 219 if (ix < 0x3fff0000)
220 220 LD_RE(ans) = atanl(ax);
221 221 else
222 222 LD_RE(ans) = half * atan2l(two, (one - ax) * (one +
223 223 one / ax));
224 224
225 225 } else if (ay > Einv * (one + ax)) {
226 226 /* INDENT OFF */
227 227 /*
228 228 * Huge y relative to 1+|x|
229 229 * |y| > Einv*(1+|x|), where Einv~2**(prec/2+3),
230 230 * 1
231 231 * A ~ --- * atan2(2x, -y*y) ~ pi/2
232 232 * 2
233 233 * y
234 234 * B ~ t*(1-2t), where t = --------------- is tiny
235 235 * (y-1)*(y-1)
236 236 */
237 237 /* INDENT ON */
238 238 LD_RE(ans) = pi_2;
239 239 t = (ay / (ay - one)) / (ay - one);
240 240 LD_IM(ans) = t * (one - (t + t));
241 241 } else if (ay == one) {
242 242 /* INDENT OFF */
243 243 /*
244 244 * y=1
245 245 * 1 1
246 246 * A = - * atan2(2x, -x*x) = --- atan2(2,-x)
247 247 * 2 2
248 248 *
249 249 * 1 [ x*x+4] 1 4 [ 0.5(log2-logx) if
250 250 * B = - log [ -----] = - log (1+ ---) = [ |x|<E, else 0.25*
251 251 * 4 [ x*x ] 4 x*x [ log1p((2/x)*(2/x))
252 252 */
253 253 /* INDENT ON */
254 254 LD_RE(ans) = half * atan2l(two, -ax);
255 255 if (ax < E)
256 256 LD_IM(ans) = half * (ln2 - logl(ax));
257 257 else {
258 258 t = two / ax;
259 259 LD_IM(ans) = 0.25L * log1pl(t * t);
260 260 }
261 261 } else if (ax > Einv * Einv) {
262 262 /* INDENT OFF */
263 263 /*
264 264 * Huge x:
265 265 * when |x| > 1/E^2,
266 266 * 1 pi
267 267 * A ~ --- * atan2(2x, -x*x-y*y) ~ ---
268 268 * 2 2
269 269 * y y/x
270 270 * B ~ t*(1-2t), where t = --------------- = (-------------- )/x
271 271 * x*x+(y-1)*(y-1) 1+((y-1)/x)^2
272 272 */
273 273 /* INDENT ON */
274 274 LD_RE(ans) = pi_2;
275 275 t = ((ay / ax) / (one + ((ay - one) / ax) * ((ay - one) /
276 276 ax))) / ax;
277 277 LD_IM(ans) = t * (one - (t + t));
278 278 } else if (ax < E * E * E * E) {
279 279 /* INDENT OFF */
280 280 /*
281 281 * Tiny x:
282 282 * when |x| < E^4, (note that y != 1)
283 283 * 1 1
284 284 * A = --- * atan2(2x, 1-x*x-y*y) ~ --- * atan2(2x,1-y*y)
285 285 * 2 2
286 286 *
287 287 * 1 [ (y+1)*(y+1) ] 1 2 1 2y
288 288 * B = - log [ ----------- ] = - log (1+ ---) or - log(1+ ----)
289 289 * 4 [ (y-1)*(y-1) ] 2 y-1 2 1-y
290 290 */
291 291 /* INDENT ON */
292 292 LD_RE(ans) = half * atan2l(ax + ax, (one - ay) * (one + ay));
293 293 if (ay > one) /* y>1 */
294 294 LD_IM(ans) = half * log1pl(two / (ay - one));
295 295 else /* y<1 */
296 296 LD_IM(ans) = half * log1pl((ay + ay) / (one - ay));
297 297 } else {
298 298 /* INDENT OFF */
299 299 /*
300 300 * normal x,y
301 301 * 1
302 302 * A = --- * atan2(2x, 1-x*x-y*y)
303 303 * 2
304 304 *
305 305 * 1 [ x*x+(y+1)*(y+1) ] 1 4y
306 306 * B = - log [ --------------- ] = - log (1+ -----------------)
307 307 * 4 [ x*x+(y-1)*(y-1) ] 4 x*x + (y-1)*(y-1)
308 308 */
309 309 /* INDENT ON */
310 310 t = one - ay;
311 311 if (iy >= 0x3ffe0000 && iy < 0x40000000) {
312 312 /* y close to 1 */
313 313 LD_RE(ans) = half * (atan2l((ax + ax), (t * (one +
314 314 ay) - ax * ax)));
315 315 } else if (ix >= 0x3ffe0000 && ix < 0x40000000) {
316 316 /* x close to 1 */
317 317 LD_RE(ans) = half * atan2l((ax + ax), ((one - ax) *
318 318 (one + ax) - ay * ay));
319 319 } else
320 320 LD_RE(ans) = half * atan2l((ax + ax), ((one - ax *
321 321 ax) - ay * ay));
322 322 LD_IM(ans) = 0.25L * log1pl((4.0L * ay) / (ax * ax + t * t));
323 323 }
324 324 if (hx < 0)
325 325 LD_RE(ans) = -LD_RE(ans);
326 326 if (hy < 0)
327 327 LD_IM(ans) = -LD_IM(ans);
328 328 return (ans);
329 329 }
|
↓ open down ↓ |
289 lines elided |
↑ open up ↑ |
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX