]>
Commit | Line | Data |
---|---|---|
f7eac6eb RM |
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 | |
ba1ffaa1 | 11 | * software is freely granted, provided that this notice |
f7eac6eb RM |
12 | * is preserved. |
13 | * ==================================================== | |
14 | */ | |
15 | ||
354691b7 | 16 | #include <errno.h> |
be254932 | 17 | #include <float.h> |
1ed0291c RH |
18 | #include <math.h> |
19 | #include <math_private.h> | |
f7eac6eb | 20 | |
f7eac6eb | 21 | static const float |
f7eac6eb RM |
22 | two = 2.0000000000e+00, /* 0x40000000 */ |
23 | one = 1.0000000000e+00; /* 0x3F800000 */ | |
24 | ||
f7eac6eb | 25 | static const float zero = 0.0000000000e+00; |
f7eac6eb | 26 | |
0ac5ae23 UD |
27 | float |
28 | __ieee754_jnf(int n, float x) | |
f7eac6eb RM |
29 | { |
30 | int32_t i,hx,ix, sgn; | |
31 | float a, b, temp, di; | |
32 | float z, w; | |
33 | ||
34 | /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x) | |
35 | * Thus, J(-n,x) = J(n,-x) | |
36 | */ | |
37 | GET_FLOAT_WORD(hx,x); | |
38 | ix = 0x7fffffff&hx; | |
39 | /* if J(n,NaN) is NaN */ | |
0ac5ae23 | 40 | if(__builtin_expect(ix>0x7f800000, 0)) return x+x; |
ba1ffaa1 | 41 | if(n<0){ |
f7eac6eb RM |
42 | n = -n; |
43 | x = -x; | |
44 | hx ^= 0x80000000; | |
45 | } | |
46 | if(n==0) return(__ieee754_j0f(x)); | |
47 | if(n==1) return(__ieee754_j1f(x)); | |
48 | sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */ | |
49 | x = fabsf(x); | |
0ac5ae23 | 50 | if(__builtin_expect(ix==0||ix>=0x7f800000, 0)) /* if x is 0 or inf */ |
f7eac6eb | 51 | b = zero; |
ba1ffaa1 | 52 | else if((float)n<=x) { |
f7eac6eb RM |
53 | /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ |
54 | a = __ieee754_j0f(x); | |
55 | b = __ieee754_j1f(x); | |
56 | for(i=1;i<n;i++){ | |
57 | temp = b; | |
b7835e32 | 58 | b = b*((double)(i+i)/x) - a; /* avoid underflow */ |
f7eac6eb RM |
59 | a = temp; |
60 | } | |
61 | } else { | |
62 | if(ix<0x30800000) { /* x < 2**-29 */ | |
ba1ffaa1 | 63 | /* x is tiny, return the first Taylor expansion of J(n,x) |
f7eac6eb RM |
64 | * J(n,x) = 1/n!*(x/2)^n - ... |
65 | */ | |
66 | if(n>33) /* underflow */ | |
67 | b = zero; | |
68 | else { | |
69 | temp = x*(float)0.5; b = temp; | |
70 | for (a=one,i=2;i<=n;i++) { | |
71 | a *= (float)i; /* a = n! */ | |
72 | b *= temp; /* b = (x/2)^n */ | |
73 | } | |
74 | b = b/a; | |
75 | } | |
76 | } else { | |
77 | /* use backward recurrence */ | |
c0ad824e | 78 | /* x x^2 x^2 |
f7eac6eb RM |
79 | * J(n,x)/J(n-1,x) = ---- ------ ------ ..... |
80 | * 2n - 2(n+1) - 2(n+2) | |
81 | * | |
c0ad824e | 82 | * 1 1 1 |
f7eac6eb RM |
83 | * (for large x) = ---- ------ ------ ..... |
84 | * 2n 2(n+1) 2(n+2) | |
ba1ffaa1 | 85 | * -- - ------ - ------ - |
f7eac6eb RM |
86 | * x x x |
87 | * | |
88 | * Let w = 2n/x and h=2/x, then the above quotient | |
89 | * is equal to the continued fraction: | |
90 | * 1 | |
91 | * = ----------------------- | |
92 | * 1 | |
93 | * w - ----------------- | |
94 | * 1 | |
0ac5ae23 | 95 | * w+h - --------- |
f7eac6eb RM |
96 | * w+2h - ... |
97 | * | |
98 | * To determine how many terms needed, let | |
99 | * Q(0) = w, Q(1) = w(w+h) - 1, | |
100 | * Q(k) = (w+k*h)*Q(k-1) - Q(k-2), | |
ba1ffaa1 UD |
101 | * When Q(k) > 1e4 good for single |
102 | * When Q(k) > 1e9 good for double | |
103 | * When Q(k) > 1e17 good for quadruple | |
f7eac6eb RM |
104 | */ |
105 | /* determine k */ | |
106 | float t,v; | |
107 | float q0,q1,h,tmp; int32_t k,m; | |
108 | w = (n+n)/(float)x; h = (float)2.0/(float)x; | |
109 | q0 = w; z = w+h; q1 = w*z - (float)1.0; k=1; | |
110 | while(q1<(float)1.0e9) { | |
111 | k += 1; z += h; | |
112 | tmp = z*q1 - q0; | |
113 | q0 = q1; | |
114 | q1 = tmp; | |
115 | } | |
116 | m = n+n; | |
117 | for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t); | |
118 | a = t; | |
119 | b = one; | |
120 | /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) | |
121 | * Hence, if n*(log(2n/x)) > ... | |
122 | * single 8.8722839355e+01 | |
123 | * double 7.09782712893383973096e+02 | |
124 | * long double 1.1356523406294143949491931077970765006170e+04 | |
ba1ffaa1 | 125 | * then recurrent value may overflow and the result is |
f7eac6eb RM |
126 | * likely underflow to zero |
127 | */ | |
128 | tmp = n; | |
129 | v = two/x; | |
130 | tmp = tmp*__ieee754_logf(fabsf(v*tmp)); | |
131 | if(tmp<(float)8.8721679688e+01) { | |
c0ad824e | 132 | for(i=n-1,di=(float)(i+i);i>0;i--){ |
0ac5ae23 | 133 | temp = b; |
f7eac6eb RM |
134 | b *= di; |
135 | b = b/x - a; | |
0ac5ae23 | 136 | a = temp; |
f7eac6eb | 137 | di -= two; |
c0ad824e | 138 | } |
f7eac6eb | 139 | } else { |
c0ad824e | 140 | for(i=n-1,di=(float)(i+i);i>0;i--){ |
0ac5ae23 | 141 | temp = b; |
f7eac6eb RM |
142 | b *= di; |
143 | b = b/x - a; | |
0ac5ae23 | 144 | a = temp; |
f7eac6eb RM |
145 | di -= two; |
146 | /* scale b to avoid spurious overflow */ | |
147 | if(b>(float)1e10) { | |
148 | a /= b; | |
149 | t /= b; | |
150 | b = one; | |
151 | } | |
c0ad824e | 152 | } |
f7eac6eb | 153 | } |
1248c1c4 PB |
154 | /* j0() and j1() suffer enormous loss of precision at and |
155 | * near zero; however, we know that their zero points never | |
156 | * coincide, so just choose the one further away from zero. | |
157 | */ | |
158 | z = __ieee754_j0f (x); | |
159 | w = __ieee754_j1f (x); | |
160 | if (fabsf (z) >= fabsf (w)) | |
161 | b = (t * z / b); | |
162 | else | |
163 | b = (t * w / a); | |
f7eac6eb RM |
164 | } |
165 | } | |
166 | if(sgn==1) return -b; else return b; | |
167 | } | |
0ac5ae23 | 168 | strong_alias (__ieee754_jnf, __jnf_finite) |
f7eac6eb | 169 | |
0ac5ae23 UD |
170 | float |
171 | __ieee754_ynf(int n, float x) | |
f7eac6eb | 172 | { |
be254932 JM |
173 | float ret; |
174 | { | |
ba1ffaa1 UD |
175 | int32_t i,hx,ix; |
176 | u_int32_t ib; | |
f7eac6eb RM |
177 | int32_t sign; |
178 | float a, b, temp; | |
179 | ||
180 | GET_FLOAT_WORD(hx,x); | |
181 | ix = 0x7fffffff&hx; | |
182 | /* if Y(n,NaN) is NaN */ | |
0ac5ae23 UD |
183 | if(__builtin_expect(ix>0x7f800000, 0)) return x+x; |
184 | if(__builtin_expect(ix==0, 0)) | |
185 | return -HUGE_VALF+x; /* -inf and overflow exception. */ | |
186 | if(__builtin_expect(hx<0, 0)) return zero/(zero*x); | |
f7eac6eb RM |
187 | sign = 1; |
188 | if(n<0){ | |
189 | n = -n; | |
190 | sign = 1 - ((n&1)<<1); | |
191 | } | |
192 | if(n==0) return(__ieee754_y0f(x)); | |
be254932 JM |
193 | SET_RESTORE_ROUNDF (FE_TONEAREST); |
194 | if(n==1) { | |
195 | ret = sign*__ieee754_y1f(x); | |
196 | goto out; | |
197 | } | |
0ac5ae23 | 198 | if(__builtin_expect(ix==0x7f800000, 0)) return zero; |
f7eac6eb RM |
199 | |
200 | a = __ieee754_y0f(x); | |
201 | b = __ieee754_y1f(x); | |
202 | /* quit if b is -inf */ | |
203 | GET_FLOAT_WORD(ib,b); | |
ba1ffaa1 | 204 | for(i=1;i<n&&ib!=0xff800000;i++){ |
f7eac6eb | 205 | temp = b; |
b7835e32 | 206 | b = ((double)(i+i)/x)*b - a; |
f7eac6eb RM |
207 | GET_FLOAT_WORD(ib,b); |
208 | a = temp; | |
209 | } | |
354691b7 | 210 | /* If B is +-Inf, set up errno accordingly. */ |
d81f90cc | 211 | if (! isfinite (b)) |
354691b7 | 212 | __set_errno (ERANGE); |
be254932 JM |
213 | if(sign>0) ret = b; else ret = -b; |
214 | } | |
215 | out: | |
d81f90cc | 216 | if (isinf (ret)) |
be254932 JM |
217 | ret = __copysignf (FLT_MAX, ret) * FLT_MAX; |
218 | return ret; | |
f7eac6eb | 219 | } |
0ac5ae23 | 220 | strong_alias (__ieee754_ynf, __ynf_finite) |