lambert_w 2.0.3

// Copyright 2024-2026 Johanna Sörngård
// SPDX-License-Identifier: MIT OR Apache-2.0

//! This module contains an implementation of the approximation of the secondary
//! branch of the Lambert W function
//! with 24 bits of accuracy from Fukushima's paper.
//! It is based on the Fortran implementation of the name "swm1" by Fukushima.

// The coefficients in these rational minimax functions all have excessive precision.
// By keeping the full precision in the source code we can ensure that there is no confusion
// when comparing with the paper.
#![allow(clippy::excessive_precision)]

use crate::{
    generic_math::{ln, rational_function, sqrt},
    INV_SQRT_E, NEG_INV_E,
};

/// The secondary branch of the Lambert W function computed to 24 bits of accuracy on 64-bit floats with Fukushima's method[^1].
///
/// # Examples
///
/// #### Basic usage
///
/// ```
/// use lambert_w::sp_lambert_wm1;
/// use approx::assert_abs_diff_eq;
///
/// let mln4 = sp_lambert_wm1(-f64::ln(2.0) / 2.0);
///
/// assert_abs_diff_eq!(mln4, -f64::ln(4.0), epsilon = f64::from(f32::EPSILON));
/// ```
///
/// #### Special cases
///
/// For inputs of -1/e and 0 the function returns exactly -1 and [`NEG_INFINITY`](f64::NEG_INFINITY) respectively:
///
/// ```
/// use lambert_w::{sp_lambert_wm1, NEG_INV_E};
///
/// assert_eq!(sp_lambert_wm1(NEG_INV_E), -1.0);
/// assert_eq!(sp_lambert_wm1(0.0), f64::NEG_INFINITY);
/// ```
///
/// Inputs smaller than -1/e or larger than 0, as well as inputs of [`NAN`](f64::NAN), result in [`NAN`](f64::NAN):
///
/// ```
/// use lambert_w::{sp_lambert_wm1, NEG_INV_E};
///
/// assert!(sp_lambert_wm1(NEG_INV_E.next_down()).is_nan());
/// assert!(sp_lambert_wm1(f64::MIN_POSITIVE).is_nan());
/// assert!(sp_lambert_wm1(f64::NAN).is_nan());
/// ```
///
/// # Reference
///
/// [^1]: Toshio Fukushima. **Precise and fast computation of Lambert W function by piecewise minimax rational function approximation with variable transformation**. DOI: [10.13140/RG.2.2.30264.37128](https://doi.org/10.13140/RG.2.2.30264.37128). November 2020.
#[must_use = "this is a pure function that only returns a value and has no side effects"]
pub fn sp_lambert_wm1(z: f64) -> f64 {
    // The critical arguments used in the if statements are related to the numbers in table 4 of the paper, column one.
    // The coefficients in the rational functions are related to the tables 8 and 9 in the paper.
    // The actual numbers are taken from Fukushima's Fortran implementation, where they have higher precision.

    if z < NEG_INV_E {
        f64::NAN
    } else if z == NEG_INV_E {
        -1.0
    } else if z <= -0.207_293_777_640_384_138_99 {
        // W >= -2.483, Y_-1

        rational_function(
            -z / (INV_SQRT_E + sqrt(z - NEG_INV_E)),
            [
                -6.383_722_822_801_905_640_2,
                -74.968_653_259_594_050_421,
                -19.714_821_552_432_482_711,
                70.677_326_667_809_236_886,
            ],
            [
                1.0,
                24.295_836_951_878_692_244,
                64.112_460_611_386_039_710,
                17.994_497_369_039_313_585,
            ],
        )
    } else if z <= -0.071_507_705_083_841_949_275 {
        // W >= -4.032, Y_-2

        rational_function(
            -z / (INV_SQRT_E + sqrt(z - NEG_INV_E)),
            [
                -7.723_328_481_229_978_044_5,
                -352.484_690_970_429_148_66,
                -1_242.008_890_368_567_715_4,
                1_171.647_596_062_049_903_3,
            ],
            [
                1.0,
                77.681_242_588_997_413_555,
                648.564_312_140_752_460_23,
                566.701_549_764_361_648_18,
            ],
        )
    } else if z <= -0.020_704_412_621_717_479_786 {
        // W >= -5.600, Y_-3

        rational_function(
            -z / (INV_SQRT_E + sqrt(z - NEG_INV_E)),
            [
                -9.137_773_141_758_154_626_2,
                -1_644.724_479_150_888_940_7,
                -28_105.096_098_779_681_470,
                3_896.079_810_390_921_391_5,
            ],
            [
                1.0,
                272.375_261_351_239_732_61,
                7_929.224_261_291_349_736_3,
                23_980.122_860_821_312_510,
            ],
        )
    } else if z <= -0.005_480_012_945_209_443_573_3 {
        // W >= -7.178, Y_-4

        rational_function(
            -z / (INV_SQRT_E + sqrt(z - NEG_INV_E)),
            [
                -10.603_388_239_566_372_960,
                -7_733.348_521_498_647_619_3,
                -575_482.407_079_644_221_74,
                -2.154_552_604_188_978_313_2e6,
            ],
            [
                1.0,
                1_021.793_856_606_681_728_7,
                111_300.229_154_865_209_18,
                1.261_425_640_008_844_218_6e6,
            ],
        )
    } else if z <= -0.001_367_466_989_250_804_268_7 {
        // W >= -8.766, Y_-5

        rational_function(
            -z / (INV_SQRT_E + sqrt(z - NEG_INV_E)),
            [
                -12.108_699_273_343_437_997,
                -36_896.535_108_166_377_289,
                -1.183_112_672_010_605_357_4e7,
                -2.756_583_081_394_092_428_6e8,
            ],
            [
                1.0,
                4_044.975_306_488_070_069_7,
                1.741_827_761_903_000_471_3e6,
                7.843_690_738_080_689_728_0e7,
            ],
        )
    } else if z <= -0.000_326_142_267_310_725_659_15 {
        // W >= -10.367, Y_-6

        rational_function(
            -z / (INV_SQRT_E + sqrt(z - NEG_INV_E)),
            [
                -13.646_761_936_746_191_398,
                -179_086.115_857_151_485_53,
                -2.508_463_493_521_464_102_0e8,
                -2.934_370_049_483_371_501_9e10,
            ],
            [
                1.0,
                16_743.826_607_737_142_476,
                2.980_965_094_601_174_304_7e7,
                5.573_951_481_695_800_120_2e9,
            ],
        )
    } else if z <= -0.000_074_906_612_036_101_435_088 {
        // W >= -11.983, Y_-7

        rational_function(
            -z / (INV_SQRT_E + sqrt(z - NEG_INV_E)),
            [
                -15.212_958_142_001_646_273,
                -884_954.687_981_689_538_07,
                -5.529_815_437_863_347_537_9e9,
                -3.093_418_743_531_466_820_3e12,
            ],
            [
                1.0,
                72_009.255_525_210_120_244,
                5.505_900_767_187_580_347_0e8,
                4.432_489_486_700_034_049_5e11,
            ],
        )
    } else if z <= -1.096_244_452_641_099_391_9e-19 {
        // W >= -47.518, V_-8

        rational_function(
            ln(-z),
            [
                -0.032_401_163_177_791_086_145,
                2.028_194_214_474_250_266_7,
                -0.527_524_312_425_927_139_92,
                0.017_340_294_772_717_582_230,
            ],
            [
                1.0,
                -0.450_042_744_438_917_372_96,
                0.017_154_705_753_566_296_216,
                -5.243_819_620_271_835_822_6e-7,
            ],
        )
    } else if z <= -2.509_609_929_994_590_145_5e-136 {
        // W >= -317.993, V_-9

        rational_function(
            ln(-z),
            [
                -1.441_124_659_581_209_802_9,
                1.281_926_963_998_047_809_5,
                -0.074_979_356_113_812_324_952,
                0.000_476_363_091_620_691_506_55,
            ],
            [
                1.0,
                -0.072_000_873_723_868_650_424,
                0.000_475_489_329_895_970_299_52,
                -4.171_497_924_754_684_061_9e-10,
            ],
        )
    } else if z < 0.0 {
        // V_-10

        rational_function(
            ln(-z),
            [
                -3.310_876_091_171_044_755_4,
                1.050_067_880_993_517_477_5,
                -0.008_236_749_582_134_320_494_0,
                5.528_956_159_491_018_844_1e-6,
            ],
            [
                1.0,
                -0.008_189_272_743_331_551_186_7,
                5.528_007_600_971_195_403_2e-6,
                -3.922_277_308_457_406_059_4e-14,
            ],
        )
    } else if z == 0.0 {
        f64::NEG_INFINITY
    } else {
        f64::NAN
    }
}