URL
https://opencores.org/ocsvn/openrisc/openrisc/trunk
Subversion Repositories openrisc
[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libgfortran/] [intrinsics/] [c99_functions.c] - Rev 742
Go to most recent revision | Compare with Previous | Blame | View Log
/* 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 #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 /* 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_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
Go to most recent revision | Compare with Previous | Blame | View Log