1 /* e_jnf.c -- float version of e_jn.c.
23 static const volatile float vone = 1, vzero = 0;
32 jnf(int n, float x)  in jnf()  argument
38     /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)  in jnf()
39      * Thus, J(-n,x) = J(n,-x)  in jnf()
43     /* if J(n,NaN) is NaN */  in jnf()
45 	if(n<0){  in jnf()
46 		n = -n;  in jnf()
47 		x = -x;  in jnf()
50 	if(n==0) return(j0f(x));  in jnf()
51 	if(n==1) return(j1f(x));  in jnf()
52 	sgn = (n&1)&(hx>>31);	/* even n -- 0, odd n -- sign(x) */  in jnf()
56 	else if((float)n<=x) {  in jnf()
57 		/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */  in jnf()
60 	    for(i=1;i<n;i++){  in jnf()
62 		b = b*((float)(i+i)/x) - a; /* avoid underflow */  in jnf()
66 	    if(ix<0x30800000) {	/* x < 2**-29 */  in jnf()
67     /* x is tiny, return the first Taylor expansion of J(n,x)  in jnf()
68      * J(n,x) = 1/n!*(x/2)^n  - ...  in jnf()
70 		if(n>33)	/* underflow */  in jnf()
74 		    for (a=one,i=2;i<=n;i++) {  in jnf()
75 			a *= (float)i;		/* a = n! */  in jnf()
76 			b *= temp;		/* b = (x/2)^n */  in jnf()
83 		 *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....  in jnf()
84 		 *			2n  - 2(n+1) - 2(n+2)  in jnf()
86 		 * 			1      1        1  in jnf()
87 		 *  (for large x)   =  ----  ------   ------   .....  in jnf()
88 		 *			2n   2(n+1)   2(n+2)  in jnf()
89 		 *			-- - ------ - ------ -  in jnf()
92 		 * Let w = 2n/x and h=2/x, then the above quotient  in jnf()
94 		 *		    1  in jnf()
95 		 *	= -----------------------  in jnf()
96 		 *		       1  in jnf()
97 		 *	   w - -----------------  in jnf()
98 		 *			  1  in jnf()
99 		 * 	        w+h - ---------  in jnf()
100 		 *		       w+2h - ...  in jnf()
103 		 * Q(0) = w, Q(1) = w(w+h) - 1,  in jnf()
104 		 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),  in jnf()
105 		 * When Q(k) > 1e4	good for single  in jnf()
106 		 * When Q(k) > 1e9	good for double  in jnf()
107 		 * When Q(k) > 1e17	good for quadruple  in jnf()
112 		w  = (n+n)/(float)x; h = (float)2.0/(float)x;  in jnf()
113 		q0 = w;  z = w+h; q1 = w*z - (float)1.0; k=1;  in jnf()
115 			k += 1; z += h;  in jnf()
116 			tmp = z*q1 - q0;  in jnf()
120 		m = n+n;  in jnf()
121 		for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);  in jnf()
124 		/*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)  in jnf()
125 		 *  Hence, if n*(log(2n/x)) > ...  in jnf()
132 		tmp = n;  in jnf()
136 	    	    for(i=n-1,di=(float)(i+i);i>0;i--){  in jnf()
139 			b  = b/x - a;  in jnf()
141 			di -= two;  in jnf()
144 	    	    for(i=n-1,di=(float)(i+i);i>0;i--){  in jnf()
147 			b  = b/x - a;  in jnf()
149 			di -= two;  in jnf()
151 			if(b>(float)1e10) {  in jnf()
166 	if(sgn==1) return -b; else return b;  in jnf()
170 ynf(int n, float x)  in ynf()  argument
179 	if(ix==0) return -one/vzero;  in ynf()
181 	sign = 1;  in ynf()
182 	if(n<0){  in ynf()
183 		n = -n;  in ynf()
184 		sign = 1 - ((n&1)<<1);  in ynf()
186 	if(n==0) return(y0f(x));  in ynf()
187 	if(n==1) return(sign*y1f(x));  in ynf()
192 	/* quit if b is -inf */  in ynf()
194 	for(i=1;i<n&&ib!=0xff800000;i++){  in ynf()
196 	    b = ((float)(i+i)/x)*b - a;  in ynf()
200 	if(sign>0) return b; else return -b;  in ynf()