]>
git.proxmox.com Git - wasi-libc.git/blob - libc-top-half/musl/src/math/jnf.c
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_jnf.c */
3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
6 * ====================================================
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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
13 * ====================================================
19 float jnf(int n
, float x
)
25 GET_FLOAT_WORD(ix
, x
);
28 if (ix
> 0x7f800000) /* nan */
31 /* J(-n,x) = J(n,-x), use |n|-1 to avoid overflow in -n */
43 sign
&= n
; /* even n: 0, odd n: signbit(x) */
45 if (ix
== 0 || ix
== 0x7f800000) /* if x is 0 or inf */
48 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
58 if (ix
< 0x35800000) { /* x < 2**-20 */
59 /* x is tiny, return the first Taylor expansion of J(n,x)
60 * J(n,x) = 1/n!*(x/2)^n - ...
62 if (nm1
> 8) /* underflow */
67 for (i
=2; i
<=nm1
+1; i
++) {
68 a
*= (float)i
; /* a = n! */
69 b
*= temp
; /* b = (x/2)^n */
73 /* use backward recurrence */
75 * J(n,x)/J(n-1,x) = ---- ------ ------ .....
76 * 2n - 2(n+1) - 2(n+2)
79 * (for large x) = ---- ------ ------ .....
81 * -- - ------ - ------ -
84 * Let w = 2n/x and h=2/x, then the above quotient
85 * is equal to the continued fraction:
87 * = -----------------------
89 * w - -----------------
94 * To determine how many terms needed, let
95 * Q(0) = w, Q(1) = w(w+h) - 1,
96 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
97 * When Q(k) > 1e4 good for single
98 * When Q(k) > 1e9 good for double
99 * When Q(k) > 1e17 good for quadruple
102 float t
,q0
,q1
,w
,h
,z
,tmp
,nf
;
112 while (q1
< 1.0e4f
) {
119 for (t
=0.0f
, i
=k
; i
>=0; i
--)
120 t
= 1.0f
/(2*(i
+nf
)/x
-t
);
123 /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
124 * Hence, if n*(log(2n/x)) > ...
125 * single 8.8722839355e+01
126 * double 7.09782712893383973096e+02
127 * long double 1.1356523406294143949491931077970765006170e+04
128 * then recurrent value may overflow and the result is
129 * likely underflow to zero
131 tmp
= nf
*logf(fabsf(w
));
132 if (tmp
< 88.721679688f
) {
133 for (i
=nm1
; i
>0; i
--) {
139 for (i
=nm1
; i
>0; i
--){
143 /* scale b to avoid spurious overflow */
153 if (fabsf(z
) >= fabsf(w
))
159 return sign
? -b
: b
;
162 float ynf(int n
, float x
)
168 GET_FLOAT_WORD(ix
, x
);
171 if (ix
> 0x7f800000) /* nan */
173 if (sign
&& ix
!= 0) /* x < 0 */
175 if (ix
== 0x7f800000)
188 return sign
? -y1f(x
) : y1f(x
);
192 /* quit if b is -inf */
193 GET_FLOAT_WORD(ib
,b
);
194 for (i
= 0; i
< nm1
&& ib
!= 0xff800000; ) {
197 b
= (2.0f
*i
/x
)*b
- a
;
198 GET_FLOAT_WORD(ib
, b
);
201 return sign
? -b
: b
;