| catrig.c (251404) | catrig.c (275819) |
|---|---|
| 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 --- 11 unchanged lines hidden (view full) --- 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 27#include <sys/cdefs.h> | 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 --- 11 unchanged lines hidden (view full) --- 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 27#include <sys/cdefs.h> |
| 28__FBSDID("$FreeBSD: head/lib/msun/src/catrig.c 251404 2013-06-05 05:33:01Z das $"); | 28__FBSDID("$FreeBSD: head/lib/msun/src/catrig.c 275819 2014-12-16 09:21:56Z ed $"); |
| 29 30#include <complex.h> 31#include <float.h> 32 33#include "math.h" 34#include "math_private.h" 35 36#undef isinf --- 244 unchanged lines hidden (view full) --- 281 x = creal(z); 282 y = cimag(z); 283 ax = fabs(x); 284 ay = fabs(y); 285 286 if (isnan(x) || isnan(y)) { 287 /* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */ 288 if (isinf(x)) | 29 30#include <complex.h> 31#include <float.h> 32 33#include "math.h" 34#include "math_private.h" 35 36#undef isinf --- 244 unchanged lines hidden (view full) --- 281 x = creal(z); 282 y = cimag(z); 283 ax = fabs(x); 284 ay = fabs(y); 285 286 if (isnan(x) || isnan(y)) { 287 /* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */ 288 if (isinf(x)) |
| 289 return (cpack(x, y + y)); | 289 return (CMPLX(x, y + y)); |
| 290 /* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */ 291 if (isinf(y)) | 290 /* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */ 291 if (isinf(y)) |
| 292 return (cpack(y, x + x)); | 292 return (CMPLX(y, x + x)); |
| 293 /* casinh(NaN + I*0) = NaN + I*0 */ 294 if (y == 0) | 293 /* casinh(NaN + I*0) = NaN + I*0 */ 294 if (y == 0) |
| 295 return (cpack(x + x, y)); | 295 return (CMPLX(x + x, y)); |
| 296 /* 297 * All other cases involving NaN return NaN + I*NaN. 298 * C99 leaves it optional whether to raise invalid if one of 299 * the arguments is not NaN, so we opt not to raise it. 300 */ | 296 /* 297 * All other cases involving NaN return NaN + I*NaN. 298 * C99 leaves it optional whether to raise invalid if one of 299 * the arguments is not NaN, so we opt not to raise it. 300 */ |
| 301 return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); | 301 return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); |
| 302 } 303 304 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { 305 /* clog...() will raise inexact unless x or y is infinite. */ 306 if (signbit(x) == 0) 307 w = clog_for_large_values(z) + m_ln2; 308 else 309 w = clog_for_large_values(-z) + m_ln2; | 302 } 303 304 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { 305 /* clog...() will raise inexact unless x or y is infinite. */ 306 if (signbit(x) == 0) 307 w = clog_for_large_values(z) + m_ln2; 308 else 309 w = clog_for_large_values(-z) + m_ln2; |
| 310 return (cpack(copysign(creal(w), x), copysign(cimag(w), y))); | 310 return (CMPLX(copysign(creal(w), x), copysign(cimag(w), y))); |
| 311 } 312 313 /* Avoid spuriously raising inexact for z = 0. */ 314 if (x == 0 && y == 0) 315 return (z); 316 317 /* All remaining cases are inexact. */ 318 raise_inexact(); 319 320 if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) 321 return (z); 322 323 do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y); 324 if (B_is_usable) 325 ry = asin(B); 326 else 327 ry = atan2(new_y, sqrt_A2my2); | 311 } 312 313 /* Avoid spuriously raising inexact for z = 0. */ 314 if (x == 0 && y == 0) 315 return (z); 316 317 /* All remaining cases are inexact. */ 318 raise_inexact(); 319 320 if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) 321 return (z); 322 323 do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y); 324 if (B_is_usable) 325 ry = asin(B); 326 else 327 ry = atan2(new_y, sqrt_A2my2); |
| 328 return (cpack(copysign(rx, x), copysign(ry, y))); | 328 return (CMPLX(copysign(rx, x), copysign(ry, y))); |
| 329} 330 331/* 332 * casin(z) = reverse(casinh(reverse(z))) 333 * where reverse(x + I*y) = y + I*x = I*conj(z). 334 */ 335double complex 336casin(double complex z) 337{ | 329} 330 331/* 332 * casin(z) = reverse(casinh(reverse(z))) 333 * where reverse(x + I*y) = y + I*x = I*conj(z). 334 */ 335double complex 336casin(double complex z) 337{ |
| 338 double complex w = casinh(cpack(cimag(z), creal(z))); | 338 double complex w = casinh(CMPLX(cimag(z), creal(z))); |
| 339 | 339 |
| 340 return (cpack(cimag(w), creal(w))); | 340 return (CMPLX(cimag(w), creal(w))); |
| 341} 342 343/* 344 * cacos(z) = PI/2 - casin(z) 345 * but do the computation carefully so cacos(z) is accurate when z is 346 * close to 1. 347 * 348 * cacos(z) = PI/2 - z + O(z^3) as z -> 0 --- 16 unchanged lines hidden (view full) --- 365 sx = signbit(x); 366 sy = signbit(y); 367 ax = fabs(x); 368 ay = fabs(y); 369 370 if (isnan(x) || isnan(y)) { 371 /* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */ 372 if (isinf(x)) | 341} 342 343/* 344 * cacos(z) = PI/2 - casin(z) 345 * but do the computation carefully so cacos(z) is accurate when z is 346 * close to 1. 347 * 348 * cacos(z) = PI/2 - z + O(z^3) as z -> 0 --- 16 unchanged lines hidden (view full) --- 365 sx = signbit(x); 366 sy = signbit(y); 367 ax = fabs(x); 368 ay = fabs(y); 369 370 if (isnan(x) || isnan(y)) { 371 /* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */ 372 if (isinf(x)) |
| 373 return (cpack(y + y, -INFINITY)); | 373 return (CMPLX(y + y, -INFINITY)); |
| 374 /* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */ 375 if (isinf(y)) | 374 /* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */ 375 if (isinf(y)) |
| 376 return (cpack(x + x, -y)); | 376 return (CMPLX(x + x, -y)); |
| 377 /* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */ 378 if (x == 0) | 377 /* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */ 378 if (x == 0) |
| 379 return (cpack(pio2_hi + pio2_lo, y + y)); | 379 return (CMPLX(pio2_hi + pio2_lo, y + y)); |
| 380 /* 381 * All other cases involving NaN return NaN + I*NaN. 382 * C99 leaves it optional whether to raise invalid if one of 383 * the arguments is not NaN, so we opt not to raise it. 384 */ | 380 /* 381 * All other cases involving NaN return NaN + I*NaN. 382 * C99 leaves it optional whether to raise invalid if one of 383 * the arguments is not NaN, so we opt not to raise it. 384 */ |
| 385 return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); | 385 return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); |
| 386 } 387 388 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { 389 /* clog...() will raise inexact unless x or y is infinite. */ 390 w = clog_for_large_values(z); 391 rx = fabs(cimag(w)); 392 ry = creal(w) + m_ln2; 393 if (sy == 0) 394 ry = -ry; | 386 } 387 388 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { 389 /* clog...() will raise inexact unless x or y is infinite. */ 390 w = clog_for_large_values(z); 391 rx = fabs(cimag(w)); 392 ry = creal(w) + m_ln2; 393 if (sy == 0) 394 ry = -ry; |
| 395 return (cpack(rx, ry)); | 395 return (CMPLX(rx, ry)); |
| 396 } 397 398 /* Avoid spuriously raising inexact for z = 1. */ 399 if (x == 1 && y == 0) | 396 } 397 398 /* Avoid spuriously raising inexact for z = 1. */ 399 if (x == 1 && y == 0) |
| 400 return (cpack(0, -y)); | 400 return (CMPLX(0, -y)); |
| 401 402 /* All remaining cases are inexact. */ 403 raise_inexact(); 404 405 if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) | 401 402 /* All remaining cases are inexact. */ 403 raise_inexact(); 404 405 if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) |
| 406 return (cpack(pio2_hi - (x - pio2_lo), -y)); | 406 return (CMPLX(pio2_hi - (x - pio2_lo), -y)); |
| 407 408 do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x); 409 if (B_is_usable) { 410 if (sx == 0) 411 rx = acos(B); 412 else 413 rx = acos(-B); 414 } else { 415 if (sx == 0) 416 rx = atan2(sqrt_A2mx2, new_x); 417 else 418 rx = atan2(sqrt_A2mx2, -new_x); 419 } 420 if (sy == 0) 421 ry = -ry; | 407 408 do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x); 409 if (B_is_usable) { 410 if (sx == 0) 411 rx = acos(B); 412 else 413 rx = acos(-B); 414 } else { 415 if (sx == 0) 416 rx = atan2(sqrt_A2mx2, new_x); 417 else 418 rx = atan2(sqrt_A2mx2, -new_x); 419 } 420 if (sy == 0) 421 ry = -ry; |
| 422 return (cpack(rx, ry)); | 422 return (CMPLX(rx, ry)); |
| 423} 424 425/* 426 * cacosh(z) = I*cacos(z) or -I*cacos(z) 427 * where the sign is chosen so Re(cacosh(z)) >= 0. 428 */ 429double complex 430cacosh(double complex z) 431{ 432 double complex w; 433 double rx, ry; 434 435 w = cacos(z); 436 rx = creal(w); 437 ry = cimag(w); 438 /* cacosh(NaN + I*NaN) = NaN + I*NaN */ 439 if (isnan(rx) && isnan(ry)) | 423} 424 425/* 426 * cacosh(z) = I*cacos(z) or -I*cacos(z) 427 * where the sign is chosen so Re(cacosh(z)) >= 0. 428 */ 429double complex 430cacosh(double complex z) 431{ 432 double complex w; 433 double rx, ry; 434 435 w = cacos(z); 436 rx = creal(w); 437 ry = cimag(w); 438 /* cacosh(NaN + I*NaN) = NaN + I*NaN */ 439 if (isnan(rx) && isnan(ry)) |
| 440 return (cpack(ry, rx)); | 440 return (CMPLX(ry, rx)); |
| 441 /* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */ 442 /* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */ 443 if (isnan(rx)) | 441 /* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */ 442 /* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */ 443 if (isnan(rx)) |
| 444 return (cpack(fabs(ry), rx)); | 444 return (CMPLX(fabs(ry), rx)); |
| 445 /* cacosh(0 + I*NaN) = NaN + I*NaN */ 446 if (isnan(ry)) | 445 /* cacosh(0 + I*NaN) = NaN + I*NaN */ 446 if (isnan(ry)) |
| 447 return (cpack(ry, ry)); 448 return (cpack(fabs(ry), copysign(rx, cimag(z)))); | 447 return (CMPLX(ry, ry)); 448 return (CMPLX(fabs(ry), copysign(rx, cimag(z)))); |
| 449} 450 451/* 452 * Optimized version of clog() for |z| finite and larger than ~RECIP_EPSILON. 453 */ 454static double complex 455clog_for_large_values(double complex z) 456{ --- 13 unchanged lines hidden (view full) --- 470 /* 471 * Avoid overflow in hypot() when x and y are both very large. 472 * Divide x and y by E, and then add 1 to the logarithm. This depends 473 * on E being larger than sqrt(2). 474 * Dividing by E causes an insignificant loss of accuracy; however 475 * this method is still poor since it is uneccessarily slow. 476 */ 477 if (ax > DBL_MAX / 2) | 449} 450 451/* 452 * Optimized version of clog() for |z| finite and larger than ~RECIP_EPSILON. 453 */ 454static double complex 455clog_for_large_values(double complex z) 456{ --- 13 unchanged lines hidden (view full) --- 470 /* 471 * Avoid overflow in hypot() when x and y are both very large. 472 * Divide x and y by E, and then add 1 to the logarithm. This depends 473 * on E being larger than sqrt(2). 474 * Dividing by E causes an insignificant loss of accuracy; however 475 * this method is still poor since it is uneccessarily slow. 476 */ 477 if (ax > DBL_MAX / 2) |
| 478 return (cpack(log(hypot(x / m_e, y / m_e)) + 1, atan2(y, x))); | 478 return (CMPLX(log(hypot(x / m_e, y / m_e)) + 1, atan2(y, x))); |
| 479 480 /* 481 * Avoid overflow when x or y is large. Avoid underflow when x or 482 * y is small. 483 */ 484 if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) | 479 480 /* 481 * Avoid overflow when x or y is large. Avoid underflow when x or 482 * y is small. 483 */ 484 if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) |
| 485 return (cpack(log(hypot(x, y)), atan2(y, x))); | 485 return (CMPLX(log(hypot(x, y)), atan2(y, x))); |
| 486 | 486 |
| 487 return (cpack(log(ax * ax + ay * ay) / 2, atan2(y, x))); | 487 return (CMPLX(log(ax * ax + ay * ay) / 2, atan2(y, x))); |
| 488} 489 490/* 491 * ================= 492 * | catanh, catan | 493 * ================= 494 */ 495 --- 74 unchanged lines hidden (view full) --- 570 571 x = creal(z); 572 y = cimag(z); 573 ax = fabs(x); 574 ay = fabs(y); 575 576 /* This helps handle many cases. */ 577 if (y == 0 && ax <= 1) | 488} 489 490/* 491 * ================= 492 * | catanh, catan | 493 * ================= 494 */ 495 --- 74 unchanged lines hidden (view full) --- 570 571 x = creal(z); 572 y = cimag(z); 573 ax = fabs(x); 574 ay = fabs(y); 575 576 /* This helps handle many cases. */ 577 if (y == 0 && ax <= 1) |
| 578 return (cpack(atanh(x), y)); | 578 return (CMPLX(atanh(x), y)); |
| 579 580 /* To ensure the same accuracy as atan(), and to filter out z = 0. */ 581 if (x == 0) | 579 580 /* To ensure the same accuracy as atan(), and to filter out z = 0. */ 581 if (x == 0) |
| 582 return (cpack(x, atan(y))); | 582 return (CMPLX(x, atan(y))); |
| 583 584 if (isnan(x) || isnan(y)) { 585 /* catanh(+-Inf + I*NaN) = +-0 + I*NaN */ 586 if (isinf(x)) | 583 584 if (isnan(x) || isnan(y)) { 585 /* catanh(+-Inf + I*NaN) = +-0 + I*NaN */ 586 if (isinf(x)) |
| 587 return (cpack(copysign(0, x), y + y)); | 587 return (CMPLX(copysign(0, x), y + y)); |
| 588 /* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */ 589 if (isinf(y)) | 588 /* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */ 589 if (isinf(y)) |
| 590 return (cpack(copysign(0, x), | 590 return (CMPLX(copysign(0, x), |
| 591 copysign(pio2_hi + pio2_lo, y))); 592 /* 593 * All other cases involving NaN return NaN + I*NaN. 594 * C99 leaves it optional whether to raise invalid if one of 595 * the arguments is not NaN, so we opt not to raise it. 596 */ | 591 copysign(pio2_hi + pio2_lo, y))); 592 /* 593 * All other cases involving NaN return NaN + I*NaN. 594 * C99 leaves it optional whether to raise invalid if one of 595 * the arguments is not NaN, so we opt not to raise it. 596 */ |
| 597 return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); | 597 return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); |
| 598 } 599 600 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) | 598 } 599 600 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) |
| 601 return (cpack(real_part_reciprocal(x, y), | 601 return (CMPLX(real_part_reciprocal(x, y), |
| 602 copysign(pio2_hi + pio2_lo, y))); 603 604 if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { 605 /* 606 * z = 0 was filtered out above. All other cases must raise 607 * inexact, but this is the only only that needs to do it 608 * explicitly. 609 */ --- 8 unchanged lines hidden (view full) --- 618 619 if (ax == 1) 620 ry = atan2(2, -ay) / 2; 621 else if (ay < DBL_EPSILON) 622 ry = atan2(2 * ay, (1 - ax) * (1 + ax)) / 2; 623 else 624 ry = atan2(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; 625 | 602 copysign(pio2_hi + pio2_lo, y))); 603 604 if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { 605 /* 606 * z = 0 was filtered out above. All other cases must raise 607 * inexact, but this is the only only that needs to do it 608 * explicitly. 609 */ --- 8 unchanged lines hidden (view full) --- 618 619 if (ax == 1) 620 ry = atan2(2, -ay) / 2; 621 else if (ay < DBL_EPSILON) 622 ry = atan2(2 * ay, (1 - ax) * (1 + ax)) / 2; 623 else 624 ry = atan2(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; 625 |
| 626 return (cpack(copysign(rx, x), copysign(ry, y))); | 626 return (CMPLX(copysign(rx, x), copysign(ry, y))); |
| 627} 628 629/* 630 * catan(z) = reverse(catanh(reverse(z))) 631 * where reverse(x + I*y) = y + I*x = I*conj(z). 632 */ 633double complex 634catan(double complex z) 635{ | 627} 628 629/* 630 * catan(z) = reverse(catanh(reverse(z))) 631 * where reverse(x + I*y) = y + I*x = I*conj(z). 632 */ 633double complex 634catan(double complex z) 635{ |
| 636 double complex w = catanh(cpack(cimag(z), creal(z))); | 636 double complex w = catanh(CMPLX(cimag(z), creal(z))); |
| 637 | 637 |
| 638 return (cpack(cimag(w), creal(w))); | 638 return (CMPLX(cimag(w), creal(w))); |
| 639} | 639} |