/*

 * // 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());

    }

}