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