/*  jn.c  */

asm ("        export Libm");

mesa double Libm←y0(), Libm←y1();

#define y0(a) (Libm←y0(a))
#define y1(a) (Libm←y1(a))

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

mesa double Libm←sin(), Libm←cos(), Libm←tan();
mesa double Libm←asin(), Libm←acos(), Libm←atan(), Libm←atan2();
mesa double Libm←sinh(), Libm←cosh(), Libm←tanh();
mesa double Libm←asinh(), Libm←acosh(), Libm←atanh();
mesa double Libm←log(), Libm←log1p(), Libm←lgamma();
mesa double Libm←exp(), Libm←frexp(), Libm←ldexp(), Libm←expm1(), Libm←pow();
mesa double Libm←fmod(), Libm←modf();
mesa double Libm←floor(), Libm←ceil(), Libm←rint();
mesa double Libm←cabs(), Libm←hypot();
mesa double Libm←sqrt(), Libm←cbrt(); 
mesa double Libm←j0(), Libm←j1();
mesa double Libm←fabs();
mesa int Libm←abs();

#define sin(x) (Libm←sin(x))
#define cos(x) (Libm←cos(x))
#define tan(x) (Libm←tan(x))
#define asin(x) (Libm←asin(x))
#define acos(x) (Libm←acos(x))
#define atan(x) (Libm←atan(x))
#define atan2(x,y) (Libm←atan2(x,y))
#define sinh(x) (Libm←sinh(x))
#define cosh(x) (Libm←cosh(x))
#define tanh(x) (Libm←tanh(x))
#define asinh(x) (Libm←asinh(x))
#define acosh(x) (Libm←acosh(x))
#define atanh(x) (Libm←atanh(x))
#define log(x) (Libm←log(x))
#define log1p(x) (Libm←log1p(x))
#define lgamma(a) (Libm←lgamma(a))
#define exp(x) (Libm←exp(x))
#define frexp(x,i) (Libm←frexp(x,i))
#define ldexp(x,e) (Libm←ldexp(x,e))
#define expm1(x) (Libm←expm1(x))
#define pow(x,y) (Libm←pow(x,y))
#define fmod(x,y) (Libm←fmod(x,y))
#define modf(d,i) (Libm←modf(d,i))
#define floor(d) (Libm←floor(d))
#define ceil(d) (Libm←ceil(d))
#define rint(x) (Libm←rint(x))
#define cabs(z) (Libm←cabs(z))
#define hypot(x,y) (Libm←hypot(x,y))
#define sqrt(x) (Libm←sqrt(x))
#define cbrt(x) (Libm←cbrt(x))
#define j0(a) (Libm←j0(a))
#define j1(a) (Libm←j1(a))
#define fabs(d) (Libm←fabs(d))
#define abs(i) (Libm←abs(i))

#define HUGE	1.701411733192644270e38


static double zero = 0.e0;

asm ("        exportproc ←jn, Libm");
double
jn(n,x) int n; double x;{
	int i;
	double a, b, temp;
	double xsq, t;
	double j0(), j1();

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


asm ("        exportproc ←yn, Libm");
double
yn(n,x) int n; double x;{
	int i;
	int sign;
	double a, b, temp;
	double y0(), y1();

	if (x <= 0) {

		return(zero/zero);	/* IEEE machines: invalid operation */
	}
	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);
}