1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
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 /*
27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 #pragma weak __catan = catan
32
33
34 /*
35 * dcomplex catan(dcomplex z);
36 *
37 * If
38 * z = x + iy,
39 *
40 * then
41 * 1 ( 2x ) 1 2 2
42 * Re w = - arctan(-----------) = - ATAN2(2x, 1 - x - y )
43 * 2 ( 2 2) 2
44 * (1 - x - y )
45 *
46 * ( 2 2)
47 * 1 (x + (y+1) ) 1 4y
48 * Im w = - log(------------) .= --- log [ 1 + ------------- ]
49 * 4 ( 2 2) 4 2 2
50 * (x + (y-1) ) x + (y-1)
51 *
52 * 2 16 3 y
53 * = t - 2t + -- t - ..., where t = -----------------
54 * 3 x*x + (y-1)*(y-1)
55 *
56 * Note that: if catan( x, y) = ( u, v), then
57 * catan(-x, y) = (-u, v)
58 * catan( x,-y) = ( u,-v)
59 *
60 * Also, catan(x,y) = -i*catanh(-y,x), or
61 * catanh(x,y) = i*catan(-y,x)
62 * So, if catanh(y,x) = (v,u), then catan(x,y) = -i*(-v,u) = (u,v), i.e.,
63 * catan(x,y) = (u,v)
64 *
65 * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)):
66 * catan( 0 , 0 ) = (0 , 0 )
67 * catan( NaN, 0 ) = (NaN , 0 )
68 * catan( 0 , 1 ) = (0 , +inf) with divide-by-zero
69 * catan( inf, y ) = (pi/2 , 0 ) for finite +y
70 * catan( NaN, y ) = (NaN , NaN ) with invalid for finite y != 0
71 * catan( x , inf ) = (pi/2 , 0 ) for finite +x
72 * catan( inf, inf ) = (pi/2 , 0 )
73 * catan( NaN, inf ) = (NaN , 0 )
74 * catan( x , NaN ) = (NaN , NaN ) with invalid for finite x
75 * catan( inf, NaN ) = (pi/2 , +-0 )
76 */
77
78 #include "libm.h" /* atan/atan2/fabs/log/log1p */
79 #include "complex_wrapper.h"
80
81 static const double pi_2 = 1.570796326794896558e+00,
82 zero = 0.0,
83 half = 0.5,
84 two = 2.0,
85 ln2 = 6.931471805599453094172321214581765680755e-0001,
86 one = 1.0;
87
88
89 dcomplex
90 catan(dcomplex z)
91 {
92 dcomplex ans;
93 double x, y, ax, ay, t;
94 int hx, hy, ix, iy;
95 unsigned lx, ly;
96
97 x = D_RE(z);
98 y = D_IM(z);
99 ax = fabs(x);
100 ay = fabs(y);
101 hx = HI_WORD(x);
102 lx = LO_WORD(x);
103 hy = HI_WORD(y);
104 ly = LO_WORD(y);
105 ix = hx & 0x7fffffff;
106 iy = hy & 0x7fffffff;
107
108 /* x is inf or NaN */
109 if (ix >= 0x7ff00000) {
110 if (ISINF(ix, lx)) {
111 D_RE(ans) = pi_2;
112 D_IM(ans) = zero;
113 } else {
114 D_RE(ans) = x + x;
115
116 if ((iy | ly) == 0 || (ISINF(iy, ly)))
117 D_IM(ans) = zero;
118 else
119 D_IM(ans) = (fabs(y) - ay) / (fabs(y) - ay);
120 }
121 } else if (iy >= 0x7ff00000) {
122 /* y is inf or NaN */
123 if (ISINF(iy, ly)) {
124 D_RE(ans) = pi_2;
125 D_IM(ans) = zero;
126 } else {
127 D_RE(ans) = (fabs(x) - ax) / (fabs(x) - ax);
128 D_IM(ans) = y;
129 }
130 } else if ((ix | lx) == 0) {
131 /* BEGIN CSTYLED */
132 /*
133 * x = 0
134 * 1 1
135 * A = --- * atan2(2x, 1-x*x-y*y) = --- atan2(0,1-|y|)
136 * 2 2
137 *
138 * 1 [ (y+1)*(y+1) ] 1 2 1 2y
139 * B = - log [ ------------ ] = - log (1+ ---) or - log(1+ ----)
140 * 4 [ (y-1)*(y-1) ] 2 y-1 2 1-y
141 */
142 /* END CSTYLED */
143 t = one - ay;
144
145 if (((iy - 0x3ff00000) | ly) == 0) {
146 /* y=1: catan(0,1)=(0,+inf) with 1/0 signal */
147 D_IM(ans) = ay / ax;
148 D_RE(ans) = zero;
149 } else if (iy >= 0x3ff00000) { /* y>1 */
150 D_IM(ans) = half * log1p(two / (-t));
151 D_RE(ans) = pi_2;
152 } else { /* y<1 */
153 D_IM(ans) = half * log1p((ay + ay) / t);
154 D_RE(ans) = zero;
155 }
156 } else if (iy < 0x3e200000 || ((ix - iy) >> 20) >= 30) {
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 */
173 if (ix < 0x3ff00000)
174 D_RE(ans) = atan(ax);
175 else
176 D_RE(ans) = half * atan2(two, (one - ax) * (one + one /
177 ax));
178
179 if ((iy | ly) == 0) {
180 D_IM(ans) = ay;
181 } else {
182 if (ix < 0x3e200000)
183 t = ay / ((ay - one) * (ay - one));
184 else if (ix > 0x41c00000)
185 t = (ay / ax) / ax;
186 else
187 t = ay / (ax * ax + (ay - one) * (ay - one));
188
189 D_IM(ans) = t * (one - (t + t));
190 }
191 } else if (iy >= 0x41c00000 && ((iy - ix) >> 20) >= 30) {
192 /* BEGIN CSTYLED */
193 /*
194 * Huge y relative to 1+|x|
195 * |y| > Einv*(1+|x|), where Einv~2**(prec/2+3),
196 * 1
197 * A ~ --- * atan2(2x, -y*y) ~ pi/2
198 * 2
199 * y
200 * B ~ t*(1-2t), where t = --------------- is tiny
201 * (y-1)*(y-1)
202 */
203 /* END CSTYLED */
204 D_RE(ans) = pi_2;
205 t = (ay / (ay - one)) / (ay - one);
206 D_IM(ans) = t * (one - (t + t));
207 } else if (((iy - 0x3ff00000) | ly) == 0) {
208 /* BEGIN CSTYLED */
209 /*
210 * y = 1
211 * 1 1
212 * A = --- * atan2(2x, -x*x) = --- atan2(2,-x)
213 * 2 2
214 *
215 * 1 [x*x + 4] 1 4 [ 0.5(log2-logx) if
216 * B = - log [-------] = - log (1+ ---) = [ |x|<E, else 0.25*
217 * 4 [ x*x ] 4 x*x [ log1p((2/x)*(2/x))
218 */
219 /* END CSTYLED */
220 D_RE(ans) = half * atan2(two, -ax);
221
222 if (ix < 0x3e200000) {
223 D_IM(ans) = half * (ln2 - log(ax));
224 } else {
225 t = two / ax;
226 D_IM(ans) = 0.25 * log1p(t * t);
227 }
228 } else if (ix >= 0x43900000) {
229 /* BEGIN CSTYLED */
230 /*
231 * Huge x:
232 * when |x| > 1/E^2,
233 * 1 pi
234 * A ~ --- * atan2(2x, -x*x-y*y) ~ ---
235 * 2 2
236 * y y/x
237 * B ~ t*(1-2t), where t = --------------- = (-------------- )/x
238 * x*x+(y-1)*(y-1) 1+((y-1)/x)^2
239 */
240 /* END CSTYLED */
241 D_RE(ans) = pi_2;
242 t = ((ay / ax) / (one + ((ay - one) / ax) * ((ay - one) /
243 ax))) / ax;
244 D_IM(ans) = t * (one - (t + t));
245 } else if (ix < 0x38b00000) {
246 /* BEGIN CSTYLED */
247 /*
248 * Tiny x:
249 * when |x| < E^4, (note that y != 1)
250 * 1 1
251 * A = --- * atan2(2x, 1-x*x-y*y) ~ --- * atan2(2x,(1-y)*(1+y))
252 * 2 2
253 *
254 * 1 [(y+1)*(y+1)] 1 2 1 2y
255 * B = - log [-----------] = - log (1+ ---) or - log(1+ ----)
256 * 4 [(y-1)*(y-1)] 2 y-1 2 1-y
257 */
258 /* END CSTYLED */
259 D_RE(ans) = half * atan2(ax + ax, (one - ay) * (one + ay));
260
261 if (iy >= 0x3ff00000)
262 D_IM(ans) = half * log1p(two / (ay - one));
263 else
264 D_IM(ans) = half * log1p((ay + ay) / (one - ay));
265 } else {
266 /* BEGIN CSTYLED */
267 /*
268 * normal x,y
269 * 1
270 * A = --- * atan2(2x, 1-x*x-y*y)
271 * 2
272 *
273 * 1 [x*x+(y+1)*(y+1)] 1 4y
274 * B = - log [---------------] = - log (1+ -----------------)
275 * 4 [x*x+(y-1)*(y-1)] 4 x*x + (y-1)*(y-1)
276 */
277 /* END CSTYLED */
278 t = one - ay;
279
280 if (iy >= 0x3fe00000 && iy < 0x40000000) {
281 /* y close to 1 */
282 D_RE(ans) = half * (atan2((ax + ax), (t * (one + ay) -
283 ax * ax)));
284 } else if (ix >= 0x3fe00000 && ix < 0x40000000) {
285 /* x close to 1 */
286 D_RE(ans) = half * atan2((ax + ax), ((one - ax) * (one +
287 ax) - ay * ay));
288 } else {
289 D_RE(ans) = half * atan2((ax + ax), ((one - ax * ax) -
290 ay * ay));
291 }
292
293 D_IM(ans) = 0.25 * log1p((4.0 * ay) / (ax * ax + t * t));
294 }
295
296 if (hx < 0)
297 D_RE(ans) = -D_RE(ans);
298
299 if (hy < 0)
300 D_IM(ans) = -D_IM(ans);
301
302 return (ans);
303 }