| 1 | /*- |
| 2 | * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG> |
| 3 | * All rights reserved. |
| 4 | * |
| 5 | * Redistribution and use in source and binary forms, with or without |
| 6 | * modification, are permitted provided that the following conditions |
| 7 | * are met: |
| 8 | * 1. Redistributions of source code must retain the above copyright |
| 9 | * notice, this list of conditions and the following disclaimer. |
| 10 | * 2. Redistributions in binary form must reproduce the above copyright |
| 11 | * notice, this list of conditions and the following disclaimer in the |
| 12 | * documentation and/or other materials provided with the distribution. |
| 13 | * |
| 14 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND |
| 15 | * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
| 16 | * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
| 17 | * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE |
| 18 | * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL |
| 19 | * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS |
| 20 | * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) |
| 21 | * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT |
| 22 | * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY |
| 23 | * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF |
| 24 | * SUCH DAMAGE. |
| 25 | * |
| 26 | * $FreeBSD: head/lib/msun/src/catrigf.c 251404 2013-06-05 05:33:01Z das $ |
| 27 | */ |
| 28 | |
| 29 | /* |
| 30 | * The algorithm is very close to that in "Implementing the complex arcsine |
| 31 | * and arccosine functions using exception handling" by T. E. Hull, Thomas F. |
| 32 | * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on |
| 33 | * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335, |
| 34 | * http://dl.acm.org/citation.cfm?id=275324. |
| 35 | * |
| 36 | * See catrig.c for complete comments. |
| 37 | * |
| 38 | * XXX comments were removed automatically, and even short ones on the right |
| 39 | * of statements were removed (all of them), contrary to normal style. Only |
| 40 | * a few comments on the right of declarations remain. |
| 41 | */ |
| 42 | |
| 43 | #include <complex.h> |
| 44 | #include <float.h> |
| 45 | |
| 46 | #include "math.h" |
| 47 | #include "math_private.h" |
| 48 | |
| 49 | #undef isinf |
| 50 | #define isinf(x) (fabsf(x) == INFINITY) |
| 51 | #undef isnan |
| 52 | #define isnan(x) ((x) != (x)) |
| 53 | #define raise_inexact() do { volatile float junk = 1 + tiny; } while(0) |
| 54 | #undef signbit |
| 55 | #define signbit(x) (__builtin_signbitf(x)) |
| 56 | |
| 57 | static const float |
| 58 | A_crossover = 10, |
| 59 | B_crossover = 0.6417, |
| 60 | FOUR_SQRT_MIN = 0x1p-61, |
| 61 | QUARTER_SQRT_MAX = 0x1p61, |
| 62 | m_e = 2.7182818285e0, /* 0xadf854.0p-22 */ |
| 63 | m_ln2 = 6.9314718056e-1, /* 0xb17218.0p-24 */ |
| 64 | pio2_hi = 1.5707962513e0, /* 0xc90fda.0p-23 */ |
| 65 | RECIP_EPSILON = 1 / FLT_EPSILON, |
| 66 | SQRT_3_EPSILON = 5.9801995673e-4, /* 0x9cc471.0p-34 */ |
| 67 | SQRT_6_EPSILON = 8.4572793338e-4, /* 0xddb3d7.0p-34 */ |
| 68 | SQRT_MIN = 0x1p-63; |
| 69 | |
| 70 | static const volatile float |
| 71 | pio2_lo = 7.5497899549e-8, /* 0xa22169.0p-47 */ |
| 72 | tiny = 0x1p-100; |
| 73 | |
| 74 | static float complex clog_for_large_values(float complex z); |
| 75 | |
| 76 | static inline float |
| 77 | f(float a, float b, float hypot_a_b) |
| 78 | { |
| 79 | if (b < 0) |
| 80 | return ((hypot_a_b - b) / 2); |
| 81 | if (b == 0) |
| 82 | return (a / 2); |
| 83 | return (a * a / (hypot_a_b + b) / 2); |
| 84 | } |
| 85 | |
| 86 | static inline void |
| 87 | do_hard_work(float x, float y, float *rx, int *B_is_usable, float *B, |
| 88 | float *sqrt_A2my2, float *new_y) |
| 89 | { |
| 90 | float R, S, A; |
| 91 | float Am1, Amy; |
| 92 | |
| 93 | R = hypotf(x, y + 1); |
| 94 | S = hypotf(x, y - 1); |
| 95 | |
| 96 | A = (R + S) / 2; |
| 97 | if (A < 1) |
| 98 | A = 1; |
| 99 | |
| 100 | if (A < A_crossover) { |
| 101 | if (y == 1 && x < FLT_EPSILON * FLT_EPSILON / 128) { |
| 102 | *rx = sqrtf(x); |
| 103 | } else if (x >= FLT_EPSILON * fabsf(y - 1)) { |
| 104 | Am1 = f(x, 1 + y, R) + f(x, 1 - y, S); |
| 105 | *rx = log1pf(Am1 + sqrtf(Am1 * (A + 1))); |
| 106 | } else if (y < 1) { |
| 107 | *rx = x / sqrtf((1 - y) * (1 + y)); |
| 108 | } else { |
| 109 | *rx = log1pf((y - 1) + sqrtf((y - 1) * (y + 1))); |
| 110 | } |
| 111 | } else { |
| 112 | *rx = logf(A + sqrtf(A * A - 1)); |
| 113 | } |
| 114 | |
| 115 | *new_y = y; |
| 116 | |
| 117 | if (y < FOUR_SQRT_MIN) { |
| 118 | *B_is_usable = 0; |
| 119 | *sqrt_A2my2 = A * (2 / FLT_EPSILON); |
| 120 | *new_y = y * (2 / FLT_EPSILON); |
| 121 | return; |
| 122 | } |
| 123 | |
| 124 | *B = y / A; |
| 125 | *B_is_usable = 1; |
| 126 | |
| 127 | if (*B > B_crossover) { |
| 128 | *B_is_usable = 0; |
| 129 | if (y == 1 && x < FLT_EPSILON / 128) { |
| 130 | *sqrt_A2my2 = sqrtf(x) * sqrtf((A + y) / 2); |
| 131 | } else if (x >= FLT_EPSILON * fabsf(y - 1)) { |
| 132 | Amy = f(x, y + 1, R) + f(x, y - 1, S); |
| 133 | *sqrt_A2my2 = sqrtf(Amy * (A + y)); |
| 134 | } else if (y > 1) { |
| 135 | *sqrt_A2my2 = x * (4 / FLT_EPSILON / FLT_EPSILON) * y / |
| 136 | sqrtf((y + 1) * (y - 1)); |
| 137 | *new_y = y * (4 / FLT_EPSILON / FLT_EPSILON); |
| 138 | } else { |
| 139 | *sqrt_A2my2 = sqrtf((1 - y) * (1 + y)); |
| 140 | } |
| 141 | } |
| 142 | } |
| 143 | |
| 144 | float complex |
| 145 | casinhf(float complex z) |
| 146 | { |
| 147 | float x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y; |
| 148 | int B_is_usable; |
| 149 | float complex w; |
| 150 | |
| 151 | x = crealf(z); |
| 152 | y = cimagf(z); |
| 153 | ax = fabsf(x); |
| 154 | ay = fabsf(y); |
| 155 | |
| 156 | if (isnan(x) || isnan(y)) { |
| 157 | if (isinf(x)) |
| 158 | return (cpackf(x, y + y)); |
| 159 | if (isinf(y)) |
| 160 | return (cpackf(y, x + x)); |
| 161 | if (y == 0) |
| 162 | return (cpackf(x + x, y)); |
| 163 | return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); |
| 164 | } |
| 165 | |
| 166 | if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { |
| 167 | if (signbit(x) == 0) |
| 168 | w = clog_for_large_values(z) + m_ln2; |
| 169 | else |
| 170 | w = clog_for_large_values(-z) + m_ln2; |
| 171 | return (cpackf(copysignf(crealf(w), x), |
| 172 | copysignf(cimagf(w), y))); |
| 173 | } |
| 174 | |
| 175 | if (x == 0 && y == 0) |
| 176 | return (z); |
| 177 | |
| 178 | raise_inexact(); |
| 179 | |
| 180 | if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) |
| 181 | return (z); |
| 182 | |
| 183 | do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y); |
| 184 | if (B_is_usable) |
| 185 | ry = asinf(B); |
| 186 | else |
| 187 | ry = atan2f(new_y, sqrt_A2my2); |
| 188 | return (cpackf(copysignf(rx, x), copysignf(ry, y))); |
| 189 | } |
| 190 | |
| 191 | float complex |
| 192 | casinf(float complex z) |
| 193 | { |
| 194 | float complex w = casinhf(cpackf(cimagf(z), crealf(z))); |
| 195 | |
| 196 | return (cpackf(cimagf(w), crealf(w))); |
| 197 | } |
| 198 | |
| 199 | float complex |
| 200 | cacosf(float complex z) |
| 201 | { |
| 202 | float x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x; |
| 203 | int sx, sy; |
| 204 | int B_is_usable; |
| 205 | float complex w; |
| 206 | |
| 207 | x = crealf(z); |
| 208 | y = cimagf(z); |
| 209 | sx = signbit(x); |
| 210 | sy = signbit(y); |
| 211 | ax = fabsf(x); |
| 212 | ay = fabsf(y); |
| 213 | |
| 214 | if (isnan(x) || isnan(y)) { |
| 215 | if (isinf(x)) |
| 216 | return (cpackf(y + y, -INFINITY)); |
| 217 | if (isinf(y)) |
| 218 | return (cpackf(x + x, -y)); |
| 219 | if (x == 0) |
| 220 | return (cpackf(pio2_hi + pio2_lo, y + y)); |
| 221 | return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); |
| 222 | } |
| 223 | |
| 224 | if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { |
| 225 | w = clog_for_large_values(z); |
| 226 | rx = fabsf(cimagf(w)); |
| 227 | ry = crealf(w) + m_ln2; |
| 228 | if (sy == 0) |
| 229 | ry = -ry; |
| 230 | return (cpackf(rx, ry)); |
| 231 | } |
| 232 | |
| 233 | if (x == 1 && y == 0) |
| 234 | return (cpackf(0, -y)); |
| 235 | |
| 236 | raise_inexact(); |
| 237 | |
| 238 | if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) |
| 239 | return (cpackf(pio2_hi - (x - pio2_lo), -y)); |
| 240 | |
| 241 | do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x); |
| 242 | if (B_is_usable) { |
| 243 | if (sx == 0) |
| 244 | rx = acosf(B); |
| 245 | else |
| 246 | rx = acosf(-B); |
| 247 | } else { |
| 248 | if (sx == 0) |
| 249 | rx = atan2f(sqrt_A2mx2, new_x); |
| 250 | else |
| 251 | rx = atan2f(sqrt_A2mx2, -new_x); |
| 252 | } |
| 253 | if (sy == 0) |
| 254 | ry = -ry; |
| 255 | return (cpackf(rx, ry)); |
| 256 | } |
| 257 | |
| 258 | float complex |
| 259 | cacoshf(float complex z) |
| 260 | { |
| 261 | float complex w; |
| 262 | float rx, ry; |
| 263 | |
| 264 | w = cacosf(z); |
| 265 | rx = crealf(w); |
| 266 | ry = cimagf(w); |
| 267 | if (isnan(rx) && isnan(ry)) |
| 268 | return (cpackf(ry, rx)); |
| 269 | if (isnan(rx)) |
| 270 | return (cpackf(fabsf(ry), rx)); |
| 271 | if (isnan(ry)) |
| 272 | return (cpackf(ry, ry)); |
| 273 | return (cpackf(fabsf(ry), copysignf(rx, cimagf(z)))); |
| 274 | } |
| 275 | |
| 276 | static float complex |
| 277 | clog_for_large_values(float complex z) |
| 278 | { |
| 279 | float x, y; |
| 280 | float ax, ay, t; |
| 281 | |
| 282 | x = crealf(z); |
| 283 | y = cimagf(z); |
| 284 | ax = fabsf(x); |
| 285 | ay = fabsf(y); |
| 286 | if (ax < ay) { |
| 287 | t = ax; |
| 288 | ax = ay; |
| 289 | ay = t; |
| 290 | } |
| 291 | |
| 292 | if (ax > FLT_MAX / 2) |
| 293 | return (cpackf(logf(hypotf(x / m_e, y / m_e)) + 1, |
| 294 | atan2f(y, x))); |
| 295 | |
| 296 | if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) |
| 297 | return (cpackf(logf(hypotf(x, y)), atan2f(y, x))); |
| 298 | |
| 299 | return (cpackf(logf(ax * ax + ay * ay) / 2, atan2f(y, x))); |
| 300 | } |
| 301 | |
| 302 | static inline float |
| 303 | sum_squares(float x, float y) |
| 304 | { |
| 305 | |
| 306 | if (y < SQRT_MIN) |
| 307 | return (x * x); |
| 308 | |
| 309 | return (x * x + y * y); |
| 310 | } |
| 311 | |
| 312 | static inline float |
| 313 | real_part_reciprocal(float x, float y) |
| 314 | { |
| 315 | float scale; |
| 316 | uint32_t hx, hy; |
| 317 | int32_t ix, iy; |
| 318 | |
| 319 | GET_FLOAT_WORD(hx, x); |
| 320 | ix = hx & 0x7f800000; |
| 321 | GET_FLOAT_WORD(hy, y); |
| 322 | iy = hy & 0x7f800000; |
| 323 | #define BIAS (FLT_MAX_EXP - 1) |
| 324 | #define CUTOFF (FLT_MANT_DIG / 2 + 1) |
| 325 | if (ix - iy >= CUTOFF << 23 || isinf(x)) |
| 326 | return (1 / x); |
| 327 | if (iy - ix >= CUTOFF << 23) |
| 328 | return (x / y / y); |
| 329 | if (ix <= (BIAS + FLT_MAX_EXP / 2 - CUTOFF) << 23) |
| 330 | return (x / (x * x + y * y)); |
| 331 | SET_FLOAT_WORD(scale, 0x7f800000 - ix); |
| 332 | x *= scale; |
| 333 | y *= scale; |
| 334 | return (x / (x * x + y * y) * scale); |
| 335 | } |
| 336 | |
| 337 | float complex |
| 338 | catanhf(float complex z) |
| 339 | { |
| 340 | float x, y, ax, ay, rx, ry; |
| 341 | |
| 342 | x = crealf(z); |
| 343 | y = cimagf(z); |
| 344 | ax = fabsf(x); |
| 345 | ay = fabsf(y); |
| 346 | |
| 347 | if (y == 0 && ax <= 1) |
| 348 | return (cpackf(atanhf(x), y)); |
| 349 | |
| 350 | if (x == 0) |
| 351 | return (cpackf(x, atanf(y))); |
| 352 | |
| 353 | if (isnan(x) || isnan(y)) { |
| 354 | if (isinf(x)) |
| 355 | return (cpackf(copysignf(0, x), y + y)); |
| 356 | if (isinf(y)) |
| 357 | return (cpackf(copysignf(0, x), |
| 358 | copysignf(pio2_hi + pio2_lo, y))); |
| 359 | return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); |
| 360 | } |
| 361 | |
| 362 | if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) |
| 363 | return (cpackf(real_part_reciprocal(x, y), |
| 364 | copysignf(pio2_hi + pio2_lo, y))); |
| 365 | |
| 366 | if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { |
| 367 | raise_inexact(); |
| 368 | return (z); |
| 369 | } |
| 370 | |
| 371 | if (ax == 1 && ay < FLT_EPSILON) |
| 372 | rx = (m_ln2 - logf(ay)) / 2; |
| 373 | else |
| 374 | rx = log1pf(4 * ax / sum_squares(ax - 1, ay)) / 4; |
| 375 | |
| 376 | if (ax == 1) |
| 377 | ry = atan2f(2, -ay) / 2; |
| 378 | else if (ay < FLT_EPSILON) |
| 379 | ry = atan2f(2 * ay, (1 - ax) * (1 + ax)) / 2; |
| 380 | else |
| 381 | ry = atan2f(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; |
| 382 | |
| 383 | return (cpackf(copysignf(rx, x), copysignf(ry, y))); |
| 384 | } |
| 385 | |
| 386 | float complex |
| 387 | catanf(float complex z) |
| 388 | { |
| 389 | float complex w = catanhf(cpackf(cimagf(z), crealf(z))); |
| 390 | |
| 391 | return (cpackf(cimagf(w), crealf(w))); |
| 392 | } |