1 /* ef_jn.c -- float version of e_jn.c.
26     __int32_t i, hx, ix, sgn;  in jnf()  local
30     /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)  in jnf()
31      * Thus, J(-n,x) = J(n,-x)  in jnf()
39         n = -n;  in jnf()
40         x = -x;  in jnf()
47     sgn = (n & 1) & (hx >> 31); /* even n -- 0, odd n -- sign(x) */  in jnf()
52         /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */  in jnf()
55         for (i = 1; i < n; i++) {  in jnf()
57             b = b * ((float)(i + i) / x) - a; /* avoid underflow */  in jnf()
61         if (ix < 0x30800000) { /* x < 2**-29 */  in jnf()
63      * J(n,x) = 1/n!*(x/2)^n  - ...  in jnf()
70                 for (a = one, i = 2; i <= n; i++) {  in jnf()
71                     a *= (float)i; /* a = n! */  in jnf()
79 		 *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....  in jnf()
80 		 *			2n  - 2(n+1) - 2(n+2)  in jnf()
83 		 *  (for large x)   =  ----  ------   ------   .....  in jnf()
85 		 *			-- - ------ - ------ -  in jnf()
91 		 *	= -----------------------  in jnf()
93 		 *	   w - -----------------  in jnf()
95 		 * 	        w+h - ---------  in jnf()
96 		 *		       w+2h - ...  in jnf()
99 		 * Q(0) = w, Q(1) = w(w+h) - 1,  in jnf()
100 		 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),  in jnf()
113             q1 = w * z - (float)1.0;  in jnf()
118                 tmp = z * q1 - q0;  in jnf()
123             for (t = zero, i = 2 * (n + k); i >= m; i -= 2)  in jnf()
124                 t = one / (i / x - t);  in jnf()
139                 for (i = n - 1, di = (float)(i + i); i > 0; i--) {  in jnf()
142                     b = b / x - a;  in jnf()
144                     di -= two;  in jnf()
147                 for (i = n - 1, di = (float)(i + i); i > 0; i--) {  in jnf()
150                     b = b / x - a;  in jnf()
152                     di -= two;  in jnf()
165         return -b;  in jnf()
173     __int32_t i, hx, ix, ib;  in ynf()  local
194         n = -n;  in ynf()
195         sign = 1 - ((n & 1) << 1);  in ynf()
204     /* quit if b is -inf */  in ynf()
206     for (i = 1; i < n && ib != (__int32_t)0xff800000; i++) {  in ynf()
208         b = ((float)(i + i) / x) * b - a;  in ynf()
215         return -b;  in ynf()