pxfm 0.1.28

/*
 * // Copyright (c) Radzivon Bartoshyk 9/2025. All rights reserved.
 * //
 * // Redistribution and use in source and binary forms, with or without modification,
 * // are permitted provided that the following conditions are met:
 * //
 * // 1.  Redistributions of source code must retain the above copyright notice, this
 * // list of conditions and the following disclaimer.
 * //
 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
 * // this list of conditions and the following disclaimer in the documentation
 * // and/or other materials provided with the distribution.
 * //
 * // 3.  Neither the name of the copyright holder nor the names of its
 * // contributors may be used to endorse or promote products derived from
 * // this software without specific prior written permission.
 * //
 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
use crate::exponents::core_expf;
use crate::polyeval::f_polyeval6;

/// Logistic function
///
/// Computes 1/(1+exp(-x))
pub fn f_logisticf(x: f32) -> f32 {
    let x_u = x.to_bits();
    let x_abs = x_u & 0x7fff_ffffu32;
    if x_abs >= 0x42b2_0000u32 || x_abs <= 0x3d00_0000u32 {
        // |x| >= 89 or |x| <= 2^-5
        let exp = ((x_u >> 23) & 0xFF) as i32;
        // |x| < 2^-25
        if exp <= 101i32 {
            return (1. / (2. + -x as f64)) as f32;
        }

        if x_abs <= 0x3400_0000u32 {
            // |x| < f32::EPSILON
            // logistic(x) ~ 0.5 + x / 4 + O(x^3)
            #[cfg(any(
                all(
                    any(target_arch = "x86", target_arch = "x86_64"),
                    target_feature = "fma"
                ),
                target_arch = "aarch64"
            ))]
            {
                use crate::common::f_fmlaf;
                return f_fmlaf(0.25, x, 0.5);
            }
            #[cfg(not(any(
                all(
                    any(target_arch = "x86", target_arch = "x86_64"),
                    target_feature = "fma"
                ),
                target_arch = "aarch64"
            )))]
            {
                use crate::common::f_fmla;
                return f_fmla(0.25, x as f64, 0.5) as f32;
            }
        }

        // When x < log(2^-150) or nan
        if x_u >= 0xc2cf_f1b5u32 {
            // logistic(-Inf) = 0
            if x.is_infinite() {
                return 0.0;
            }
            // logistic(nan) = nan
            if x.is_nan() {
                return x;
            }
            return if x.is_sign_positive() {
                f32::INFINITY
            } else {
                0.
            };
        }
        if x.is_sign_positive() && x_u >= 0x42d0_0000u32 {
            // x >= 104
            if x.is_nan() {
                return f32::NAN;
            }
            return 1.;
        }
        if x.is_sign_negative() && x_u >= 0xc2cb_fe00u32 {
            // x <= -101.99609
            return 0.;
        }

        if x_abs <= 0x3d00_0000u32 {
            // |x| <= 2^-5

            // Polynomial for logistics function generated by Sollya:
            // d = [-2^-5, 2^-5];
            // f_logistic = 1/(1+exp(-x));
            // Q = fpminimax(f_logistic, 5, [|D...|], d);

            let dx = x as f64;
            let p = f_polyeval6(
                dx,
                f64::from_bits(0x3fe0000000000000),
                f64::from_bits(0x3fcffffffffffcf5),
                f64::from_bits(0x3cfbce9c01a4987c),
                f64::from_bits(0xbf955555524c34a4),
                f64::from_bits(0xbda6322b5fac64d4),
                f64::from_bits(0x3f61104f3a0dedc3),
            );
            return p as f32;
        }
    }

    let v_exp = core_expf(-x);
    let logistic = 1. / (1. + v_exp);
    logistic as f32
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_logisticf() {
        assert_eq!(f_logisticf(-89.), 2.227364e-39);
        assert_eq!(f_logisticf(-104.), 0.);
        assert_eq!(f_logisticf(-103.), 0.);
        assert_eq!(f_logisticf(-88.9), 2.461613e-39);
        assert_eq!(f_logisticf(-88.), 6.054601e-39);
        assert_eq!(f_logisticf(-80.), 1.8048513e-35);
        assert_eq!(f_logisticf(-60.), 8.756511e-27);
        assert_eq!(f_logisticf(-40.), 4.248354e-18);
        assert_eq!(f_logisticf(-20.), 2.0611537e-9);
        assert_eq!(f_logisticf(-1.591388e29), 0.);
        assert_eq!(f_logisticf(-3.), 0.047425874);
        assert_eq!(f_logisticf(3.), 0.95257413);
        assert_eq!(f_logisticf(20.), 1.);
        assert_eq!(f_logisticf(55.), 1.);
        assert_eq!(f_logisticf(-104.), 0.);
        assert_eq!(f_logisticf(-90.), 8.19401e-40);
        assert_eq!(f_logisticf(0.00000000000524323), 0.5);
        assert_eq!(f_logisticf(0.00000000524323), 0.5);
        assert_eq!(f_logisticf(0.02), 0.5049998);
        assert_eq!(f_logisticf(-0.02), 0.49500015);
        assert_eq!(f_logisticf(90.), 1.);
        assert_eq!(f_logisticf(105.), 1.);
        assert_eq!(f_logisticf(150.), 1.);
        assert_eq!(f_logisticf(f32::INFINITY), 1.0);
        assert_eq!(f_logisticf(f32::NEG_INFINITY), 0.);
        assert!(f_logisticf(f32::NAN).is_nan());
    }
}