Annotation of src/lib/libm/complex/catrig.c, Revision 1.2
1.2 ! christos 1: /* $NetBSD: catrig.c,v 1.1 2016/09/19 22:05:05 christos Exp $ */
1.1 christos 2: /*-
3: * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG>
4: * All rights reserved.
5: *
6: * Redistribution and use in source and binary forms, with or without
7: * modification, are permitted provided that the following conditions
8: * are met:
9: * 1. Redistributions of source code must retain the above copyright
10: * notice, this list of conditions and the following disclaimer.
11: * 2. Redistributions in binary form must reproduce the above copyright
12: * notice, this list of conditions and the following disclaimer in the
13: * documentation and/or other materials provided with the distribution.
14: *
15: * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
16: * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17: * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
18: * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
19: * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20: * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
21: * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
22: * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
23: * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
24: * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
25: * SUCH DAMAGE.
26: */
27:
28: #include <sys/cdefs.h>
29: #if 0
30: __FBSDID("$FreeBSD: head/lib/msun/src/catrig.c 275819 2014-12-16 09:21:56Z ed $");
31: #endif
1.2 ! christos 32: __RCSID("$NetBSD: catrig.c,v 1.1 2016/09/19 22:05:05 christos Exp $");
1.1 christos 33:
34: #include "namespace.h"
35: #ifdef __weak_alias
36: __weak_alias(casin, _casin)
37: #endif
38: #ifdef __weak_alias
39: __weak_alias(catan, _catan)
40: #endif
41:
42: #include <complex.h>
43: #include <float.h>
44:
45: #include "math.h"
46: #include "math_private.h"
47:
48:
49:
50: #undef isinf
51: #define isinf(x) (fabs(x) == INFINITY)
52: #undef isnan
53: #define isnan(x) ((x) != (x))
54: #define raise_inexact() do { volatile float junk __unused = /*LINTED*/1 + tiny; } while(/*CONSTCOND*/0)
55: #undef signbit
56: #define signbit(x) (__builtin_signbit(x))
57:
58: /* We need that DBL_EPSILON^2/128 is larger than FOUR_SQRT_MIN. */
59: static const double
60: A_crossover = 10, /* Hull et al suggest 1.5, but 10 works better */
61: B_crossover = 0.6417, /* suggested by Hull et al */
62: m_e = 2.7182818284590452e0, /* 0x15bf0a8b145769.0p-51 */
63: m_ln2 = 6.9314718055994531e-1, /* 0x162e42fefa39ef.0p-53 */
64: pio2_hi = 1.5707963267948966e0, /* 0x1921fb54442d18.0p-52 */
65: RECIP_EPSILON = 1 / DBL_EPSILON,
66: SQRT_3_EPSILON = 2.5809568279517849e-8, /* 0x1bb67ae8584caa.0p-78 */
67: SQRT_6_EPSILON = 3.6500241499888571e-8, /* 0x13988e1409212e.0p-77 */
1.2 ! christos 68: #if DBL_MAX_EXP == 1024 /* IEEE */
! 69: FOUR_SQRT_MIN = 0x1p-509, /* >= 4 * sqrt(DBL_MIN) */
! 70: QUARTER_SQRT_MAX = 0x1p509, /* <= sqrt(DBL_MAX) / 4 */
1.1 christos 71: SQRT_MIN = 0x1p-511; /* >= sqrt(DBL_MIN) */
1.2 ! christos 72: #elif DBL_MAX_EXP == 127 /* VAX */
! 73: FOUR_SQRT_MIN = 0x1p-62, /* >= 4 * sqrt(DBL_MIN) */
! 74: QUARTER_SQRT_MAX = 0x1p62, /* <= sqrt(DBL_MAX) / 4 */
! 75: SQRT_MIN = 0x1p-64; /* >= sqrt(DBL_MIN) */
! 76: #else
! 77: #error "unsupported floating point format"
! 78: #endif
! 79:
1.1 christos 80:
81: static const volatile double
82: pio2_lo = 6.1232339957367659e-17; /* 0x11a62633145c07.0p-106 */
83: static const volatile float
84: tiny = 0x1p-100;
85:
86: static double complex clog_for_large_values(double complex z);
87:
88: /*
89: * Testing indicates that all these functions are accurate up to 4 ULP.
90: * The functions casin(h) and cacos(h) are about 2.5 times slower than asinh.
91: * The functions catan(h) are a little under 2 times slower than atanh.
92: *
93: * The code for casinh, casin, cacos, and cacosh comes first. The code is
94: * rather complicated, and the four functions are highly interdependent.
95: *
96: * The code for catanh and catan comes at the end. It is much simpler than
97: * the other functions, and the code for these can be disconnected from the
98: * rest of the code.
99: */
100:
101: /*
102: * ================================
103: * | casinh, casin, cacos, cacosh |
104: * ================================
105: */
106:
107: /*
108: * The algorithm is very close to that in "Implementing the complex arcsine
109: * and arccosine functions using exception handling" by T. E. Hull, Thomas F.
110: * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on
111: * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335,
112: * http://dl.acm.org/citation.cfm?id=275324.
113: *
114: * Throughout we use the convention z = x + I*y.
115: *
116: * casinh(z) = sign(x)*log(A+sqrt(A*A-1)) + I*asin(B)
117: * where
118: * A = (|z+I| + |z-I|) / 2
119: * B = (|z+I| - |z-I|) / 2 = y/A
120: *
121: * These formulas become numerically unstable:
122: * (a) for Re(casinh(z)) when z is close to the line segment [-I, I] (that
123: * is, Re(casinh(z)) is close to 0);
124: * (b) for Im(casinh(z)) when z is close to either of the intervals
125: * [I, I*infinity) or (-I*infinity, -I] (that is, |Im(casinh(z))| is
126: * close to PI/2).
127: *
128: * These numerical problems are overcome by defining
129: * f(a, b) = (hypot(a, b) - b) / 2 = a*a / (hypot(a, b) + b) / 2
130: * Then if A < A_crossover, we use
131: * log(A + sqrt(A*A-1)) = log1p((A-1) + sqrt((A-1)*(A+1)))
132: * A-1 = f(x, 1+y) + f(x, 1-y)
133: * and if B > B_crossover, we use
134: * asin(B) = atan2(y, sqrt(A*A - y*y)) = atan2(y, sqrt((A+y)*(A-y)))
135: * A-y = f(x, y+1) + f(x, y-1)
136: * where without loss of generality we have assumed that x and y are
137: * non-negative.
138: *
139: * Much of the difficulty comes because the intermediate computations may
140: * produce overflows or underflows. This is dealt with in the paper by Hull
141: * et al by using exception handling. We do this by detecting when
142: * computations risk underflow or overflow. The hardest part is handling the
143: * underflows when computing f(a, b).
144: *
145: * Note that the function f(a, b) does not appear explicitly in the paper by
146: * Hull et al, but the idea may be found on pages 308 and 309. Introducing the
147: * function f(a, b) allows us to concentrate many of the clever tricks in this
148: * paper into one function.
149: */
150:
151: /*
152: * Function f(a, b, hypot_a_b) = (hypot(a, b) - b) / 2.
153: * Pass hypot(a, b) as the third argument.
154: */
155: static inline double
156: f(double a, double b, double hypot_a_b)
157: {
158: if (b < 0)
159: return ((hypot_a_b - b) / 2);
160: if (b == 0)
161: return (a / 2);
162: return (a * a / (hypot_a_b + b) / 2);
163: }
164:
165: /*
166: * All the hard work is contained in this function.
167: * x and y are assumed positive or zero, and less than RECIP_EPSILON.
168: * Upon return:
169: * rx = Re(casinh(z)) = -Im(cacos(y + I*x)).
170: * B_is_usable is set to 1 if the value of B is usable.
171: * If B_is_usable is set to 0, sqrt_A2my2 = sqrt(A*A - y*y), and new_y = y.
172: * If returning sqrt_A2my2 has potential to result in an underflow, it is
173: * rescaled, and new_y is similarly rescaled.
174: */
175: static inline void
176: do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B,
177: double *sqrt_A2my2, double *new_y)
178: {
179: double R, S, A; /* A, B, R, and S are as in Hull et al. */
180: double Am1, Amy; /* A-1, A-y. */
181:
182: R = hypot(x, y + 1); /* |z+I| */
183: S = hypot(x, y - 1); /* |z-I| */
184:
185: /* A = (|z+I| + |z-I|) / 2 */
186: A = (R + S) / 2;
187: /*
188: * Mathematically A >= 1. There is a small chance that this will not
189: * be so because of rounding errors. So we will make certain it is
190: * so.
191: */
192: if (A < 1)
193: A = 1;
194:
195: if (A < A_crossover) {
196: /*
197: * Am1 = fp + fm, where fp = f(x, 1+y), and fm = f(x, 1-y).
198: * rx = log1p(Am1 + sqrt(Am1*(A+1)))
199: */
200: if (y == 1 && x < DBL_EPSILON * DBL_EPSILON / 128) {
201: /*
202: * fp is of order x^2, and fm = x/2.
203: * A = 1 (inexactly).
204: */
205: *rx = sqrt(x);
206: } else if (x >= DBL_EPSILON * fabs(y - 1)) {
207: /*
208: * Underflow will not occur because
209: * x >= DBL_EPSILON^2/128 >= FOUR_SQRT_MIN
210: */
211: Am1 = f(x, 1 + y, R) + f(x, 1 - y, S);
212: *rx = log1p(Am1 + sqrt(Am1 * (A + 1)));
213: } else if (y < 1) {
214: /*
215: * fp = x*x/(1+y)/4, fm = x*x/(1-y)/4, and
216: * A = 1 (inexactly).
217: */
218: *rx = x / sqrt((1 - y) * (1 + y));
219: } else { /* if (y > 1) */
220: /*
221: * A-1 = y-1 (inexactly).
222: */
223: *rx = log1p((y - 1) + sqrt((y - 1) * (y + 1)));
224: }
225: } else {
226: *rx = log(A + sqrt(A * A - 1));
227: }
228:
229: *new_y = y;
230:
231: if (y < FOUR_SQRT_MIN) {
232: /*
233: * Avoid a possible underflow caused by y/A. For casinh this
234: * would be legitimate, but will be picked up by invoking atan2
235: * later on. For cacos this would not be legitimate.
236: */
237: *B_is_usable = 0;
238: *sqrt_A2my2 = A * (2 / DBL_EPSILON);
239: *new_y = y * (2 / DBL_EPSILON);
240: return;
241: }
242:
243: /* B = (|z+I| - |z-I|) / 2 = y/A */
244: *B = y / A;
245: *B_is_usable = 1;
246:
247: if (*B > B_crossover) {
248: *B_is_usable = 0;
249: /*
250: * Amy = fp + fm, where fp = f(x, y+1), and fm = f(x, y-1).
251: * sqrt_A2my2 = sqrt(Amy*(A+y))
252: */
253: if (y == 1 && x < DBL_EPSILON / 128) {
254: /*
255: * fp is of order x^2, and fm = x/2.
256: * A = 1 (inexactly).
257: */
258: *sqrt_A2my2 = sqrt(x) * sqrt((A + y) / 2);
259: } else if (x >= DBL_EPSILON * fabs(y - 1)) {
260: /*
261: * Underflow will not occur because
262: * x >= DBL_EPSILON/128 >= FOUR_SQRT_MIN
263: * and
264: * x >= DBL_EPSILON^2 >= FOUR_SQRT_MIN
265: */
266: Amy = f(x, y + 1, R) + f(x, y - 1, S);
267: *sqrt_A2my2 = sqrt(Amy * (A + y));
268: } else if (y > 1) {
269: /*
270: * fp = x*x/(y+1)/4, fm = x*x/(y-1)/4, and
271: * A = y (inexactly).
272: *
273: * y < RECIP_EPSILON. So the following
274: * scaling should avoid any underflow problems.
275: */
276: *sqrt_A2my2 = x * (4 / DBL_EPSILON / DBL_EPSILON) * y /
277: sqrt((y + 1) * (y - 1));
278: *new_y = y * (4 / DBL_EPSILON / DBL_EPSILON);
279: } else { /* if (y < 1) */
280: /*
281: * fm = 1-y >= DBL_EPSILON, fp is of order x^2, and
282: * A = 1 (inexactly).
283: */
284: *sqrt_A2my2 = sqrt((1 - y) * (1 + y));
285: }
286: }
287: }
288:
289: /*
290: * casinh(z) = z + O(z^3) as z -> 0
291: *
292: * casinh(z) = sign(x)*clog(sign(x)*z) + O(1/z^2) as z -> infinity
293: * The above formula works for the imaginary part as well, because
294: * Im(casinh(z)) = sign(x)*atan2(sign(x)*y, fabs(x)) + O(y/z^3)
295: * as z -> infinity, uniformly in y
296: */
297: double complex
298: casinh(double complex z)
299: {
300: double x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y;
301: int B_is_usable;
302: double complex w;
303:
304: x = creal(z);
305: y = cimag(z);
306: ax = fabs(x);
307: ay = fabs(y);
308:
309: if (isnan(x) || isnan(y)) {
310: /* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */
311: if (isinf(x))
312: return (CMPLX(x, y + y));
313: /* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */
314: if (isinf(y))
315: return (CMPLX(y, x + x));
316: /* casinh(NaN + I*0) = NaN + I*0 */
317: if (y == 0)
318: return (CMPLX(x + x, y));
319: /*
320: * All other cases involving NaN return NaN + I*NaN.
321: * C99 leaves it optional whether to raise invalid if one of
322: * the arguments is not NaN, so we opt not to raise it.
323: */
324: return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
325: }
326:
327: if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
328: /* clog...() will raise inexact unless x or y is infinite. */
329: if (signbit(x) == 0)
330: w = clog_for_large_values(z) + m_ln2;
331: else
332: w = clog_for_large_values(-z) + m_ln2;
333: return (CMPLX(copysign(creal(w), x), copysign(cimag(w), y)));
334: }
335:
336: /* Avoid spuriously raising inexact for z = 0. */
337: if (x == 0 && y == 0)
338: return (z);
339:
340: /* All remaining cases are inexact. */
341: raise_inexact();
342:
343: if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
344: return (z);
345:
346: do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y);
347: if (B_is_usable)
348: ry = asin(B);
349: else
350: ry = atan2(new_y, sqrt_A2my2);
351: return (CMPLX(copysign(rx, x), copysign(ry, y)));
352: }
353:
354: /*
355: * casin(z) = reverse(casinh(reverse(z)))
356: * where reverse(x + I*y) = y + I*x = I*conj(z).
357: */
358: double complex
359: casin(double complex z)
360: {
361: double complex w = casinh(CMPLX(cimag(z), creal(z)));
362:
363: return (CMPLX(cimag(w), creal(w)));
364: }
365:
366: /*
367: * cacos(z) = PI/2 - casin(z)
368: * but do the computation carefully so cacos(z) is accurate when z is
369: * close to 1.
370: *
371: * cacos(z) = PI/2 - z + O(z^3) as z -> 0
372: *
373: * cacos(z) = -sign(y)*I*clog(z) + O(1/z^2) as z -> infinity
374: * The above formula works for the real part as well, because
375: * Re(cacos(z)) = atan2(fabs(y), x) + O(y/z^3)
376: * as z -> infinity, uniformly in y
377: */
378: double complex
379: cacos(double complex z)
380: {
381: double x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x;
382: int sx, sy;
383: int B_is_usable;
384: double complex w;
385:
386: x = creal(z);
387: y = cimag(z);
388: sx = signbit(x);
389: sy = signbit(y);
390: ax = fabs(x);
391: ay = fabs(y);
392:
393: if (isnan(x) || isnan(y)) {
394: /* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */
395: if (isinf(x))
396: return (CMPLX(y + y, -INFINITY));
397: /* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */
398: if (isinf(y))
399: return (CMPLX(x + x, -y));
400: /* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */
401: if (x == 0)
402: return (CMPLX(pio2_hi + pio2_lo, y + y));
403: /*
404: * All other cases involving NaN return NaN + I*NaN.
405: * C99 leaves it optional whether to raise invalid if one of
406: * the arguments is not NaN, so we opt not to raise it.
407: */
408: return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
409: }
410:
411: if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
412: /* clog...() will raise inexact unless x or y is infinite. */
413: w = clog_for_large_values(z);
414: rx = fabs(cimag(w));
415: ry = creal(w) + m_ln2;
416: if (sy == 0)
417: ry = -ry;
418: return (CMPLX(rx, ry));
419: }
420:
421: /* Avoid spuriously raising inexact for z = 1. */
422: if (x == 1 && y == 0)
423: return (CMPLX(0, -y));
424:
425: /* All remaining cases are inexact. */
426: raise_inexact();
427:
428: if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
429: return (CMPLX(pio2_hi - (x - pio2_lo), -y));
430:
431: do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
432: if (B_is_usable) {
433: if (sx == 0)
434: rx = acos(B);
435: else
436: rx = acos(-B);
437: } else {
438: if (sx == 0)
439: rx = atan2(sqrt_A2mx2, new_x);
440: else
441: rx = atan2(sqrt_A2mx2, -new_x);
442: }
443: if (sy == 0)
444: ry = -ry;
445: return (CMPLX(rx, ry));
446: }
447:
448: /*
449: * cacosh(z) = I*cacos(z) or -I*cacos(z)
450: * where the sign is chosen so Re(cacosh(z)) >= 0.
451: */
452: double complex
453: cacosh(double complex z)
454: {
455: double complex w;
456: double rx, ry;
457:
458: w = cacos(z);
459: rx = creal(w);
460: ry = cimag(w);
461: /* cacosh(NaN + I*NaN) = NaN + I*NaN */
462: if (isnan(rx) && isnan(ry))
463: return (CMPLX(ry, rx));
464: /* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */
465: /* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */
466: if (isnan(rx))
467: return (CMPLX(fabs(ry), rx));
468: /* cacosh(0 + I*NaN) = NaN + I*NaN */
469: if (isnan(ry))
470: return (CMPLX(ry, ry));
471: return (CMPLX(fabs(ry), copysign(rx, cimag(z))));
472: }
473:
474: /*
475: * Optimized version of clog() for |z| finite and larger than ~RECIP_EPSILON.
476: */
477: static double complex
478: clog_for_large_values(double complex z)
479: {
480: double x, y;
481: double ax, ay, t;
482:
483: x = creal(z);
484: y = cimag(z);
485: ax = fabs(x);
486: ay = fabs(y);
487: if (ax < ay) {
488: t = ax;
489: ax = ay;
490: ay = t;
491: }
492:
493: /*
494: * Avoid overflow in hypot() when x and y are both very large.
495: * Divide x and y by E, and then add 1 to the logarithm. This depends
496: * on E being larger than sqrt(2).
497: * Dividing by E causes an insignificant loss of accuracy; however
498: * this method is still poor since it is uneccessarily slow.
499: */
500: if (ax > DBL_MAX / 2)
501: return (CMPLX(log(hypot(x / m_e, y / m_e)) + 1, atan2(y, x)));
502:
503: /*
504: * Avoid overflow when x or y is large. Avoid underflow when x or
505: * y is small.
506: */
507: if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
508: return (CMPLX(log(hypot(x, y)), atan2(y, x)));
509:
510: return (CMPLX(log(ax * ax + ay * ay) / 2, atan2(y, x)));
511: }
512:
513: /*
514: * =================
515: * | catanh, catan |
516: * =================
517: */
518:
519: /*
520: * sum_squares(x,y) = x*x + y*y (or just x*x if y*y would underflow).
521: * Assumes x*x and y*y will not overflow.
522: * Assumes x and y are finite.
523: * Assumes y is non-negative.
524: * Assumes fabs(x) >= DBL_EPSILON.
525: */
526: static inline double
527: sum_squares(double x, double y)
528: {
529:
530: /* Avoid underflow when y is small. */
531: if (y < SQRT_MIN)
532: return (x * x);
533:
534: return (x * x + y * y);
535: }
536:
537: /*
538: * real_part_reciprocal(x, y) = Re(1/(x+I*y)) = x/(x*x + y*y).
539: * Assumes x and y are not NaN, and one of x and y is larger than
540: * RECIP_EPSILON. We avoid unwarranted underflow. It is important to not use
541: * the code creal(1/z), because the imaginary part may produce an unwanted
542: * underflow.
543: * This is only called in a context where inexact is always raised before
544: * the call, so no effort is made to avoid or force inexact.
545: */
546: static inline double
547: real_part_reciprocal(double x, double y)
548: {
549: double scale;
550: uint32_t hx, hy;
551: int32_t ix, iy;
552:
553: /*
554: * This code is inspired by the C99 document n1124.pdf, Section G.5.1,
555: * example 2.
556: */
557: GET_HIGH_WORD(hx, x);
558: ix = hx & 0x7ff00000;
559: GET_HIGH_WORD(hy, y);
560: iy = hy & 0x7ff00000;
561: #define BIAS (DBL_MAX_EXP - 1)
562: /* XXX more guard digits are useful iff there is extra precision. */
563: #define CUTOFF (DBL_MANT_DIG / 2 + 1) /* just half or 1 guard digit */
564: if (ix - iy >= CUTOFF << 20 || isinf(x))
565: return (1 / x); /* +-Inf -> +-0 is special */
566: if (iy - ix >= CUTOFF << 20)
567: return (x / y / y); /* should avoid double div, but hard */
568: if (ix <= (BIAS + DBL_MAX_EXP / 2 - CUTOFF) << 20)
569: return (x / (x * x + y * y));
570: scale = 1;
571: SET_HIGH_WORD(scale, 0x7ff00000 - ix); /* 2**(1-ilogb(x)) */
572: x *= scale;
573: y *= scale;
574: return (x / (x * x + y * y) * scale);
575: }
576:
577: /*
578: * catanh(z) = log((1+z)/(1-z)) / 2
579: * = log1p(4*x / |z-1|^2) / 4
580: * + I * atan2(2*y, (1-x)*(1+x)-y*y) / 2
581: *
582: * catanh(z) = z + O(z^3) as z -> 0
583: *
584: * catanh(z) = 1/z + sign(y)*I*PI/2 + O(1/z^3) as z -> infinity
585: * The above formula works for the real part as well, because
586: * Re(catanh(z)) = x/|z|^2 + O(x/z^4)
587: * as z -> infinity, uniformly in x
588: */
589: double complex
590: catanh(double complex z)
591: {
592: double x, y, ax, ay, rx, ry;
593:
594: x = creal(z);
595: y = cimag(z);
596: ax = fabs(x);
597: ay = fabs(y);
598:
599: /* This helps handle many cases. */
600: if (y == 0 && ax <= 1)
601: return (CMPLX(atanh(x), y));
602:
603: /* To ensure the same accuracy as atan(), and to filter out z = 0. */
604: if (x == 0)
605: return (CMPLX(x, atan(y)));
606:
607: if (isnan(x) || isnan(y)) {
608: /* catanh(+-Inf + I*NaN) = +-0 + I*NaN */
609: if (isinf(x))
610: return (CMPLX(copysign(0, x), y + y));
611: /* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */
612: if (isinf(y))
613: return (CMPLX(copysign(0, x),
614: copysign(pio2_hi + pio2_lo, y)));
615: /*
616: * All other cases involving NaN return NaN + I*NaN.
617: * C99 leaves it optional whether to raise invalid if one of
618: * the arguments is not NaN, so we opt not to raise it.
619: */
620: return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
621: }
622:
623: if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
624: return (CMPLX(real_part_reciprocal(x, y),
625: copysign(pio2_hi + pio2_lo, y)));
626:
627: if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
628: /*
629: * z = 0 was filtered out above. All other cases must raise
630: * inexact, but this is the only only that needs to do it
631: * explicitly.
632: */
633: raise_inexact();
634: return (z);
635: }
636:
637: if (ax == 1 && ay < DBL_EPSILON)
638: rx = (m_ln2 - log(ay)) / 2;
639: else
640: rx = log1p(4 * ax / sum_squares(ax - 1, ay)) / 4;
641:
642: if (ax == 1)
643: ry = atan2(2, -ay) / 2;
644: else if (ay < DBL_EPSILON)
645: ry = atan2(2 * ay, (1 - ax) * (1 + ax)) / 2;
646: else
647: ry = atan2(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
648:
649: return (CMPLX(copysign(rx, x), copysign(ry, y)));
650: }
651:
652: /*
653: * catan(z) = reverse(catanh(reverse(z)))
654: * where reverse(x + I*y) = y + I*x = I*conj(z).
655: */
656: double complex
657: catan(double complex z)
658: {
659: double complex w = catanh(CMPLX(cimag(z), creal(z)));
660:
661: return (CMPLX(cimag(w), creal(w)));
662: }
CVSweb <webmaster@jp.NetBSD.org>
![[BACK]](/icons/back.gif){kind=icon}
Return to catrig.c CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [cvs.NetBSD.org] / src / lib / libm / complex