diff options
Diffstat (limited to 'libgfortran/intrinsics/c99_functions.c')
| -rw-r--r-- | libgfortran/intrinsics/c99_functions.c | 2135 |
1 files changed, 2135 insertions, 0 deletions
diff --git a/libgfortran/intrinsics/c99_functions.c b/libgfortran/intrinsics/c99_functions.c new file mode 100644 index 000000000..9ba5544a0 --- /dev/null +++ b/libgfortran/intrinsics/c99_functions.c @@ -0,0 +1,2135 @@ +/* Implementation of various C99 functions + Copyright (C) 2004, 2009, 2010 Free Software Foundation, Inc. + +This file is part of the GNU Fortran 95 runtime library (libgfortran). + +Libgfortran is free software; you can redistribute it and/or +modify it under the terms of the GNU General Public +License as published by the Free Software Foundation; either +version 3 of the License, or (at your option) any later version. + +Libgfortran is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +Under Section 7 of GPL version 3, you are granted additional +permissions described in the GCC Runtime Library Exception, version +3.1, as published by the Free Software Foundation. + +You should have received a copy of the GNU General Public License and +a copy of the GCC Runtime Library Exception along with this program; +see the files COPYING3 and COPYING.RUNTIME respectively. If not, see +<http://www.gnu.org/licenses/>. */ + +#include "config.h" + +#define C99_PROTOS_H WE_DONT_WANT_PROTOS_NOW +#include "libgfortran.h" + +/* IRIX's <math.h> declares a non-C99 compliant implementation of cabs, + which takes two floating point arguments instead of a single complex. + If <complex.h> is missing this prevents building of c99_functions.c. + To work around this we redirect cabs{,f,l} calls to __gfc_cabs{,f,l}. */ + +#if defined(__sgi__) && !defined(HAVE_COMPLEX_H) +#undef HAVE_CABS +#undef HAVE_CABSF +#undef HAVE_CABSL +#define cabs __gfc_cabs +#define cabsf __gfc_cabsf +#define cabsl __gfc_cabsl +#endif + +/* Tru64's <math.h> declares a non-C99 compliant implementation of cabs, + which takes two floating point arguments instead of a single complex. + To work around this we redirect cabs{,f,l} calls to __gfc_cabs{,f,l}. */ + +#ifdef __osf__ +#undef HAVE_CABS +#undef HAVE_CABSF +#undef HAVE_CABSL +#define cabs __gfc_cabs +#define cabsf __gfc_cabsf +#define cabsl __gfc_cabsl +#endif + +/* On a C99 system "I" (with I*I = -1) should be defined in complex.h; + if not, we define a fallback version here. */ +#ifndef I +# if defined(_Imaginary_I) +# define I _Imaginary_I +# elif defined(_Complex_I) +# define I _Complex_I +# else +# define I (1.0fi) +# endif +#endif + +/* Prototypes are included to silence -Wstrict-prototypes + -Wmissing-prototypes. */ + + +/* Wrappers for systems without the various C99 single precision Bessel + functions. */ + +#if defined(HAVE_J0) && ! defined(HAVE_J0F) +#define HAVE_J0F 1 +float j0f (float); + +float +j0f (float x) +{ + return (float) j0 ((double) x); +} +#endif + +#if defined(HAVE_J1) && !defined(HAVE_J1F) +#define HAVE_J1F 1 +float j1f (float); + +float j1f (float x) +{ + return (float) j1 ((double) x); +} +#endif + +#if defined(HAVE_JN) && !defined(HAVE_JNF) +#define HAVE_JNF 1 +float jnf (int, float); + +float +jnf (int n, float x) +{ + return (float) jn (n, (double) x); +} +#endif + +#if defined(HAVE_Y0) && !defined(HAVE_Y0F) +#define HAVE_Y0F 1 +float y0f (float); + +float +y0f (float x) +{ + return (float) y0 ((double) x); +} +#endif + +#if defined(HAVE_Y1) && !defined(HAVE_Y1F) +#define HAVE_Y1F 1 +float y1f (float); + +float +y1f (float x) +{ + return (float) y1 ((double) x); +} +#endif + +#if defined(HAVE_YN) && !defined(HAVE_YNF) +#define HAVE_YNF 1 +float ynf (int, float); + +float +ynf (int n, float x) +{ + return (float) yn (n, (double) x); +} +#endif + + +/* Wrappers for systems without the C99 erff() and erfcf() functions. */ + +#if defined(HAVE_ERF) && !defined(HAVE_ERFF) +#define HAVE_ERFF 1 +float erff (float); + +float +erff (float x) +{ + return (float) erf ((double) x); +} +#endif + +#if defined(HAVE_ERFC) && !defined(HAVE_ERFCF) +#define HAVE_ERFCF 1 +float erfcf (float); + +float +erfcf (float x) +{ + return (float) erfc ((double) x); +} +#endif + + +#ifndef HAVE_ACOSF +#define HAVE_ACOSF 1 +float acosf (float x); + +float +acosf (float x) +{ + return (float) acos (x); +} +#endif + +#if HAVE_ACOSH && !HAVE_ACOSHF +float acoshf (float x); + +float +acoshf (float x) +{ + return (float) acosh ((double) x); +} +#endif + +#ifndef HAVE_ASINF +#define HAVE_ASINF 1 +float asinf (float x); + +float +asinf (float x) +{ + return (float) asin (x); +} +#endif + +#if HAVE_ASINH && !HAVE_ASINHF +float asinhf (float x); + +float +asinhf (float x) +{ + return (float) asinh ((double) x); +} +#endif + +#ifndef HAVE_ATAN2F +#define HAVE_ATAN2F 1 +float atan2f (float y, float x); + +float +atan2f (float y, float x) +{ + return (float) atan2 (y, x); +} +#endif + +#ifndef HAVE_ATANF +#define HAVE_ATANF 1 +float atanf (float x); + +float +atanf (float x) +{ + return (float) atan (x); +} +#endif + +#if HAVE_ATANH && !HAVE_ATANHF +float atanhf (float x); + +float +atanhf (float x) +{ + return (float) atanh ((double) x); +} +#endif + +#ifndef HAVE_CEILF +#define HAVE_CEILF 1 +float ceilf (float x); + +float +ceilf (float x) +{ + return (float) ceil (x); +} +#endif + +#ifndef HAVE_COPYSIGNF +#define HAVE_COPYSIGNF 1 +float copysignf (float x, float y); + +float +copysignf (float x, float y) +{ + return (float) copysign (x, y); +} +#endif + +#ifndef HAVE_COSF +#define HAVE_COSF 1 +float cosf (float x); + +float +cosf (float x) +{ + return (float) cos (x); +} +#endif + +#ifndef HAVE_COSHF +#define HAVE_COSHF 1 +float coshf (float x); + +float +coshf (float x) +{ + return (float) cosh (x); +} +#endif + +#ifndef HAVE_EXPF +#define HAVE_EXPF 1 +float expf (float x); + +float +expf (float x) +{ + return (float) exp (x); +} +#endif + +#ifndef HAVE_FABSF +#define HAVE_FABSF 1 +float fabsf (float x); + +float +fabsf (float x) +{ + return (float) fabs (x); +} +#endif + +#ifndef HAVE_FLOORF +#define HAVE_FLOORF 1 +float floorf (float x); + +float +floorf (float x) +{ + return (float) floor (x); +} +#endif + +#ifndef HAVE_FMODF +#define HAVE_FMODF 1 +float fmodf (float x, float y); + +float +fmodf (float x, float y) +{ + return (float) fmod (x, y); +} +#endif + +#ifndef HAVE_FREXPF +#define HAVE_FREXPF 1 +float frexpf (float x, int *exp); + +float +frexpf (float x, int *exp) +{ + return (float) frexp (x, exp); +} +#endif + +#ifndef HAVE_HYPOTF +#define HAVE_HYPOTF 1 +float hypotf (float x, float y); + +float +hypotf (float x, float y) +{ + return (float) hypot (x, y); +} +#endif + +#ifndef HAVE_LOGF +#define HAVE_LOGF 1 +float logf (float x); + +float +logf (float x) +{ + return (float) log (x); +} +#endif + +#ifndef HAVE_LOG10F +#define HAVE_LOG10F 1 +float log10f (float x); + +float +log10f (float x) +{ + return (float) log10 (x); +} +#endif + +#ifndef HAVE_SCALBN +#define HAVE_SCALBN 1 +double scalbn (double x, int y); + +double +scalbn (double x, int y) +{ +#if (FLT_RADIX == 2) && defined(HAVE_LDEXP) + return ldexp (x, y); +#else + return x * pow (FLT_RADIX, y); +#endif +} +#endif + +#ifndef HAVE_SCALBNF +#define HAVE_SCALBNF 1 +float scalbnf (float x, int y); + +float +scalbnf (float x, int y) +{ + return (float) scalbn (x, y); +} +#endif + +#ifndef HAVE_SINF +#define HAVE_SINF 1 +float sinf (float x); + +float +sinf (float x) +{ + return (float) sin (x); +} +#endif + +#ifndef HAVE_SINHF +#define HAVE_SINHF 1 +float sinhf (float x); + +float +sinhf (float x) +{ + return (float) sinh (x); +} +#endif + +#ifndef HAVE_SQRTF +#define HAVE_SQRTF 1 +float sqrtf (float x); + +float +sqrtf (float x) +{ + return (float) sqrt (x); +} +#endif + +#ifndef HAVE_TANF +#define HAVE_TANF 1 +float tanf (float x); + +float +tanf (float x) +{ + return (float) tan (x); +} +#endif + +#ifndef HAVE_TANHF +#define HAVE_TANHF 1 +float tanhf (float x); + +float +tanhf (float x) +{ + return (float) tanh (x); +} +#endif + +#ifndef HAVE_TRUNC +#define HAVE_TRUNC 1 +double trunc (double x); + +double +trunc (double x) +{ + if (!isfinite (x)) + return x; + + if (x < 0.0) + return - floor (-x); + else + return floor (x); +} +#endif + +#ifndef HAVE_TRUNCF +#define HAVE_TRUNCF 1 +float truncf (float x); + +float +truncf (float x) +{ + return (float) trunc (x); +} +#endif + +#ifndef HAVE_NEXTAFTERF +#define HAVE_NEXTAFTERF 1 +/* This is a portable implementation of nextafterf that is intended to be + independent of the floating point format or its in memory representation. + This implementation works correctly with denormalized values. */ +float nextafterf (float x, float y); + +float +nextafterf (float x, float y) +{ + /* This variable is marked volatile to avoid excess precision problems + on some platforms, including IA-32. */ + volatile float delta; + float absx, denorm_min; + + if (isnan (x) || isnan (y)) + return x + y; + if (x == y) + return x; + if (!isfinite (x)) + return x > 0 ? __FLT_MAX__ : - __FLT_MAX__; + + /* absx = fabsf (x); */ + absx = (x < 0.0) ? -x : x; + + /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals. */ + if (__FLT_DENORM_MIN__ == 0.0f) + denorm_min = __FLT_MIN__; + else + denorm_min = __FLT_DENORM_MIN__; + + if (absx < __FLT_MIN__) + delta = denorm_min; + else + { + float frac; + int exp; + + /* Discard the fraction from x. */ + frac = frexpf (absx, &exp); + delta = scalbnf (0.5f, exp); + + /* Scale x by the epsilon of the representation. By rights we should + have been able to combine this with scalbnf, but some targets don't + get that correct with denormals. */ + delta *= __FLT_EPSILON__; + + /* If we're going to be reducing the absolute value of X, and doing so + would reduce the exponent of X, then the delta to be applied is + one exponent smaller. */ + if (frac == 0.5f && (y < x) == (x > 0)) + delta *= 0.5f; + + /* If that underflows to zero, then we're back to the minimum. */ + if (delta == 0.0f) + delta = denorm_min; + } + + if (y < x) + delta = -delta; + + return x + delta; +} +#endif + + +#if !defined(HAVE_POWF) || defined(HAVE_BROKEN_POWF) +#ifndef HAVE_POWF +#define HAVE_POWF 1 +#endif +float powf (float x, float y); + +float +powf (float x, float y) +{ + return (float) pow (x, y); +} +#endif + + +/* Algorithm by Steven G. Kargl. */ + +#if !defined(HAVE_ROUNDL) +#define HAVE_ROUNDL 1 +long double roundl (long double x); + +#if defined(HAVE_CEILL) +/* Round to nearest integral value. If the argument is halfway between two + integral values then round away from zero. */ + +long double +roundl (long double x) +{ + long double t; + if (!isfinite (x)) + return (x); + + if (x >= 0.0) + { + t = ceill (x); + if (t - x > 0.5) + t -= 1.0; + return (t); + } + else + { + t = ceill (-x); + if (t + x > 0.5) + t -= 1.0; + return (-t); + } +} +#else + +/* Poor version of roundl for system that don't have ceill. */ +long double +roundl (long double x) +{ + if (x > DBL_MAX || x < -DBL_MAX) + { +#ifdef HAVE_NEXTAFTERL + long double prechalf = nextafterl (0.5L, LDBL_MAX); +#else + static long double prechalf = 0.5L; +#endif + return (GFC_INTEGER_LARGEST) (x + (x > 0 ? prechalf : -prechalf)); + } + else + /* Use round(). */ + return round ((double) x); +} + +#endif +#endif + +#ifndef HAVE_ROUND +#define HAVE_ROUND 1 +/* Round to nearest integral value. If the argument is halfway between two + integral values then round away from zero. */ +double round (double x); + +double +round (double x) +{ + double t; + if (!isfinite (x)) + return (x); + + if (x >= 0.0) + { + t = floor (x); + if (t - x <= -0.5) + t += 1.0; + return (t); + } + else + { + t = floor (-x); + if (t + x <= -0.5) + t += 1.0; + return (-t); + } +} +#endif + +#ifndef HAVE_ROUNDF +#define HAVE_ROUNDF 1 +/* Round to nearest integral value. If the argument is halfway between two + integral values then round away from zero. */ +float roundf (float x); + +float +roundf (float x) +{ + float t; + if (!isfinite (x)) + return (x); + + if (x >= 0.0) + { + t = floorf (x); + if (t - x <= -0.5) + t += 1.0; + return (t); + } + else + { + t = floorf (-x); + if (t + x <= -0.5) + t += 1.0; + return (-t); + } +} +#endif + + +/* lround{f,,l} and llround{f,,l} functions. */ + +#if !defined(HAVE_LROUNDF) && defined(HAVE_ROUNDF) +#define HAVE_LROUNDF 1 +long int lroundf (float x); + +long int +lroundf (float x) +{ + return (long int) roundf (x); +} +#endif + +#if !defined(HAVE_LROUND) && defined(HAVE_ROUND) +#define HAVE_LROUND 1 +long int lround (double x); + +long int +lround (double x) +{ + return (long int) round (x); +} +#endif + +#if !defined(HAVE_LROUNDL) && defined(HAVE_ROUNDL) +#define HAVE_LROUNDL 1 +long int lroundl (long double x); + +long int +lroundl (long double x) +{ + return (long long int) roundl (x); +} +#endif + +#if !defined(HAVE_LLROUNDF) && defined(HAVE_ROUNDF) +#define HAVE_LLROUNDF 1 +long long int llroundf (float x); + +long long int +llroundf (float x) +{ + return (long long int) roundf (x); +} +#endif + +#if !defined(HAVE_LLROUND) && defined(HAVE_ROUND) +#define HAVE_LLROUND 1 +long long int llround (double x); + +long long int +llround (double x) +{ + return (long long int) round (x); +} +#endif + +#if !defined(HAVE_LLROUNDL) && defined(HAVE_ROUNDL) +#define HAVE_LLROUNDL 1 +long long int llroundl (long double x); + +long long int +llroundl (long double x) +{ + return (long long int) roundl (x); +} +#endif + + +#ifndef HAVE_LOG10L +#define HAVE_LOG10L 1 +/* log10 function for long double variables. The version provided here + reduces the argument until it fits into a double, then use log10. */ +long double log10l (long double x); + +long double +log10l (long double x) +{ +#if LDBL_MAX_EXP > DBL_MAX_EXP + if (x > DBL_MAX) + { + double val; + int p2_result = 0; + if (x > 0x1p16383L) { p2_result += 16383; x /= 0x1p16383L; } + if (x > 0x1p8191L) { p2_result += 8191; x /= 0x1p8191L; } + if (x > 0x1p4095L) { p2_result += 4095; x /= 0x1p4095L; } + if (x > 0x1p2047L) { p2_result += 2047; x /= 0x1p2047L; } + if (x > 0x1p1023L) { p2_result += 1023; x /= 0x1p1023L; } + val = log10 ((double) x); + return (val + p2_result * .30102999566398119521373889472449302L); + } +#endif +#if LDBL_MIN_EXP < DBL_MIN_EXP + if (x < DBL_MIN) + { + double val; + int p2_result = 0; + if (x < 0x1p-16380L) { p2_result += 16380; x /= 0x1p-16380L; } + if (x < 0x1p-8189L) { p2_result += 8189; x /= 0x1p-8189L; } + if (x < 0x1p-4093L) { p2_result += 4093; x /= 0x1p-4093L; } + if (x < 0x1p-2045L) { p2_result += 2045; x /= 0x1p-2045L; } + if (x < 0x1p-1021L) { p2_result += 1021; x /= 0x1p-1021L; } + val = fabs (log10 ((double) x)); + return (- val - p2_result * .30102999566398119521373889472449302L); + } +#endif + return log10 (x); +} +#endif + + +#ifndef HAVE_FLOORL +#define HAVE_FLOORL 1 +long double floorl (long double x); + +long double +floorl (long double x) +{ + /* Zero, possibly signed. */ + if (x == 0) + return x; + + /* Large magnitude. */ + if (x > DBL_MAX || x < (-DBL_MAX)) + return x; + + /* Small positive values. */ + if (x >= 0 && x < DBL_MIN) + return 0; + + /* Small negative values. */ + if (x < 0 && x > (-DBL_MIN)) + return -1; + + return floor (x); +} +#endif + + +#ifndef HAVE_FMODL +#define HAVE_FMODL 1 +long double fmodl (long double x, long double y); + +long double +fmodl (long double x, long double y) +{ + if (y == 0.0L) + return 0.0L; + + /* Need to check that the result has the same sign as x and magnitude + less than the magnitude of y. */ + return x - floorl (x / y) * y; +} +#endif + + +#if !defined(HAVE_CABSF) +#define HAVE_CABSF 1 +float cabsf (float complex z); + +float +cabsf (float complex z) +{ + return hypotf (REALPART (z), IMAGPART (z)); +} +#endif + +#if !defined(HAVE_CABS) +#define HAVE_CABS 1 +double cabs (double complex z); + +double +cabs (double complex z) +{ + return hypot (REALPART (z), IMAGPART (z)); +} +#endif + +#if !defined(HAVE_CABSL) && defined(HAVE_HYPOTL) +#define HAVE_CABSL 1 +long double cabsl (long double complex z); + +long double +cabsl (long double complex z) +{ + return hypotl (REALPART (z), IMAGPART (z)); +} +#endif + + +#if !defined(HAVE_CARGF) +#define HAVE_CARGF 1 +float cargf (float complex z); + +float +cargf (float complex z) +{ + return atan2f (IMAGPART (z), REALPART (z)); +} +#endif + +#if !defined(HAVE_CARG) +#define HAVE_CARG 1 +double carg (double complex z); + +double +carg (double complex z) +{ + return atan2 (IMAGPART (z), REALPART (z)); +} +#endif + +#if !defined(HAVE_CARGL) && defined(HAVE_ATAN2L) +#define HAVE_CARGL 1 +long double cargl (long double complex z); + +long double +cargl (long double complex z) +{ + return atan2l (IMAGPART (z), REALPART (z)); +} +#endif + + +/* exp(z) = exp(a)*(cos(b) + i sin(b)) */ +#if !defined(HAVE_CEXPF) +#define HAVE_CEXPF 1 +float complex cexpf (float complex z); + +float complex +cexpf (float complex z) +{ + float a, b; + float complex v; + + a = REALPART (z); + b = IMAGPART (z); + COMPLEX_ASSIGN (v, cosf (b), sinf (b)); + return expf (a) * v; +} +#endif + +#if !defined(HAVE_CEXP) +#define HAVE_CEXP 1 +double complex cexp (double complex z); + +double complex +cexp (double complex z) +{ + double a, b; + double complex v; + + a = REALPART (z); + b = IMAGPART (z); + COMPLEX_ASSIGN (v, cos (b), sin (b)); + return exp (a) * v; +} +#endif + +#if !defined(HAVE_CEXPL) && defined(HAVE_COSL) && defined(HAVE_SINL) && defined(EXPL) +#define HAVE_CEXPL 1 +long double complex cexpl (long double complex z); + +long double complex +cexpl (long double complex z) +{ + long double a, b; + long double complex v; + + a = REALPART (z); + b = IMAGPART (z); + COMPLEX_ASSIGN (v, cosl (b), sinl (b)); + return expl (a) * v; +} +#endif + + +/* log(z) = log (cabs(z)) + i*carg(z) */ +#if !defined(HAVE_CLOGF) +#define HAVE_CLOGF 1 +float complex clogf (float complex z); + +float complex +clogf (float complex z) +{ + float complex v; + + COMPLEX_ASSIGN (v, logf (cabsf (z)), cargf (z)); + return v; +} +#endif + +#if !defined(HAVE_CLOG) +#define HAVE_CLOG 1 +double complex clog (double complex z); + +double complex +clog (double complex z) +{ + double complex v; + + COMPLEX_ASSIGN (v, log (cabs (z)), carg (z)); + return v; +} +#endif + +#if !defined(HAVE_CLOGL) && defined(HAVE_LOGL) && defined(HAVE_CABSL) && defined(HAVE_CARGL) +#define HAVE_CLOGL 1 +long double complex clogl (long double complex z); + +long double complex +clogl (long double complex z) +{ + long double complex v; + + COMPLEX_ASSIGN (v, logl (cabsl (z)), cargl (z)); + return v; +} +#endif + + +/* log10(z) = log10 (cabs(z)) + i*carg(z) */ +#if !defined(HAVE_CLOG10F) +#define HAVE_CLOG10F 1 +float complex clog10f (float complex z); + +float complex +clog10f (float complex z) +{ + float complex v; + + COMPLEX_ASSIGN (v, log10f (cabsf (z)), cargf (z)); + return v; +} +#endif + +#if !defined(HAVE_CLOG10) +#define HAVE_CLOG10 1 +double complex clog10 (double complex z); + +double complex +clog10 (double complex z) +{ + double complex v; + + COMPLEX_ASSIGN (v, log10 (cabs (z)), carg (z)); + return v; +} +#endif + +#if !defined(HAVE_CLOG10L) && defined(HAVE_LOG10L) && defined(HAVE_CABSL) && defined(HAVE_CARGL) +#define HAVE_CLOG10L 1 +long double complex clog10l (long double complex z); + +long double complex +clog10l (long double complex z) +{ + long double complex v; + + COMPLEX_ASSIGN (v, log10l (cabsl (z)), cargl (z)); + return v; +} +#endif + + +/* pow(base, power) = cexp (power * clog (base)) */ +#if !defined(HAVE_CPOWF) +#define HAVE_CPOWF 1 +float complex cpowf (float complex base, float complex power); + +float complex +cpowf (float complex base, float complex power) +{ + return cexpf (power * clogf (base)); +} +#endif + +#if !defined(HAVE_CPOW) +#define HAVE_CPOW 1 +double complex cpow (double complex base, double complex power); + +double complex +cpow (double complex base, double complex power) +{ + return cexp (power * clog (base)); +} +#endif + +#if !defined(HAVE_CPOWL) && defined(HAVE_CEXPL) && defined(HAVE_CLOGL) +#define HAVE_CPOWL 1 +long double complex cpowl (long double complex base, long double complex power); + +long double complex +cpowl (long double complex base, long double complex power) +{ + return cexpl (power * clogl (base)); +} +#endif + + +/* sqrt(z). Algorithm pulled from glibc. */ +#if !defined(HAVE_CSQRTF) +#define HAVE_CSQRTF 1 +float complex csqrtf (float complex z); + +float complex +csqrtf (float complex z) +{ + float re, im; + float complex v; + + re = REALPART (z); + im = IMAGPART (z); + if (im == 0) + { + if (re < 0) + { + COMPLEX_ASSIGN (v, 0, copysignf (sqrtf (-re), im)); + } + else + { + COMPLEX_ASSIGN (v, fabsf (sqrtf (re)), copysignf (0, im)); + } + } + else if (re == 0) + { + float r; + + r = sqrtf (0.5 * fabsf (im)); + + COMPLEX_ASSIGN (v, r, copysignf (r, im)); + } + else + { + float d, r, s; + + d = hypotf (re, im); + /* Use the identity 2 Re res Im res = Im x + to avoid cancellation error in d +/- Re x. */ + if (re > 0) + { + r = sqrtf (0.5 * d + 0.5 * re); + s = (0.5 * im) / r; + } + else + { + s = sqrtf (0.5 * d - 0.5 * re); + r = fabsf ((0.5 * im) / s); + } + + COMPLEX_ASSIGN (v, r, copysignf (s, im)); + } + return v; +} +#endif + +#if !defined(HAVE_CSQRT) +#define HAVE_CSQRT 1 +double complex csqrt (double complex z); + +double complex +csqrt (double complex z) +{ + double re, im; + double complex v; + + re = REALPART (z); + im = IMAGPART (z); + if (im == 0) + { + if (re < 0) + { + COMPLEX_ASSIGN (v, 0, copysign (sqrt (-re), im)); + } + else + { + COMPLEX_ASSIGN (v, fabs (sqrt (re)), copysign (0, im)); + } + } + else if (re == 0) + { + double r; + + r = sqrt (0.5 * fabs (im)); + + COMPLEX_ASSIGN (v, r, copysign (r, im)); + } + else + { + double d, r, s; + + d = hypot (re, im); + /* Use the identity 2 Re res Im res = Im x + to avoid cancellation error in d +/- Re x. */ + if (re > 0) + { + r = sqrt (0.5 * d + 0.5 * re); + s = (0.5 * im) / r; + } + else + { + s = sqrt (0.5 * d - 0.5 * re); + r = fabs ((0.5 * im) / s); + } + + COMPLEX_ASSIGN (v, r, copysign (s, im)); + } + return v; +} +#endif + +#if !defined(HAVE_CSQRTL) && defined(HAVE_COPYSIGNL) && defined(HAVE_SQRTL) && defined(HAVE_FABSL) && defined(HAVE_HYPOTL) +#define HAVE_CSQRTL 1 +long double complex csqrtl (long double complex z); + +long double complex +csqrtl (long double complex z) +{ + long double re, im; + long double complex v; + + re = REALPART (z); + im = IMAGPART (z); + if (im == 0) + { + if (re < 0) + { + COMPLEX_ASSIGN (v, 0, copysignl (sqrtl (-re), im)); + } + else + { + COMPLEX_ASSIGN (v, fabsl (sqrtl (re)), copysignl (0, im)); + } + } + else if (re == 0) + { + long double r; + + r = sqrtl (0.5 * fabsl (im)); + + COMPLEX_ASSIGN (v, copysignl (r, im), r); + } + else + { + long double d, r, s; + + d = hypotl (re, im); + /* Use the identity 2 Re res Im res = Im x + to avoid cancellation error in d +/- Re x. */ + if (re > 0) + { + r = sqrtl (0.5 * d + 0.5 * re); + s = (0.5 * im) / r; + } + else + { + s = sqrtl (0.5 * d - 0.5 * re); + r = fabsl ((0.5 * im) / s); + } + + COMPLEX_ASSIGN (v, r, copysignl (s, im)); + } + return v; +} +#endif + + +/* sinh(a + i b) = sinh(a) cos(b) + i cosh(a) sin(b) */ +#if !defined(HAVE_CSINHF) +#define HAVE_CSINHF 1 +float complex csinhf (float complex a); + +float complex +csinhf (float complex a) +{ + float r, i; + float complex v; + + r = REALPART (a); + i = IMAGPART (a); + COMPLEX_ASSIGN (v, sinhf (r) * cosf (i), coshf (r) * sinf (i)); + return v; +} +#endif + +#if !defined(HAVE_CSINH) +#define HAVE_CSINH 1 +double complex csinh (double complex a); + +double complex +csinh (double complex a) +{ + double r, i; + double complex v; + + r = REALPART (a); + i = IMAGPART (a); + COMPLEX_ASSIGN (v, sinh (r) * cos (i), cosh (r) * sin (i)); + return v; +} +#endif + +#if !defined(HAVE_CSINHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL) +#define HAVE_CSINHL 1 +long double complex csinhl (long double complex a); + +long double complex +csinhl (long double complex a) +{ + long double r, i; + long double complex v; + + r = REALPART (a); + i = IMAGPART (a); + COMPLEX_ASSIGN (v, sinhl (r) * cosl (i), coshl (r) * sinl (i)); + return v; +} +#endif + + +/* cosh(a + i b) = cosh(a) cos(b) + i sinh(a) sin(b) */ +#if !defined(HAVE_CCOSHF) +#define HAVE_CCOSHF 1 +float complex ccoshf (float complex a); + +float complex +ccoshf (float complex a) +{ + float r, i; + float complex v; + + r = REALPART (a); + i = IMAGPART (a); + COMPLEX_ASSIGN (v, coshf (r) * cosf (i), sinhf (r) * sinf (i)); + return v; +} +#endif + +#if !defined(HAVE_CCOSH) +#define HAVE_CCOSH 1 +double complex ccosh (double complex a); + +double complex +ccosh (double complex a) +{ + double r, i; + double complex v; + + r = REALPART (a); + i = IMAGPART (a); + COMPLEX_ASSIGN (v, cosh (r) * cos (i), sinh (r) * sin (i)); + return v; +} +#endif + +#if !defined(HAVE_CCOSHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL) +#define HAVE_CCOSHL 1 +long double complex ccoshl (long double complex a); + +long double complex +ccoshl (long double complex a) +{ + long double r, i; + long double complex v; + + r = REALPART (a); + i = IMAGPART (a); + COMPLEX_ASSIGN (v, coshl (r) * cosl (i), sinhl (r) * sinl (i)); + return v; +} +#endif + + +/* tanh(a + i b) = (tanh(a) + i tan(b)) / (1 + i tanh(a) tan(b)) */ +#if !defined(HAVE_CTANHF) +#define HAVE_CTANHF 1 +float complex ctanhf (float complex a); + +float complex +ctanhf (float complex a) +{ + float rt, it; + float complex n, d; + + rt = tanhf (REALPART (a)); + it = tanf (IMAGPART (a)); + COMPLEX_ASSIGN (n, rt, it); + COMPLEX_ASSIGN (d, 1, rt * it); + + return n / d; +} +#endif + +#if !defined(HAVE_CTANH) +#define HAVE_CTANH 1 +double complex ctanh (double complex a); +double complex +ctanh (double complex a) +{ + double rt, it; + double complex n, d; + + rt = tanh (REALPART (a)); + it = tan (IMAGPART (a)); + COMPLEX_ASSIGN (n, rt, it); + COMPLEX_ASSIGN (d, 1, rt * it); + + return n / d; +} +#endif + +#if !defined(HAVE_CTANHL) && defined(HAVE_TANL) && defined(HAVE_TANHL) +#define HAVE_CTANHL 1 +long double complex ctanhl (long double complex a); + +long double complex +ctanhl (long double complex a) +{ + long double rt, it; + long double complex n, d; + + rt = tanhl (REALPART (a)); + it = tanl (IMAGPART (a)); + COMPLEX_ASSIGN (n, rt, it); + COMPLEX_ASSIGN (d, 1, rt * it); + + return n / d; +} +#endif + + +/* sin(a + i b) = sin(a) cosh(b) + i cos(a) sinh(b) */ +#if !defined(HAVE_CSINF) +#define HAVE_CSINF 1 +float complex csinf (float complex a); + +float complex +csinf (float complex a) +{ + float r, i; + float complex v; + + r = REALPART (a); + i = IMAGPART (a); + COMPLEX_ASSIGN (v, sinf (r) * coshf (i), cosf (r) * sinhf (i)); + return v; +} +#endif + +#if !defined(HAVE_CSIN) +#define HAVE_CSIN 1 +double complex csin (double complex a); + +double complex +csin (double complex a) +{ + double r, i; + double complex v; + + r = REALPART (a); + i = IMAGPART (a); + COMPLEX_ASSIGN (v, sin (r) * cosh (i), cos (r) * sinh (i)); + return v; +} +#endif + +#if !defined(HAVE_CSINL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL) +#define HAVE_CSINL 1 +long double complex csinl (long double complex a); + +long double complex +csinl (long double complex a) +{ + long double r, i; + long double complex v; + + r = REALPART (a); + i = IMAGPART (a); + COMPLEX_ASSIGN (v, sinl (r) * coshl (i), cosl (r) * sinhl (i)); + return v; +} +#endif + + +/* cos(a + i b) = cos(a) cosh(b) - i sin(a) sinh(b) */ +#if !defined(HAVE_CCOSF) +#define HAVE_CCOSF 1 +float complex ccosf (float complex a); + +float complex +ccosf (float complex a) +{ + float r, i; + float complex v; + + r = REALPART (a); + i = IMAGPART (a); + COMPLEX_ASSIGN (v, cosf (r) * coshf (i), - (sinf (r) * sinhf (i))); + return v; +} +#endif + +#if !defined(HAVE_CCOS) +#define HAVE_CCOS 1 +double complex ccos (double complex a); + +double complex +ccos (double complex a) +{ + double r, i; + double complex v; + + r = REALPART (a); + i = IMAGPART (a); + COMPLEX_ASSIGN (v, cos (r) * cosh (i), - (sin (r) * sinh (i))); + return v; +} +#endif + +#if !defined(HAVE_CCOSL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL) +#define HAVE_CCOSL 1 +long double complex ccosl (long double complex a); + +long double complex +ccosl (long double complex a) +{ + long double r, i; + long double complex v; + + r = REALPART (a); + i = IMAGPART (a); + COMPLEX_ASSIGN (v, cosl (r) * coshl (i), - (sinl (r) * sinhl (i))); + return v; +} +#endif + + +/* tan(a + i b) = (tan(a) + i tanh(b)) / (1 - i tan(a) tanh(b)) */ +#if !defined(HAVE_CTANF) +#define HAVE_CTANF 1 +float complex ctanf (float complex a); + +float complex +ctanf (float complex a) +{ + float rt, it; + float complex n, d; + + rt = tanf (REALPART (a)); + it = tanhf (IMAGPART (a)); + COMPLEX_ASSIGN (n, rt, it); + COMPLEX_ASSIGN (d, 1, - (rt * it)); + + return n / d; +} +#endif + +#if !defined(HAVE_CTAN) +#define HAVE_CTAN 1 +double complex ctan (double complex a); + +double complex +ctan (double complex a) +{ + double rt, it; + double complex n, d; + + rt = tan (REALPART (a)); + it = tanh (IMAGPART (a)); + COMPLEX_ASSIGN (n, rt, it); + COMPLEX_ASSIGN (d, 1, - (rt * it)); + + return n / d; +} +#endif + +#if !defined(HAVE_CTANL) && defined(HAVE_TANL) && defined(HAVE_TANHL) +#define HAVE_CTANL 1 +long double complex ctanl (long double complex a); + +long double complex +ctanl (long double complex a) +{ + long double rt, it; + long double complex n, d; + + rt = tanl (REALPART (a)); + it = tanhl (IMAGPART (a)); + COMPLEX_ASSIGN (n, rt, it); + COMPLEX_ASSIGN (d, 1, - (rt * it)); + + return n / d; +} +#endif + + +/* Complex ASIN. Returns wrongly NaN for infinite arguments. + Algorithm taken from Abramowitz & Stegun. */ + +#if !defined(HAVE_CASINF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF) +#define HAVE_CASINF 1 +complex float casinf (complex float z); + +complex float +casinf (complex float z) +{ + return -I*clogf (I*z + csqrtf (1.0f-z*z)); +} +#endif + + +#if !defined(HAVE_CASIN) && defined(HAVE_CLOG) && defined(HAVE_CSQRT) +#define HAVE_CASIN 1 +complex double casin (complex double z); + +complex double +casin (complex double z) +{ + return -I*clog (I*z + csqrt (1.0-z*z)); +} +#endif + + +#if !defined(HAVE_CASINL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL) +#define HAVE_CASINL 1 +complex long double casinl (complex long double z); + +complex long double +casinl (complex long double z) +{ + return -I*clogl (I*z + csqrtl (1.0L-z*z)); +} +#endif + + +/* Complex ACOS. Returns wrongly NaN for infinite arguments. + Algorithm taken from Abramowitz & Stegun. */ + +#if !defined(HAVE_CACOSF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF) +#define HAVE_CACOSF 1 +complex float cacosf (complex float z); + +complex float +cacosf (complex float z) +{ + return -I*clogf (z + I*csqrtf (1.0f-z*z)); +} +#endif + + +#if !defined(HAVE_CACOS) && defined(HAVE_CLOG) && defined(HAVE_CSQRT) +#define HAVE_CACOS 1 +complex double cacos (complex double z); + +complex double +cacos (complex double z) +{ + return -I*clog (z + I*csqrt (1.0-z*z)); +} +#endif + + +#if !defined(HAVE_CACOSL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL) +#define HAVE_CACOSL 1 +complex long double cacosl (complex long double z); + +complex long double +cacosl (complex long double z) +{ + return -I*clogl (z + I*csqrtl (1.0L-z*z)); +} +#endif + + +/* Complex ATAN. Returns wrongly NaN for infinite arguments. + Algorithm taken from Abramowitz & Stegun. */ + +#if !defined(HAVE_CATANF) && defined(HAVE_CLOGF) +#define HAVE_CACOSF 1 +complex float catanf (complex float z); + +complex float +catanf (complex float z) +{ + return I*clogf ((I+z)/(I-z))/2.0f; +} +#endif + + +#if !defined(HAVE_CATAN) && defined(HAVE_CLOG) +#define HAVE_CACOS 1 +complex double catan (complex double z); + +complex double +catan (complex double z) +{ + return I*clog ((I+z)/(I-z))/2.0; +} +#endif + + +#if !defined(HAVE_CATANL) && defined(HAVE_CLOGL) +#define HAVE_CACOSL 1 +complex long double catanl (complex long double z); + +complex long double +catanl (complex long double z) +{ + return I*clogl ((I+z)/(I-z))/2.0L; +} +#endif + + +/* Complex ASINH. Returns wrongly NaN for infinite arguments. + Algorithm taken from Abramowitz & Stegun. */ + +#if !defined(HAVE_CASINHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF) +#define HAVE_CASINHF 1 +complex float casinhf (complex float z); + +complex float +casinhf (complex float z) +{ + return clogf (z + csqrtf (z*z+1.0f)); +} +#endif + + +#if !defined(HAVE_CASINH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT) +#define HAVE_CASINH 1 +complex double casinh (complex double z); + +complex double +casinh (complex double z) +{ + return clog (z + csqrt (z*z+1.0)); +} +#endif + + +#if !defined(HAVE_CASINHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL) +#define HAVE_CASINHL 1 +complex long double casinhl (complex long double z); + +complex long double +casinhl (complex long double z) +{ + return clogl (z + csqrtl (z*z+1.0L)); +} +#endif + + +/* Complex ACOSH. Returns wrongly NaN for infinite arguments. + Algorithm taken from Abramowitz & Stegun. */ + +#if !defined(HAVE_CACOSHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF) +#define HAVE_CACOSHF 1 +complex float cacoshf (complex float z); + +complex float +cacoshf (complex float z) +{ + return clogf (z + csqrtf (z-1.0f) * csqrtf (z+1.0f)); +} +#endif + + +#if !defined(HAVE_CACOSH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT) +#define HAVE_CACOSH 1 +complex double cacosh (complex double z); + +complex double +cacosh (complex double z) +{ + return clog (z + csqrt (z-1.0) * csqrt (z+1.0)); +} +#endif + + +#if !defined(HAVE_CACOSHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL) +#define HAVE_CACOSHL 1 +complex long double cacoshl (complex long double z); + +complex long double +cacoshl (complex long double z) +{ + return clogl (z + csqrtl (z-1.0L) * csqrtl (z+1.0L)); +} +#endif + + +/* Complex ATANH. Returns wrongly NaN for infinite arguments. + Algorithm taken from Abramowitz & Stegun. */ + +#if !defined(HAVE_CATANHF) && defined(HAVE_CLOGF) +#define HAVE_CATANHF 1 +complex float catanhf (complex float z); + +complex float +catanhf (complex float z) +{ + return clogf ((1.0f+z)/(1.0f-z))/2.0f; +} +#endif + + +#if !defined(HAVE_CATANH) && defined(HAVE_CLOG) +#define HAVE_CATANH 1 +complex double catanh (complex double z); + +complex double +catanh (complex double z) +{ + return clog ((1.0+z)/(1.0-z))/2.0; +} +#endif + +#if !defined(HAVE_CATANHL) && defined(HAVE_CLOGL) +#define HAVE_CATANHL 1 +complex long double catanhl (complex long double z); + +complex long double +catanhl (complex long double z) +{ + return clogl ((1.0L+z)/(1.0L-z))/2.0L; +} +#endif + + +#if !defined(HAVE_TGAMMA) +#define HAVE_TGAMMA 1 +double tgamma (double); + +/* Fallback tgamma() function. Uses the algorithm from + http://www.netlib.org/specfun/gamma and references therein. */ + +#undef SQRTPI +#define SQRTPI 0.9189385332046727417803297 + +#undef PI +#define PI 3.1415926535897932384626434 + +double +tgamma (double x) +{ + int i, n, parity; + double fact, res, sum, xden, xnum, y, y1, ysq, z; + + static double p[8] = { + -1.71618513886549492533811e0, 2.47656508055759199108314e1, + -3.79804256470945635097577e2, 6.29331155312818442661052e2, + 8.66966202790413211295064e2, -3.14512729688483675254357e4, + -3.61444134186911729807069e4, 6.64561438202405440627855e4 }; + + static double q[8] = { + -3.08402300119738975254353e1, 3.15350626979604161529144e2, + -1.01515636749021914166146e3, -3.10777167157231109440444e3, + 2.25381184209801510330112e4, 4.75584627752788110767815e3, + -1.34659959864969306392456e5, -1.15132259675553483497211e5 }; + + static double c[7] = { -1.910444077728e-03, + 8.4171387781295e-04, -5.952379913043012e-04, + 7.93650793500350248e-04, -2.777777777777681622553e-03, + 8.333333333333333331554247e-02, 5.7083835261e-03 }; + + static const double xminin = 2.23e-308; + static const double xbig = 171.624; + static const double xnan = __builtin_nan ("0x0"), xinf = __builtin_inf (); + static double eps = 0; + + if (eps == 0) + eps = nextafter (1., 2.) - 1.; + + parity = 0; + fact = 1; + n = 0; + y = x; + + if (isnan (x)) + return x; + + if (y <= 0) + { + y = -x; + y1 = trunc (y); + res = y - y1; + + if (res != 0) + { + if (y1 != trunc (y1*0.5l)*2) + parity = 1; + fact = -PI / sin (PI*res); + y = y + 1; + } + else + return x == 0 ? copysign (xinf, x) : xnan; + } + + if (y < eps) + { + if (y >= xminin) + res = 1 / y; + else + return xinf; + } + else if (y < 13) + { + y1 = y; + if (y < 1) + { + z = y; + y = y + 1; + } + else + { + n = (int)y - 1; + y = y - n; + z = y - 1; + } + + xnum = 0; + xden = 1; + for (i = 0; i < 8; i++) + { + xnum = (xnum + p[i]) * z; + xden = xden * z + q[i]; + } + + res = xnum / xden + 1; + + if (y1 < y) + res = res / y1; + else if (y1 > y) + for (i = 1; i <= n; i++) + { + res = res * y; + y = y + 1; + } + } + else + { + if (y < xbig) + { + ysq = y * y; + sum = c[6]; + for (i = 0; i < 6; i++) + sum = sum / ysq + c[i]; + + sum = sum/y - y + SQRTPI; + sum = sum + (y - 0.5) * log (y); + res = exp (sum); + } + else + return x < 0 ? xnan : xinf; + } + + if (parity) + res = -res; + if (fact != 1) + res = fact / res; + + return res; +} +#endif + + + +#if !defined(HAVE_LGAMMA) +#define HAVE_LGAMMA 1 +double lgamma (double); + +/* Fallback lgamma() function. Uses the algorithm from + http://www.netlib.org/specfun/algama and references therein, + except for negative arguments (where netlib would return +Inf) + where we use the following identity: + lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y) + */ + +double +lgamma (double y) +{ + +#undef SQRTPI +#define SQRTPI 0.9189385332046727417803297 + +#undef PI +#define PI 3.1415926535897932384626434 + +#define PNT68 0.6796875 +#define D1 -0.5772156649015328605195174 +#define D2 0.4227843350984671393993777 +#define D4 1.791759469228055000094023 + + static double p1[8] = { + 4.945235359296727046734888e0, 2.018112620856775083915565e2, + 2.290838373831346393026739e3, 1.131967205903380828685045e4, + 2.855724635671635335736389e4, 3.848496228443793359990269e4, + 2.637748787624195437963534e4, 7.225813979700288197698961e3 }; + static double q1[8] = { + 6.748212550303777196073036e1, 1.113332393857199323513008e3, + 7.738757056935398733233834e3, 2.763987074403340708898585e4, + 5.499310206226157329794414e4, 6.161122180066002127833352e4, + 3.635127591501940507276287e4, 8.785536302431013170870835e3 }; + static double p2[8] = { + 4.974607845568932035012064e0, 5.424138599891070494101986e2, + 1.550693864978364947665077e4, 1.847932904445632425417223e5, + 1.088204769468828767498470e6, 3.338152967987029735917223e6, + 5.106661678927352456275255e6, 3.074109054850539556250927e6 }; + static double q2[8] = { + 1.830328399370592604055942e2, 7.765049321445005871323047e3, + 1.331903827966074194402448e5, 1.136705821321969608938755e6, + 5.267964117437946917577538e6, 1.346701454311101692290052e7, + 1.782736530353274213975932e7, 9.533095591844353613395747e6 }; + static double p4[8] = { + 1.474502166059939948905062e4, 2.426813369486704502836312e6, + 1.214755574045093227939592e8, 2.663432449630976949898078e9, + 2.940378956634553899906876e10, 1.702665737765398868392998e11, + 4.926125793377430887588120e11, 5.606251856223951465078242e11 }; + static double q4[8] = { + 2.690530175870899333379843e3, 6.393885654300092398984238e5, + 4.135599930241388052042842e7, 1.120872109616147941376570e9, + 1.488613728678813811542398e10, 1.016803586272438228077304e11, + 3.417476345507377132798597e11, 4.463158187419713286462081e11 }; + static double c[7] = { + -1.910444077728e-03, 8.4171387781295e-04, + -5.952379913043012e-04, 7.93650793500350248e-04, + -2.777777777777681622553e-03, 8.333333333333333331554247e-02, + 5.7083835261e-03 }; + + static double xbig = 2.55e305, xinf = __builtin_inf (), eps = 0, + frtbig = 2.25e76; + + int i; + double corr, res, xden, xm1, xm2, xm4, xnum, ysq; + + if (eps == 0) + eps = __builtin_nextafter (1., 2.) - 1.; + + if ((y > 0) && (y <= xbig)) + { + if (y <= eps) + res = -log (y); + else if (y <= 1.5) + { + if (y < PNT68) + { + corr = -log (y); + xm1 = y; + } + else + { + corr = 0; + xm1 = (y - 0.5) - 0.5; + } + + if ((y <= 0.5) || (y >= PNT68)) + { + xden = 1; + xnum = 0; + for (i = 0; i < 8; i++) + { + xnum = xnum*xm1 + p1[i]; + xden = xden*xm1 + q1[i]; + } + res = corr + (xm1 * (D1 + xm1*(xnum/xden))); + } + else + { + xm2 = (y - 0.5) - 0.5; + xden = 1; + xnum = 0; + for (i = 0; i < 8; i++) + { + xnum = xnum*xm2 + p2[i]; + xden = xden*xm2 + q2[i]; + } + res = corr + xm2 * (D2 + xm2*(xnum/xden)); + } + } + else if (y <= 4) + { + xm2 = y - 2; + xden = 1; + xnum = 0; + for (i = 0; i < 8; i++) + { + xnum = xnum*xm2 + p2[i]; + xden = xden*xm2 + q2[i]; + } + res = xm2 * (D2 + xm2*(xnum/xden)); + } + else if (y <= 12) + { + xm4 = y - 4; + xden = -1; + xnum = 0; + for (i = 0; i < 8; i++) + { + xnum = xnum*xm4 + p4[i]; + xden = xden*xm4 + q4[i]; + } + res = D4 + xm4*(xnum/xden); + } + else + { + res = 0; + if (y <= frtbig) + { + res = c[6]; + ysq = y * y; + for (i = 0; i < 6; i++) + res = res / ysq + c[i]; + } + res = res/y; + corr = log (y); + res = res + SQRTPI - 0.5*corr; + res = res + y*(corr-1); + } + } + else if (y < 0 && __builtin_floor (y) != y) + { + /* lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y) + For abs(y) very close to zero, we use a series expansion to + the first order in y to avoid overflow. */ + if (y > -1.e-100) + res = -2 * log (fabs (y)) - lgamma (-y); + else + res = log (PI / fabs (y * sin (PI * y))) - lgamma (-y); + } + else + res = xinf; + + return res; +} +#endif + + +#if defined(HAVE_TGAMMA) && !defined(HAVE_TGAMMAF) +#define HAVE_TGAMMAF 1 +float tgammaf (float); + +float +tgammaf (float x) +{ + return (float) tgamma ((double) x); +} +#endif + +#if defined(HAVE_LGAMMA) && !defined(HAVE_LGAMMAF) +#define HAVE_LGAMMAF 1 +float lgammaf (float); + +float +lgammaf (float x) +{ + return (float) lgamma ((double) x); +} +#endif |