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