/* Implementation of various C99 functions
- Copyright (C) 2004 Free Software Foundation, Inc.
+ 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 Lesser General Public
+modify it under the terms of the GNU General Public
License as published by the Free Software Foundation; either
-version 2.1 of the License, or (at your option) any later version.
+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 Lesser General Public License for more details.
+GNU General Public License for more details.
-You should have received a copy of the GNU Lesser General Public
-License along with libgfortran; see the file COPYING.LIB. If not,
-write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-Boston, MA 02111-1307, USA. */
+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"
-#include <sys/types.h>
-#include <float.h>
-#include <math.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
-acosf(float x)
+acoshf (float x)
{
- return (float) acos(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
-asinf(float x)
+asinhf (float x)
{
- return (float) asin(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)
+atan2f (float y, float x)
{
- return (float) atan2(y, 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
-atanf(float x)
+atanhf (float x)
{
- return (float) atan(x);
+ return (float) atanh ((double) x);
}
#endif
#ifndef HAVE_CEILF
+#define HAVE_CEILF 1
+float ceilf (float x);
+
float
-ceilf(float x)
+ceilf (float x)
{
- return (float) ceil(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)
+copysignf (float x, float y)
{
- return (float) copysign(x, y);
+ return (float) copysign (x, y);
}
#endif
#ifndef HAVE_COSF
+#define HAVE_COSF 1
+float cosf (float x);
+
float
-cosf(float x)
+cosf (float x)
{
- return (float) cos(x);
+ return (float) cos (x);
}
#endif
#ifndef HAVE_COSHF
+#define HAVE_COSHF 1
+float coshf (float x);
+
float
-coshf(float x)
+coshf (float x)
{
- return (float) cosh(x);
+ return (float) cosh (x);
}
#endif
#ifndef HAVE_EXPF
+#define HAVE_EXPF 1
+float expf (float x);
+
float
-expf(float x)
+expf (float x)
{
- return (float) exp(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
-floorf(float x)
+fmodf (float x, float y)
{
- return (float) floor(x);
+ 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)
+frexpf (float x, int *exp)
{
- return (float) frexp(x, 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)
+hypotf (float x, float y)
{
- return (float) hypot(x, y);
+ return (float) hypot (x, y);
}
#endif
#ifndef HAVE_LOGF
+#define HAVE_LOGF 1
+float logf (float x);
+
float
-logf(float x)
+logf (float x)
{
- return (float) log(x);
+ return (float) log (x);
}
#endif
#ifndef HAVE_LOG10F
+#define HAVE_LOG10F 1
+float log10f (float x);
+
float
-log10f(float x)
+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)
{
- return (float) log10(x);
+#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)
+scalbnf (float x, int y)
{
- return (float) scalbn(x, y);
+ return (float) scalbn (x, y);
}
#endif
#ifndef HAVE_SINF
+#define HAVE_SINF 1
+float sinf (float x);
+
float
-sinf(float x)
+sinf (float x)
{
- return (float) sin(x);
+ return (float) sin (x);
}
#endif
#ifndef HAVE_SINHF
+#define HAVE_SINHF 1
+float sinhf (float x);
+
float
-sinhf(float x)
+sinhf (float x)
{
- return (float) sinh(x);
+ return (float) sinh (x);
}
#endif
#ifndef HAVE_SQRTF
+#define HAVE_SQRTF 1
+float sqrtf (float x);
+
float
-sqrtf(float x)
+sqrtf (float x)
{
- return (float) sqrt(x);
+ return (float) sqrt (x);
}
#endif
#ifndef HAVE_TANF
+#define HAVE_TANF 1
+float tanf (float x);
+
float
-tanf(float x)
+tanf (float x)
{
- return (float) tan(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
-tanhf(float x)
+truncf (float x)
{
- return (float) tanh(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)
+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))
+ 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;
#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)
+powf (float x, float y)
{
- return (float) pow(x, y);
+ return (float) pow (x, y);
}
#endif
-/* Note that if fpclassify is not defined, then NaN is not handled */
-
-/* Algorithm by Steven G. Kargl. */
#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)
+round (double x)
{
double t;
-#if defined(fpclassify)
- int i;
- i = fpclassify(x);
- if (i == FP_INFINITE || i == FP_NAN)
+ if (!isfinite (x))
return (x);
-#endif
if (x >= 0.0)
{
- t = ceil(x);
+ 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 = ceil(-x);
+ 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)
+roundf (float x)
{
float t;
-#if defined(fpclassify)
- int i;
-
- i = fpclassify(x);
- if (i == FP_INFINITE || i == FP_NAN)
+ if (!isfinite (x))
return (x);
-#endif
if (x >= 0.0)
{
- t = ceilf(x);
- if (t - x > 0.5)
- t -= 1.0;
+ t = floorf (x);
+ if (t - x <= -0.5)
+ t += 1.0;
return (t);
}
else
{
- t = ceilf(-x);
- if (t + x > 0.5)
- t -= 1.0;
+ 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
OSDN Git Service
