Import sources from 2011-03-30 iso image

1 files changed, 116 insertions, 0 deletions

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;
+}