+
+/* 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