lambert_w 2.0.1

// 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 50 bits of accuracy from Fukushima's paper.
//! It is based on the Fortran implementation of the name "dwm1c" 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 50 bits of accuracy on 64-bit floats with Fukushima's method[^1].
///
/// # Examples
///
/// #### Basic usage
///
/// ```
/// use lambert_w::lambert_wm1;
/// use approx::assert_abs_diff_eq;
///
/// let mln4 = lambert_wm1(-f64::ln(2.0) / 2.0);
///
/// assert_abs_diff_eq!(mln4, -f64::ln(4.0));
/// ```
///
/// #### Special cases
///
/// For inputs of -1/e and 0 the function returns exactly -1 and [`NEG_INFINITY`](f64::NEG_INFINITY) respectively:
///
/// ```
/// use lambert_w::{lambert_wm1, NEG_INV_E};
///
/// assert_eq!(lambert_wm1(NEG_INV_E), -1.0);
/// assert_eq!(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::{lambert_wm1, NEG_INV_E};
///
/// assert!(lambert_wm1(NEG_INV_E.next_down()).is_nan());
/// assert!(lambert_wm1(f64::MIN_POSITIVE).is_nan());
/// assert!(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 lambert_wm1(z: f64) -> f64 {
    // The critical arguments used in the if statements are the numbers in table 4 of the paper, column two, with 1/e added, as well as equation 20.
    // The coefficients in the rational functions are from tables 15 through 18 in the paper.

    if z < NEG_INV_E {
        f64::NAN
    } else if z == NEG_INV_E {
        -1.0
    } else if z <= -0.354_291_330_944_216_4 {
        // W >= -1.3, X_-1

        rational_function(
            sqrt(z - NEG_INV_E),
            [
                -1.000_000_000_000_000_111_0,
                4.296_301_617_877_712_700_9,
                -4.099_140_792_400_745_761_2,
                -6.844_284_220_083_330_972_4,
                17.084_773_793_345_271_001,
                -13.015_133_123_886_661_124,
                3.930_360_862_953_985_104_9,
                -0.346_367_465_122_474_573_19,
            ],
            [
                1.0,
                -6.627_945_599_474_762_405_9,
                17.740_962_374_121_397_994,
                -24.446_872_319_343_475_890,
                18.249_006_287_190_617_068,
                -7.058_075_875_662_479_055_0,
                1.197_878_676_279_400_354_5,
                -0.053_875_778_140_352_599_789,
            ],
        )
    } else if z <= -0.188_726_882_822_894_340_49 {
        // W >= -2.637, Y_-1

        rational_function(
            -z / (INV_SQRT_E + sqrt(z - NEG_INV_E)),
            [
                -8.225_315_526_444_684_485_4,
                -813.207_067_320_014_871_78,
                -15_270.113_237_678_509_000,
                -79_971.585_089_674_149_237,
                -103_667.542_158_083_765_11,
                42_284.755_505_061_257_427,
                74_953.525_397_605_484_884,
                10_554.369_146_366_736_811,
            ],
            [
                1.0,
                146.363_151_616_695_676_59,
                3_912.476_137_253_924_071_2,
                31_912.693_749_754_847_460,
                92_441.293_717_108_619_527,
                94_918.733_120_470_346_165,
                29_531.165_406_571_745_340,
                1_641.680_896_033_037_098_7,
            ],
        )
    } else if z <= -0.060_497_597_226_958_343_647 {
        // W >= -4.253, Y_-2

        rational_function(
            -z / (INV_SQRT_E + sqrt(z - NEG_INV_E)),
            [
                -9.618_412_744_335_402_429_5,
                -3_557.856_904_301_800_412_1,
                -254_015.593_112_843_810_43,
                -5.392_389_363_067_063_939_1e6,
                -3.663_825_741_753_689_679_8e7,
                -6.148_431_948_622_696_621_3e7,
                3.042_169_037_744_613_445_1e7,
                3.972_813_905_487_932_045_2e7,
            ],
            [
                1.0,
                507.405_256_285_233_008_01,
                46_852.747_159_777_876_192,
                1.316_830_464_009_143_629_7e6,
                1.311_169_069_371_241_524_2e7,
                4.614_211_644_525_801_519_5e7,
                4.898_226_895_620_883_087_6e7,
                9.195_910_098_798_385_512_2e6,
            ],
        )
    } else if z <= -0.017_105_334_740_676_008_194 {
        // W >= -5.832, Y_-3

        rational_function(
            -z / (INV_SQRT_E + sqrt(z - NEG_INV_E)),
            [
                -11.038_489_462_297_466_388,
                -15_575.812_882_656_619_195,
                -4.249_294_730_489_777_343_3e6,
                -3.517_024_593_880_342_376_8e8,
                -9.865_916_303_661_136_464_0e9,
                -8.619_537_230_330_500_390_8e10,
                -1.328_633_557_402_761_600_0e11,
                1.598_954_643_442_066_046_2e11,
            ],
            [
                1.0,
                1_837.077_069_301_716_681_8,
                612_840.975_855_950_927_61,
                6.214_918_139_846_548_303_7e7,
                2.230_401_131_444_308_396_9e9,
                2.825_423_248_527_369_802_1e10,
                1.077_086_663_954_315_616_5e11,
                7.196_469_887_604_913_199_2e10,
            ],
        )
    } else if z <= -0.004_595_496_212_794_370_643_3 {
        // W >= -7.382, Y_-4

        rational_function(
            -z / (INV_SQRT_E + sqrt(z - NEG_INV_E)),
            [
                -12.474_405_916_395_746_052,
                -68_180.335_575_543_773_385,
                -7.184_659_984_562_009_327_8e7,
                -2.314_268_822_175_918_115_1e10,
                -2.580_137_833_794_529_513_0e12,
                -9.518_274_816_138_631_461_6e13,
                -8.607_325_098_621_032_176_6e14,
                1.404_194_185_333_996_143_9e14,
            ],
            [
                1.0,
                6_852.581_373_443_110_097_1,
                8.515_300_102_546_654_437_9e6,
                3.214_602_823_968_569_465_5e9,
                4.292_980_741_745_319_611_3e11,
                2.023_438_116_163_808_435_9e13,
                2.869_993_326_823_392_384_2e14,
                7.121_013_665_152_547_709_6e14,
            ],
        )
    } else if z <= -0.001_200_161_067_219_772_417_3 {
        // W >= -8.913, Y_-5

        rational_function(
            -z / (INV_SQRT_E + sqrt(z - NEG_INV_E)),
            [
                -13.921_651_376_890_072_595,
                -298_789.564_823_880_655_26,
                -1.231_301_993_732_209_233_4e9,
                -1.555_614_908_189_950_897_0e12,
                -6.868_534_110_677_270_873_4e14,
                -1.029_061_627_593_326_683_5e17,
                -4.140_468_370_161_964_847_1e18,
                -1.442_330_999_800_636_839_7e19,
            ],
            [
                1.0,
                26_154.955_236_499_142_433,
                1.239_308_727_744_204_149_4e8,
                1.783_292_270_247_076_111_3e11,
                9.077_260_816_381_085_044_6e13,
                1.631_473_474_005_425_274_1e16,
                8.837_132_386_123_350_453_3e17,
                8.416_662_064_338_501_338_4e18,
            ],
        )
    } else if z <= -0.000_307_288_059_321_914_998_44 {
        // W >= -10.433, Y_-6

        rational_function(
            -z / (INV_SQRT_E + sqrt(z - NEG_INV_E)),
            [
                -15.377_894_224_591_557_534,
                -1.312_231_200_509_697_995_2e6,
                -2.140_815_702_211_173_788_8e10,
                -1.071_828_743_155_781_180_8e14,
                -1.884_935_352_402_773_445_6e17,
                -1.139_485_860_730_931_199_5e20,
                -1.926_155_508_872_914_159_0e22,
                -3.997_845_208_667_690_129_6e23,
            ],
            [
                1.0,
                101_712.867_717_606_200_46,
                1.872_854_594_505_038_118_8e9,
                1.046_961_741_666_440_275_7e13,
                2.070_434_906_012_044_304_9e16,
                1.446_490_790_238_607_449_6e19,
                3.051_043_220_560_890_094_9e21,
                1.139_758_913_979_073_971_7e23,
            ],
        )
    } else if z <= -0.000_077_447_159_838_062_184_354 {
        // W >= -11.946, Y_-7

        rational_function(
            -z / (INV_SQRT_E + sqrt(z - NEG_INV_E)),
            [
                -16.841_701_411_264_981_596,
                -5.779_082_325_757_713_841_6e6,
                -3.775_723_079_125_640_411_6e11,
                -7.571_213_374_258_986_094_1e15,
                -5.347_933_891_601_146_568_5e19,
                -1.308_271_173_229_786_547_6e23,
                -9.146_277_700_452_142_744_0e25,
                -8.960_276_811_926_362_934_0e27,
            ],
            [
                1.,
                401_820.466_662_307_253_28,
                2.921_151_813_690_049_204_6e10,
                6.445_613_537_341_028_907_9e14,
                5.031_180_957_649_953_028_1e18,
                1.387_904_123_971_628_947_8e22,
                1.157_514_616_751_351_622_5e25,
                1.719_922_018_594_775_665_4e27,
            ],
        )
    } else if z <= -4.580_811_969_815_817_317_4e-17 {
        // W >= -41.344, V_-8

        rational_function(
            ln(-z),
            [
                -2.083_626_038_401_643_926_5,
                1.612_243_624_227_149_571_0,
                5.446_426_495_963_720_761_9,
                -3.088_633_112_831_716_010_5,
                0.461_078_291_553_701_378_80,
                -0.023_553_839_118_456_381_330,
                0.000_405_389_041_702_534_047_80,
                -1.794_815_692_251_682_545_8e-6,
            ],
            [
                1.,
                2.369_964_891_270_301_561_0,
                -2.124_944_970_740_481_284_7,
                0.384_809_800_985_884_839_13,
                -0.021_720_009_380_176_605_969,
                0.000_394_058_628_906_086_368_76,
                -1.790_931_206_686_595_790_5e-6,
                3.115_367_330_813_367_145_2e-12,
            ],
        )
    } else if z <= -6.107_367_223_659_479_298_2e-79 {
        // W >= -185.316, V_-9

        rational_function(
            ln(-z),
            [
                0.160_453_837_665_705_414_09,
                2.221_418_252_446_151_402_9,
                -0.941_196_624_920_508_929_71,
                0.091_921_523_818_747_869_300,
                -0.002_906_976_053_317_166_322_4,
                0.000_032_707_247_990_255_961_149,
                -1.248_667_233_688_989_301_8e-7,
                1.224_743_827_986_178_529_1e-10,
            ],
            [
                1.0,
                -0.702_549_960_878_703_322_89,
                0.080_974_347_786_703_195_026,
                -0.002_746_985_002_956_315_393_9,
                0.000_031_943_362_385_183_657_062,
                -1.239_062_068_732_166_643_9e-7,
                1.224_163_611_516_820_199_9e-10,
                -1.027_571_802_054_676_540_0e-17,
            ],
        )
    } else if z < 0.0 {
        // V_-10

        rational_function(
            ln(-z),
            [
                -1.274_217_970_307_544_056_4,
                1.369_665_880_542_138_376_5,
                -0.125_193_453_875_587_832_23,
                0.002_515_572_246_076_384_473_7,
                -0.000_015_748_033_750_499_977_208,
                3.431_608_538_691_378_641_0e-8,
                -2.502_524_288_534_043_853_3e-11,
                4.642_388_501_409_958_335_1e-15,
            ],
            [
                1.,
                -0.114_200_064_741_524_656_94,
                0.002_428_523_383_212_259_594_2,
                -0.000_015_520_907_512_751_723_152,
                3.412_053_476_039_600_226_0e-8,
                -2.498_105_618_645_027_458_7e-11,
                4.641_976_809_305_970_607_9e-15,
                -1.360_871_393_694_260_298_5e-23,
            ],
        )
    } else if z == 0.0 {
        f64::NEG_INFINITY
    } else {
        f64::NAN
    }
}