1251121Sdas/*- 2251121Sdas * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG> 3251121Sdas * All rights reserved. 4251121Sdas * 5251121Sdas * Redistribution and use in source and binary forms, with or without 6251121Sdas * modification, are permitted provided that the following conditions 7251121Sdas * are met: 8251121Sdas * 1. Redistributions of source code must retain the above copyright 9251121Sdas * notice, this list of conditions and the following disclaimer. 10251121Sdas * 2. Redistributions in binary form must reproduce the above copyright 11251121Sdas * notice, this list of conditions and the following disclaimer in the 12251121Sdas * documentation and/or other materials provided with the distribution. 13251121Sdas * 14251121Sdas * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 15251121Sdas * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 16251121Sdas * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 17251121Sdas * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 18251121Sdas * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 19251121Sdas * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 20251121Sdas * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 21251121Sdas * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 22251121Sdas * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 23251121Sdas * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 24251121Sdas * SUCH DAMAGE. 25251121Sdas */ 26251121Sdas 27251121Sdas/* 28251121Sdas * The algorithm is very close to that in "Implementing the complex arcsine 29251121Sdas * and arccosine functions using exception handling" by T. E. Hull, Thomas F. 30251121Sdas * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on 31251121Sdas * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335, 32251121Sdas * http://dl.acm.org/citation.cfm?id=275324. 33251121Sdas * 34251404Sdas * See catrig.c for complete comments. 35251404Sdas * 36251404Sdas * XXX comments were removed automatically, and even short ones on the right 37251404Sdas * of statements were removed (all of them), contrary to normal style. Only 38251404Sdas * a few comments on the right of declarations remain. 39251121Sdas */ 40251121Sdas 41251121Sdas#include <sys/cdefs.h> 42251121Sdas__FBSDID("$FreeBSD$"); 43251121Sdas 44251121Sdas#include <complex.h> 45251121Sdas#include <float.h> 46251121Sdas 47251121Sdas#include "math.h" 48251121Sdas#include "math_private.h" 49251121Sdas 50251121Sdas#undef isinf 51251121Sdas#define isinf(x) (fabsf(x) == INFINITY) 52251121Sdas#undef isnan 53251121Sdas#define isnan(x) ((x) != (x)) 54251121Sdas#define raise_inexact() do { volatile float junk = 1 + tiny; } while(0) 55251121Sdas#undef signbit 56251121Sdas#define signbit(x) (__builtin_signbitf(x)) 57251121Sdas 58251121Sdasstatic const float 59251121SdasA_crossover = 10, 60251121SdasB_crossover = 0.6417, 61251121SdasFOUR_SQRT_MIN = 0x1p-61, 62251121SdasQUARTER_SQRT_MAX = 0x1p61, 63251121Sdasm_e = 2.7182818285e0, /* 0xadf854.0p-22 */ 64251121Sdasm_ln2 = 6.9314718056e-1, /* 0xb17218.0p-24 */ 65251121Sdaspio2_hi = 1.5707962513e0, /* 0xc90fda.0p-23 */ 66251121SdasRECIP_EPSILON = 1 / FLT_EPSILON, 67251121SdasSQRT_3_EPSILON = 5.9801995673e-4, /* 0x9cc471.0p-34 */ 68251121SdasSQRT_6_EPSILON = 8.4572793338e-4, /* 0xddb3d7.0p-34 */ 69251121SdasSQRT_MIN = 0x1p-63; 70251121Sdas 71251121Sdasstatic const volatile float 72251121Sdaspio2_lo = 7.5497899549e-8, /* 0xa22169.0p-47 */ 73251121Sdastiny = 0x1p-100; 74251121Sdas 75251121Sdasstatic float complex clog_for_large_values(float complex z); 76251121Sdas 77251121Sdasstatic inline float 78251121Sdasf(float a, float b, float hypot_a_b) 79251121Sdas{ 80251121Sdas if (b < 0) 81251121Sdas return ((hypot_a_b - b) / 2); 82251121Sdas if (b == 0) 83251121Sdas return (a / 2); 84251121Sdas return (a * a / (hypot_a_b + b) / 2); 85251121Sdas} 86251121Sdas 87251121Sdasstatic inline void 88251121Sdasdo_hard_work(float x, float y, float *rx, int *B_is_usable, float *B, 89251404Sdas float *sqrt_A2my2, float *new_y) 90251121Sdas{ 91251121Sdas float R, S, A; 92251121Sdas float Am1, Amy; 93251121Sdas 94251121Sdas R = hypotf(x, y + 1); 95251121Sdas S = hypotf(x, y - 1); 96251121Sdas 97251121Sdas A = (R + S) / 2; 98251121Sdas if (A < 1) 99251121Sdas A = 1; 100251121Sdas 101251121Sdas if (A < A_crossover) { 102251121Sdas if (y == 1 && x < FLT_EPSILON * FLT_EPSILON / 128) { 103251121Sdas *rx = sqrtf(x); 104251121Sdas } else if (x >= FLT_EPSILON * fabsf(y - 1)) { 105251121Sdas Am1 = f(x, 1 + y, R) + f(x, 1 - y, S); 106251121Sdas *rx = log1pf(Am1 + sqrtf(Am1 * (A + 1))); 107251121Sdas } else if (y < 1) { 108251404Sdas *rx = x / sqrtf((1 - y) * (1 + y)); 109251121Sdas } else { 110251121Sdas *rx = log1pf((y - 1) + sqrtf((y - 1) * (y + 1))); 111251121Sdas } 112251121Sdas } else { 113251121Sdas *rx = logf(A + sqrtf(A * A - 1)); 114251121Sdas } 115251121Sdas 116251121Sdas *new_y = y; 117251121Sdas 118251121Sdas if (y < FOUR_SQRT_MIN) { 119251121Sdas *B_is_usable = 0; 120251121Sdas *sqrt_A2my2 = A * (2 / FLT_EPSILON); 121251121Sdas *new_y = y * (2 / FLT_EPSILON); 122251121Sdas return; 123251121Sdas } 124251121Sdas 125251121Sdas *B = y / A; 126251121Sdas *B_is_usable = 1; 127251121Sdas 128251121Sdas if (*B > B_crossover) { 129251121Sdas *B_is_usable = 0; 130251121Sdas if (y == 1 && x < FLT_EPSILON / 128) { 131251121Sdas *sqrt_A2my2 = sqrtf(x) * sqrtf((A + y) / 2); 132251121Sdas } else if (x >= FLT_EPSILON * fabsf(y - 1)) { 133251121Sdas Amy = f(x, y + 1, R) + f(x, y - 1, S); 134251121Sdas *sqrt_A2my2 = sqrtf(Amy * (A + y)); 135251121Sdas } else if (y > 1) { 136251121Sdas *sqrt_A2my2 = x * (4 / FLT_EPSILON / FLT_EPSILON) * y / 137251404Sdas sqrtf((y + 1) * (y - 1)); 138251121Sdas *new_y = y * (4 / FLT_EPSILON / FLT_EPSILON); 139251121Sdas } else { 140251121Sdas *sqrt_A2my2 = sqrtf((1 - y) * (1 + y)); 141251121Sdas } 142251121Sdas } 143251121Sdas} 144251121Sdas 145251121Sdasfloat complex 146251121Sdascasinhf(float complex z) 147251121Sdas{ 148251121Sdas float x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y; 149251121Sdas int B_is_usable; 150251121Sdas float complex w; 151251121Sdas 152251121Sdas x = crealf(z); 153251121Sdas y = cimagf(z); 154251121Sdas ax = fabsf(x); 155251121Sdas ay = fabsf(y); 156251121Sdas 157251121Sdas if (isnan(x) || isnan(y)) { 158251121Sdas if (isinf(x)) 159251121Sdas return (cpackf(x, y + y)); 160251121Sdas if (isinf(y)) 161251121Sdas return (cpackf(y, x + x)); 162251121Sdas if (y == 0) 163251121Sdas return (cpackf(x + x, y)); 164251121Sdas return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); 165251121Sdas } 166251121Sdas 167251121Sdas if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { 168251121Sdas if (signbit(x) == 0) 169251121Sdas w = clog_for_large_values(z) + m_ln2; 170251121Sdas else 171251121Sdas w = clog_for_large_values(-z) + m_ln2; 172251121Sdas return (cpackf(copysignf(crealf(w), x), 173251404Sdas copysignf(cimagf(w), y))); 174251121Sdas } 175251121Sdas 176251121Sdas if (x == 0 && y == 0) 177251121Sdas return (z); 178251121Sdas 179251121Sdas raise_inexact(); 180251121Sdas 181251121Sdas if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) 182251121Sdas return (z); 183251121Sdas 184251121Sdas do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y); 185251121Sdas if (B_is_usable) 186251121Sdas ry = asinf(B); 187251121Sdas else 188251121Sdas ry = atan2f(new_y, sqrt_A2my2); 189251121Sdas return (cpackf(copysignf(rx, x), copysignf(ry, y))); 190251121Sdas} 191251121Sdas 192251121Sdasfloat complex 193251121Sdascasinf(float complex z) 194251121Sdas{ 195251121Sdas float complex w = casinhf(cpackf(cimagf(z), crealf(z))); 196251404Sdas 197251121Sdas return (cpackf(cimagf(w), crealf(w))); 198251121Sdas} 199251121Sdas 200251121Sdasfloat complex 201251121Sdascacosf(float complex z) 202251121Sdas{ 203251121Sdas float x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x; 204251121Sdas int sx, sy; 205251121Sdas int B_is_usable; 206251121Sdas float complex w; 207251121Sdas 208251121Sdas x = crealf(z); 209251121Sdas y = cimagf(z); 210251121Sdas sx = signbit(x); 211251121Sdas sy = signbit(y); 212251121Sdas ax = fabsf(x); 213251121Sdas ay = fabsf(y); 214251121Sdas 215251121Sdas if (isnan(x) || isnan(y)) { 216251121Sdas if (isinf(x)) 217251121Sdas return (cpackf(y + y, -INFINITY)); 218251121Sdas if (isinf(y)) 219251121Sdas return (cpackf(x + x, -y)); 220251404Sdas if (x == 0) 221251404Sdas return (cpackf(pio2_hi + pio2_lo, y + y)); 222251121Sdas return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); 223251121Sdas } 224251121Sdas 225251121Sdas if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { 226251121Sdas w = clog_for_large_values(z); 227251121Sdas rx = fabsf(cimagf(w)); 228251121Sdas ry = crealf(w) + m_ln2; 229251121Sdas if (sy == 0) 230251121Sdas ry = -ry; 231251121Sdas return (cpackf(rx, ry)); 232251121Sdas } 233251121Sdas 234251121Sdas if (x == 1 && y == 0) 235251121Sdas return (cpackf(0, -y)); 236251121Sdas 237251121Sdas raise_inexact(); 238251121Sdas 239251121Sdas if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) 240251121Sdas return (cpackf(pio2_hi - (x - pio2_lo), -y)); 241251121Sdas 242251121Sdas do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x); 243251121Sdas if (B_is_usable) { 244251404Sdas if (sx == 0) 245251121Sdas rx = acosf(B); 246251121Sdas else 247251121Sdas rx = acosf(-B); 248251121Sdas } else { 249251404Sdas if (sx == 0) 250251121Sdas rx = atan2f(sqrt_A2mx2, new_x); 251251121Sdas else 252251121Sdas rx = atan2f(sqrt_A2mx2, -new_x); 253251121Sdas } 254251404Sdas if (sy == 0) 255251121Sdas ry = -ry; 256251121Sdas return (cpackf(rx, ry)); 257251121Sdas} 258251121Sdas 259251121Sdasfloat complex 260251121Sdascacoshf(float complex z) 261251121Sdas{ 262251121Sdas float complex w; 263251121Sdas float rx, ry; 264251121Sdas 265251121Sdas w = cacosf(z); 266251121Sdas rx = crealf(w); 267251121Sdas ry = cimagf(w); 268251121Sdas if (isnan(rx) && isnan(ry)) 269251121Sdas return (cpackf(ry, rx)); 270251121Sdas if (isnan(rx)) 271251121Sdas return (cpackf(fabsf(ry), rx)); 272251121Sdas if (isnan(ry)) 273251121Sdas return (cpackf(ry, ry)); 274251121Sdas return (cpackf(fabsf(ry), copysignf(rx, cimagf(z)))); 275251121Sdas} 276251121Sdas 277251121Sdasstatic float complex 278251121Sdasclog_for_large_values(float complex z) 279251121Sdas{ 280251121Sdas float x, y; 281251121Sdas float ax, ay, t; 282251121Sdas 283251121Sdas x = crealf(z); 284251121Sdas y = cimagf(z); 285251121Sdas ax = fabsf(x); 286251121Sdas ay = fabsf(y); 287251121Sdas if (ax < ay) { 288251121Sdas t = ax; 289251121Sdas ax = ay; 290251121Sdas ay = t; 291251121Sdas } 292251121Sdas 293251404Sdas if (ax > FLT_MAX / 2) 294251121Sdas return (cpackf(logf(hypotf(x / m_e, y / m_e)) + 1, 295251404Sdas atan2f(y, x))); 296251121Sdas 297251121Sdas if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) 298251121Sdas return (cpackf(logf(hypotf(x, y)), atan2f(y, x))); 299251121Sdas 300251121Sdas return (cpackf(logf(ax * ax + ay * ay) / 2, atan2f(y, x))); 301251121Sdas} 302251121Sdas 303251121Sdasstatic inline float 304251121Sdassum_squares(float x, float y) 305251121Sdas{ 306251121Sdas 307251121Sdas if (y < SQRT_MIN) 308251404Sdas return (x * x); 309251404Sdas 310251404Sdas return (x * x + y * y); 311251121Sdas} 312251121Sdas 313251121Sdasstatic inline float 314251121Sdasreal_part_reciprocal(float x, float y) 315251121Sdas{ 316251121Sdas float scale; 317251121Sdas uint32_t hx, hy; 318251121Sdas int32_t ix, iy; 319251121Sdas 320251121Sdas GET_FLOAT_WORD(hx, x); 321251121Sdas ix = hx & 0x7f800000; 322251121Sdas GET_FLOAT_WORD(hy, y); 323251121Sdas iy = hy & 0x7f800000; 324251121Sdas#define BIAS (FLT_MAX_EXP - 1) 325251121Sdas#define CUTOFF (FLT_MANT_DIG / 2 + 1) 326251121Sdas if (ix - iy >= CUTOFF << 23 || isinf(x)) 327251404Sdas return (1 / x); 328251121Sdas if (iy - ix >= CUTOFF << 23) 329251404Sdas return (x / y / y); 330251121Sdas if (ix <= (BIAS + FLT_MAX_EXP / 2 - CUTOFF) << 23) 331251121Sdas return (x / (x * x + y * y)); 332251121Sdas SET_FLOAT_WORD(scale, 0x7f800000 - ix); 333251121Sdas x *= scale; 334251121Sdas y *= scale; 335251121Sdas return (x / (x * x + y * y) * scale); 336251121Sdas} 337251121Sdas 338251121Sdasfloat complex 339251121Sdascatanhf(float complex z) 340251121Sdas{ 341251121Sdas float x, y, ax, ay, rx, ry; 342251121Sdas 343251121Sdas x = crealf(z); 344251121Sdas y = cimagf(z); 345251121Sdas ax = fabsf(x); 346251121Sdas ay = fabsf(y); 347251121Sdas 348251121Sdas if (y == 0 && ax <= 1) 349251404Sdas return (cpackf(atanhf(x), y)); 350251121Sdas 351251121Sdas if (x == 0) 352251121Sdas return (cpackf(x, atanf(y))); 353251121Sdas 354251121Sdas if (isnan(x) || isnan(y)) { 355251121Sdas if (isinf(x)) 356251404Sdas return (cpackf(copysignf(0, x), y + y)); 357251404Sdas if (isinf(y)) 358251121Sdas return (cpackf(copysignf(0, x), 359251404Sdas copysignf(pio2_hi + pio2_lo, y))); 360251121Sdas return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); 361251121Sdas } 362251121Sdas 363251404Sdas if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) 364251121Sdas return (cpackf(real_part_reciprocal(x, y), 365251404Sdas copysignf(pio2_hi + pio2_lo, y))); 366251121Sdas 367251121Sdas if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { 368251121Sdas raise_inexact(); 369251121Sdas return (z); 370251121Sdas } 371251121Sdas 372251121Sdas if (ax == 1 && ay < FLT_EPSILON) 373251404Sdas rx = (m_ln2 - logf(ay)) / 2; 374251121Sdas else 375251121Sdas rx = log1pf(4 * ax / sum_squares(ax - 1, ay)) / 4; 376251121Sdas 377251121Sdas if (ax == 1) 378251121Sdas ry = atan2f(2, -ay) / 2; 379251121Sdas else if (ay < FLT_EPSILON) 380251121Sdas ry = atan2f(2 * ay, (1 - ax) * (1 + ax)) / 2; 381251121Sdas else 382251121Sdas ry = atan2f(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; 383251121Sdas 384251121Sdas return (cpackf(copysignf(rx, x), copysignf(ry, y))); 385251121Sdas} 386251121Sdas 387251121Sdasfloat complex 388251121Sdascatanf(float complex z) 389251121Sdas{ 390251121Sdas float complex w = catanhf(cpackf(cimagf(z), crealf(z))); 391251404Sdas 392251121Sdas return (cpackf(cimagf(w), crealf(w))); 393251121Sdas} 394