catrigf.c revision 275819
155714Skris/*- 255714Skris * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG> 355714Skris * All rights reserved. 455714Skris * 555714Skris * Redistribution and use in source and binary forms, with or without 655714Skris * modification, are permitted provided that the following conditions 755714Skris * are met: 8296341Sdelphij * 1. Redistributions of source code must retain the above copyright 955714Skris * notice, this list of conditions and the following disclaimer. 1055714Skris * 2. Redistributions in binary form must reproduce the above copyright 1155714Skris * notice, this list of conditions and the following disclaimer in the 1255714Skris * documentation and/or other materials provided with the distribution. 1355714Skris * 1455714Skris * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 15296341Sdelphij * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 1655714Skris * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 1755714Skris * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 1855714Skris * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 1955714Skris * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 2055714Skris * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 2155714Skris * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 22296341Sdelphij * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 2355714Skris * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 2455714Skris * SUCH DAMAGE. 2555714Skris */ 2655714Skris 2755714Skris/* 2855714Skris * The algorithm is very close to that in "Implementing the complex arcsine 2955714Skris * and arccosine functions using exception handling" by T. E. Hull, Thomas F. 3055714Skris * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on 3155714Skris * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335, 3255714Skris * http://dl.acm.org/citation.cfm?id=275324. 3355714Skris * 3455714Skris * See catrig.c for complete comments. 3555714Skris * 3655714Skris * XXX comments were removed automatically, and even short ones on the right 37296341Sdelphij * of statements were removed (all of them), contrary to normal style. Only 3855714Skris * a few comments on the right of declarations remain. 3955714Skris */ 40296341Sdelphij 4155714Skris#include <sys/cdefs.h> 4255714Skris__FBSDID("$FreeBSD: head/lib/msun/src/catrigf.c 275819 2014-12-16 09:21:56Z ed $"); 4355714Skris 4455714Skris#include <complex.h> 4555714Skris#include <float.h> 4655714Skris 4755714Skris#include "math.h" 4855714Skris#include "math_private.h" 4955714Skris 5055714Skris#undef isinf 5155714Skris#define isinf(x) (fabsf(x) == INFINITY) 52296341Sdelphij#undef isnan 5355714Skris#define isnan(x) ((x) != (x)) 5455714Skris#define raise_inexact() do { volatile float junk = 1 + tiny; } while(0) 5555714Skris#undef signbit 5655714Skris#define signbit(x) (__builtin_signbitf(x)) 5755714Skris 58160814Ssimonstatic const float 59238405SjkimA_crossover = 10, 60238405SjkimB_crossover = 0.6417, 61238405SjkimFOUR_SQRT_MIN = 0x1p-61, 62238405SjkimQUARTER_SQRT_MAX = 0x1p61, 63238405Sjkimm_e = 2.7182818285e0, /* 0xadf854.0p-22 */ 64238405Sjkimm_ln2 = 6.9314718056e-1, /* 0xb17218.0p-24 */ 65238405Sjkimpio2_hi = 1.5707962513e0, /* 0xc90fda.0p-23 */ 66296341SdelphijRECIP_EPSILON = 1 / FLT_EPSILON, 67238405SjkimSQRT_3_EPSILON = 5.9801995673e-4, /* 0x9cc471.0p-34 */ 68238405SjkimSQRT_6_EPSILON = 8.4572793338e-4, /* 0xddb3d7.0p-34 */ 69238405SjkimSQRT_MIN = 0x1p-63; 70238405Sjkim 71238405Sjkimstatic const volatile float 72238405Sjkimpio2_lo = 7.5497899549e-8, /* 0xa22169.0p-47 */ 73238405Sjkimtiny = 0x1p-100; 74238405Sjkim 75238405Sjkimstatic float complex clog_for_large_values(float complex z); 76238405Sjkim 77238405Sjkimstatic inline float 78238405Sjkimf(float a, float b, float hypot_a_b) 79238405Sjkim{ 80238405Sjkim if (b < 0) 81238405Sjkim return ((hypot_a_b - b) / 2); 82238405Sjkim if (b == 0) 83238405Sjkim return (a / 2); 84238405Sjkim return (a * a / (hypot_a_b + b) / 2); 85238405Sjkim} 86238405Sjkim 87238405Sjkimstatic inline void 88238405Sjkimdo_hard_work(float x, float y, float *rx, int *B_is_usable, float *B, 89238405Sjkim float *sqrt_A2my2, float *new_y) 90238405Sjkim{ 91238405Sjkim float R, S, A; 92238405Sjkim float Am1, Amy; 93238405Sjkim 94238405Sjkim R = hypotf(x, y + 1); 95238405Sjkim S = hypotf(x, y - 1); 96238405Sjkim 97238405Sjkim A = (R + S) / 2; 98238405Sjkim if (A < 1) 99238405Sjkim A = 1; 100238405Sjkim 101238405Sjkim if (A < A_crossover) { 102238405Sjkim if (y == 1 && x < FLT_EPSILON * FLT_EPSILON / 128) { 103238405Sjkim *rx = sqrtf(x); 104238405Sjkim } else if (x >= FLT_EPSILON * fabsf(y - 1)) { 105238405Sjkim Am1 = f(x, 1 + y, R) + f(x, 1 - y, S); 106238405Sjkim *rx = log1pf(Am1 + sqrtf(Am1 * (A + 1))); 107238405Sjkim } else if (y < 1) { 108238405Sjkim *rx = x / sqrtf((1 - y) * (1 + y)); 109238405Sjkim } else { 110238405Sjkim *rx = log1pf((y - 1) + sqrtf((y - 1) * (y + 1))); 111238405Sjkim } 112160814Ssimon } else { 113160814Ssimon *rx = logf(A + sqrtf(A * A - 1)); 114296341Sdelphij } 115160814Ssimon 116160814Ssimon *new_y = y; 117160814Ssimon 118160814Ssimon if (y < FOUR_SQRT_MIN) { 119160814Ssimon *B_is_usable = 0; 120160814Ssimon *sqrt_A2my2 = A * (2 / FLT_EPSILON); 121160814Ssimon *new_y = y * (2 / FLT_EPSILON); 122160814Ssimon return; 123160814Ssimon } 124238405Sjkim 125238405Sjkim *B = y / A; 126238405Sjkim *B_is_usable = 1; 127238405Sjkim 128238405Sjkim if (*B > B_crossover) { 129238405Sjkim *B_is_usable = 0; 130238405Sjkim if (y == 1 && x < FLT_EPSILON / 128) { 131238405Sjkim *sqrt_A2my2 = sqrtf(x) * sqrtf((A + y) / 2); 132238405Sjkim } else if (x >= FLT_EPSILON * fabsf(y - 1)) { 133238405Sjkim Amy = f(x, y + 1, R) + f(x, y - 1, S); 134238405Sjkim *sqrt_A2my2 = sqrtf(Amy * (A + y)); 135238405Sjkim } else if (y > 1) { 136238405Sjkim *sqrt_A2my2 = x * (4 / FLT_EPSILON / FLT_EPSILON) * y / 137238405Sjkim sqrtf((y + 1) * (y - 1)); 138238405Sjkim *new_y = y * (4 / FLT_EPSILON / FLT_EPSILON); 139238405Sjkim } else { 140238405Sjkim *sqrt_A2my2 = sqrtf((1 - y) * (1 + y)); 141238405Sjkim } 142238405Sjkim } 143238405Sjkim} 144238405Sjkim 145238405Sjkimfloat complex 146238405Sjkimcasinhf(float complex z) 147238405Sjkim{ 148238405Sjkim float x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y; 149238405Sjkim int B_is_usable; 15055714Skris float complex w; 151296341Sdelphij 152296341Sdelphij x = crealf(z); 15355714Skris y = cimagf(z); 154296341Sdelphij ax = fabsf(x); 15555714Skris ay = fabsf(y); 15655714Skris 15755714Skris if (isnan(x) || isnan(y)) { 15855714Skris if (isinf(x)) 15955714Skris return (CMPLXF(x, y + y)); 160296341Sdelphij if (isinf(y)) 16155714Skris return (CMPLXF(y, x + x)); 162296341Sdelphij if (y == 0) 163296341Sdelphij return (CMPLXF(x + x, y)); 164296341Sdelphij return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); 165296341Sdelphij } 166238405Sjkim 167296341Sdelphij if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { 168296341Sdelphij if (signbit(x) == 0) 169273399Sdelphij w = clog_for_large_values(z) + m_ln2; 170296341Sdelphij else 171296341Sdelphij w = clog_for_large_values(-z) + m_ln2; 172238405Sjkim return (CMPLXF(copysignf(crealf(w), x), 173296341Sdelphij copysignf(cimagf(w), y))); 174296341Sdelphij } 17555714Skris 176296341Sdelphij if (x == 0 && y == 0) 177296341Sdelphij return (z); 178238405Sjkim 179296341Sdelphij raise_inexact(); 180296341Sdelphij 181238405Sjkim if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) 182296341Sdelphij return (z); 183296341Sdelphij 184296341Sdelphij do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y); 185296341Sdelphij if (B_is_usable) 186296341Sdelphij ry = asinf(B); 187296341Sdelphij else 188296341Sdelphij ry = atan2f(new_y, sqrt_A2my2); 189296341Sdelphij return (CMPLXF(copysignf(rx, x), copysignf(ry, y))); 190296341Sdelphij} 191296341Sdelphij 192296341Sdelphijfloat complex 193296341Sdelphijcasinf(float complex z) 194296341Sdelphij{ 195194206Ssimon float complex w = casinhf(CMPLXF(cimagf(z), crealf(z))); 196296341Sdelphij 197296341Sdelphij return (CMPLXF(cimagf(w), crealf(w))); 198296341Sdelphij} 199296341Sdelphij 200296341Sdelphijfloat complex 201296341Sdelphijcacosf(float complex z) 20255714Skris{ 203238405Sjkim float x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x; 204296341Sdelphij int sx, sy; 205296341Sdelphij int B_is_usable; 206296341Sdelphij float complex w; 207296341Sdelphij 208296341Sdelphij x = crealf(z); 209296341Sdelphij y = cimagf(z); 210238405Sjkim sx = signbit(x); 211296341Sdelphij sy = signbit(y); 212238405Sjkim ax = fabsf(x); 213238405Sjkim ay = fabsf(y); 214296341Sdelphij 215296341Sdelphij if (isnan(x) || isnan(y)) { 216238405Sjkim if (isinf(x)) 217238405Sjkim return (CMPLXF(y + y, -INFINITY)); 218296341Sdelphij if (isinf(y)) 219238405Sjkim return (CMPLXF(x + x, -y)); 220238405Sjkim if (x == 0) 221296341Sdelphij return (CMPLXF(pio2_hi + pio2_lo, y + y)); 222296341Sdelphij return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); 223238405Sjkim } 224238405Sjkim 225296341Sdelphij if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { 226238405Sjkim w = clog_for_large_values(z); 227238405Sjkim rx = fabsf(cimagf(w)); 228296341Sdelphij ry = crealf(w) + m_ln2; 229238405Sjkim if (sy == 0) 230238405Sjkim ry = -ry; 231296341Sdelphij return (CMPLXF(rx, ry)); 232238405Sjkim } 233238405Sjkim 234296341Sdelphij if (x == 1 && y == 0) 235238405Sjkim return (CMPLXF(0, -y)); 236296341Sdelphij 237296341Sdelphij raise_inexact(); 238264331Sjkim 239264331Sjkim if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) 240264331Sjkim return (CMPLXF(pio2_hi - (x - pio2_lo), -y)); 241296341Sdelphij 242264331Sjkim do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x); 243238405Sjkim if (B_is_usable) { 244296341Sdelphij if (sx == 0) 245194206Ssimon rx = acosf(B); 246238405Sjkim else 247296341Sdelphij rx = acosf(-B); 248296341Sdelphij } else { 249296341Sdelphij if (sx == 0) 250296341Sdelphij rx = atan2f(sqrt_A2mx2, new_x); 251296341Sdelphij else 252296341Sdelphij rx = atan2f(sqrt_A2mx2, -new_x); 253296341Sdelphij } 254296341Sdelphij if (sy == 0) 255238405Sjkim ry = -ry; 256205128Ssimon return (CMPLXF(rx, ry)); 257296341Sdelphij} 258205128Ssimon 259296341Sdelphijfloat complex 260238405Sjkimcacoshf(float complex z) 261296341Sdelphij{ 262296341Sdelphij float complex w; 263238405Sjkim float rx, ry; 264194206Ssimon 265296341Sdelphij w = cacosf(z); 266194206Ssimon rx = crealf(w); 267296341Sdelphij ry = cimagf(w); 268194206Ssimon if (isnan(rx) && isnan(ry)) 269238405Sjkim return (CMPLXF(ry, rx)); 270296341Sdelphij if (isnan(rx)) 271296341Sdelphij return (CMPLXF(fabsf(ry), rx)); 272296341Sdelphij if (isnan(ry)) 273296341Sdelphij return (CMPLXF(ry, ry)); 274296341Sdelphij return (CMPLXF(fabsf(ry), copysignf(rx, cimagf(z)))); 275238405Sjkim} 276238405Sjkim 277238405Sjkimstatic float complex 278296341Sdelphijclog_for_large_values(float complex z) 279296341Sdelphij{ 280296341Sdelphij float x, y; 281296341Sdelphij float ax, ay, t; 282238405Sjkim 283296341Sdelphij x = crealf(z); 284296341Sdelphij y = cimagf(z); 285296341Sdelphij ax = fabsf(x); 286296341Sdelphij ay = fabsf(y); 287296341Sdelphij if (ax < ay) { 288296341Sdelphij t = ax; 289296341Sdelphij ax = ay; 290238405Sjkim ay = t; 291296341Sdelphij } 292194206Ssimon 293296341Sdelphij if (ax > FLT_MAX / 2) 294194206Ssimon return (CMPLXF(logf(hypotf(x / m_e, y / m_e)) + 1, 295238405Sjkim atan2f(y, x))); 296238405Sjkim 297296341Sdelphij if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) 298296341Sdelphij return (CMPLXF(logf(hypotf(x, y)), atan2f(y, x))); 299238405Sjkim 300238405Sjkim return (CMPLXF(logf(ax * ax + ay * ay) / 2, atan2f(y, x))); 301296341Sdelphij} 302296341Sdelphij 303238405Sjkimstatic inline float 304238405Sjkimsum_squares(float x, float y) 305296341Sdelphij{ 306296341Sdelphij 307296341Sdelphij if (y < SQRT_MIN) 308194206Ssimon return (x * x); 309296341Sdelphij 310194206Ssimon return (x * x + y * y); 311194206Ssimon} 312296341Sdelphij 313194206Ssimonstatic inline float 314194206Ssimonreal_part_reciprocal(float x, float y) 315296341Sdelphij{ 316194206Ssimon float scale; 317194206Ssimon uint32_t hx, hy; 318296341Sdelphij int32_t ix, iy; 319194206Ssimon 320194206Ssimon GET_FLOAT_WORD(hx, x); 321296341Sdelphij ix = hx & 0x7f800000; 322194206Ssimon GET_FLOAT_WORD(hy, y); 323194206Ssimon iy = hy & 0x7f800000; 324296341Sdelphij#define BIAS (FLT_MAX_EXP - 1) 325194206Ssimon#define CUTOFF (FLT_MANT_DIG / 2 + 1) 326194206Ssimon if (ix - iy >= CUTOFF << 23 || isinf(x)) 327296341Sdelphij return (1 / x); 328194206Ssimon if (iy - ix >= CUTOFF << 23) 329194206Ssimon return (x / y / y); 330296341Sdelphij if (ix <= (BIAS + FLT_MAX_EXP / 2 - CUTOFF) << 23) 331194206Ssimon return (x / (x * x + y * y)); 332194206Ssimon SET_FLOAT_WORD(scale, 0x7f800000 - ix); 333296341Sdelphij x *= scale; 334194206Ssimon y *= scale; 335194206Ssimon return (x / (x * x + y * y) * scale); 336296341Sdelphij} 337194206Ssimon 338194206Ssimonfloat complex 339296341Sdelphijcatanhf(float complex z) 340194206Ssimon{ 341194206Ssimon float x, y, ax, ay, rx, ry; 342296341Sdelphij 343296341Sdelphij x = crealf(z); 344296341Sdelphij y = cimagf(z); 345296341Sdelphij ax = fabsf(x); 346194206Ssimon ay = fabsf(y); 347296341Sdelphij 348194206Ssimon if (y == 0 && ax <= 1) 349194206Ssimon return (CMPLXF(atanhf(x), y)); 350296341Sdelphij 351296341Sdelphij if (x == 0) 352296341Sdelphij return (CMPLXF(x, atanf(y))); 353296341Sdelphij 354194206Ssimon if (isnan(x) || isnan(y)) { 355296341Sdelphij if (isinf(x)) 356194206Ssimon return (CMPLXF(copysignf(0, x), y + y)); 357194206Ssimon if (isinf(y)) 358296341Sdelphij return (CMPLXF(copysignf(0, x), 359194206Ssimon copysignf(pio2_hi + pio2_lo, y))); 360194206Ssimon return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); 361296341Sdelphij } 362238405Sjkim 363296341Sdelphij if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) 364238405Sjkim return (CMPLXF(real_part_reciprocal(x, y), 365296341Sdelphij copysignf(pio2_hi + pio2_lo, y))); 366238405Sjkim 367238405Sjkim if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { 368296341Sdelphij raise_inexact(); 369194206Ssimon return (z); 370194206Ssimon } 371296341Sdelphij 372296341Sdelphij if (ax == 1 && ay < FLT_EPSILON) 373296341Sdelphij rx = (m_ln2 - logf(ay)) / 2; 374296341Sdelphij else 375238405Sjkim rx = log1pf(4 * ax / sum_squares(ax - 1, ay)) / 4; 376296341Sdelphij 377238405Sjkim if (ax == 1) 378296341Sdelphij ry = atan2f(2, -ay) / 2; 379238405Sjkim else if (ay < FLT_EPSILON) 380296341Sdelphij ry = atan2f(2 * ay, (1 - ax) * (1 + ax)) / 2; 381296341Sdelphij else 382194206Ssimon ry = atan2f(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; 383238405Sjkim 384296341Sdelphij return (CMPLXF(copysignf(rx, x), copysignf(ry, y))); 385296341Sdelphij} 386296341Sdelphij 387296341Sdelphijfloat complex 388238405Sjkimcatanf(float complex z) 389296341Sdelphij{ 390296341Sdelphij float complex w = catanhf(CMPLXF(cimagf(z), crealf(z))); 391296341Sdelphij 392296341Sdelphij return (CMPLXF(cimagf(w), crealf(w))); 393296341Sdelphij} 394296341Sdelphij