1 /* origin: FreeBSD /usr/src/lib/msun/src/e_lgammaf_r.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 * ====================================================
16 use super::{floorf, k_cosf, k_sinf, logf}
;
18 const PI
: f32 = 3.1415927410e+00; /* 0x40490fdb */
19 const A0
: f32 = 7.7215664089e-02; /* 0x3d9e233f */
20 const A1
: f32 = 3.2246702909e-01; /* 0x3ea51a66 */
21 const A2
: f32 = 6.7352302372e-02; /* 0x3d89f001 */
22 const A3
: f32 = 2.0580807701e-02; /* 0x3ca89915 */
23 const A4
: f32 = 7.3855509982e-03; /* 0x3bf2027e */
24 const A5
: f32 = 2.8905137442e-03; /* 0x3b3d6ec6 */
25 const A6
: f32 = 1.1927076848e-03; /* 0x3a9c54a1 */
26 const A7
: f32 = 5.1006977446e-04; /* 0x3a05b634 */
27 const A8
: f32 = 2.2086278477e-04; /* 0x39679767 */
28 const A9
: f32 = 1.0801156895e-04; /* 0x38e28445 */
29 const A10
: f32 = 2.5214456400e-05; /* 0x37d383a2 */
30 const A11
: f32 = 4.4864096708e-05; /* 0x383c2c75 */
31 const TC
: f32 = 1.4616321325e+00; /* 0x3fbb16c3 */
32 const TF
: f32 = -1.2148628384e-01; /* 0xbdf8cdcd */
33 /* TT = -(tail of TF) */
34 const TT
: f32 = 6.6971006518e-09; /* 0x31e61c52 */
35 const T0
: f32 = 4.8383611441e-01; /* 0x3ef7b95e */
36 const T1
: f32 = -1.4758771658e-01; /* 0xbe17213c */
37 const T2
: f32 = 6.4624942839e-02; /* 0x3d845a15 */
38 const T3
: f32 = -3.2788541168e-02; /* 0xbd064d47 */
39 const T4
: f32 = 1.7970675603e-02; /* 0x3c93373d */
40 const T5
: f32 = -1.0314224288e-02; /* 0xbc28fcfe */
41 const T6
: f32 = 6.1005386524e-03; /* 0x3bc7e707 */
42 const T7
: f32 = -3.6845202558e-03; /* 0xbb7177fe */
43 const T8
: f32 = 2.2596477065e-03; /* 0x3b141699 */
44 const T9
: f32 = -1.4034647029e-03; /* 0xbab7f476 */
45 const T10
: f32 = 8.8108185446e-04; /* 0x3a66f867 */
46 const T11
: f32 = -5.3859531181e-04; /* 0xba0d3085 */
47 const T12
: f32 = 3.1563205994e-04; /* 0x39a57b6b */
48 const T13
: f32 = -3.1275415677e-04; /* 0xb9a3f927 */
49 const T14
: f32 = 3.3552918467e-04; /* 0x39afe9f7 */
50 const U0
: f32 = -7.7215664089e-02; /* 0xbd9e233f */
51 const U1
: f32 = 6.3282704353e-01; /* 0x3f2200f4 */
52 const U2
: f32 = 1.4549225569e+00; /* 0x3fba3ae7 */
53 const U3
: f32 = 9.7771751881e-01; /* 0x3f7a4bb2 */
54 const U4
: f32 = 2.2896373272e-01; /* 0x3e6a7578 */
55 const U5
: f32 = 1.3381091878e-02; /* 0x3c5b3c5e */
56 const V1
: f32 = 2.4559779167e+00; /* 0x401d2ebe */
57 const V2
: f32 = 2.1284897327e+00; /* 0x4008392d */
58 const V3
: f32 = 7.6928514242e-01; /* 0x3f44efdf */
59 const V4
: f32 = 1.0422264785e-01; /* 0x3dd572af */
60 const V5
: f32 = 3.2170924824e-03; /* 0x3b52d5db */
61 const S0
: f32 = -7.7215664089e-02; /* 0xbd9e233f */
62 const S1
: f32 = 2.1498242021e-01; /* 0x3e5c245a */
63 const S2
: f32 = 3.2577878237e-01; /* 0x3ea6cc7a */
64 const S3
: f32 = 1.4635047317e-01; /* 0x3e15dce6 */
65 const S4
: f32 = 2.6642270386e-02; /* 0x3cda40e4 */
66 const S5
: f32 = 1.8402845599e-03; /* 0x3af135b4 */
67 const S6
: f32 = 3.1947532989e-05; /* 0x3805ff67 */
68 const R1
: f32 = 1.3920053244e+00; /* 0x3fb22d3b */
69 const R2
: f32 = 7.2193557024e-01; /* 0x3f38d0c5 */
70 const R3
: f32 = 1.7193385959e-01; /* 0x3e300f6e */
71 const R4
: f32 = 1.8645919859e-02; /* 0x3c98bf54 */
72 const R5
: f32 = 7.7794247773e-04; /* 0x3a4beed6 */
73 const R6
: f32 = 7.3266842264e-06; /* 0x36f5d7bd */
74 const W0
: f32 = 4.1893854737e-01; /* 0x3ed67f1d */
75 const W1
: f32 = 8.3333335817e-02; /* 0x3daaaaab */
76 const W2
: f32 = -2.7777778450e-03; /* 0xbb360b61 */
77 const W3
: f32 = 7.9365057172e-04; /* 0x3a500cfd */
78 const W4
: f32 = -5.9518753551e-04; /* 0xba1c065c */
79 const W5
: f32 = 8.3633989561e-04; /* 0x3a5b3dd2 */
80 const W6
: f32 = -1.6309292987e-03; /* 0xbad5c4e8 */
82 /* sin(PI*x) assuming x > 2^-100, if sin(PI*x)==0 the sign is arbitrary */
83 fn sin_pi(mut x
: f32) -> f32 {
87 /* spurious inexact if odd int */
88 x
= 2.0 * (x
* 0.5 - floorf(x
* 0.5)); /* x mod 2.0 */
90 n
= (x
* 4.0) as isize;
92 y
= (x
as f64) - (n
as f64) * 0.5;
93 y
*= 3.14159265358979323846;
102 pub fn lgammaf_r(mut x
: f32) -> (f32, i32) {
118 let mut signgam
: i32;
120 /* purge off +-inf, NaN, +-0, tiny and negative arguments */
122 sign
= (u
>> 31) != 0;
124 if ix
>= 0x7f800000 {
125 return (x
* x
, signgam
);
128 /* |x| < 2**-21, return -log(|x|) */
133 return (-logf(x
), signgam
);
140 return (1.0 / (x
- x
), signgam
);
147 nadj
= logf(PI
/ (t
* x
));
152 /* purge off 1 and 2 */
153 if ix
== 0x3f800000 || ix
== 0x40000000 {
157 else if ix
< 0x40000000 {
158 if ix
<= 0x3f666666 {
159 /* lgamma(x) = lgamma(x+1)-log(x) */
161 if ix
>= 0x3f3b4a20 {
164 } else if ix
>= 0x3e6d3308 {
173 if ix
>= 0x3fdda618 {
177 } else if ix
>= 0x3F9da620 {
189 p1
= A0
+ z
* (A2
+ z
* (A4
+ z
* (A6
+ z
* (A8
+ z
* A10
))));
190 p2
= z
* (A1
+ z
* (A3
+ z
* (A5
+ z
* (A7
+ z
* (A9
+ z
* A11
)))));
197 p1
= T0
+ w
* (T3
+ w
* (T6
+ w
* (T9
+ w
* T12
))); /* parallel comp */
198 p2
= T1
+ w
* (T4
+ w
* (T7
+ w
* (T10
+ w
* T13
)));
199 p3
= T2
+ w
* (T5
+ w
* (T8
+ w
* (T11
+ w
* T14
)));
200 p
= z
* p1
- (TT
- w
* (p2
+ y
* p3
));
204 p1
= y
* (U0
+ y
* (U1
+ y
* (U2
+ y
* (U3
+ y
* (U4
+ y
* U5
)))));
205 p2
= 1.0 + y
* (V1
+ y
* (V2
+ y
* (V3
+ y
* (V4
+ y
* V5
))));
206 r
+= -0.5 * y
+ p1
/ p2
;
208 #[cfg(debug_assertions)]
210 #[cfg(not(debug_assertions))]
213 } else if ix
< 0x41000000 {
217 p
= y
* (S0
+ y
* (S1
+ y
* (S2
+ y
* (S3
+ y
* (S4
+ y
* (S5
+ y
* S6
))))));
218 q
= 1.0 + y
* (R1
+ y
* (R2
+ y
* (R3
+ y
* (R4
+ y
* (R5
+ y
* R6
)))));
220 z
= 1.0; /* lgamma(1+s) = log(s) + lgamma(s) */
221 // TODO: In C, this was implemented using switch jumps with fallthrough.
222 // Does this implementation have performance problems?
239 } else if ix
< 0x5c800000 {
240 /* 8.0 <= x < 2**58 */
244 w
= W0
+ z
* (W1
+ y
* (W2
+ y
* (W3
+ y
* (W4
+ y
* (W5
+ y
* W6
)))));
245 r
= (x
- 0.5) * (t
- 1.0) + w
;
247 /* 2**58 <= x <= inf */
248 r
= x
* (logf(x
) - 1.0);