Import sources from 2011-03-30 iso image

21 files changed, 1481 insertions, 0 deletions

diff --git a/sys/src/ape/lib/ap/math/asin.c b/sys/src/ape/lib/ap/math/asin.c
new file mode 100755
index 000000000..efd3a8cc8
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/asin.c
@@ -0,0 +1,47 @@
+/*
+	asin(arg) and acos(arg) return the arcsin, arccos,
+	respectively of their arguments.
+
+	Arctan is called after appropriate range reduction.
+ */
+
+#include <math.h>
+#include <errno.h>
+
+static double pio2	= 1.570796326794896619231e0;
+
+double
+asin(double arg)
+{
+	double temp;
+	int sign;
+
+	sign = 0;
+	if(arg < 0) {
+		arg = -arg;
+		sign++;
+	}
+	if(arg > 1) {
+		errno = EDOM;
+		return 0;
+	}
+	temp = sqrt(1 - arg*arg);
+	if(arg > 0.7)
+		temp = pio2 - atan(temp/arg);
+	else
+		temp = atan(arg/temp);
+
+	if(sign)
+		temp = -temp;
+	return temp;
+}
+
+double
+acos(double arg)
+{
+	if(arg > 1 || arg < -1) {
+		errno = EDOM;
+		return 0;
+	}
+	return pio2 - asin(arg);
+}
diff --git a/sys/src/ape/lib/ap/math/atan.c b/sys/src/ape/lib/ap/math/atan.c
new file mode 100755
index 000000000..a5674d848
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/atan.c
@@ -0,0 +1,86 @@
+/*
+	floating-point arctangent
+
+	atan returns the value of the arctangent of its
+	argument in the range [-pi/2,pi/2].
+
+	atan2 returns the arctangent of arg1/arg2
+	in the range [-pi,pi].
+
+	there are no error returns.
+
+	coefficients are #5077 from Hart & Cheney. (19.56D)
+*/
+
+#include <math.h>
+
+#define sq2p1 2.414213562373095048802e0
+#define sq2m1  .414213562373095048802e0
+#define pio2 1.570796326794896619231e0
+#define pio4  .785398163397448309615e0
+#define p4  .161536412982230228262e2
+#define p3  .26842548195503973794141e3
+#define p2  .11530293515404850115428136e4
+#define p1  .178040631643319697105464587e4
+#define p0  .89678597403663861959987488e3
+#define q4  .5895697050844462222791e2
+#define q3  .536265374031215315104235e3
+#define q2  .16667838148816337184521798e4
+#define q1  .207933497444540981287275926e4
+#define q0  .89678597403663861962481162e3
+
+
+/*
+	xatan evaluates a series valid in the
+	range [-0.414...,+0.414...].
+ */
+
+static
+double
+xatan(double arg)
+{
+	double argsq, value;
+
+	/* get denormalized add in following if range arg**10 is much smaller
+	    than q1, so check for that case
+	*/
+	if(-.01 < arg && arg < .01)
+		value = p0/q0;
+	else {
+		argsq = arg*arg;
+		value = ((((p4*argsq + p3)*argsq + p2)*argsq + p1)*argsq + p0);
+		value = value/(((((argsq + q4)*argsq + q3)*argsq + q2)*argsq + q1)*argsq + q0);
+	}
+	return value*arg;
+}
+
+/*
+	satan reduces its argument (known to be positive)
+	to the range [0,0.414...] and calls xatan.
+ */
+
+static
+double
+satan(double arg)
+{
+
+	if(arg < sq2m1)
+		return xatan(arg);
+	if(arg > sq2p1)
+		return pio2 - xatan(1.0/arg);
+	return pio4 + xatan((arg-1.0)/(arg+1.0));
+}
+
+/*
+	atan makes its argument positive and
+	calls the inner routine satan.
+ */
+
+double
+atan(double arg)
+{
+
+	if(arg > 0)
+		return satan(arg);
+	return -satan(-arg);
+}
diff --git a/sys/src/ape/lib/ap/math/atan2.c b/sys/src/ape/lib/ap/math/atan2.c
new file mode 100755
index 000000000..51c535329
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/atan2.c
@@ -0,0 +1,30 @@
+#include <math.h>
+#include <errno.h>
+/*
+	atan2 discovers what quadrant the angle
+	is in and calls atan.
+*/
+#define pio2 1.5707963267948966192313217
+#define pi   3.1415926535897932384626434;
+
+double
+atan2(double arg1, double arg2)
+{
+
+	if(arg1 == 0.0 && arg2 == 0.0){
+		errno = EDOM;
+		return 0.0;
+	}
+	if(arg1+arg2 == arg1) {
+		if(arg1 >= 0)
+			return pio2;
+		return -pio2;
+	}
+	arg1 = atan(arg1/arg2);
+	if(arg2 < 0) {
+		if(arg1 <= 0)
+			return arg1 + pi;
+		return arg1 - pi;
+	}
+	return arg1;
+}
diff --git a/sys/src/ape/lib/ap/math/erf.c b/sys/src/ape/lib/ap/math/erf.c
new file mode 100755
index 000000000..60ba4ee05
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/erf.c
@@ -0,0 +1,116 @@
+#include <math.h>
+#include <errno.h>
+/*
+	C program for floating point error function
+
+	erf(x) returns the error function of its argument
+	erfc(x) returns 1 - erf(x)
+
+	erf(x) is defined by
+	${2 over sqrt(pi)} int from 0 to x e sup {-t sup 2} dt$
+
+	the entry for erfc is provided because of the
+	extreme loss of relative accuracy if erf(x) is
+	called for large x and the result subtracted
+	from 1. (e.g. for x= 10, 12 places are lost).
+
+	There are no error returns.
+
+	Calls exp.
+
+	Coefficients for large x are #5667 from Hart & Cheney (18.72D).
+*/
+
+#define M 7
+#define N 9
+static double torp = 1.1283791670955125738961589031;
+static double p1[] = {
+	0.804373630960840172832162e5,
+	0.740407142710151470082064e4,
+	0.301782788536507577809226e4,
+	0.380140318123903008244444e2,
+	0.143383842191748205576712e2,
+	-.288805137207594084924010e0,
+	0.007547728033418631287834e0,
+};
+static double q1[]  = {
+	0.804373630960840172826266e5,
+	0.342165257924628539769006e5,
+	0.637960017324428279487120e4,
+	0.658070155459240506326937e3,
+	0.380190713951939403753468e2,
+	0.100000000000000000000000e1,
+	0.0,
+};
+static double p2[]  = {
+	0.18263348842295112592168999e4,
+	0.28980293292167655611275846e4,
+	0.2320439590251635247384768711e4,
+	0.1143262070703886173606073338e4,
+	0.3685196154710010637133875746e3,
+	0.7708161730368428609781633646e2,
+	0.9675807882987265400604202961e1,
+	0.5641877825507397413087057563e0,
+	0.0,
+};
+static double q2[]  = {
+	0.18263348842295112595576438e4,
+	0.495882756472114071495438422e4,
+	0.60895424232724435504633068e4,
+	0.4429612803883682726711528526e4,
+	0.2094384367789539593790281779e4,
+	0.6617361207107653469211984771e3,
+	0.1371255960500622202878443578e3,
+	0.1714980943627607849376131193e2,
+	1.0,
+};
+
+double erfc(double);
+
+double
+erf(double arg)
+{
+	int sign;
+	double argsq;
+	double d, n;
+	int i;
+
+	errno = 0;
+	sign = 1;
+	if(arg < 0) {
+		arg = -arg;
+		sign = -1;
+	}
+	if(arg < 0.5) {
+		argsq = arg*arg;
+		for(n=0,d=0,i=M-1; i>=0; i--) {
+			n = n*argsq + p1[i];
+			d = d*argsq + q1[i];
+		}
+		return sign*torp*arg*n/d;
+	}
+	if(arg >= 10)
+		return sign;
+	return sign*(1 - erfc(arg));
+}
+
+double
+erfc(double arg)
+{
+	double n, d;
+	int i;
+
+	errno = 0;
+	if(arg < 0)
+		return 2 - erfc(-arg);
+	if(arg < 0.5)
+		return 1 - erf(arg);
+	if(arg >= 10)
+		return 0;
+
+	for(n=0,d=0,i=N-1; i>=0; i--) {
+		n = n*arg + p2[i];
+		d = d*arg + q2[i];
+	}
+	return exp(-arg*arg)*n/d;
+}
diff --git a/sys/src/ape/lib/ap/math/exp.c b/sys/src/ape/lib/ap/math/exp.c
new file mode 100755
index 000000000..06d21d70f
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/exp.c
@@ -0,0 +1,44 @@
+/*
+	exp returns the exponential function of its
+	floating-point argument.
+
+	The coefficients are #1069 from Hart and Cheney. (22.35D)
+*/
+
+#include <math.h>
+#include <errno.h>
+
+#define	p0  .2080384346694663001443843411e7
+#define	p1  .3028697169744036299076048876e5
+#define	p2  .6061485330061080841615584556e2
+#define	q0  .6002720360238832528230907598e7
+#define	q1  .3277251518082914423057964422e6
+#define	q2  .1749287689093076403844945335e4
+#define	log2e  1.4426950408889634073599247
+#define	sqrt2  1.4142135623730950488016887
+#define	maxf  10000
+
+double
+exp(double arg)
+{
+	double fract, temp1, temp2, xsq;
+	int ent;
+
+	if(arg == 0)
+		return 1;
+	if(arg < -maxf) {
+		errno = ERANGE;
+		return 0;
+	}
+	if(arg > maxf) {
+		errno = ERANGE;
+		return HUGE_VAL;
+	}
+	arg *= log2e;
+	ent = floor(arg);
+	fract = (arg-ent) - 0.5;
+	xsq = fract*fract;
+	temp1 = ((p2*xsq+p1)*xsq+p0)*fract;
+	temp2 = ((xsq+q2)*xsq+q1)*xsq + q0;
+	return ldexp(sqrt2*(temp2+temp1)/(temp2-temp1), ent);
+}
diff --git a/sys/src/ape/lib/ap/math/fabs.c b/sys/src/ape/lib/ap/math/fabs.c
new file mode 100755
index 000000000..5566e165b
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/fabs.c
@@ -0,0 +1,10 @@
+#include <math.h>
+
+double
+fabs(double arg)
+{
+
+	if(arg < 0)
+		return -arg;
+	return arg;
+}
diff --git a/sys/src/ape/lib/ap/math/floor.c b/sys/src/ape/lib/ap/math/floor.c
new file mode 100755
index 000000000..dfab000c1
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/floor.c
@@ -0,0 +1,26 @@
+#include <math.h>
+/*
+ * floor and ceil-- greatest integer <= arg
+ * (resp least >=)
+ */
+
+double
+floor(double d)
+{
+	double fract;
+
+	if(d < 0) {
+		fract = modf(-d, &d);
+		if(fract != 0.0)
+			d += 1;
+		d = -d;
+	} else
+		modf(d, &d);
+	return d;
+}
+
+double
+ceil(double d)
+{
+	return -floor(-d);
+}
diff --git a/sys/src/ape/lib/ap/math/fmod.c b/sys/src/ape/lib/ap/math/fmod.c
new file mode 100755
index 000000000..876c64c39
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/fmod.c
@@ -0,0 +1,27 @@
+/* floating-point mod function without infinity or NaN checking */
+#include <math.h>
+double
+fmod (double x, double y)
+{
+	int sign = 0, yexp;
+	double r, yfr;
+
+	if (y == 0)
+		return 0;
+	if (y < 0)
+		y = -y;
+	yfr = frexp (y, &yexp);
+	if (x < 0) {
+		sign = 1;
+		r = -x;
+	} else
+		r = x;
+	while (r >= y) {
+		int rexp;
+		double rfr = frexp (r, &rexp);
+		r -= ldexp (y, rexp - yexp - (rfr < yfr));
+	}
+	if (sign)
+		r = -r;
+	return r;
+}
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));
+}
diff --git a/sys/src/ape/lib/ap/math/hypot.c b/sys/src/ape/lib/ap/math/hypot.c
new file mode 100755
index 000000000..25457d8a9
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/hypot.c
@@ -0,0 +1,29 @@
+#include <math.h>
+/*
+ * sqrt(a^2 + b^2)
+ *	(but carefully)
+ */
+
+double
+hypot(double a, double b)
+{
+	double t;
+
+	if(a < 0)
+		a = -a;
+	if(b < 0)
+		b = -b;
+	if(a > b) {
+		t = a;
+		a = b;
+		b = t;
+	}
+	if(b == 0) 
+		return 0;
+	a /= b;
+	/*
+	 * pathological overflow possible
+	 * in the next line.
+	 */
+	return b * sqrt(1 + a*a);
+}
diff --git a/sys/src/ape/lib/ap/math/j0.c b/sys/src/ape/lib/ap/math/j0.c
new file mode 100755
index 000000000..7dea3c9d1
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/j0.c
@@ -0,0 +1,188 @@
+#include <math.h>
+#include <errno.h>
+/*
+	floating point Bessel's function
+	of the first and second kinds
+	of order zero
+
+	j0(x) returns the value of J0(x)
+	for all real values of x.
+
+	There are no error returns.
+	Calls sin, cos, sqrt.
+
+	There is a niggling bug in J0 which
+	causes errors up to 2e-16 for x in the
+	interval [-8,8].
+	The bug is caused by an inappropriate order
+	of summation of the series.  rhm will fix it
+	someday.
+
+	Coefficients are from Hart & Cheney.
+	#5849 (19.22D)
+	#6549 (19.25D)
+	#6949 (19.41D)
+
+	y0(x) returns the value of Y0(x)
+	for positive real values of x.
+	For x<=0, error number EDOM is set and a
+	large negative value is returned.
+
+	Calls sin, cos, sqrt, log, j0.
+
+	The values of Y0 have not been checked
+	to more than ten places.
+
+	Coefficients are from Hart & Cheney.
+	#6245 (18.78D)
+	#6549 (19.25D)
+	#6949 (19.41D)
+*/
+
+static double pzero, qzero;
+static double tpi	= .6366197723675813430755350535e0;
+static double pio4	= .7853981633974483096156608458e0;
+static double p1[] = {
+	0.4933787251794133561816813446e21,
+	-.1179157629107610536038440800e21,
+	0.6382059341072356562289432465e19,
+	-.1367620353088171386865416609e18,
+	0.1434354939140344111664316553e16,
+	-.8085222034853793871199468171e13,
+	0.2507158285536881945555156435e11,
+	-.4050412371833132706360663322e8,
+	0.2685786856980014981415848441e5,
+};
+static double q1[] = {
+	0.4933787251794133562113278438e21,
+	0.5428918384092285160200195092e19,
+	0.3024635616709462698627330784e17,
+	0.1127756739679798507056031594e15,
+	0.3123043114941213172572469442e12,
+	0.6699987672982239671814028660e9,
+	0.1114636098462985378182402543e7,
+	0.1363063652328970604442810507e4,
+	1.0
+};
+static double p2[] = {
+	0.5393485083869438325262122897e7,
+	0.1233238476817638145232406055e8,
+	0.8413041456550439208464315611e7,
+	0.2016135283049983642487182349e7,
+	0.1539826532623911470917825993e6,
+	0.2485271928957404011288128951e4,
+	0.0,
+};
+static double q2[] = {
+	0.5393485083869438325560444960e7,
+	0.1233831022786324960844856182e8,
+	0.8426449050629797331554404810e7,
+	0.2025066801570134013891035236e7,
+	0.1560017276940030940592769933e6,
+	0.2615700736920839685159081813e4,
+	1.0,
+};
+static double p3[] = {
+	-.3984617357595222463506790588e4,
+	-.1038141698748464093880530341e5,
+	-.8239066313485606568803548860e4,
+	-.2365956170779108192723612816e4,
+	-.2262630641933704113967255053e3,
+	-.4887199395841261531199129300e1,
+	0.0,
+};
+static double q3[] = {
+	0.2550155108860942382983170882e6,
+	0.6667454239319826986004038103e6,
+	0.5332913634216897168722255057e6,
+	0.1560213206679291652539287109e6,
+	0.1570489191515395519392882766e5,
+	0.4087714673983499223402830260e3,
+	1.0,
+};
+static double p4[] = {
+	-.2750286678629109583701933175e20,
+	0.6587473275719554925999402049e20,
+	-.5247065581112764941297350814e19,
+	0.1375624316399344078571335453e18,
+	-.1648605817185729473122082537e16,
+	0.1025520859686394284509167421e14,
+	-.3436371222979040378171030138e11,
+	0.5915213465686889654273830069e8,
+	-.4137035497933148554125235152e5,
+};
+static double q4[] = {
+	0.3726458838986165881989980e21,
+	0.4192417043410839973904769661e19,
+	0.2392883043499781857439356652e17,
+	0.9162038034075185262489147968e14,
+	0.2613065755041081249568482092e12,
+	0.5795122640700729537480087915e9,
+	0.1001702641288906265666651753e7,
+	0.1282452772478993804176329391e4,
+	1.0,
+};
+
+static
+asympt(double arg)
+{
+	double zsq, n, d;
+	int i;
+
+	zsq = 64 / (arg*arg);
+	for(n=0,d=0,i=6;i>=0;i--) {
+		n = n*zsq + p2[i];
+		d = d*zsq + q2[i];
+	}
+	pzero = n/d;
+	for(n=0,d=0,i=6;i>=0;i--) {
+		n = n*zsq + p3[i];
+		d = d*zsq + q3[i];
+	}
+	qzero = (8/arg)*(n/d);
+}
+
+double
+j0(double arg)
+{
+	double argsq, n, d;
+	int i;
+
+	if(arg < 0)
+		arg = -arg;
+	if(arg > 8) {
+		asympt(arg);
+		n = arg - pio4;
+		return sqrt(tpi/arg)*(pzero*cos(n) - qzero*sin(n));
+	}
+	argsq = arg*arg;
+	for(n=0,d=0,i=8;i>=0;i--) {
+		n = n*argsq + p1[i];
+		d = d*argsq + q1[i];
+	}
+	return n/d;
+}
+
+double
+y0(double arg)
+{
+	double argsq, n, d;
+	int i;
+
+	errno = 0;
+	if(arg <= 0) {
+		errno = EDOM;
+		return(-HUGE_VAL);
+	}
+	if(arg > 8) {
+		asympt(arg);
+		n = arg - pio4;
+		return sqrt(tpi/arg)*(pzero*sin(n) + qzero*cos(n));
+	}
+	argsq = arg*arg;
+	for(n=0,d=0,i=8;i>=0;i--) {
+		n = n*argsq + p4[i];
+		d = d*argsq + q4[i];
+	}
+	return n/d + tpi*j0(arg)*log(arg);
+}
diff --git a/sys/src/ape/lib/ap/math/j1.c b/sys/src/ape/lib/ap/math/j1.c
new file mode 100755
index 000000000..38c73fe2b
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/j1.c
@@ -0,0 +1,195 @@
+#include <math.h>
+#include <errno.h>
+/*
+	floating point Bessel's function
+	of the first and second kinds
+	of order one
+
+	j1(x) returns the value of J1(x)
+	for all real values of x.
+
+	There are no error returns.
+	Calls sin, cos, sqrt.
+
+	There is a niggling bug in J1 which
+	causes errors up to 2e-16 for x in the
+	interval [-8,8].
+	The bug is caused by an inappropriate order
+	of summation of the series.  rhm will fix it
+	someday.
+
+	Coefficients are from Hart & Cheney.
+	#6050 (20.98D)
+	#6750 (19.19D)
+	#7150 (19.35D)
+
+	y1(x) returns the value of Y1(x)
+	for positive real values of x.
+	For x<=0, error number EDOM is set and a
+	large negative value is returned.
+
+	Calls sin, cos, sqrt, log, j1.
+
+	The values of Y1 have not been checked
+	to more than ten places.
+
+	Coefficients are from Hart & Cheney.
+	#6447 (22.18D)
+	#6750 (19.19D)
+	#7150 (19.35D)
+*/
+
+static double pzero, qzero;
+static double tpi	= .6366197723675813430755350535e0;
+static double pio4	= .7853981633974483096156608458e0;
+static double p1[] = {
+	0.581199354001606143928050809e21,
+	-.6672106568924916298020941484e20,
+	0.2316433580634002297931815435e19,
+	-.3588817569910106050743641413e17,
+	0.2908795263834775409737601689e15,
+	-.1322983480332126453125473247e13,
+	0.3413234182301700539091292655e10,
+	-.4695753530642995859767162166e7,
+	0.2701122710892323414856790990e4,
+};
+static double q1[] = {
+	0.1162398708003212287858529400e22,
+	0.1185770712190320999837113348e20,
+	0.6092061398917521746105196863e17,
+	0.2081661221307607351240184229e15,
+	0.5243710262167649715406728642e12,
+	0.1013863514358673989967045588e10,
+	0.1501793594998585505921097578e7,
+	0.1606931573481487801970916749e4,
+	1.0,
+};
+static double p2[] = {
+	-.4435757816794127857114720794e7,
+	-.9942246505077641195658377899e7,
+	-.6603373248364939109255245434e7,
+	-.1523529351181137383255105722e7,
+	-.1098240554345934672737413139e6,
+	-.1611616644324610116477412898e4,
+	0.0,
+};
+static double q2[] = {
+	-.4435757816794127856828016962e7,
+	-.9934124389934585658967556309e7,
+	-.6585339479723087072826915069e7,
+	-.1511809506634160881644546358e7,
+	-.1072638599110382011903063867e6,
+	-.1455009440190496182453565068e4,
+	1.0,
+};
+static double p3[] = {
+	0.3322091340985722351859704442e5,
+	0.8514516067533570196555001171e5,
+	0.6617883658127083517939992166e5,
+	0.1849426287322386679652009819e5,
+	0.1706375429020768002061283546e4,
+	0.3526513384663603218592175580e2,
+	0.0,
+};
+static double q3[] = {
+	0.7087128194102874357377502472e6,
+	0.1819458042243997298924553839e7,
+	0.1419460669603720892855755253e7,
+	0.4002944358226697511708610813e6,
+	0.3789022974577220264142952256e5,
+	0.8638367769604990967475517183e3,
+	1.0,
+};
+static double p4[] = {
+	-.9963753424306922225996744354e23,
+	0.2655473831434854326894248968e23,
+	-.1212297555414509577913561535e22,
+	0.2193107339917797592111427556e20,
+	-.1965887462722140658820322248e18,
+	0.9569930239921683481121552788e15,
+	-.2580681702194450950541426399e13,
+	0.3639488548124002058278999428e10,
+	-.2108847540133123652824139923e7,
+	0.0,
+};
+static double q4[] = {
+	0.5082067366941243245314424152e24,
+	0.5435310377188854170800653097e22,
+	0.2954987935897148674290758119e20,
+	0.1082258259408819552553850180e18,
+	0.2976632125647276729292742282e15,
+	0.6465340881265275571961681500e12,
+	0.1128686837169442121732366891e10,
+	0.1563282754899580604737366452e7,
+	0.1612361029677000859332072312e4,
+	1.0,
+};
+
+static
+asympt(double arg)
+{
+	double zsq, n, d;
+	int i;
+
+	zsq = 64/(arg*arg);
+	for(n=0,d=0,i=6;i>=0;i--) {
+		n = n*zsq + p2[i];
+		d = d*zsq + q2[i];
+	}
+	pzero = n/d;
+	for(n=0,d=0,i=6;i>=0;i--) {
+		n = n*zsq + p3[i];
+		d = d*zsq + q3[i];
+	}
+	qzero = (8/arg)*(n/d);
+}
+
+double
+j1(double arg)
+{
+	double xsq, n, d, x;
+	int i;
+
+	x = arg;
+	if(x < 0)
+		x = -x;
+	if(x > 8) {
+		asympt(x);
+		n = x - 3*pio4;
+		n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n));
+		if(arg < 0)
+			n = -n;
+		return n;
+	}
+	xsq = x*x;
+	for(n=0,d=0,i=8;i>=0;i--) {
+		n = n*xsq + p1[i];
+		d = d*xsq + q1[i];
+	}
+	return arg*n/d;
+}
+
+double
+y1(double arg)
+{
+	double xsq, n, d, x;
+	int i;
+
+	errno = 0;
+	x = arg;
+	if(x <= 0) {
+		errno = EDOM;
+		return -HUGE_VAL;
+	}
+	if(x > 8) {
+		asympt(x);
+		n = x - 3*pio4;
+		return sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n));
+	}
+	xsq = x*x;
+	for(n=0,d=0,i=9;i>=0;i--) {
+		n = n*xsq + p4[i];
+		d = d*xsq + q4[i];
+	}
+	return x*n/d + tpi*(j1(x)*log(x)-1/x);
+}
diff --git a/sys/src/ape/lib/ap/math/jn.c b/sys/src/ape/lib/ap/math/jn.c
new file mode 100755
index 000000000..55e985d48
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/jn.c
@@ -0,0 +1,116 @@
+#include <math.h>
+#include <errno.h>
+
+/*
+	floating point Bessel's function of
+	the first and second kinds and of
+	integer order.
+
+	int n;
+	double x;
+	jn(n,x);
+
+	returns the value of Jn(x) for all
+	integer values of n and all real values
+	of x.
+
+	There are no error returns.
+	Calls j0, j1.
+
+	For n=0, j0(x) is called,
+	for n=1, j1(x) is called,
+	for n<x, forward recursion us used starting
+	from values of j0(x) and j1(x).
+	for n>x, a continued fraction approximation to
+	j(n,x)/j(n-1,x) is evaluated and then backward
+	recursion is used starting from a supposed value
+	for j(n,x). The resulting value of j(0,x) is
+	compared with the actual value to correct the
+	supposed value of j(n,x).
+
+	yn(n,x) is similar in all respects, except
+	that forward recursion is used for all
+	values of n>1.
+*/
+
+double	j0(double);
+double	j1(double);
+double	y0(double);
+double	y1(double);
+
+double
+jn(int n, double x)
+{
+	int i;
+	double a, b, temp;
+	double xsq, t;
+
+	if(n < 0) {
+		n = -n;
+		x = -x;
+	}
+	if(n == 0)
+		return j0(x);
+	if(n == 1)
+		return j1(x);
+	if(x == 0)
+		return 0;
+	if(n > x)
+		goto recurs;
+
+	a = j0(x);
+	b = j1(x);
+	for(i=1; i<n; i++) {
+		temp = b;
+		b = (2*i/x)*b - a;
+		a = temp;
+	}
+	return b;
+
+recurs:
+	xsq = x*x;
+	for(t=0,i=n+16; i>n; i--)
+		t = xsq/(2*i - t);
+	t = x/(2*n-t);
+
+	a = t;
+	b = 1;
+	for(i=n-1; i>0; i--) {
+		temp = b;
+		b = (2*i/x)*b - a;
+		a = temp;
+	}
+	return t*j0(x)/b;
+}
+
+double
+yn(int n, double x)
+{
+	int i;
+	int sign;
+	double a, b, temp;
+
+	if (x <= 0) {
+		errno = EDOM;
+		return -HUGE_VAL;
+	}
+	sign = 1;
+	if(n < 0) {
+		n = -n;
+		if(n%2 == 1)
+			sign = -1;
+	}
+	if(n == 0)
+		return y0(x);
+	if(n == 1)
+		return sign*y1(x);
+
+	a = y0(x);
+	b = y1(x);
+	for(i=1; i<n; i++) {
+		temp = b;
+		b = (2*i/x)*b - a;
+		a = temp;
+	}
+	return sign*b;
+}
diff --git a/sys/src/ape/lib/ap/math/log.c b/sys/src/ape/lib/ap/math/log.c
new file mode 100755
index 000000000..bc1c2e1cb
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/log.c
@@ -0,0 +1,62 @@
+/*
+	log returns the natural logarithm of its floating
+	point argument.
+
+	The coefficients are #2705 from Hart & Cheney. (19.38D)
+
+	It calls frexp.
+*/
+
+#include <math.h>
+#include <errno.h>
+
+#define	log2    0.693147180559945309e0
+#define	ln10o1  .4342944819032518276511
+#define	sqrto2  0.707106781186547524e0
+#define	p0      -.240139179559210510e2
+#define	p1      0.309572928215376501e2
+#define	p2      -.963769093377840513e1
+#define	p3      0.421087371217979714e0
+#define	q0      -.120069589779605255e2
+#define	q1      0.194809660700889731e2
+#define	q2      -.891110902798312337e1
+
+double
+log(double arg)
+{
+	double x, z, zsq, temp;
+	int exp;
+
+	if(arg <= 0) {
+		errno = (arg==0)? ERANGE : EDOM;
+		return -HUGE_VAL;
+	}
+	x = frexp(arg, &exp);
+	while(x < 0.5) {
+		x *= 2;
+		exp--;
+	}
+	if(x < sqrto2) {
+		x *= 2;
+		exp--;
+	}
+
+	z = (x-1) / (x+1);
+	zsq = z*z;
+
+	temp = ((p3*zsq + p2)*zsq + p1)*zsq + p0;
+	temp = temp/(((zsq + q2)*zsq + q1)*zsq + q0);
+	temp = temp*z + exp*log2;
+	return temp;
+}
+
+double
+log10(double arg)
+{
+
+	if(arg <= 0) {
+		errno = (arg==0)? ERANGE : EDOM;
+		return -HUGE_VAL;
+	}
+	return log(arg) * ln10o1;
+}
diff --git a/sys/src/ape/lib/ap/math/mkfile b/sys/src/ape/lib/ap/math/mkfile
new file mode 100755
index 000000000..5e2d7483f
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/mkfile
@@ -0,0 +1,27 @@
+APE=/sys/src/ape
+<$APE/config
+LIB=/$objtype/lib/ape/libap.a
+OFILES=\
+	asin.$O\
+	atan.$O\
+	atan2.$O\
+	erf.$O\
+	exp.$O\
+	fabs.$O\
+	floor.$O\
+	fmod.$O\
+	gamma.$O\
+	hypot.$O\
+	j0.$O\
+	j1.$O\
+	jn.$O\
+	log.$O\
+	pow.$O\
+	sin.$O\
+	sinh.$O\
+	tan.$O\
+	tanh.$O\
+
+</sys/src/cmd/mksyslib
+
+CFLAGS=-c -D_POSIX_SOURCE
diff --git a/sys/src/ape/lib/ap/math/modf.c b/sys/src/ape/lib/ap/math/modf.c
new file mode 100755
index 000000000..4326cf44a
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/modf.c
@@ -0,0 +1,53 @@
+#include <math.h>
+#include <errno.h>
+
+/* modf suitable for IEEE double-precision */
+
+#define	MASK	0x7ffL
+#define SIGN	0x80000000
+#define	SHIFT	20
+#define	BIAS	1022L
+
+typedef	union
+{
+	double	d;
+	struct
+	{
+		long	ms;
+		long	ls;
+	} i;
+} Cheat;
+
+double
+modf(double d, double *ip)
+{
+	Cheat x;
+	int e;
+
+	if(-1 < d && d < 1) {
+		*ip = 0;
+		return d;
+	}
+	x.d = d;
+	x.i.ms &= ~SIGN;
+	e = (x.i.ms >> SHIFT) & MASK;
+	if(e == MASK || e == 0){
+		errno = EDOM;
+		*ip = (d > 0)? HUGE_VAL : -HUGE_VAL;
+		return 0;
+	}
+	e -= BIAS;
+	if(e <= SHIFT+1) {
+		x.i.ms &= ~(0x1fffffL >> e);
+		x.i.ls = 0;
+	} else
+	if(e <= SHIFT+33)
+		x.i.ls &= ~(0x7fffffffL >> (e-SHIFT-2));
+	if(d > 0){
+		*ip = x.d;
+		return d - x.d;
+	}else{
+		*ip = -x.d;
+		return d + x.d;
+	}
+}
diff --git a/sys/src/ape/lib/ap/math/pow.c b/sys/src/ape/lib/ap/math/pow.c
new file mode 100755
index 000000000..9512a649b
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/pow.c
@@ -0,0 +1,76 @@
+#include <math.h>
+#include <errno.h>
+#include <limits.h>
+
+double
+pow(double x, double y) /* return x ^ y (exponentiation) */
+{
+	double xy, y1, ye;
+	long i;
+	int ex, ey, flip;
+
+	if(y == 0.0)
+		return 1.0;
+
+	flip = 0;
+	if(y < 0.){
+		y = -y;
+		flip = 1;
+	}
+	y1 = modf(y, &ye);
+	if(y1 != 0.0){
+		if(x <= 0.)
+			goto zreturn;
+		if(y1 > 0.5) {
+			y1 -= 1.;
+			ye += 1.;
+		}
+		xy = exp(y1 * log(x));
+	}else
+		xy = 1.0;
+	if(ye > LONG_MAX){
+		if(x <= 0.){
+ zreturn:
+			if(x || flip)
+				errno = EDOM;
+			return 0.;
+		}
+		if(flip){
+			if(y == .5)
+				return 1./sqrt(x);
+			y = -y;
+		}else if(y == .5)
+			return sqrt(x);
+		return exp(y * log(x));
+	}
+	x = frexp(x, &ex);
+	ey = 0;
+	i = ye;
+	if(i)
+		for(;;){
+			if(i & 1){
+				xy *= x;
+				ey += ex;
+			}
+			i >>= 1;
+			if(i == 0)
+				break;
+			x *= x;
+			ex <<= 1;
+			if(x < .5){
+				x += x;
+				ex -= 1;
+			}
+		}
+	if(flip){
+		xy = 1. / xy;
+		ey = -ey;
+	}
+	errno = 0;
+	x = ldexp(xy, ey);
+	if(errno && ey < 0) {
+		errno = 0;
+		x = 0.;
+	}
+	return x;
+}
diff --git a/sys/src/ape/lib/ap/math/sin.c b/sys/src/ape/lib/ap/math/sin.c
new file mode 100755
index 000000000..694a223f4
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/sin.c
@@ -0,0 +1,68 @@
+/*
+	C program for floating point sin/cos.
+	Calls modf.
+	There are no error exits.
+	Coefficients are #3370 from Hart & Cheney (18.80D).
+*/
+
+#include <math.h>
+
+#define	PIO2	1.570796326794896619231e0
+#define p0      .1357884097877375669092680e8
+#define p1     -.4942908100902844161158627e7
+#define p2      .4401030535375266501944918e6
+#define p3     -.1384727249982452873054457e5
+#define p4      .1459688406665768722226959e3
+#define q0      .8644558652922534429915149e7
+#define q1      .4081792252343299749395779e6
+#define q2      .9463096101538208180571257e4
+#define q3      .1326534908786136358911494e3
+
+static
+double
+sinus(double arg, int quad)
+{
+	double e, f, ysq, x, y, temp1, temp2;
+	int k;
+
+	x = arg;
+	if(x < 0) {
+		x = -x;
+		quad += 2;
+	}
+	x *= 1/PIO2;	/* underflow? */
+	if(x > 32764) {
+		y = modf(x, &e);
+		e += quad;
+		modf(0.25*e, &f);
+		quad = e - 4*f;
+	} else {
+		k = x;
+		y = x - k;
+		quad += k;
+		quad &= 3;
+	}
+	if(quad & 1)
+		y = 1-y;
+	if(quad > 1)
+		y = -y;
+
+	ysq = y*y;
+	temp1 = ((((p4*ysq+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y;
+	temp2 = ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0);
+	return temp1/temp2;
+}
+
+double
+cos(double arg)
+{
+	if(arg < 0)
+		arg = -arg;
+	return sinus(arg, 1);
+}
+
+double
+sin(double arg)
+{
+	return sinus(arg, 0);
+}
diff --git a/sys/src/ape/lib/ap/math/sinh.c b/sys/src/ape/lib/ap/math/sinh.c
new file mode 100755
index 000000000..58c261b76
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/sinh.c
@@ -0,0 +1,71 @@
+/*
+	sinh(arg) returns the hyperbolic sine of its floating-
+	point argument.
+
+	The exponential function is called for arguments
+	greater in magnitude than 0.5.
+
+	A series is used for arguments smaller in magnitude than 0.5.
+	The coefficients are #2029 from Hart & Cheney. (20.36D)
+
+	cosh(arg) is computed from the exponential function for
+	all arguments.
+*/
+
+#include <math.h>
+#include <errno.h>
+
+static double p0  = -0.6307673640497716991184787251e+6;
+static double p1  = -0.8991272022039509355398013511e+5;
+static double p2  = -0.2894211355989563807284660366e+4;
+static double p3  = -0.2630563213397497062819489e+2;
+static double q0  = -0.6307673640497716991212077277e+6;
+static double q1   = 0.1521517378790019070696485176e+5;
+static double q2  = -0.173678953558233699533450911e+3;
+
+double
+sinh(double arg)
+{
+	double temp, argsq;
+	int sign;
+
+	sign = 1;
+	if(arg < 0) {
+		arg = - arg;
+		sign = -1;
+	}
+	if(arg > 21) {
+		if(arg >= HUGE_VAL){
+			errno = ERANGE;
+			temp = HUGE_VAL;
+		} else
+			temp = exp(arg)/2;
+		if(sign > 0)
+			return temp;
+		else
+			return -temp;
+	}
+
+	if(arg > 0.5)
+		return sign*(exp(arg) - exp(-arg))/2;
+
+	argsq = arg*arg;
+	temp = (((p3*argsq+p2)*argsq+p1)*argsq+p0)*arg;
+	temp /= (((argsq+q2)*argsq+q1)*argsq+q0);
+	return sign*temp;
+}
+
+double
+cosh(double arg)
+{
+	if(arg < 0)
+		arg = - arg;
+	if(arg > 21) {
+		if(arg >= HUGE_VAL){
+			errno = ERANGE;
+			return HUGE_VAL;
+		} else
+			return(exp(arg)/2);
+	}
+	return (exp(arg) + exp(-arg))/2;
+}
diff --git a/sys/src/ape/lib/ap/math/tan.c b/sys/src/ape/lib/ap/math/tan.c
new file mode 100755
index 000000000..5147f8bc1
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/tan.c
@@ -0,0 +1,70 @@
+/*
+	floating point tangent
+
+	A series is used after range reduction.
+	Coefficients are #4285 from Hart & Cheney. (19.74D)
+ */
+
+#include <math.h>
+#include <errno.h>
+
+static double invpi	  = 1.27323954473516268;
+static double p0	 = -0.1306820264754825668269611177e+5;
+static double p1	  = 0.1055970901714953193602353981e+4;
+static double p2	 = -0.1550685653483266376941705728e+2;
+static double p3	  = 0.3422554387241003435328470489e-1;
+static double p4	  = 0.3386638642677172096076369e-4;
+static double q0	 = -0.1663895238947119001851464661e+5;
+static double q1	  = 0.4765751362916483698926655581e+4;
+static double q2	 = -0.1555033164031709966900124574e+3;
+
+double
+tan(double arg)
+{
+	double sign, temp, e, x, xsq;
+	int flag, i;
+
+	flag = 0;
+	sign = 1;
+	if(arg < 0){
+		arg = -arg;
+		sign = -1;
+	}
+	arg = arg*invpi;   /* overflow? */
+	x = modf(arg, &e);
+	i = e;
+	switch(i%4) {
+	case 1:
+		x = 1 - x;
+		flag = 1;
+		break;
+
+	case 2:
+		sign = - sign;
+		flag = 1;
+		break;
+
+	case 3:
+		x = 1 - x;
+		sign = - sign;
+		break;
+
+	case 0:
+		break;
+	}
+
+	xsq = x*x;
+	temp = ((((p4*xsq+p3)*xsq+p2)*xsq+p1)*xsq+p0)*x;
+	temp = temp/(((xsq+q2)*xsq+q1)*xsq+q0);
+
+	if(flag == 1) {
+		if(temp == 0) {
+			errno = EDOM;
+			if (sign > 0)
+				return HUGE_VAL;
+			return -HUGE_VAL;
+		}
+		temp = 1/temp;
+	}
+	return sign*temp;
+}
diff --git a/sys/src/ape/lib/ap/math/tanh.c b/sys/src/ape/lib/ap/math/tanh.c
new file mode 100755
index 000000000..cad68362b
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/tanh.c
@@ -0,0 +1,24 @@
+#include <math.h>
+
+/*
+	tanh(arg) computes the hyperbolic tangent of its floating
+	point argument.
+
+	sinh and cosh are called except for large arguments, which
+	would cause overflow improperly.
+ */
+
+double
+tanh(double arg)
+{
+
+	if(arg < 0) {
+		arg = -arg;
+		if(arg > 21)
+			return -1;
+		return -sinh(arg)/cosh(arg);
+	}
+	if(arg > 21)
+		return 1;
+	return sinh(arg)/cosh(arg);
+}