xref: /titanic_50/usr/src/lib/libm/common/complex/catan.c (revision ddc0e0b53c661f6e439e3b7072b3ef353eadb4af)
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
catan(dcomplex z)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 }
293