| /* Implementation of various C99 functions |
| Copyright (C) 2004, 2009 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 |
| |
| /* Prototypes to silence -Wstrict-prototypes -Wmissing-prototypes. */ |
| |
| float cabsf(float complex); |
| double cabs(double complex); |
| long double cabsl(long double complex); |
| |
| float cargf(float complex); |
| double carg(double complex); |
| long double cargl(long double complex); |
| |
| float complex clog10f(float complex); |
| double complex clog10(double complex); |
| long double complex clog10l(long double complex); |
| |
| |
| /* Wrappers for systems without the various C99 single precision Bessel |
| functions. */ |
| |
| #if defined(HAVE_J0) && ! defined(HAVE_J0F) |
| #define HAVE_J0F 1 |
| extern float j0f (float); |
| |
| float |
| j0f (float x) |
| { |
| return (float) j0 ((double) x); |
| } |
| #endif |
| |
| #if defined(HAVE_J1) && !defined(HAVE_J1F) |
| #define HAVE_J1F 1 |
| extern float j1f (float); |
| |
| float j1f (float x) |
| { |
| return (float) j1 ((double) x); |
| } |
| #endif |
| |
| #if defined(HAVE_JN) && !defined(HAVE_JNF) |
| #define HAVE_JNF 1 |
| extern 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 |
| extern float y0f (float); |
| |
| float |
| y0f (float x) |
| { |
| return (float) y0 ((double) x); |
| } |
| #endif |
| |
| #if defined(HAVE_Y1) && !defined(HAVE_Y1F) |
| #define HAVE_Y1F 1 |
| extern float y1f (float); |
| |
| float |
| y1f (float x) |
| { |
| return (float) y1 ((double) x); |
| } |
| #endif |
| |
| #if defined(HAVE_YN) && !defined(HAVE_YNF) |
| #define HAVE_YNF 1 |
| extern 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 |
| extern float erff (float); |
| |
| float |
| erff (float x) |
| { |
| return (float) erf ((double) x); |
| } |
| #endif |
| |
| #if defined(HAVE_ERFC) && !defined(HAVE_ERFCF) |
| #define HAVE_ERFCF 1 |
| extern float erfcf (float); |
| |
| float |
| erfcf (float x) |
| { |
| return (float) erfc ((double) x); |
| } |
| #endif |
| |
| |
| #ifndef HAVE_ACOSF |
| #define HAVE_ACOSF 1 |
| float |
| acosf(float x) |
| { |
| return (float) acos(x); |
| } |
| #endif |
| |
| #if HAVE_ACOSH && !HAVE_ACOSHF |
| float |
| acoshf (float x) |
| { |
| return (float) acosh ((double) x); |
| } |
| #endif |
| |
| #ifndef HAVE_ASINF |
| #define HAVE_ASINF 1 |
| float |
| asinf(float x) |
| { |
| return (float) asin(x); |
| } |
| #endif |
| |
| #if HAVE_ASINH && !HAVE_ASINHF |
| float |
| asinhf (float x) |
| { |
| return (float) asinh ((double) x); |
| } |
| #endif |
| |
| #ifndef HAVE_ATAN2F |
| #define HAVE_ATAN2F 1 |
| float |
| atan2f(float y, float x) |
| { |
| return (float) atan2(y, x); |
| } |
| #endif |
| |
| #ifndef HAVE_ATANF |
| #define HAVE_ATANF 1 |
| float |
| atanf(float x) |
| { |
| return (float) atan(x); |
| } |
| #endif |
| |
| #if HAVE_ATANH && !HAVE_ATANHF |
| float |
| atanhf (float x) |
| { |
| return (float) atanh ((double) x); |
| } |
| #endif |
| |
| #ifndef HAVE_CEILF |
| #define HAVE_CEILF 1 |
| float |
| ceilf(float x) |
| { |
| return (float) ceil(x); |
| } |
| #endif |
| |
| #ifndef HAVE_COPYSIGNF |
| #define HAVE_COPYSIGNF 1 |
| float |
| copysignf(float x, float y) |
| { |
| return (float) copysign(x, y); |
| } |
| #endif |
| |
| #ifndef HAVE_COSF |
| #define HAVE_COSF 1 |
| float |
| cosf(float x) |
| { |
| return (float) cos(x); |
| } |
| #endif |
| |
| #ifndef HAVE_COSHF |
| #define HAVE_COSHF 1 |
| float |
| coshf(float x) |
| { |
| return (float) cosh(x); |
| } |
| #endif |
| |
| #ifndef HAVE_EXPF |
| #define HAVE_EXPF 1 |
| float |
| expf(float x) |
| { |
| return (float) exp(x); |
| } |
| #endif |
| |
| #ifndef HAVE_FABSF |
| #define HAVE_FABSF 1 |
| float |
| fabsf(float x) |
| { |
| return (float) fabs(x); |
| } |
| #endif |
| |
| #ifndef HAVE_FLOORF |
| #define HAVE_FLOORF 1 |
| float |
| floorf(float x) |
| { |
| return (float) floor(x); |
| } |
| #endif |
| |
| #ifndef HAVE_FMODF |
| #define HAVE_FMODF 1 |
| 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) |
| { |
| return (float) frexp(x, exp); |
| } |
| #endif |
| |
| #ifndef HAVE_HYPOTF |
| #define HAVE_HYPOTF 1 |
| float |
| hypotf(float x, float y) |
| { |
| return (float) hypot(x, y); |
| } |
| #endif |
| |
| #ifndef HAVE_LOGF |
| #define HAVE_LOGF 1 |
| float |
| logf(float x) |
| { |
| return (float) log(x); |
| } |
| #endif |
| |
| #ifndef HAVE_LOG10F |
| #define HAVE_LOG10F 1 |
| float |
| log10f(float x) |
| { |
| return (float) log10(x); |
| } |
| #endif |
| |
| #ifndef HAVE_SCALBN |
| #define HAVE_SCALBN 1 |
| 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) |
| { |
| return (float) scalbn(x, y); |
| } |
| #endif |
| |
| #ifndef HAVE_SINF |
| #define HAVE_SINF 1 |
| float |
| sinf(float x) |
| { |
| return (float) sin(x); |
| } |
| #endif |
| |
| #ifndef HAVE_SINHF |
| #define HAVE_SINHF 1 |
| float |
| sinhf(float x) |
| { |
| return (float) sinh(x); |
| } |
| #endif |
| |
| #ifndef HAVE_SQRTF |
| #define HAVE_SQRTF 1 |
| float |
| sqrtf(float x) |
| { |
| return (float) sqrt(x); |
| } |
| #endif |
| |
| #ifndef HAVE_TANF |
| #define HAVE_TANF 1 |
| float |
| tanf(float x) |
| { |
| return (float) tan(x); |
| } |
| #endif |
| |
| #ifndef HAVE_TANHF |
| #define HAVE_TANHF 1 |
| float |
| tanhf(float x) |
| { |
| return (float) tanh(x); |
| } |
| #endif |
| |
| #ifndef HAVE_TRUNC |
| #define HAVE_TRUNC 1 |
| 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) |
| { |
| 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) |
| { |
| /* 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) |
| { |
| return (float) pow(x, y); |
| } |
| #endif |
| |
| /* Note that if fpclassify is not defined, then NaN is not handled */ |
| |
| /* Algorithm by Steven G. Kargl. */ |
| |
| #if !defined(HAVE_ROUNDL) |
| #define HAVE_ROUNDL 1 |
| #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 |
| static long double prechalf = nexafterl (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 t; |
| if (!isfinite (x)) |
| return (x); |
| |
| if (x >= 0.0) |
| { |
| t = ceil(x); |
| if (t - x > 0.5) |
| t -= 1.0; |
| return (t); |
| } |
| else |
| { |
| t = ceil(-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 t; |
| if (!isfinite (x)) |
| return (x); |
| |
| if (x >= 0.0) |
| { |
| t = ceilf(x); |
| if (t - x > 0.5) |
| t -= 1.0; |
| return (t); |
| } |
| else |
| { |
| t = ceilf(-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) |
| { |
| return (long int) roundf (x); |
| } |
| #endif |
| |
| #if !defined(HAVE_LROUND) && defined(HAVE_ROUND) |
| #define HAVE_LROUND 1 |
| 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) |
| { |
| return (long long int) roundl (x); |
| } |
| #endif |
| |
| #if !defined(HAVE_LLROUNDF) && defined(HAVE_ROUNDF) |
| #define HAVE_LLROUNDF 1 |
| 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) |
| { |
| 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) |
| { |
| 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) |
| { |
| #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) |
| { |
| /* 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) |
| { |
| 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) |
| { |
| return hypotf (REALPART (z), IMAGPART (z)); |
| } |
| #endif |
| |
| #if !defined(HAVE_CABS) |
| #define HAVE_CABS 1 |
| 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) |
| { |
| return hypotl (REALPART (z), IMAGPART (z)); |
| } |
| #endif |
| |
| |
| #if !defined(HAVE_CARGF) |
| #define HAVE_CARGF 1 |
| 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) |
| { |
| 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) |
| { |
| 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 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 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 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 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 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 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 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 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 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) |
| { |
| return cexpf (power * clogf (base)); |
| } |
| #endif |
| |
| #if !defined(HAVE_CPOW) |
| #define HAVE_CPOW 1 |
| 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) |
| { |
| 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 |
| |
| |
| #if !defined(HAVE_TGAMMA) |
| #define HAVE_TGAMMA 1 |
| |
| extern 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 (__builtin_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 |
| |
| extern 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 |
| extern float tgammaf (float); |
| |
| float |
| tgammaf (float x) |
| { |
| return (float) tgamma ((double) x); |
| } |
| #endif |
| |
| #if defined(HAVE_LGAMMA) && !defined(HAVE_LGAMMAF) |
| #define HAVE_LGAMMAF 1 |
| extern float lgammaf (float); |
| |
| float |
| lgammaf (float x) |
| { |
| return (float) lgamma ((double) x); |
| } |
| #endif |
