1 /* e_jnf.c -- float version of e_jn.c.
2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
3 */
4
5 /*
6 * ====================================================
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8 *
9 * Developed at SunPro, a Sun Microsystems, Inc. business.
10 * Permission to use, copy, modify, and distribute this
11 * software is freely granted, provided that this notice
12 * is preserved.
13 * ====================================================
14 */
15
16 #include "cdefs-compat.h"
17 //__FBSDID("$FreeBSD: src/lib/msun/src/e_jnf.c,v 1.11 2010/11/13 10:54:10 uqs Exp $");
18
19 #include <openlibm_math.h>
20
21 #include "math_private.h"
22
23 static const float
24 two = 2.0000000000e+00, /* 0x40000000 */
25 one = 1.0000000000e+00; /* 0x3F800000 */
26
27 static const float zero = 0.0000000000e+00;
28
29 OLM_DLLEXPORT float
__ieee754_jnf(int n,float x)30 __ieee754_jnf(int n, float x)
31 {
32 int32_t i,hx,ix, sgn;
33 float a, b, temp, di;
34 float z, w;
35
36 /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
37 * Thus, J(-n,x) = J(n,-x)
38 */
39 GET_FLOAT_WORD(hx,x);
40 ix = 0x7fffffff&hx;
41 /* if J(n,NaN) is NaN */
42 if(ix>0x7f800000) return x+x;
43 if(n<0){
44 n = -n;
45 x = -x;
46 hx ^= 0x80000000;
47 }
48 if(n==0) return(__ieee754_j0f(x));
49 if(n==1) return(__ieee754_j1f(x));
50 sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
51 x = fabsf(x);
52 if(ix==0||ix>=0x7f800000) /* if x is 0 or inf */
53 b = zero;
54 else if((float)n<=x) {
55 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
56 a = __ieee754_j0f(x);
57 b = __ieee754_j1f(x);
58 for(i=1;i<n;i++){
59 temp = b;
60 b = b*((float)(i+i)/x) - a; /* avoid underflow */
61 a = temp;
62 }
63 } else {
64 if(ix<0x30800000) { /* x < 2**-29 */
65 /* x is tiny, return the first Taylor expansion of J(n,x)
66 * J(n,x) = 1/n!*(x/2)^n - ...
67 */
68 if(n>33) /* underflow */
69 b = zero;
70 else {
71 temp = x*(float)0.5; b = temp;
72 for (a=one,i=2;i<=n;i++) {
73 a *= (float)i; /* a = n! */
74 b *= temp; /* b = (x/2)^n */
75 }
76 b = b/a;
77 }
78 } else {
79 /* use backward recurrence */
80 /* x x^2 x^2
81 * J(n,x)/J(n-1,x) = ---- ------ ------ .....
82 * 2n - 2(n+1) - 2(n+2)
83 *
84 * 1 1 1
85 * (for large x) = ---- ------ ------ .....
86 * 2n 2(n+1) 2(n+2)
87 * -- - ------ - ------ -
88 * x x x
89 *
90 * Let w = 2n/x and h=2/x, then the above quotient
91 * is equal to the continued fraction:
92 * 1
93 * = -----------------------
94 * 1
95 * w - -----------------
96 * 1
97 * w+h - ---------
98 * w+2h - ...
99 *
100 * To determine how many terms needed, let
101 * Q(0) = w, Q(1) = w(w+h) - 1,
102 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
103 * When Q(k) > 1e4 good for single
104 * When Q(k) > 1e9 good for double
105 * When Q(k) > 1e17 good for quadruple
106 */
107 /* determine k */
108 float t,v;
109 float q0,q1,h,tmp; int32_t k,m;
110 w = (n+n)/(float)x; h = (float)2.0/(float)x;
111 q0 = w; z = w+h; q1 = w*z - (float)1.0; k=1;
112 while(q1<(float)1.0e9) {
113 k += 1; z += h;
114 tmp = z*q1 - q0;
115 q0 = q1;
116 q1 = tmp;
117 }
118 m = n+n;
119 for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);
120 a = t;
121 b = one;
122 /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
123 * Hence, if n*(log(2n/x)) > ...
124 * single 8.8722839355e+01
125 * double 7.09782712893383973096e+02
126 * long double 1.1356523406294143949491931077970765006170e+04
127 * then recurrent value may overflow and the result is
128 * likely underflow to zero
129 */
130 tmp = n;
131 v = two/x;
132 tmp = tmp*__ieee754_logf(fabsf(v*tmp));
133 if(tmp<(float)8.8721679688e+01) {
134 for(i=n-1,di=(float)(i+i);i>0;i--){
135 temp = b;
136 b *= di;
137 b = b/x - a;
138 a = temp;
139 di -= two;
140 }
141 } else {
142 for(i=n-1,di=(float)(i+i);i>0;i--){
143 temp = b;
144 b *= di;
145 b = b/x - a;
146 a = temp;
147 di -= two;
148 /* scale b to avoid spurious overflow */
149 if(b>(float)1e10) {
150 a /= b;
151 t /= b;
152 b = one;
153 }
154 }
155 }
156 z = __ieee754_j0f(x);
157 w = __ieee754_j1f(x);
158 if (fabsf(z) >= fabsf(w))
159 b = (t*z/b);
160 else
161 b = (t*w/a);
162 }
163 }
164 if(sgn==1) return -b; else return b;
165 }
166
167 OLM_DLLEXPORT float
__ieee754_ynf(int n,float x)168 __ieee754_ynf(int n, float x)
169 {
170 int32_t i,hx,ix,ib;
171 int32_t sign;
172 float a, b, temp;
173
174 GET_FLOAT_WORD(hx,x);
175 ix = 0x7fffffff&hx;
176 /* if Y(n,NaN) is NaN */
177 if(ix>0x7f800000) return x+x;
178 if(ix==0) return -one/zero;
179 if(hx<0) return zero/zero;
180 sign = 1;
181 if(n<0){
182 n = -n;
183 sign = 1 - ((n&1)<<1);
184 }
185 if(n==0) return(__ieee754_y0f(x));
186 if(n==1) return(sign*__ieee754_y1f(x));
187 if(ix==0x7f800000) return zero;
188
189 a = __ieee754_y0f(x);
190 b = __ieee754_y1f(x);
191 /* quit if b is -inf */
192 GET_FLOAT_WORD(ib,b);
193 for(i=1;i<n&&ib!=0xff800000;i++){
194 temp = b;
195 b = ((float)(i+i)/x)*b - a;
196 GET_FLOAT_WORD(ib,b);
197 a = temp;
198 }
199 if(sign>0) return b; else return -b;
200 }
201