Import sources from 2011-03-30 iso image

1 files changed, 116 insertions, 0 deletions

diff --git a/sys/src/ape/lib/ap/math/gamma.c b/sys/src/ape/lib/ap/math/gamma.c
new file mode 100755
index 000000000..7b541042b
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/gamma.c
@@ -0,0 +1,116 @@
+#include <math.h>
+#include <errno.h>
+
+/*
+	C program for floating point log gamma function
+
+	gamma(x) computes the log of the absolute
+	value of the gamma function.
+	The sign of the gamma function is returned in the
+	external quantity signgam.
+
+	The coefficients for expansion around zero
+	are #5243 from Hart & Cheney; for expansion
+	around infinity they are #5404.
+
+	Calls log and sin.
+*/
+
+int	signgam;
+static double goobie	= 0.9189385332046727417803297;
+static double pi	= 3.1415926535897932384626434;
+
+#define M 6
+#define N 8
+static double p1[] = {
+	0.83333333333333101837e-1,
+	-.277777777735865004e-2,
+	0.793650576493454e-3,
+	-.5951896861197e-3,
+	0.83645878922e-3,
+	-.1633436431e-2,
+};
+static double p2[] = {
+	-.42353689509744089647e5,
+	-.20886861789269887364e5,
+	-.87627102978521489560e4,
+	-.20085274013072791214e4,
+	-.43933044406002567613e3,
+	-.50108693752970953015e2,
+	-.67449507245925289918e1,
+	0.0,
+};
+static double q2[] = {
+	-.42353689509744090010e5,
+	-.29803853309256649932e4,
+	0.99403074150827709015e4,
+	-.15286072737795220248e4,
+	-.49902852662143904834e3,
+	0.18949823415702801641e3,
+	-.23081551524580124562e2,
+	0.10000000000000000000e1,
+};
+
+static double
+asym(double arg)
+{
+	double n, argsq;
+	int i;
+
+	argsq = 1 / (arg*arg);
+	n = 0;
+	for(i=M-1; i>=0; i--)
+		n = n*argsq + p1[i];
+	return (arg-.5)*log(arg) - arg + goobie + n/arg;
+}
+
+static double
+pos(double arg)
+{
+	double n, d, s;
+	int i;
+
+	if(arg < 2)
+		return pos(arg+1)/arg;
+	if(arg > 3)
+		return (arg-1)*pos(arg-1);
+
+	s = arg - 2;
+	n = 0;
+	d = 0;
+	for(i=N-1; i>=0; i--){
+		n = n*s + p2[i];
+		d = d*s + q2[i];
+	}
+	return n/d;
+}
+
+static double
+neg(double arg)
+{
+	double temp;
+
+	arg = -arg;
+	temp = sin(pi*arg);
+	if(temp == 0) {
+		errno = EDOM;
+		return HUGE_VAL;
+	}
+	if(temp < 0)
+		temp = -temp;
+	else
+		signgam = -1;
+	return -log(arg*pos(arg)*temp/pi);
+}
+
+double
+gamma(double arg)
+{
+
+	signgam = 1;
+	if(arg <= 0)
+		return neg(arg);
+	if(arg > 8)
+		return asym(arg);
+	return log(pos(arg));
+}