| | 1 | /* Implementation of various C99 functions |
| | 2 | Copyright (C) 2004-2020 Free Software Foundation, Inc. |
| | 3 | |
| | 4 | This file is part of the GNU Fortran 95 runtime library (libgfortran). |
| | 5 | |
| | 6 | Libgfortran is free software; you can redistribute it and/or |
| | 7 | modify it under the terms of the GNU General Public |
| | 8 | License as published by the Free Software Foundation; either |
| | 9 | version 3 of the License, or (at your option) any later version. |
| | 10 | |
| | 11 | Libgfortran is distributed in the hope that it will be useful, |
| | 12 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
| | 13 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| | 14 | GNU General Public License for more details. |
| | 15 | |
| | 16 | Under Section 7 of GPL version 3, you are granted additional |
| | 17 | permissions described in the GCC Runtime Library Exception, version |
| | 18 | 3.1, as published by the Free Software Foundation. |
| | 19 | |
| | 20 | You should have received a copy of the GNU General Public License and |
| | 21 | a copy of the GCC Runtime Library Exception along with this program; |
| | 22 | see the files COPYING3 and COPYING.RUNTIME respectively. If not, see |
| | 23 | <http://www.gnu.org/licenses/>. */ |
| | 24 | |
| | 25 | #include "config.h" |
| | 26 | |
| | 27 | #define C99_PROTOS_H WE_DONT_WANT_PROTOS_NOW |
| | 28 | #include "libgfortran.h" |
| | 29 | |
| | 30 | /* On a C99 system "I" (with I*I = -1) should be defined in complex.h; |
| | 31 | if not, we define a fallback version here. */ |
| | 32 | #ifndef I |
| | 33 | # if defined(_Imaginary_I) |
| | 34 | # define I _Imaginary_I |
| | 35 | # elif defined(_Complex_I) |
| | 36 | # define I _Complex_I |
| | 37 | # else |
| | 38 | # define I (1.0fi) |
| | 39 | # endif |
| | 40 | #endif |
| | 41 | |
| | 42 | /* Macros to get real and imaginary parts of a complex, and set |
| | 43 | a complex value. */ |
| | 44 | #define REALPART(z) (__real__(z)) |
| | 45 | #define IMAGPART(z) (__imag__(z)) |
| | 46 | #define COMPLEX_ASSIGN(z_, r_, i_) {__real__(z_) = (r_); __imag__(z_) = (i_);} |
| | 47 | |
| | 48 | |
| | 49 | /* Prototypes are included to silence -Wstrict-prototypes |
| | 50 | -Wmissing-prototypes. */ |
| | 51 | |
| | 52 | |
| | 53 | /* Wrappers for systems without the various C99 single precision Bessel |
| | 54 | functions. */ |
| | 55 | |
| | 56 | #if defined(HAVE_J0) && ! defined(HAVE_J0F) |
| | 57 | #define HAVE_J0F 1 |
| | 58 | float j0f (float); |
| | 59 | |
| | 60 | float |
| | 61 | j0f (float x) |
| | 62 | { |
| | 63 | return (float) j0 ((double) x); |
| | 64 | } |
| | 65 | #endif |
| | 66 | |
| | 67 | #if defined(HAVE_J1) && !defined(HAVE_J1F) |
| | 68 | #define HAVE_J1F 1 |
| | 69 | float j1f (float); |
| | 70 | |
| | 71 | float j1f (float x) |
| | 72 | { |
| | 73 | return (float) j1 ((double) x); |
| | 74 | } |
| | 75 | #endif |
| | 76 | |
| | 77 | #if defined(HAVE_JN) && !defined(HAVE_JNF) |
| | 78 | #define HAVE_JNF 1 |
| | 79 | float jnf (int, float); |
| | 80 | |
| | 81 | float |
| | 82 | jnf (int n, float x) |
| | 83 | { |
| | 84 | return (float) jn (n, (double) x); |
| | 85 | } |
| | 86 | #endif |
| | 87 | |
| | 88 | #if defined(HAVE_Y0) && !defined(HAVE_Y0F) |
| | 89 | #define HAVE_Y0F 1 |
| | 90 | float y0f (float); |
| | 91 | |
| | 92 | float |
| | 93 | y0f (float x) |
| | 94 | { |
| | 95 | return (float) y0 ((double) x); |
| | 96 | } |
| | 97 | #endif |
| | 98 | |
| | 99 | #if defined(HAVE_Y1) && !defined(HAVE_Y1F) |
| | 100 | #define HAVE_Y1F 1 |
| | 101 | float y1f (float); |
| | 102 | |
| | 103 | float |
| | 104 | y1f (float x) |
| | 105 | { |
| | 106 | return (float) y1 ((double) x); |
| | 107 | } |
| | 108 | #endif |
| | 109 | |
| | 110 | #if defined(HAVE_YN) && !defined(HAVE_YNF) |
| | 111 | #define HAVE_YNF 1 |
| | 112 | float ynf (int, float); |
| | 113 | |
| | 114 | float |
| | 115 | ynf (int n, float x) |
| | 116 | { |
| | 117 | return (float) yn (n, (double) x); |
| | 118 | } |
| | 119 | #endif |
| | 120 | |
| | 121 | |
| | 122 | /* Wrappers for systems without the C99 erff() and erfcf() functions. */ |
| | 123 | |
| | 124 | #if defined(HAVE_ERF) && !defined(HAVE_ERFF) |
| | 125 | #define HAVE_ERFF 1 |
| | 126 | float erff (float); |
| | 127 | |
| | 128 | float |
| | 129 | erff (float x) |
| | 130 | { |
| | 131 | return (float) erf ((double) x); |
| | 132 | } |
| | 133 | #endif |
| | 134 | |
| | 135 | #if defined(HAVE_ERFC) && !defined(HAVE_ERFCF) |
| | 136 | #define HAVE_ERFCF 1 |
| | 137 | float erfcf (float); |
| | 138 | |
| | 139 | float |
| | 140 | erfcf (float x) |
| | 141 | { |
| | 142 | return (float) erfc ((double) x); |
| | 143 | } |
| | 144 | #endif |
| | 145 | |
| | 146 | |
| | 147 | #ifndef HAVE_ACOSF |
| | 148 | #define HAVE_ACOSF 1 |
| | 149 | float acosf (float x); |
| | 150 | |
| | 151 | float |
| | 152 | acosf (float x) |
| | 153 | { |
| | 154 | return (float) acos (x); |
| | 155 | } |
| | 156 | #endif |
| | 157 | |
| | 158 | #if HAVE_ACOSH && !HAVE_ACOSHF |
| | 159 | float acoshf (float x); |
| | 160 | |
| | 161 | float |
| | 162 | acoshf (float x) |
| | 163 | { |
| | 164 | return (float) acosh ((double) x); |
| | 165 | } |
| | 166 | #endif |
| | 167 | |
| | 168 | #ifndef HAVE_ASINF |
| | 169 | #define HAVE_ASINF 1 |
| | 170 | float asinf (float x); |
| | 171 | |
| | 172 | float |
| | 173 | asinf (float x) |
| | 174 | { |
| | 175 | return (float) asin (x); |
| | 176 | } |
| | 177 | #endif |
| | 178 | |
| | 179 | #if HAVE_ASINH && !HAVE_ASINHF |
| | 180 | float asinhf (float x); |
| | 181 | |
| | 182 | float |
| | 183 | asinhf (float x) |
| | 184 | { |
| | 185 | return (float) asinh ((double) x); |
| | 186 | } |
| | 187 | #endif |
| | 188 | |
| | 189 | #ifndef HAVE_ATAN2F |
| | 190 | #define HAVE_ATAN2F 1 |
| | 191 | float atan2f (float y, float x); |
| | 192 | |
| | 193 | float |
| | 194 | atan2f (float y, float x) |
| | 195 | { |
| | 196 | return (float) atan2 (y, x); |
| | 197 | } |
| | 198 | #endif |
| | 199 | |
| | 200 | #ifndef HAVE_ATANF |
| | 201 | #define HAVE_ATANF 1 |
| | 202 | float atanf (float x); |
| | 203 | |
| | 204 | float |
| | 205 | atanf (float x) |
| | 206 | { |
| | 207 | return (float) atan (x); |
| | 208 | } |
| | 209 | #endif |
| | 210 | |
| | 211 | #if HAVE_ATANH && !HAVE_ATANHF |
| | 212 | float atanhf (float x); |
| | 213 | |
| | 214 | float |
| | 215 | atanhf (float x) |
| | 216 | { |
| | 217 | return (float) atanh ((double) x); |
| | 218 | } |
| | 219 | #endif |
| | 220 | |
| | 221 | #ifndef HAVE_CEILF |
| | 222 | #define HAVE_CEILF 1 |
| | 223 | float ceilf (float x); |
| | 224 | |
| | 225 | float |
| | 226 | ceilf (float x) |
| | 227 | { |
| | 228 | return (float) ceil (x); |
| | 229 | } |
| | 230 | #endif |
| | 231 | |
| | 232 | #ifndef HAVE_COPYSIGNF |
| | 233 | #define HAVE_COPYSIGNF 1 |
| | 234 | float copysignf (float x, float y); |
| | 235 | |
| | 236 | float |
| | 237 | copysignf (float x, float y) |
| | 238 | { |
| | 239 | return (float) copysign (x, y); |
| | 240 | } |
| | 241 | #endif |
| | 242 | |
| | 243 | #ifndef HAVE_COSF |
| | 244 | #define HAVE_COSF 1 |
| | 245 | float cosf (float x); |
| | 246 | |
| | 247 | float |
| | 248 | cosf (float x) |
| | 249 | { |
| | 250 | return (float) cos (x); |
| | 251 | } |
| | 252 | #endif |
| | 253 | |
| | 254 | #ifndef HAVE_COSHF |
| | 255 | #define HAVE_COSHF 1 |
| | 256 | float coshf (float x); |
| | 257 | |
| | 258 | float |
| | 259 | coshf (float x) |
| | 260 | { |
| | 261 | return (float) cosh (x); |
| | 262 | } |
| | 263 | #endif |
| | 264 | |
| | 265 | #ifndef HAVE_EXPF |
| | 266 | #define HAVE_EXPF 1 |
| | 267 | float expf (float x); |
| | 268 | |
| | 269 | float |
| | 270 | expf (float x) |
| | 271 | { |
| | 272 | return (float) exp (x); |
| | 273 | } |
| | 274 | #endif |
| | 275 | |
| | 276 | #ifndef HAVE_FABSF |
| | 277 | #define HAVE_FABSF 1 |
| | 278 | float fabsf (float x); |
| | 279 | |
| | 280 | float |
| | 281 | fabsf (float x) |
| | 282 | { |
| | 283 | return (float) fabs (x); |
| | 284 | } |
| | 285 | #endif |
| | 286 | |
| | 287 | #ifndef HAVE_FLOORF |
| | 288 | #define HAVE_FLOORF 1 |
| | 289 | float floorf (float x); |
| | 290 | |
| | 291 | float |
| | 292 | floorf (float x) |
| | 293 | { |
| | 294 | return (float) floor (x); |
| | 295 | } |
| | 296 | #endif |
| | 297 | |
| | 298 | #ifndef HAVE_FMODF |
| | 299 | #define HAVE_FMODF 1 |
| | 300 | float fmodf (float x, float y); |
| | 301 | |
| | 302 | float |
| | 303 | fmodf (float x, float y) |
| | 304 | { |
| | 305 | return (float) fmod (x, y); |
| | 306 | } |
| | 307 | #endif |
| | 308 | |
| | 309 | #ifndef HAVE_FREXPF |
| | 310 | #define HAVE_FREXPF 1 |
| | 311 | float frexpf (float x, int *exp); |
| | 312 | |
| | 313 | float |
| | 314 | frexpf (float x, int *exp) |
| | 315 | { |
| | 316 | return (float) frexp (x, exp); |
| | 317 | } |
| | 318 | #endif |
| | 319 | |
| | 320 | #ifndef HAVE_HYPOTF |
| | 321 | #define HAVE_HYPOTF 1 |
| | 322 | float hypotf (float x, float y); |
| | 323 | |
| | 324 | float |
| | 325 | hypotf (float x, float y) |
| | 326 | { |
| | 327 | return (float) hypot (x, y); |
| | 328 | } |
| | 329 | #endif |
| | 330 | |
| | 331 | #ifndef HAVE_LOGF |
| | 332 | #define HAVE_LOGF 1 |
| | 333 | float logf (float x); |
| | 334 | |
| | 335 | float |
| | 336 | logf (float x) |
| | 337 | { |
| | 338 | return (float) log (x); |
| | 339 | } |
| | 340 | #endif |
| | 341 | |
| | 342 | #ifndef HAVE_LOG10F |
| | 343 | #define HAVE_LOG10F 1 |
| | 344 | float log10f (float x); |
| | 345 | |
| | 346 | float |
| | 347 | log10f (float x) |
| | 348 | { |
| | 349 | return (float) log10 (x); |
| | 350 | } |
| | 351 | #endif |
| | 352 | |
| | 353 | #ifndef HAVE_SCALBN |
| | 354 | #define HAVE_SCALBN 1 |
| | 355 | double scalbn (double x, int y); |
| | 356 | |
| | 357 | double |
| | 358 | scalbn (double x, int y) |
| | 359 | { |
| | 360 | #if (FLT_RADIX == 2) && defined(HAVE_LDEXP) |
| | 361 | return ldexp (x, y); |
| | 362 | #else |
| | 363 | return x * pow (FLT_RADIX, y); |
| | 364 | #endif |
| | 365 | } |
| | 366 | #endif |
| | 367 | |
| | 368 | #ifndef HAVE_SCALBNF |
| | 369 | #define HAVE_SCALBNF 1 |
| | 370 | float scalbnf (float x, int y); |
| | 371 | |
| | 372 | float |
| | 373 | scalbnf (float x, int y) |
| | 374 | { |
| | 375 | return (float) scalbn (x, y); |
| | 376 | } |
| | 377 | #endif |
| | 378 | |
| | 379 | #ifndef HAVE_SINF |
| | 380 | #define HAVE_SINF 1 |
| | 381 | float sinf (float x); |
| | 382 | |
| | 383 | float |
| | 384 | sinf (float x) |
| | 385 | { |
| | 386 | return (float) sin (x); |
| | 387 | } |
| | 388 | #endif |
| | 389 | |
| | 390 | #ifndef HAVE_SINHF |
| | 391 | #define HAVE_SINHF 1 |
| | 392 | float sinhf (float x); |
| | 393 | |
| | 394 | float |
| | 395 | sinhf (float x) |
| | 396 | { |
| | 397 | return (float) sinh (x); |
| | 398 | } |
| | 399 | #endif |
| | 400 | |
| | 401 | #ifndef HAVE_SQRTF |
| | 402 | #define HAVE_SQRTF 1 |
| | 403 | float sqrtf (float x); |
| | 404 | |
| | 405 | float |
| | 406 | sqrtf (float x) |
| | 407 | { |
| | 408 | return (float) sqrt (x); |
| | 409 | } |
| | 410 | #endif |
| | 411 | |
| | 412 | #ifndef HAVE_TANF |
| | 413 | #define HAVE_TANF 1 |
| | 414 | float tanf (float x); |
| | 415 | |
| | 416 | float |
| | 417 | tanf (float x) |
| | 418 | { |
| | 419 | return (float) tan (x); |
| | 420 | } |
| | 421 | #endif |
| | 422 | |
| | 423 | #ifndef HAVE_TANHF |
| | 424 | #define HAVE_TANHF 1 |
| | 425 | float tanhf (float x); |
| | 426 | |
| | 427 | float |
| | 428 | tanhf (float x) |
| | 429 | { |
| | 430 | return (float) tanh (x); |
| | 431 | } |
| | 432 | #endif |
| | 433 | |
| | 434 | #ifndef HAVE_TRUNC |
| | 435 | #define HAVE_TRUNC 1 |
| | 436 | double trunc (double x); |
| | 437 | |
| | 438 | double |
| | 439 | trunc (double x) |
| | 440 | { |
| | 441 | if (!isfinite (x)) |
| | 442 | return x; |
| | 443 | |
| | 444 | if (x < 0.0) |
| | 445 | return - floor (-x); |
| | 446 | else |
| | 447 | return floor (x); |
| | 448 | } |
| | 449 | #endif |
| | 450 | |
| | 451 | #ifndef HAVE_TRUNCF |
| | 452 | #define HAVE_TRUNCF 1 |
| | 453 | float truncf (float x); |
| | 454 | |
| | 455 | float |
| | 456 | truncf (float x) |
| | 457 | { |
| | 458 | return (float) trunc (x); |
| | 459 | } |
| | 460 | #endif |
| | 461 | |
| | 462 | #ifndef HAVE_NEXTAFTERF |
| | 463 | #define HAVE_NEXTAFTERF 1 |
| | 464 | /* This is a portable implementation of nextafterf that is intended to be |
| | 465 | independent of the floating point format or its in memory representation. |
| | 466 | This implementation works correctly with denormalized values. */ |
| | 467 | float nextafterf (float x, float y); |
| | 468 | |
| | 469 | float |
| | 470 | nextafterf (float x, float y) |
| | 471 | { |
| | 472 | /* This variable is marked volatile to avoid excess precision problems |
| | 473 | on some platforms, including IA-32. */ |
| | 474 | volatile float delta; |
| | 475 | float absx, denorm_min; |
| | 476 | |
| | 477 | if (isnan (x) || isnan (y)) |
| | 478 | return x + y; |
| | 479 | if (x == y) |
| | 480 | return x; |
| | 481 | if (!isfinite (x)) |
| | 482 | return x > 0 ? __FLT_MAX__ : - __FLT_MAX__; |
| | 483 | |
| | 484 | /* absx = fabsf (x); */ |
| | 485 | absx = (x < 0.0) ? -x : x; |
| | 486 | |
| | 487 | /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals. */ |
| | 488 | if (__FLT_DENORM_MIN__ == 0.0f) |
| | 489 | denorm_min = __FLT_MIN__; |
| | 490 | else |
| | 491 | denorm_min = __FLT_DENORM_MIN__; |
| | 492 | |
| | 493 | if (absx < __FLT_MIN__) |
| | 494 | delta = denorm_min; |
| | 495 | else |
| | 496 | { |
| | 497 | float frac; |
| | 498 | int exp; |
| | 499 | |
| | 500 | /* Discard the fraction from x. */ |
| | 501 | frac = frexpf (absx, &exp); |
| | 502 | delta = scalbnf (0.5f, exp); |
| | 503 | |
| | 504 | /* Scale x by the epsilon of the representation. By rights we should |
| | 505 | have been able to combine this with scalbnf, but some targets don't |
| | 506 | get that correct with denormals. */ |
| | 507 | delta *= __FLT_EPSILON__; |
| | 508 | |
| | 509 | /* If we're going to be reducing the absolute value of X, and doing so |
| | 510 | would reduce the exponent of X, then the delta to be applied is |
| | 511 | one exponent smaller. */ |
| | 512 | if (frac == 0.5f && (y < x) == (x > 0)) |
| | 513 | delta *= 0.5f; |
| | 514 | |
| | 515 | /* If that underflows to zero, then we're back to the minimum. */ |
| | 516 | if (delta == 0.0f) |
| | 517 | delta = denorm_min; |
| | 518 | } |
| | 519 | |
| | 520 | if (y < x) |
| | 521 | delta = -delta; |
| | 522 | |
| | 523 | return x + delta; |
| | 524 | } |
| | 525 | #endif |
| | 526 | |
| | 527 | |
| | 528 | #ifndef HAVE_POWF |
| | 529 | #define HAVE_POWF 1 |
| | 530 | float powf (float x, float y); |
| | 531 | |
| | 532 | float |
| | 533 | powf (float x, float y) |
| | 534 | { |
| | 535 | return (float) pow (x, y); |
| | 536 | } |
| | 537 | #endif |
| | 538 | |
| | 539 | |
| | 540 | #ifndef HAVE_ROUND |
| | 541 | #define HAVE_ROUND 1 |
| | 542 | /* Round to nearest integral value. If the argument is halfway between two |
| | 543 | integral values then round away from zero. */ |
| | 544 | double round (double x); |
| | 545 | |
| | 546 | double |
| | 547 | round (double x) |
| | 548 | { |
| | 549 | double t; |
| | 550 | if (!isfinite (x)) |
| | 551 | return (x); |
| | 552 | |
| | 553 | if (x >= 0.0) |
| | 554 | { |
| | 555 | t = floor (x); |
| | 556 | if (t - x <= -0.5) |
| | 557 | t += 1.0; |
| | 558 | return (t); |
| | 559 | } |
| | 560 | else |
| | 561 | { |
| | 562 | t = floor (-x); |
| | 563 | if (t + x <= -0.5) |
| | 564 | t += 1.0; |
| | 565 | return (-t); |
| | 566 | } |
| | 567 | } |
| | 568 | #endif |
| | 569 | |
| | 570 | |
| | 571 | /* Algorithm by Steven G. Kargl. */ |
| | 572 | |
| | 573 | #if !defined(HAVE_ROUNDL) |
| | 574 | #define HAVE_ROUNDL 1 |
| | 575 | long double roundl (long double x); |
| | 576 | |
| | 577 | #if defined(HAVE_CEILL) |
| | 578 | /* Round to nearest integral value. If the argument is halfway between two |
| | 579 | integral values then round away from zero. */ |
| | 580 | |
| | 581 | long double |
| | 582 | roundl (long double x) |
| | 583 | { |
| | 584 | long double t; |
| | 585 | if (!isfinite (x)) |
| | 586 | return (x); |
| | 587 | |
| | 588 | if (x >= 0.0) |
| | 589 | { |
| | 590 | t = ceill (x); |
| | 591 | if (t - x > 0.5) |
| | 592 | t -= 1.0; |
| | 593 | return (t); |
| | 594 | } |
| | 595 | else |
| | 596 | { |
| | 597 | t = ceill (-x); |
| | 598 | if (t + x > 0.5) |
| | 599 | t -= 1.0; |
| | 600 | return (-t); |
| | 601 | } |
| | 602 | } |
| | 603 | #else |
| | 604 | |
| | 605 | /* Poor version of roundl for system that don't have ceill. */ |
| | 606 | long double |
| | 607 | roundl (long double x) |
| | 608 | { |
| | 609 | if (x > DBL_MAX || x < -DBL_MAX) |
| | 610 | { |
| | 611 | #ifdef HAVE_NEXTAFTERL |
| | 612 | long double prechalf = nextafterl (0.5L, LDBL_MAX); |
| | 613 | #else |
| | 614 | static long double prechalf = 0.5L; |
| | 615 | #endif |
| | 616 | return (GFC_INTEGER_LARGEST) (x + (x > 0 ? prechalf : -prechalf)); |
| | 617 | } |
| | 618 | else |
| | 619 | /* Use round(). */ |
| | 620 | return round ((double) x); |
| | 621 | } |
| | 622 | |
| | 623 | #endif |
| | 624 | #endif |
| | 625 | |
| | 626 | #ifndef HAVE_ROUNDF |
| | 627 | #define HAVE_ROUNDF 1 |
| | 628 | /* Round to nearest integral value. If the argument is halfway between two |
| | 629 | integral values then round away from zero. */ |
| | 630 | float roundf (float x); |
| | 631 | |
| | 632 | float |
| | 633 | roundf (float x) |
| | 634 | { |
| | 635 | float t; |
| | 636 | if (!isfinite (x)) |
| | 637 | return (x); |
| | 638 | |
| | 639 | if (x >= 0.0) |
| | 640 | { |
| | 641 | t = floorf (x); |
| | 642 | if (t - x <= -0.5) |
| | 643 | t += 1.0; |
| | 644 | return (t); |
| | 645 | } |
| | 646 | else |
| | 647 | { |
| | 648 | t = floorf (-x); |
| | 649 | if (t + x <= -0.5) |
| | 650 | t += 1.0; |
| | 651 | return (-t); |
| | 652 | } |
| | 653 | } |
| | 654 | #endif |
| | 655 | |
| | 656 | |
| | 657 | /* lround{f,,l} and llround{f,,l} functions. */ |
| | 658 | |
| | 659 | #if !defined(HAVE_LROUNDF) && defined(HAVE_ROUNDF) |
| | 660 | #define HAVE_LROUNDF 1 |
| | 661 | long int lroundf (float x); |
| | 662 | |
| | 663 | long int |
| | 664 | lroundf (float x) |
| | 665 | { |
| | 666 | return (long int) roundf (x); |
| | 667 | } |
| | 668 | #endif |
| | 669 | |
| | 670 | #if !defined(HAVE_LROUND) && defined(HAVE_ROUND) |
| | 671 | #define HAVE_LROUND 1 |
| | 672 | long int lround (double x); |
| | 673 | |
| | 674 | long int |
| | 675 | lround (double x) |
| | 676 | { |
| | 677 | return (long int) round (x); |
| | 678 | } |
| | 679 | #endif |
| | 680 | |
| | 681 | #if !defined(HAVE_LROUNDL) && defined(HAVE_ROUNDL) |
| | 682 | #define HAVE_LROUNDL 1 |
| | 683 | long int lroundl (long double x); |
| | 684 | |
| | 685 | long int |
| | 686 | lroundl (long double x) |
| | 687 | { |
| | 688 | return (long long int) roundl (x); |
| | 689 | } |
| | 690 | #endif |
| | 691 | |
| | 692 | #if !defined(HAVE_LLROUNDF) && defined(HAVE_ROUNDF) |
| | 693 | #define HAVE_LLROUNDF 1 |
| | 694 | long long int llroundf (float x); |
| | 695 | |
| | 696 | long long int |
| | 697 | llroundf (float x) |
| | 698 | { |
| | 699 | return (long long int) roundf (x); |
| | 700 | } |
| | 701 | #endif |
| | 702 | |
| | 703 | #if !defined(HAVE_LLROUND) && defined(HAVE_ROUND) |
| | 704 | #define HAVE_LLROUND 1 |
| | 705 | long long int llround (double x); |
| | 706 | |
| | 707 | long long int |
| | 708 | llround (double x) |
| | 709 | { |
| | 710 | return (long long int) round (x); |
| | 711 | } |
| | 712 | #endif |
| | 713 | |
| | 714 | #if !defined(HAVE_LLROUNDL) && defined(HAVE_ROUNDL) |
| | 715 | #define HAVE_LLROUNDL 1 |
| | 716 | long long int llroundl (long double x); |
| | 717 | |
| | 718 | long long int |
| | 719 | llroundl (long double x) |
| | 720 | { |
| | 721 | return (long long int) roundl (x); |
| | 722 | } |
| | 723 | #endif |
| | 724 | |
| | 725 | |
| | 726 | #ifndef HAVE_LOG10L |
| | 727 | #define HAVE_LOG10L 1 |
| | 728 | /* log10 function for long double variables. The version provided here |
| | 729 | reduces the argument until it fits into a double, then use log10. */ |
| | 730 | long double log10l (long double x); |
| | 731 | |
| | 732 | long double |
| | 733 | log10l (long double x) |
| | 734 | { |
| | 735 | #if LDBL_MAX_EXP > DBL_MAX_EXP |
| | 736 | if (x > DBL_MAX) |
| | 737 | { |
| | 738 | double val; |
| | 739 | int p2_result = 0; |
| | 740 | if (x > 0x1p16383L) { p2_result += 16383; x /= 0x1p16383L; } |
| | 741 | if (x > 0x1p8191L) { p2_result += 8191; x /= 0x1p8191L; } |
| | 742 | if (x > 0x1p4095L) { p2_result += 4095; x /= 0x1p4095L; } |
| | 743 | if (x > 0x1p2047L) { p2_result += 2047; x /= 0x1p2047L; } |
| | 744 | if (x > 0x1p1023L) { p2_result += 1023; x /= 0x1p1023L; } |
| | 745 | val = log10 ((double) x); |
| | 746 | return (val + p2_result * .30102999566398119521373889472449302L); |
| | 747 | } |
| | 748 | #endif |
| | 749 | #if LDBL_MIN_EXP < DBL_MIN_EXP |
| | 750 | if (x < DBL_MIN) |
| | 751 | { |
| | 752 | double val; |
| | 753 | int p2_result = 0; |
| | 754 | if (x < 0x1p-16380L) { p2_result += 16380; x /= 0x1p-16380L; } |
| | 755 | if (x < 0x1p-8189L) { p2_result += 8189; x /= 0x1p-8189L; } |
| | 756 | if (x < 0x1p-4093L) { p2_result += 4093; x /= 0x1p-4093L; } |
| | 757 | if (x < 0x1p-2045L) { p2_result += 2045; x /= 0x1p-2045L; } |
| | 758 | if (x < 0x1p-1021L) { p2_result += 1021; x /= 0x1p-1021L; } |
| | 759 | val = fabs (log10 ((double) x)); |
| | 760 | return (- val - p2_result * .30102999566398119521373889472449302L); |
| | 761 | } |
| | 762 | #endif |
| | 763 | return log10 (x); |
| | 764 | } |
| | 765 | #endif |
| | 766 | |
| | 767 | |
| | 768 | #ifndef HAVE_FLOORL |
| | 769 | #define HAVE_FLOORL 1 |
| | 770 | long double floorl (long double x); |
| | 771 | |
| | 772 | long double |
| | 773 | floorl (long double x) |
| | 774 | { |
| | 775 | /* Zero, possibly signed. */ |
| | 776 | if (x == 0) |
| | 777 | return x; |
| | 778 | |
| | 779 | /* Large magnitude. */ |
| | 780 | if (x > DBL_MAX || x < (-DBL_MAX)) |
| | 781 | return x; |
| | 782 | |
| | 783 | /* Small positive values. */ |
| | 784 | if (x >= 0 && x < DBL_MIN) |
| | 785 | return 0; |
| | 786 | |
| | 787 | /* Small negative values. */ |
| | 788 | if (x < 0 && x > (-DBL_MIN)) |
| | 789 | return -1; |
| | 790 | |
| | 791 | return floor (x); |
| | 792 | } |
| | 793 | #endif |
| | 794 | |
| | 795 | |
| | 796 | #ifndef HAVE_FMODL |
| | 797 | #define HAVE_FMODL 1 |
| | 798 | long double fmodl (long double x, long double y); |
| | 799 | |
| | 800 | long double |
| | 801 | fmodl (long double x, long double y) |
| | 802 | { |
| | 803 | if (y == 0.0L) |
| | 804 | return 0.0L; |
| | 805 | |
| | 806 | /* Need to check that the result has the same sign as x and magnitude |
| | 807 | less than the magnitude of y. */ |
| | 808 | return x - floorl (x / y) * y; |
| | 809 | } |
| | 810 | #endif |
| | 811 | |
| | 812 | |
| | 813 | #if !defined(HAVE_CABSF) |
| | 814 | #define HAVE_CABSF 1 |
| | 815 | float cabsf (float complex z); |
| | 816 | |
| | 817 | float |
| | 818 | cabsf (float complex z) |
| | 819 | { |
| | 820 | return hypotf (REALPART (z), IMAGPART (z)); |
| | 821 | } |
| | 822 | #endif |
| | 823 | |
| | 824 | #if !defined(HAVE_CABS) |
| | 825 | #define HAVE_CABS 1 |
| | 826 | double cabs (double complex z); |
| | 827 | |
| | 828 | double |
| | 829 | cabs (double complex z) |
| | 830 | { |
| | 831 | return hypot (REALPART (z), IMAGPART (z)); |
| | 832 | } |
| | 833 | #endif |
| | 834 | |
| | 835 | #if !defined(HAVE_CABSL) && defined(HAVE_HYPOTL) |
| | 836 | #define HAVE_CABSL 1 |
| | 837 | long double cabsl (long double complex z); |
| | 838 | |
| | 839 | long double |
| | 840 | cabsl (long double complex z) |
| | 841 | { |
| | 842 | return hypotl (REALPART (z), IMAGPART (z)); |
| | 843 | } |
| | 844 | #endif |
| | 845 | |
| | 846 | |
| | 847 | #if !defined(HAVE_CARGF) |
| | 848 | #define HAVE_CARGF 1 |
| | 849 | float cargf (float complex z); |
| | 850 | |
| | 851 | float |
| | 852 | cargf (float complex z) |
| | 853 | { |
| | 854 | return atan2f (IMAGPART (z), REALPART (z)); |
| | 855 | } |
| | 856 | #endif |
| | 857 | |
| | 858 | #if !defined(HAVE_CARG) |
| | 859 | #define HAVE_CARG 1 |
| | 860 | double carg (double complex z); |
| | 861 | |
| | 862 | double |
| | 863 | carg (double complex z) |
| | 864 | { |
| | 865 | return atan2 (IMAGPART (z), REALPART (z)); |
| | 866 | } |
| | 867 | #endif |
| | 868 | |
| | 869 | #if !defined(HAVE_CARGL) && defined(HAVE_ATAN2L) |
| | 870 | #define HAVE_CARGL 1 |
| | 871 | long double cargl (long double complex z); |
| | 872 | |
| | 873 | long double |
| | 874 | cargl (long double complex z) |
| | 875 | { |
| | 876 | return atan2l (IMAGPART (z), REALPART (z)); |
| | 877 | } |
| | 878 | #endif |
| | 879 | |
| | 880 | |
| | 881 | /* exp(z) = exp(a)*(cos(b) + i sin(b)) */ |
| | 882 | #if !defined(HAVE_CEXPF) |
| | 883 | #define HAVE_CEXPF 1 |
| | 884 | float complex cexpf (float complex z); |
| | 885 | |
| | 886 | float complex |
| | 887 | cexpf (float complex z) |
| | 888 | { |
| | 889 | float a, b; |
| | 890 | float complex v; |
| | 891 | |
| | 892 | a = REALPART (z); |
| | 893 | b = IMAGPART (z); |
| | 894 | COMPLEX_ASSIGN (v, cosf (b), sinf (b)); |
| | 895 | return expf (a) * v; |
| | 896 | } |
| | 897 | #endif |
| | 898 | |
| | 899 | #if !defined(HAVE_CEXP) |
| | 900 | #define HAVE_CEXP 1 |
| | 901 | double complex cexp (double complex z); |
| | 902 | |
| | 903 | double complex |
| | 904 | cexp (double complex z) |
| | 905 | { |
| | 906 | double a, b; |
| | 907 | double complex v; |
| | 908 | |
| | 909 | a = REALPART (z); |
| | 910 | b = IMAGPART (z); |
| | 911 | COMPLEX_ASSIGN (v, cos (b), sin (b)); |
| | 912 | return exp (a) * v; |
| | 913 | } |
| | 914 | #endif |
| | 915 | |
| | 916 | #if !defined(HAVE_CEXPL) && defined(HAVE_COSL) && defined(HAVE_SINL) && defined(HAVE_EXPL) |
| | 917 | #define HAVE_CEXPL 1 |
| | 918 | long double complex cexpl (long double complex z); |
| | 919 | |
| | 920 | long double complex |
| | 921 | cexpl (long double complex z) |
| | 922 | { |
| | 923 | long double a, b; |
| | 924 | long double complex v; |
| | 925 | |
| | 926 | a = REALPART (z); |
| | 927 | b = IMAGPART (z); |
| | 928 | COMPLEX_ASSIGN (v, cosl (b), sinl (b)); |
| | 929 | return expl (a) * v; |
| | 930 | } |
| | 931 | #endif |
| | 932 | |
| | 933 | |
| | 934 | /* log(z) = log (cabs(z)) + i*carg(z) */ |
| | 935 | #if !defined(HAVE_CLOGF) |
| | 936 | #define HAVE_CLOGF 1 |
| | 937 | float complex clogf (float complex z); |
| | 938 | |
| | 939 | float complex |
| | 940 | clogf (float complex z) |
| | 941 | { |
| | 942 | float complex v; |
| | 943 | |
| | 944 | COMPLEX_ASSIGN (v, logf (cabsf (z)), cargf (z)); |
| | 945 | return v; |
| | 946 | } |
| | 947 | #endif |
| | 948 | |
| | 949 | #if !defined(HAVE_CLOG) |
| | 950 | #define HAVE_CLOG 1 |
| | 951 | double complex clog (double complex z); |
| | 952 | |
| | 953 | double complex |
| | 954 | clog (double complex z) |
| | 955 | { |
| | 956 | double complex v; |
| | 957 | |
| | 958 | COMPLEX_ASSIGN (v, log (cabs (z)), carg (z)); |
| | 959 | return v; |
| | 960 | } |
| | 961 | #endif |
| | 962 | |
| | 963 | #if !defined(HAVE_CLOGL) && defined(HAVE_LOGL) && defined(HAVE_CABSL) && defined(HAVE_CARGL) |
| | 964 | #define HAVE_CLOGL 1 |
| | 965 | long double complex clogl (long double complex z); |
| | 966 | |
| | 967 | long double complex |
| | 968 | clogl (long double complex z) |
| | 969 | { |
| | 970 | long double complex v; |
| | 971 | |
| | 972 | COMPLEX_ASSIGN (v, logl (cabsl (z)), cargl (z)); |
| | 973 | return v; |
| | 974 | } |
| | 975 | #endif |
| | 976 | |
| | 977 | |
| | 978 | /* log10(z) = log10 (cabs(z)) + i*carg(z) */ |
| | 979 | #if !defined(HAVE_CLOG10F) |
| | 980 | #define HAVE_CLOG10F 1 |
| | 981 | float complex clog10f (float complex z); |
| | 982 | |
| | 983 | float complex |
| | 984 | clog10f (float complex z) |
| | 985 | { |
| | 986 | float complex v; |
| | 987 | |
| | 988 | COMPLEX_ASSIGN (v, log10f (cabsf (z)), cargf (z)); |
| | 989 | return v; |
| | 990 | } |
| | 991 | #endif |
| | 992 | |
| | 993 | #if !defined(HAVE_CLOG10) |
| | 994 | #define HAVE_CLOG10 1 |
| | 995 | double complex clog10 (double complex z); |
| | 996 | |
| | 997 | double complex |
| | 998 | clog10 (double complex z) |
| | 999 | { |
| | 1000 | double complex v; |
| | 1001 | |
| | 1002 | COMPLEX_ASSIGN (v, log10 (cabs (z)), carg (z)); |
| | 1003 | return v; |
| | 1004 | } |
| | 1005 | #endif |
| | 1006 | |
| | 1007 | #if !defined(HAVE_CLOG10L) && defined(HAVE_LOG10L) && defined(HAVE_CABSL) && defined(HAVE_CARGL) |
| | 1008 | #define HAVE_CLOG10L 1 |
| | 1009 | long double complex clog10l (long double complex z); |
| | 1010 | |
| | 1011 | long double complex |
| | 1012 | clog10l (long double complex z) |
| | 1013 | { |
| | 1014 | long double complex v; |
| | 1015 | |
| | 1016 | COMPLEX_ASSIGN (v, log10l (cabsl (z)), cargl (z)); |
| | 1017 | return v; |
| | 1018 | } |
| | 1019 | #endif |
| | 1020 | |
| | 1021 | |
| | 1022 | /* pow(base, power) = cexp (power * clog (base)) */ |
| | 1023 | #if !defined(HAVE_CPOWF) |
| | 1024 | #define HAVE_CPOWF 1 |
| | 1025 | float complex cpowf (float complex base, float complex power); |
| | 1026 | |
| | 1027 | float complex |
| | 1028 | cpowf (float complex base, float complex power) |
| | 1029 | { |
| | 1030 | return cexpf (power * clogf (base)); |
| | 1031 | } |
| | 1032 | #endif |
| | 1033 | |
| | 1034 | #if !defined(HAVE_CPOW) |
| | 1035 | #define HAVE_CPOW 1 |
| | 1036 | double complex cpow (double complex base, double complex power); |
| | 1037 | |
| | 1038 | double complex |
| | 1039 | cpow (double complex base, double complex power) |
| | 1040 | { |
| | 1041 | return cexp (power * clog (base)); |
| | 1042 | } |
| | 1043 | #endif |
| | 1044 | |
| | 1045 | #if !defined(HAVE_CPOWL) && defined(HAVE_CEXPL) && defined(HAVE_CLOGL) |
| | 1046 | #define HAVE_CPOWL 1 |
| | 1047 | long double complex cpowl (long double complex base, long double complex power); |
| | 1048 | |
| | 1049 | long double complex |
| | 1050 | cpowl (long double complex base, long double complex power) |
| | 1051 | { |
| | 1052 | return cexpl (power * clogl (base)); |
| | 1053 | } |
| | 1054 | #endif |
| | 1055 | |
| | 1056 | |
| | 1057 | /* sqrt(z). Algorithm pulled from glibc. */ |
| | 1058 | #if !defined(HAVE_CSQRTF) |
| | 1059 | #define HAVE_CSQRTF 1 |
| | 1060 | float complex csqrtf (float complex z); |
| | 1061 | |
| | 1062 | float complex |
| | 1063 | csqrtf (float complex z) |
| | 1064 | { |
| | 1065 | float re, im; |
| | 1066 | float complex v; |
| | 1067 | |
| | 1068 | re = REALPART (z); |
| | 1069 | im = IMAGPART (z); |
| | 1070 | if (im == 0) |
| | 1071 | { |
| | 1072 | if (re < 0) |
| | 1073 | { |
| | 1074 | COMPLEX_ASSIGN (v, 0, copysignf (sqrtf (-re), im)); |
| | 1075 | } |
| | 1076 | else |
| | 1077 | { |
| | 1078 | COMPLEX_ASSIGN (v, fabsf (sqrtf (re)), copysignf (0, im)); |
| | 1079 | } |
| | 1080 | } |
| | 1081 | else if (re == 0) |
| | 1082 | { |
| | 1083 | float r; |
| | 1084 | |
| | 1085 | r = sqrtf (0.5 * fabsf (im)); |
| | 1086 | |
| | 1087 | COMPLEX_ASSIGN (v, r, copysignf (r, im)); |
| | 1088 | } |
| | 1089 | else |
| | 1090 | { |
| | 1091 | float d, r, s; |
| | 1092 | |
| | 1093 | d = hypotf (re, im); |
| | 1094 | /* Use the identity 2 Re res Im res = Im x |
| | 1095 | to avoid cancellation error in d +/- Re x. */ |
| | 1096 | if (re > 0) |
| | 1097 | { |
| | 1098 | r = sqrtf (0.5 * d + 0.5 * re); |
| | 1099 | s = (0.5 * im) / r; |
| | 1100 | } |
| | 1101 | else |
| | 1102 | { |
| | 1103 | s = sqrtf (0.5 * d - 0.5 * re); |
| | 1104 | r = fabsf ((0.5 * im) / s); |
| | 1105 | } |
| | 1106 | |
| | 1107 | COMPLEX_ASSIGN (v, r, copysignf (s, im)); |
| | 1108 | } |
| | 1109 | return v; |
| | 1110 | } |
| | 1111 | #endif |
| | 1112 | |
| | 1113 | #if !defined(HAVE_CSQRT) |
| | 1114 | #define HAVE_CSQRT 1 |
| | 1115 | double complex csqrt (double complex z); |
| | 1116 | |
| | 1117 | double complex |
| | 1118 | csqrt (double complex z) |
| | 1119 | { |
| | 1120 | double re, im; |
| | 1121 | double complex v; |
| | 1122 | |
| | 1123 | re = REALPART (z); |
| | 1124 | im = IMAGPART (z); |
| | 1125 | if (im == 0) |
| | 1126 | { |
| | 1127 | if (re < 0) |
| | 1128 | { |
| | 1129 | COMPLEX_ASSIGN (v, 0, copysign (sqrt (-re), im)); |
| | 1130 | } |
| | 1131 | else |
| | 1132 | { |
| | 1133 | COMPLEX_ASSIGN (v, fabs (sqrt (re)), copysign (0, im)); |
| | 1134 | } |
| | 1135 | } |
| | 1136 | else if (re == 0) |
| | 1137 | { |
| | 1138 | double r; |
| | 1139 | |
| | 1140 | r = sqrt (0.5 * fabs (im)); |
| | 1141 | |
| | 1142 | COMPLEX_ASSIGN (v, r, copysign (r, im)); |
| | 1143 | } |
| | 1144 | else |
| | 1145 | { |
| | 1146 | double d, r, s; |
| | 1147 | |
| | 1148 | d = hypot (re, im); |
| | 1149 | /* Use the identity 2 Re res Im res = Im x |
| | 1150 | to avoid cancellation error in d +/- Re x. */ |
| | 1151 | if (re > 0) |
| | 1152 | { |
| | 1153 | r = sqrt (0.5 * d + 0.5 * re); |
| | 1154 | s = (0.5 * im) / r; |
| | 1155 | } |
| | 1156 | else |
| | 1157 | { |
| | 1158 | s = sqrt (0.5 * d - 0.5 * re); |
| | 1159 | r = fabs ((0.5 * im) / s); |
| | 1160 | } |
| | 1161 | |
| | 1162 | COMPLEX_ASSIGN (v, r, copysign (s, im)); |
| | 1163 | } |
| | 1164 | return v; |
| | 1165 | } |
| | 1166 | #endif |
| | 1167 | |
| | 1168 | #if !defined(HAVE_CSQRTL) && defined(HAVE_COPYSIGNL) && defined(HAVE_SQRTL) && defined(HAVE_FABSL) && defined(HAVE_HYPOTL) |
| | 1169 | #define HAVE_CSQRTL 1 |
| | 1170 | long double complex csqrtl (long double complex z); |
| | 1171 | |
| | 1172 | long double complex |
| | 1173 | csqrtl (long double complex z) |
| | 1174 | { |
| | 1175 | long double re, im; |
| | 1176 | long double complex v; |
| | 1177 | |
| | 1178 | re = REALPART (z); |
| | 1179 | im = IMAGPART (z); |
| | 1180 | if (im == 0) |
| | 1181 | { |
| | 1182 | if (re < 0) |
| | 1183 | { |
| | 1184 | COMPLEX_ASSIGN (v, 0, copysignl (sqrtl (-re), im)); |
| | 1185 | } |
| | 1186 | else |
| | 1187 | { |
| | 1188 | COMPLEX_ASSIGN (v, fabsl (sqrtl (re)), copysignl (0, im)); |
| | 1189 | } |
| | 1190 | } |
| | 1191 | else if (re == 0) |
| | 1192 | { |
| | 1193 | long double r; |
| | 1194 | |
| | 1195 | r = sqrtl (0.5 * fabsl (im)); |
| | 1196 | |
| | 1197 | COMPLEX_ASSIGN (v, copysignl (r, im), r); |
| | 1198 | } |
| | 1199 | else |
| | 1200 | { |
| | 1201 | long double d, r, s; |
| | 1202 | |
| | 1203 | d = hypotl (re, im); |
| | 1204 | /* Use the identity 2 Re res Im res = Im x |
| | 1205 | to avoid cancellation error in d +/- Re x. */ |
| | 1206 | if (re > 0) |
| | 1207 | { |
| | 1208 | r = sqrtl (0.5 * d + 0.5 * re); |
| | 1209 | s = (0.5 * im) / r; |
| | 1210 | } |
| | 1211 | else |
| | 1212 | { |
| | 1213 | s = sqrtl (0.5 * d - 0.5 * re); |
| | 1214 | r = fabsl ((0.5 * im) / s); |
| | 1215 | } |
| | 1216 | |
| | 1217 | COMPLEX_ASSIGN (v, r, copysignl (s, im)); |
| | 1218 | } |
| | 1219 | return v; |
| | 1220 | } |
| | 1221 | #endif |
| | 1222 | |
| | 1223 | |
| | 1224 | /* sinh(a + i b) = sinh(a) cos(b) + i cosh(a) sin(b) */ |
| | 1225 | #if !defined(HAVE_CSINHF) |
| | 1226 | #define HAVE_CSINHF 1 |
| | 1227 | float complex csinhf (float complex a); |
| | 1228 | |
| | 1229 | float complex |
| | 1230 | csinhf (float complex a) |
| | 1231 | { |
| | 1232 | float r, i; |
| | 1233 | float complex v; |
| | 1234 | |
| | 1235 | r = REALPART (a); |
| | 1236 | i = IMAGPART (a); |
| | 1237 | COMPLEX_ASSIGN (v, sinhf (r) * cosf (i), coshf (r) * sinf (i)); |
| | 1238 | return v; |
| | 1239 | } |
| | 1240 | #endif |
| | 1241 | |
| | 1242 | #if !defined(HAVE_CSINH) |
| | 1243 | #define HAVE_CSINH 1 |
| | 1244 | double complex csinh (double complex a); |
| | 1245 | |
| | 1246 | double complex |
| | 1247 | csinh (double complex a) |
| | 1248 | { |
| | 1249 | double r, i; |
| | 1250 | double complex v; |
| | 1251 | |
| | 1252 | r = REALPART (a); |
| | 1253 | i = IMAGPART (a); |
| | 1254 | COMPLEX_ASSIGN (v, sinh (r) * cos (i), cosh (r) * sin (i)); |
| | 1255 | return v; |
| | 1256 | } |
| | 1257 | #endif |
| | 1258 | |
| | 1259 | #if !defined(HAVE_CSINHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL) |
| | 1260 | #define HAVE_CSINHL 1 |
| | 1261 | long double complex csinhl (long double complex a); |
| | 1262 | |
| | 1263 | long double complex |
| | 1264 | csinhl (long double complex a) |
| | 1265 | { |
| | 1266 | long double r, i; |
| | 1267 | long double complex v; |
| | 1268 | |
| | 1269 | r = REALPART (a); |
| | 1270 | i = IMAGPART (a); |
| | 1271 | COMPLEX_ASSIGN (v, sinhl (r) * cosl (i), coshl (r) * sinl (i)); |
| | 1272 | return v; |
| | 1273 | } |
| | 1274 | #endif |
| | 1275 | |
| | 1276 | |
| | 1277 | /* cosh(a + i b) = cosh(a) cos(b) + i sinh(a) sin(b) */ |
| | 1278 | #if !defined(HAVE_CCOSHF) |
| | 1279 | #define HAVE_CCOSHF 1 |
| | 1280 | float complex ccoshf (float complex a); |
| | 1281 | |
| | 1282 | float complex |
| | 1283 | ccoshf (float complex a) |
| | 1284 | { |
| | 1285 | float r, i; |
| | 1286 | float complex v; |
| | 1287 | |
| | 1288 | r = REALPART (a); |
| | 1289 | i = IMAGPART (a); |
| | 1290 | COMPLEX_ASSIGN (v, coshf (r) * cosf (i), sinhf (r) * sinf (i)); |
| | 1291 | return v; |
| | 1292 | } |
| | 1293 | #endif |
| | 1294 | |
| | 1295 | #if !defined(HAVE_CCOSH) |
| | 1296 | #define HAVE_CCOSH 1 |
| | 1297 | double complex ccosh (double complex a); |
| | 1298 | |
| | 1299 | double complex |
| | 1300 | ccosh (double complex a) |
| | 1301 | { |
| | 1302 | double r, i; |
| | 1303 | double complex v; |
| | 1304 | |
| | 1305 | r = REALPART (a); |
| | 1306 | i = IMAGPART (a); |
| | 1307 | COMPLEX_ASSIGN (v, cosh (r) * cos (i), sinh (r) * sin (i)); |
| | 1308 | return v; |
| | 1309 | } |
| | 1310 | #endif |
| | 1311 | |
| | 1312 | #if !defined(HAVE_CCOSHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL) |
| | 1313 | #define HAVE_CCOSHL 1 |
| | 1314 | long double complex ccoshl (long double complex a); |
| | 1315 | |
| | 1316 | long double complex |
| | 1317 | ccoshl (long double complex a) |
| | 1318 | { |
| | 1319 | long double r, i; |
| | 1320 | long double complex v; |
| | 1321 | |
| | 1322 | r = REALPART (a); |
| | 1323 | i = IMAGPART (a); |
| | 1324 | COMPLEX_ASSIGN (v, coshl (r) * cosl (i), sinhl (r) * sinl (i)); |
| | 1325 | return v; |
| | 1326 | } |
| | 1327 | #endif |
| | 1328 | |
| | 1329 | |
| | 1330 | /* tanh(a + i b) = (tanh(a) + i tan(b)) / (1 + i tanh(a) tan(b)) */ |
| | 1331 | #if !defined(HAVE_CTANHF) |
| | 1332 | #define HAVE_CTANHF 1 |
| | 1333 | float complex ctanhf (float complex a); |
| | 1334 | |
| | 1335 | float complex |
| | 1336 | ctanhf (float complex a) |
| | 1337 | { |
| | 1338 | float rt, it; |
| | 1339 | float complex n, d; |
| | 1340 | |
| | 1341 | rt = tanhf (REALPART (a)); |
| | 1342 | it = tanf (IMAGPART (a)); |
| | 1343 | COMPLEX_ASSIGN (n, rt, it); |
| | 1344 | COMPLEX_ASSIGN (d, 1, rt * it); |
| | 1345 | |
| | 1346 | return n / d; |
| | 1347 | } |
| | 1348 | #endif |
| | 1349 | |
| | 1350 | #if !defined(HAVE_CTANH) |
| | 1351 | #define HAVE_CTANH 1 |
| | 1352 | double complex ctanh (double complex a); |
| | 1353 | double complex |
| | 1354 | ctanh (double complex a) |
| | 1355 | { |
| | 1356 | double rt, it; |
| | 1357 | double complex n, d; |
| | 1358 | |
| | 1359 | rt = tanh (REALPART (a)); |
| | 1360 | it = tan (IMAGPART (a)); |
| | 1361 | COMPLEX_ASSIGN (n, rt, it); |
| | 1362 | COMPLEX_ASSIGN (d, 1, rt * it); |
| | 1363 | |
| | 1364 | return n / d; |
| | 1365 | } |
| | 1366 | #endif |
| | 1367 | |
| | 1368 | #if !defined(HAVE_CTANHL) && defined(HAVE_TANL) && defined(HAVE_TANHL) |
| | 1369 | #define HAVE_CTANHL 1 |
| | 1370 | long double complex ctanhl (long double complex a); |
| | 1371 | |
| | 1372 | long double complex |
| | 1373 | ctanhl (long double complex a) |
| | 1374 | { |
| | 1375 | long double rt, it; |
| | 1376 | long double complex n, d; |
| | 1377 | |
| | 1378 | rt = tanhl (REALPART (a)); |
| | 1379 | it = tanl (IMAGPART (a)); |
| | 1380 | COMPLEX_ASSIGN (n, rt, it); |
| | 1381 | COMPLEX_ASSIGN (d, 1, rt * it); |
| | 1382 | |
| | 1383 | return n / d; |
| | 1384 | } |
| | 1385 | #endif |
| | 1386 | |
| | 1387 | |
| | 1388 | /* sin(a + i b) = sin(a) cosh(b) + i cos(a) sinh(b) */ |
| | 1389 | #if !defined(HAVE_CSINF) |
| | 1390 | #define HAVE_CSINF 1 |
| | 1391 | float complex csinf (float complex a); |
| | 1392 | |
| | 1393 | float complex |
| | 1394 | csinf (float complex a) |
| | 1395 | { |
| | 1396 | float r, i; |
| | 1397 | float complex v; |
| | 1398 | |
| | 1399 | r = REALPART (a); |
| | 1400 | i = IMAGPART (a); |
| | 1401 | COMPLEX_ASSIGN (v, sinf (r) * coshf (i), cosf (r) * sinhf (i)); |
| | 1402 | return v; |
| | 1403 | } |
| | 1404 | #endif |
| | 1405 | |
| | 1406 | #if !defined(HAVE_CSIN) |
| | 1407 | #define HAVE_CSIN 1 |
| | 1408 | double complex csin (double complex a); |
| | 1409 | |
| | 1410 | double complex |
| | 1411 | csin (double complex a) |
| | 1412 | { |
| | 1413 | double r, i; |
| | 1414 | double complex v; |
| | 1415 | |
| | 1416 | r = REALPART (a); |
| | 1417 | i = IMAGPART (a); |
| | 1418 | COMPLEX_ASSIGN (v, sin (r) * cosh (i), cos (r) * sinh (i)); |
| | 1419 | return v; |
| | 1420 | } |
| | 1421 | #endif |
| | 1422 | |
| | 1423 | #if !defined(HAVE_CSINL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL) |
| | 1424 | #define HAVE_CSINL 1 |
| | 1425 | long double complex csinl (long double complex a); |
| | 1426 | |
| | 1427 | long double complex |
| | 1428 | csinl (long double complex a) |
| | 1429 | { |
| | 1430 | long double r, i; |
| | 1431 | long double complex v; |
| | 1432 | |
| | 1433 | r = REALPART (a); |
| | 1434 | i = IMAGPART (a); |
| | 1435 | COMPLEX_ASSIGN (v, sinl (r) * coshl (i), cosl (r) * sinhl (i)); |
| | 1436 | return v; |
| | 1437 | } |
| | 1438 | #endif |
| | 1439 | |
| | 1440 | |
| | 1441 | /* cos(a + i b) = cos(a) cosh(b) - i sin(a) sinh(b) */ |
| | 1442 | #if !defined(HAVE_CCOSF) |
| | 1443 | #define HAVE_CCOSF 1 |
| | 1444 | float complex ccosf (float complex a); |
| | 1445 | |
| | 1446 | float complex |
| | 1447 | ccosf (float complex a) |
| | 1448 | { |
| | 1449 | float r, i; |
| | 1450 | float complex v; |
| | 1451 | |
| | 1452 | r = REALPART (a); |
| | 1453 | i = IMAGPART (a); |
| | 1454 | COMPLEX_ASSIGN (v, cosf (r) * coshf (i), - (sinf (r) * sinhf (i))); |
| | 1455 | return v; |
| | 1456 | } |
| | 1457 | #endif |
| | 1458 | |
| | 1459 | #if !defined(HAVE_CCOS) |
| | 1460 | #define HAVE_CCOS 1 |
| | 1461 | double complex ccos (double complex a); |
| | 1462 | |
| | 1463 | double complex |
| | 1464 | ccos (double complex a) |
| | 1465 | { |
| | 1466 | double r, i; |
| | 1467 | double complex v; |
| | 1468 | |
| | 1469 | r = REALPART (a); |
| | 1470 | i = IMAGPART (a); |
| | 1471 | COMPLEX_ASSIGN (v, cos (r) * cosh (i), - (sin (r) * sinh (i))); |
| | 1472 | return v; |
| | 1473 | } |
| | 1474 | #endif |
| | 1475 | |
| | 1476 | #if !defined(HAVE_CCOSL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL) |
| | 1477 | #define HAVE_CCOSL 1 |
| | 1478 | long double complex ccosl (long double complex a); |
| | 1479 | |
| | 1480 | long double complex |
| | 1481 | ccosl (long double complex a) |
| | 1482 | { |
| | 1483 | long double r, i; |
| | 1484 | long double complex v; |
| | 1485 | |
| | 1486 | r = REALPART (a); |
| | 1487 | i = IMAGPART (a); |
| | 1488 | COMPLEX_ASSIGN (v, cosl (r) * coshl (i), - (sinl (r) * sinhl (i))); |
| | 1489 | return v; |
| | 1490 | } |
| | 1491 | #endif |
| | 1492 | |
| | 1493 | |
| | 1494 | /* tan(a + i b) = (tan(a) + i tanh(b)) / (1 - i tan(a) tanh(b)) */ |
| | 1495 | #if !defined(HAVE_CTANF) |
| | 1496 | #define HAVE_CTANF 1 |
| | 1497 | float complex ctanf (float complex a); |
| | 1498 | |
| | 1499 | float complex |
| | 1500 | ctanf (float complex a) |
| | 1501 | { |
| | 1502 | float rt, it; |
| | 1503 | float complex n, d; |
| | 1504 | |
| | 1505 | rt = tanf (REALPART (a)); |
| | 1506 | it = tanhf (IMAGPART (a)); |
| | 1507 | COMPLEX_ASSIGN (n, rt, it); |
| | 1508 | COMPLEX_ASSIGN (d, 1, - (rt * it)); |
| | 1509 | |
| | 1510 | return n / d; |
| | 1511 | } |
| | 1512 | #endif |
| | 1513 | |
| | 1514 | #if !defined(HAVE_CTAN) |
| | 1515 | #define HAVE_CTAN 1 |
| | 1516 | double complex ctan (double complex a); |
| | 1517 | |
| | 1518 | double complex |
| | 1519 | ctan (double complex a) |
| | 1520 | { |
| | 1521 | double rt, it; |
| | 1522 | double complex n, d; |
| | 1523 | |
| | 1524 | rt = tan (REALPART (a)); |
| | 1525 | it = tanh (IMAGPART (a)); |
| | 1526 | COMPLEX_ASSIGN (n, rt, it); |
| | 1527 | COMPLEX_ASSIGN (d, 1, - (rt * it)); |
| | 1528 | |
| | 1529 | return n / d; |
| | 1530 | } |
| | 1531 | #endif |
| | 1532 | |
| | 1533 | #if !defined(HAVE_CTANL) && defined(HAVE_TANL) && defined(HAVE_TANHL) |
| | 1534 | #define HAVE_CTANL 1 |
| | 1535 | long double complex ctanl (long double complex a); |
| | 1536 | |
| | 1537 | long double complex |
| | 1538 | ctanl (long double complex a) |
| | 1539 | { |
| | 1540 | long double rt, it; |
| | 1541 | long double complex n, d; |
| | 1542 | |
| | 1543 | rt = tanl (REALPART (a)); |
| | 1544 | it = tanhl (IMAGPART (a)); |
| | 1545 | COMPLEX_ASSIGN (n, rt, it); |
| | 1546 | COMPLEX_ASSIGN (d, 1, - (rt * it)); |
| | 1547 | |
| | 1548 | return n / d; |
| | 1549 | } |
| | 1550 | #endif |
| | 1551 | |
| | 1552 | |
| | 1553 | /* Complex ASIN. Returns wrongly NaN for infinite arguments. |
| | 1554 | Algorithm taken from Abramowitz & Stegun. */ |
| | 1555 | |
| | 1556 | #if !defined(HAVE_CASINF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF) |
| | 1557 | #define HAVE_CASINF 1 |
| | 1558 | complex float casinf (complex float z); |
| | 1559 | |
| | 1560 | complex float |
| | 1561 | casinf (complex float z) |
| | 1562 | { |
| | 1563 | return -I*clogf (I*z + csqrtf (1.0f-z*z)); |
| | 1564 | } |
| | 1565 | #endif |
| | 1566 | |
| | 1567 | |
| | 1568 | #if !defined(HAVE_CASIN) && defined(HAVE_CLOG) && defined(HAVE_CSQRT) |
| | 1569 | #define HAVE_CASIN 1 |
| | 1570 | complex double casin (complex double z); |
| | 1571 | |
| | 1572 | complex double |
| | 1573 | casin (complex double z) |
| | 1574 | { |
| | 1575 | return -I*clog (I*z + csqrt (1.0-z*z)); |
| | 1576 | } |
| | 1577 | #endif |
| | 1578 | |
| | 1579 | |
| | 1580 | #if !defined(HAVE_CASINL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL) |
| | 1581 | #define HAVE_CASINL 1 |
| | 1582 | complex long double casinl (complex long double z); |
| | 1583 | |
| | 1584 | complex long double |
| | 1585 | casinl (complex long double z) |
| | 1586 | { |
| | 1587 | return -I*clogl (I*z + csqrtl (1.0L-z*z)); |
| | 1588 | } |
| | 1589 | #endif |
| | 1590 | |
| | 1591 | |
| | 1592 | /* Complex ACOS. Returns wrongly NaN for infinite arguments. |
| | 1593 | Algorithm taken from Abramowitz & Stegun. */ |
| | 1594 | |
| | 1595 | #if !defined(HAVE_CACOSF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF) |
| | 1596 | #define HAVE_CACOSF 1 |
| | 1597 | complex float cacosf (complex float z); |
| | 1598 | |
| | 1599 | complex float |
| | 1600 | cacosf (complex float z) |
| | 1601 | { |
| | 1602 | return -I*clogf (z + I*csqrtf (1.0f-z*z)); |
| | 1603 | } |
| | 1604 | #endif |
| | 1605 | |
| | 1606 | |
| | 1607 | #if !defined(HAVE_CACOS) && defined(HAVE_CLOG) && defined(HAVE_CSQRT) |
| | 1608 | #define HAVE_CACOS 1 |
| | 1609 | complex double cacos (complex double z); |
| | 1610 | |
| | 1611 | complex double |
| | 1612 | cacos (complex double z) |
| | 1613 | { |
| | 1614 | return -I*clog (z + I*csqrt (1.0-z*z)); |
| | 1615 | } |
| | 1616 | #endif |
| | 1617 | |
| | 1618 | |
| | 1619 | #if !defined(HAVE_CACOSL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL) |
| | 1620 | #define HAVE_CACOSL 1 |
| | 1621 | complex long double cacosl (complex long double z); |
| | 1622 | |
| | 1623 | complex long double |
| | 1624 | cacosl (complex long double z) |
| | 1625 | { |
| | 1626 | return -I*clogl (z + I*csqrtl (1.0L-z*z)); |
| | 1627 | } |
| | 1628 | #endif |
| | 1629 | |
| | 1630 | |
| | 1631 | /* Complex ATAN. Returns wrongly NaN for infinite arguments. |
| | 1632 | Algorithm taken from Abramowitz & Stegun. */ |
| | 1633 | |
| | 1634 | #if !defined(HAVE_CATANF) && defined(HAVE_CLOGF) |
| | 1635 | #define HAVE_CACOSF 1 |
| | 1636 | complex float catanf (complex float z); |
| | 1637 | |
| | 1638 | complex float |
| | 1639 | catanf (complex float z) |
| | 1640 | { |
| | 1641 | return I*clogf ((I+z)/(I-z))/2.0f; |
| | 1642 | } |
| | 1643 | #endif |
| | 1644 | |
| | 1645 | |
| | 1646 | #if !defined(HAVE_CATAN) && defined(HAVE_CLOG) |
| | 1647 | #define HAVE_CACOS 1 |
| | 1648 | complex double catan (complex double z); |
| | 1649 | |
| | 1650 | complex double |
| | 1651 | catan (complex double z) |
| | 1652 | { |
| | 1653 | return I*clog ((I+z)/(I-z))/2.0; |
| | 1654 | } |
| | 1655 | #endif |
| | 1656 | |
| | 1657 | |
| | 1658 | #if !defined(HAVE_CATANL) && defined(HAVE_CLOGL) |
| | 1659 | #define HAVE_CACOSL 1 |
| | 1660 | complex long double catanl (complex long double z); |
| | 1661 | |
| | 1662 | complex long double |
| | 1663 | catanl (complex long double z) |
| | 1664 | { |
| | 1665 | return I*clogl ((I+z)/(I-z))/2.0L; |
| | 1666 | } |
| | 1667 | #endif |
| | 1668 | |
| | 1669 | |
| | 1670 | /* Complex ASINH. Returns wrongly NaN for infinite arguments. |
| | 1671 | Algorithm taken from Abramowitz & Stegun. */ |
| | 1672 | |
| | 1673 | #if !defined(HAVE_CASINHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF) |
| | 1674 | #define HAVE_CASINHF 1 |
| | 1675 | complex float casinhf (complex float z); |
| | 1676 | |
| | 1677 | complex float |
| | 1678 | casinhf (complex float z) |
| | 1679 | { |
| | 1680 | return clogf (z + csqrtf (z*z+1.0f)); |
| | 1681 | } |
| | 1682 | #endif |
| | 1683 | |
| | 1684 | |
| | 1685 | #if !defined(HAVE_CASINH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT) |
| | 1686 | #define HAVE_CASINH 1 |
| | 1687 | complex double casinh (complex double z); |
| | 1688 | |
| | 1689 | complex double |
| | 1690 | casinh (complex double z) |
| | 1691 | { |
| | 1692 | return clog (z + csqrt (z*z+1.0)); |
| | 1693 | } |
| | 1694 | #endif |
| | 1695 | |
| | 1696 | |
| | 1697 | #if !defined(HAVE_CASINHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL) |
| | 1698 | #define HAVE_CASINHL 1 |
| | 1699 | complex long double casinhl (complex long double z); |
| | 1700 | |
| | 1701 | complex long double |
| | 1702 | casinhl (complex long double z) |
| | 1703 | { |
| | 1704 | return clogl (z + csqrtl (z*z+1.0L)); |
| | 1705 | } |
| | 1706 | #endif |
| | 1707 | |
| | 1708 | |
| | 1709 | /* Complex ACOSH. Returns wrongly NaN for infinite arguments. |
| | 1710 | Algorithm taken from Abramowitz & Stegun. */ |
| | 1711 | |
| | 1712 | #if !defined(HAVE_CACOSHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF) |
| | 1713 | #define HAVE_CACOSHF 1 |
| | 1714 | complex float cacoshf (complex float z); |
| | 1715 | |
| | 1716 | complex float |
| | 1717 | cacoshf (complex float z) |
| | 1718 | { |
| | 1719 | return clogf (z + csqrtf (z-1.0f) * csqrtf (z+1.0f)); |
| | 1720 | } |
| | 1721 | #endif |
| | 1722 | |
| | 1723 | |
| | 1724 | #if !defined(HAVE_CACOSH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT) |
| | 1725 | #define HAVE_CACOSH 1 |
| | 1726 | complex double cacosh (complex double z); |
| | 1727 | |
| | 1728 | complex double |
| | 1729 | cacosh (complex double z) |
| | 1730 | { |
| | 1731 | return clog (z + csqrt (z-1.0) * csqrt (z+1.0)); |
| | 1732 | } |
| | 1733 | #endif |
| | 1734 | |
| | 1735 | |
| | 1736 | #if !defined(HAVE_CACOSHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL) |
| | 1737 | #define HAVE_CACOSHL 1 |
| | 1738 | complex long double cacoshl (complex long double z); |
| | 1739 | |
| | 1740 | complex long double |
| | 1741 | cacoshl (complex long double z) |
| | 1742 | { |
| | 1743 | return clogl (z + csqrtl (z-1.0L) * csqrtl (z+1.0L)); |
| | 1744 | } |
| | 1745 | #endif |
| | 1746 | |
| | 1747 | |
| | 1748 | /* Complex ATANH. Returns wrongly NaN for infinite arguments. |
| | 1749 | Algorithm taken from Abramowitz & Stegun. */ |
| | 1750 | |
| | 1751 | #if !defined(HAVE_CATANHF) && defined(HAVE_CLOGF) |
| | 1752 | #define HAVE_CATANHF 1 |
| | 1753 | complex float catanhf (complex float z); |
| | 1754 | |
| | 1755 | complex float |
| | 1756 | catanhf (complex float z) |
| | 1757 | { |
| | 1758 | return clogf ((1.0f+z)/(1.0f-z))/2.0f; |
| | 1759 | } |
| | 1760 | #endif |
| | 1761 | |
| | 1762 | |
| | 1763 | #if !defined(HAVE_CATANH) && defined(HAVE_CLOG) |
| | 1764 | #define HAVE_CATANH 1 |
| | 1765 | complex double catanh (complex double z); |
| | 1766 | |
| | 1767 | complex double |
| | 1768 | catanh (complex double z) |
| | 1769 | { |
| | 1770 | return clog ((1.0+z)/(1.0-z))/2.0; |
| | 1771 | } |
| | 1772 | #endif |
| | 1773 | |
| | 1774 | #if !defined(HAVE_CATANHL) && defined(HAVE_CLOGL) |
| | 1775 | #define HAVE_CATANHL 1 |
| | 1776 | complex long double catanhl (complex long double z); |
| | 1777 | |
| | 1778 | complex long double |
| | 1779 | catanhl (complex long double z) |
| | 1780 | { |
| | 1781 | return clogl ((1.0L+z)/(1.0L-z))/2.0L; |
| | 1782 | } |
| | 1783 | #endif |
| | 1784 | |
| | 1785 | |
| | 1786 | #if !defined(HAVE_TGAMMA) |
| | 1787 | #define HAVE_TGAMMA 1 |
| | 1788 | double tgamma (double); |
| | 1789 | |
| | 1790 | /* Fallback tgamma() function. Uses the algorithm from |
| | 1791 | http://www.netlib.org/specfun/gamma and references therein. */ |
| | 1792 | |
| | 1793 | #undef SQRTPI |
| | 1794 | #define SQRTPI 0.9189385332046727417803297 |
| | 1795 | |
| | 1796 | #undef PI |
| | 1797 | #define PI 3.1415926535897932384626434 |
| | 1798 | |
| | 1799 | double |
| | 1800 | tgamma (double x) |
| | 1801 | { |
| | 1802 | int i, n, parity; |
| | 1803 | double fact, res, sum, xden, xnum, y, y1, ysq, z; |
| | 1804 | |
| | 1805 | static double p[8] = { |
| | 1806 | -1.71618513886549492533811e0, 2.47656508055759199108314e1, |
| | 1807 | -3.79804256470945635097577e2, 6.29331155312818442661052e2, |
| | 1808 | 8.66966202790413211295064e2, -3.14512729688483675254357e4, |
| | 1809 | -3.61444134186911729807069e4, 6.64561438202405440627855e4 }; |
| | 1810 | |
| | 1811 | static double q[8] = { |
| | 1812 | -3.08402300119738975254353e1, 3.15350626979604161529144e2, |
| | 1813 | -1.01515636749021914166146e3, -3.10777167157231109440444e3, |
| | 1814 | 2.25381184209801510330112e4, 4.75584627752788110767815e3, |
| | 1815 | -1.34659959864969306392456e5, -1.15132259675553483497211e5 }; |
| | 1816 | |
| | 1817 | static double c[7] = { -1.910444077728e-03, |
| | 1818 | 8.4171387781295e-04, -5.952379913043012e-04, |
| | 1819 | 7.93650793500350248e-04, -2.777777777777681622553e-03, |
| | 1820 | 8.333333333333333331554247e-02, 5.7083835261e-03 }; |
| | 1821 | |
| | 1822 | static const double xminin = 2.23e-308; |
| | 1823 | static const double xbig = 171.624; |
| | 1824 | static const double xnan = __builtin_nan ("0x0"), xinf = __builtin_inf (); |
| | 1825 | static double eps = 0; |
| | 1826 | |
| | 1827 | if (eps == 0) |
| | 1828 | eps = nextafter (1., 2.) - 1.; |
| | 1829 | |
| | 1830 | parity = 0; |
| | 1831 | fact = 1; |
| | 1832 | n = 0; |
| | 1833 | y = x; |
| | 1834 | |
| | 1835 | if (isnan (x)) |
| | 1836 | return x; |
| | 1837 | |
| | 1838 | if (y <= 0) |
| | 1839 | { |
| | 1840 | y = -x; |
| | 1841 | y1 = trunc (y); |
| | 1842 | res = y - y1; |
| | 1843 | |
| | 1844 | if (res != 0) |
| | 1845 | { |
| | 1846 | if (y1 != trunc (y1*0.5l)*2) |
| | 1847 | parity = 1; |
| | 1848 | fact = -PI / sin (PI*res); |
| | 1849 | y = y + 1; |
| | 1850 | } |
| | 1851 | else |
| | 1852 | return x == 0 ? copysign (xinf, x) : xnan; |
| | 1853 | } |
| | 1854 | |
| | 1855 | if (y < eps) |
| | 1856 | { |
| | 1857 | if (y >= xminin) |
| | 1858 | res = 1 / y; |
| | 1859 | else |
| | 1860 | return xinf; |
| | 1861 | } |
| | 1862 | else if (y < 13) |
| | 1863 | { |
| | 1864 | y1 = y; |
| | 1865 | if (y < 1) |
| | 1866 | { |
| | 1867 | z = y; |
| | 1868 | y = y + 1; |
| | 1869 | } |
| | 1870 | else |
| | 1871 | { |
| | 1872 | n = (int)y - 1; |
| | 1873 | y = y - n; |
| | 1874 | z = y - 1; |
| | 1875 | } |
| | 1876 | |
| | 1877 | xnum = 0; |
| | 1878 | xden = 1; |
| | 1879 | for (i = 0; i < 8; i++) |
| | 1880 | { |
| | 1881 | xnum = (xnum + p[i]) * z; |
| | 1882 | xden = xden * z + q[i]; |
| | 1883 | } |
| | 1884 | |
| | 1885 | res = xnum / xden + 1; |
| | 1886 | |
| | 1887 | if (y1 < y) |
| | 1888 | res = res / y1; |
| | 1889 | else if (y1 > y) |
| | 1890 | for (i = 1; i <= n; i++) |
| | 1891 | { |
| | 1892 | res = res * y; |
| | 1893 | y = y + 1; |
| | 1894 | } |
| | 1895 | } |
| | 1896 | else |
| | 1897 | { |
| | 1898 | if (y < xbig) |
| | 1899 | { |
| | 1900 | ysq = y * y; |
| | 1901 | sum = c[6]; |
| | 1902 | for (i = 0; i < 6; i++) |
| | 1903 | sum = sum / ysq + c[i]; |
| | 1904 | |
| | 1905 | sum = sum/y - y + SQRTPI; |
| | 1906 | sum = sum + (y - 0.5) * log (y); |
| | 1907 | res = exp (sum); |
| | 1908 | } |
| | 1909 | else |
| | 1910 | return x < 0 ? xnan : xinf; |
| | 1911 | } |
| | 1912 | |
| | 1913 | if (parity) |
| | 1914 | res = -res; |
| | 1915 | if (fact != 1) |
| | 1916 | res = fact / res; |
| | 1917 | |
| | 1918 | return res; |
| | 1919 | } |
| | 1920 | #endif |
| | 1921 | |
| | 1922 | |
| | 1923 | |
| | 1924 | #if !defined(HAVE_LGAMMA) |
| | 1925 | #define HAVE_LGAMMA 1 |
| | 1926 | double lgamma (double); |
| | 1927 | |
| | 1928 | /* Fallback lgamma() function. Uses the algorithm from |
| | 1929 | http://www.netlib.org/specfun/algama and references therein, |
| | 1930 | except for negative arguments (where netlib would return +Inf) |
| | 1931 | where we use the following identity: |
| | 1932 | lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y) |
| | 1933 | */ |
| | 1934 | |
| | 1935 | double |
| | 1936 | lgamma (double y) |
| | 1937 | { |
| | 1938 | |
| | 1939 | #undef SQRTPI |
| | 1940 | #define SQRTPI 0.9189385332046727417803297 |
| | 1941 | |
| | 1942 | #undef PI |
| | 1943 | #define PI 3.1415926535897932384626434 |
| | 1944 | |
| | 1945 | #define PNT68 0.6796875 |
| | 1946 | #define D1 -0.5772156649015328605195174 |
| | 1947 | #define D2 0.4227843350984671393993777 |
| | 1948 | #define D4 1.791759469228055000094023 |
| | 1949 | |
| | 1950 | static double p1[8] = { |
| | 1951 | 4.945235359296727046734888e0, 2.018112620856775083915565e2, |
| | 1952 | 2.290838373831346393026739e3, 1.131967205903380828685045e4, |
| | 1953 | 2.855724635671635335736389e4, 3.848496228443793359990269e4, |
| | 1954 | 2.637748787624195437963534e4, 7.225813979700288197698961e3 }; |
| | 1955 | static double q1[8] = { |
| | 1956 | 6.748212550303777196073036e1, 1.113332393857199323513008e3, |
| | 1957 | 7.738757056935398733233834e3, 2.763987074403340708898585e4, |
| | 1958 | 5.499310206226157329794414e4, 6.161122180066002127833352e4, |
| | 1959 | 3.635127591501940507276287e4, 8.785536302431013170870835e3 }; |
| | 1960 | static double p2[8] = { |
| | 1961 | 4.974607845568932035012064e0, 5.424138599891070494101986e2, |
| | 1962 | 1.550693864978364947665077e4, 1.847932904445632425417223e5, |
| | 1963 | 1.088204769468828767498470e6, 3.338152967987029735917223e6, |
| | 1964 | 5.106661678927352456275255e6, 3.074109054850539556250927e6 }; |
| | 1965 | static double q2[8] = { |
| | 1966 | 1.830328399370592604055942e2, 7.765049321445005871323047e3, |
| | 1967 | 1.331903827966074194402448e5, 1.136705821321969608938755e6, |
| | 1968 | 5.267964117437946917577538e6, 1.346701454311101692290052e7, |
| | 1969 | 1.782736530353274213975932e7, 9.533095591844353613395747e6 }; |
| | 1970 | static double p4[8] = { |
| | 1971 | 1.474502166059939948905062e4, 2.426813369486704502836312e6, |
| | 1972 | 1.214755574045093227939592e8, 2.663432449630976949898078e9, |
| | 1973 | 2.940378956634553899906876e10, 1.702665737765398868392998e11, |
| | 1974 | 4.926125793377430887588120e11, 5.606251856223951465078242e11 }; |
| | 1975 | static double q4[8] = { |
| | 1976 | 2.690530175870899333379843e3, 6.393885654300092398984238e5, |
| | 1977 | 4.135599930241388052042842e7, 1.120872109616147941376570e9, |
| | 1978 | 1.488613728678813811542398e10, 1.016803586272438228077304e11, |
| | 1979 | 3.417476345507377132798597e11, 4.463158187419713286462081e11 }; |
| | 1980 | static double c[7] = { |
| | 1981 | -1.910444077728e-03, 8.4171387781295e-04, |
| | 1982 | -5.952379913043012e-04, 7.93650793500350248e-04, |
| | 1983 | -2.777777777777681622553e-03, 8.333333333333333331554247e-02, |
| | 1984 | 5.7083835261e-03 }; |
| | 1985 | |
| | 1986 | static double xbig = 2.55e305, xinf = __builtin_inf (), eps = 0, |
| | 1987 | frtbig = 2.25e76; |
| | 1988 | |
| | 1989 | int i; |
| | 1990 | double corr, res, xden, xm1, xm2, xm4, xnum, ysq; |
| | 1991 | |
| | 1992 | if (eps == 0) |
| | 1993 | eps = __builtin_nextafter (1., 2.) - 1.; |
| | 1994 | |
| | 1995 | if ((y > 0) && (y <= xbig)) |
| | 1996 | { |
| | 1997 | if (y <= eps) |
| | 1998 | res = -log (y); |
| | 1999 | else if (y <= 1.5) |
| | 2000 | { |
| | 2001 | if (y < PNT68) |
| | 2002 | { |
| | 2003 | corr = -log (y); |
| | 2004 | xm1 = y; |
| | 2005 | } |
| | 2006 | else |
| | 2007 | { |
| | 2008 | corr = 0; |
| | 2009 | xm1 = (y - 0.5) - 0.5; |
| | 2010 | } |
| | 2011 | |
| | 2012 | if ((y <= 0.5) || (y >= PNT68)) |
| | 2013 | { |
| | 2014 | xden = 1; |
| | 2015 | xnum = 0; |
| | 2016 | for (i = 0; i < 8; i++) |
| | 2017 | { |
| | 2018 | xnum = xnum*xm1 + p1[i]; |
| | 2019 | xden = xden*xm1 + q1[i]; |
| | 2020 | } |
| | 2021 | res = corr + (xm1 * (D1 + xm1*(xnum/xden))); |
| | 2022 | } |
| | 2023 | else |
| | 2024 | { |
| | 2025 | xm2 = (y - 0.5) - 0.5; |
| | 2026 | xden = 1; |
| | 2027 | xnum = 0; |
| | 2028 | for (i = 0; i < 8; i++) |
| | 2029 | { |
| | 2030 | xnum = xnum*xm2 + p2[i]; |
| | 2031 | xden = xden*xm2 + q2[i]; |
| | 2032 | } |
| | 2033 | res = corr + xm2 * (D2 + xm2*(xnum/xden)); |
| | 2034 | } |
| | 2035 | } |
| | 2036 | else if (y <= 4) |
| | 2037 | { |
| | 2038 | xm2 = y - 2; |
| | 2039 | xden = 1; |
| | 2040 | xnum = 0; |
| | 2041 | for (i = 0; i < 8; i++) |
| | 2042 | { |
| | 2043 | xnum = xnum*xm2 + p2[i]; |
| | 2044 | xden = xden*xm2 + q2[i]; |
| | 2045 | } |
| | 2046 | res = xm2 * (D2 + xm2*(xnum/xden)); |
| | 2047 | } |
| | 2048 | else if (y <= 12) |
| | 2049 | { |
| | 2050 | xm4 = y - 4; |
| | 2051 | xden = -1; |
| | 2052 | xnum = 0; |
| | 2053 | for (i = 0; i < 8; i++) |
| | 2054 | { |
| | 2055 | xnum = xnum*xm4 + p4[i]; |
| | 2056 | xden = xden*xm4 + q4[i]; |
| | 2057 | } |
| | 2058 | res = D4 + xm4*(xnum/xden); |
| | 2059 | } |
| | 2060 | else |
| | 2061 | { |
| | 2062 | res = 0; |
| | 2063 | if (y <= frtbig) |
| | 2064 | { |
| | 2065 | res = c[6]; |
| | 2066 | ysq = y * y; |
| | 2067 | for (i = 0; i < 6; i++) |
| | 2068 | res = res / ysq + c[i]; |
| | 2069 | } |
| | 2070 | res = res/y; |
| | 2071 | corr = log (y); |
| | 2072 | res = res + SQRTPI - 0.5*corr; |
| | 2073 | res = res + y*(corr-1); |
| | 2074 | } |
| | 2075 | } |
| | 2076 | else if (y < 0 && __builtin_floor (y) != y) |
| | 2077 | { |
| | 2078 | /* lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y) |
| | 2079 | For abs(y) very close to zero, we use a series expansion to |
| | 2080 | the first order in y to avoid overflow. */ |
| | 2081 | if (y > -1.e-100) |
| | 2082 | res = -2 * log (fabs (y)) - lgamma (-y); |
| | 2083 | else |
| | 2084 | res = log (PI / fabs (y * sin (PI * y))) - lgamma (-y); |
| | 2085 | } |
| | 2086 | else |
| | 2087 | res = xinf; |
| | 2088 | |
| | 2089 | return res; |
| | 2090 | } |
| | 2091 | #endif |
| | 2092 | |
| | 2093 | |
| | 2094 | #if defined(HAVE_TGAMMA) && !defined(HAVE_TGAMMAF) |
| | 2095 | #define HAVE_TGAMMAF 1 |
| | 2096 | float tgammaf (float); |
| | 2097 | |
| | 2098 | float |
| | 2099 | tgammaf (float x) |
| | 2100 | { |
| | 2101 | return (float) tgamma ((double) x); |
| | 2102 | } |
| | 2103 | #endif |
| | 2104 | |
| | 2105 | #if defined(HAVE_LGAMMA) && !defined(HAVE_LGAMMAF) |
| | 2106 | #define HAVE_LGAMMAF 1 |
| | 2107 | float lgammaf (float); |
| | 2108 | |
| | 2109 | float |
| | 2110 | lgammaf (float x) |
| | 2111 | { |
| | 2112 | return (float) lgamma ((double) x); |
| | 2113 | } |
| | 2114 | #endif |
projects / thirdparty / gcc.git / blame_incremental