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