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