pxfm/gamma/
betaincf.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 9/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1.  Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3.  Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::common::f_fmla;
30use crate::gamma::lnbetaf::lnbetaf_core;
31use crate::{f_exp, f_log, f_log1p, f_pow, f_powf};
32
33/// Regularized incomplete beta
34pub fn f_betainc_regf(a: f32, b: f32, x: f32) -> f32 {
35    let aa = a.to_bits();
36    let ab = b.to_bits();
37    let ax = x.to_bits();
38
39    if aa >= 0xffu32 << 23
40        || aa == 0
41        || ab >= 0xffu32 << 23
42        || ab == 0
43        || ax == 0
44        || ax >= 0x3f800000
45    {
46        if (aa >> 31) != 0 || (ab >> 31) != 0 || (ax >> 31) != 0 {
47            // |a| < 0 or |b| < 0
48            return f32::NAN;
49        }
50        if ax >= 0x3f800000 {
51            // |x| > 1
52            if ax == 0x3f800000 {
53                // |x| == 1
54                return 1.;
55            }
56            return f32::NAN;
57        }
58        if aa.wrapping_shl(1) == 0 {
59            // |a| == 0
60            return 1.0;
61        }
62        if ab.wrapping_shl(1) == 0 {
63            // |b| == 0
64            return 0.;
65        }
66        if ax.wrapping_shl(1) == 0 {
67            // |x| == 0
68            return 0.;
69        }
70        if a.is_infinite() {
71            // |a| == inf
72            return 0.;
73        }
74        if b.is_infinite() {
75            // |b| == inf
76            return 1.;
77        }
78        return a + f32::NAN; // nan
79    }
80
81    if aa == 0x3f800000 {
82        // a == 1
83        return (1. - f_pow(1. - x as f64, b as f64)) as f32;
84    }
85
86    if ab == 0x3f800000 {
87        // b == 1
88        return f_powf(x, a);
89    }
90
91    /*The continued fraction converges nicely for x < (a+1)/(a+b+2)*/
92    /*Use the fact that beta is symmetrical.*/
93    let mut return_inverse = false;
94    let mut dx = x as f64;
95    let mut a = a;
96    let mut b = b;
97    if x > (a + 1.0) / (a + b + 2.0) {
98        std::mem::swap(&mut a, &mut b);
99        dx = 1.0 - dx;
100        return_inverse = true;
101    }
102    /*Find the first part before the continued fraction.*/
103    let da = a as f64;
104    let db = b as f64;
105    let w0 = f_fmla(f_log(dx), da, f_fmla(f_log1p(-dx), db, -lnbetaf_core(a, b)));
106    let front = f_exp(w0) / da;
107
108    /*Use Lentz's algorithm to evaluate the continued fraction.*/
109    let mut f = 1.0;
110    let mut c = 1.0;
111    let mut d = 0.0;
112
113    const TINY: f64 = 1.0e-30;
114    const STOP: f64 = 1.0e-8;
115
116    for i in 0..200 {
117        let m = i / 2;
118        let numerator: f64 = if i == 0 {
119            1.0 /*First numerator is 1.0.*/
120        } else if i % 2 == 0 {
121            let m = m as f64;
122            let c0 = f_fmla(2.0, m, da);
123            (m * (db - m) * dx) / ((c0 - 1.0) * c0) /*Even term.*/
124        } else {
125            let m = m as f64;
126            let c0 = f_fmla(2.0, m, da);
127            -((da + m) * (da + db + m) * dx) / (c0 * (c0 + 1.)) /*Odd term.*/
128        };
129
130        /*Do an iteration of Lentz's algorithm.*/
131        d = f_fmla(numerator, d, 1.0);
132        if d.abs() < TINY {
133            d = TINY;
134        }
135        d = 1.0 / d;
136
137        c = 1.0 + numerator / c;
138        if c.abs() < TINY {
139            c = TINY;
140        }
141
142        let cd = c * d;
143        f *= cd;
144
145        /*Check for stop.*/
146        if (1.0 - cd).abs() < STOP {
147            return if return_inverse {
148                f_fmla(-front, f - 1.0, 1.) as f32
149            } else {
150                (front * (f - 1.0)) as f32
151            };
152        }
153    }
154
155    if return_inverse {
156        f_fmla(-front, f - 1.0, 1.) as f32
157    } else {
158        (front * (f - 1.0)) as f32
159    }
160}
161
162#[cfg(test)]
163mod tests {
164    use super::*;
165
166    #[test]
167    fn test_betaincf() {
168        assert_eq!(f_betainc_regf(0.5, 2.5, 0.5), 0.9244131815783875);
169        assert_eq!(f_betainc_regf(2.5, 1.0, 0.5), 0.1767766952966368811);
170        assert_eq!(f_betainc_regf(0.5, 0.5, 0.5), 0.5);
171        assert_eq!(f_betainc_regf(54221., 23124., 0.64534), 0.0);
172        assert_eq!(f_betainc_regf(5., 1.4324, 0.1312), 8.872578e-5);
173        assert_eq!(f_betainc_regf(7., 42., 0.4324), 0.99999547);
174        assert_eq!(f_betainc_regf(0.5, 0., 1.), 1.);
175        assert_eq!(f_betainc_regf(5., 2., 1.), 1.);
176        assert_eq!(f_betainc_regf(5., 2., 0.), 0.);
177        assert_eq!(f_betainc_regf(5., 2., 0.5), 0.109375);
178        assert!(f_betainc_regf(5., 2., -1.).is_nan());
179        assert!(f_betainc_regf(5., 2., 1.1).is_nan());
180        assert!(f_betainc_regf(5., 2., f32::INFINITY).is_nan());
181        assert!(f_betainc_regf(5., 2., f32::NEG_INFINITY).is_nan());
182        assert!(f_betainc_regf(5., 2., f32::NAN).is_nan());
183        assert!(f_betainc_regf(-5., 2., 0.432).is_nan());
184        assert!(f_betainc_regf(5., -2., 0.432).is_nan());
185        assert!(f_betainc_regf(5., 2., -0.432).is_nan());
186    }
187}
pxfm/gamma/betaincf.rs

pxfm/gamma/
betaincf.rs
