diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml index c0a4818..c45c21b 100644 --- a/.github/workflows/rust.yml +++ b/.github/workflows/rust.yml @@ -41,7 +41,7 @@ jobs: - uses: actions/checkout@v4 - uses: dtolnay/rust-toolchain@stable - name: test - run: cargo test --all-features --verbose && cargo test --no-default-features -F 24 --verbose && cargo test --no-default-features -F 50 --verbose + run: cargo test --all-features --verbose && cargo test --no-default-features -F 24bits --verbose && cargo test --no-default-features -F 50bits --verbose doc: runs-on: ubuntu-latest diff --git a/Cargo.lock b/Cargo.lock index 188c485..87726c2 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -19,7 +19,7 @@ checksum = "0c4b4d0bd25bd0b74681c0ad21497610ce1b7c91b1022cd21c80c6fbdd9476b0" [[package]] name = "lambert_w" -version = "0.1.1" +version = "0.2.0" dependencies = [ "approx", ] diff --git a/Cargo.toml b/Cargo.toml index d049bc3..ecca71a 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "lambert_w" -version = "0.1.1" +version = "0.2.0" edition = "2021" authors = ["Johanna Sörngård "] categories = ["mathematics"] @@ -16,11 +16,11 @@ repository = "https://github.com/JSorngard/lambert_w" approx = "0.5.1" [features] -default = ["24", "50"] -# Enables the function versions with 50 bits of accuracy. -50 = [] -# Enables the function versions with 24 bits of accuracy. -24 = [] +default = ["24bits", "50bits"] +# Enables the more accurate function versions with 50 bits of accuracy. +50bits = [] +# Enables the faster function versions with 24 bits of accuracy. +24bits = [] # docs.rs-specific configuration. Taken from . [package.metadata.docs.rs] diff --git a/README.md b/README.md index 03ec0c6..304e1fd 100644 --- a/README.md +++ b/README.md @@ -8,8 +8,6 @@ Fast evaluation of the principal and secondary branches of the [Lambert W functi This method uses a piecewise minimax rational approximation of the function. -This crate is a Rust port of the original Fortran 90 code by T. Fukushima. - ## Examples Evaluate the principal branch of the Lambert W function to 50 bits of accuracy: @@ -32,4 +30,9 @@ use approx::assert_abs_diff_eq; let w = lambert_w_0(PI)?; assert_abs_diff_eq!(w, 1.0736581947961492, epsilon = 1e-7); -``` \ No newline at end of file +``` + +## Speed-accuracy trade-off + +The 50-bit accurate versions in the `accurate` module are more accurate, but slightly slower, than the 24-bit accurate versions in the `fast` module. +`fast::lambert_w_0` is around 15% faster than `accurate::lambert_w_0` and `fast::lambert_w_m1` is around 41% faster than `accurate::lambert_w_m1`. \ No newline at end of file diff --git a/src/accurate/dw0c.rs b/src/accurate/dw0c.rs index d44f2a3..41b41a9 100644 --- a/src/accurate/dw0c.rs +++ b/src/accurate/dw0c.rs @@ -1,7 +1,5 @@ //! The original implementation of the principal branch of the Lambert W function by Toshio Fukushima, accurate to 50 bits, ported to Rust. -use crate::LambertW0Error; - /// 50-bit accuracy computation of principal branch of the Lambert W function, W_0(z), /// by piecewise minimax rational function approximation /// @@ -15,13 +13,13 @@ use crate::LambertW0Error; /// rational function approximation with variable transformation" // Formatting this function takes a lot of time, so I have ran `cargo fmt` on it once, and now no one else has to / Johanna. #[rustfmt::skip] -pub fn dw0c(zc: f64) -> Result { +pub fn dw0c(zc: f64) -> Option { if zc < 0.0 { - Err(LambertW0Error::new()) + None } else if zc <= 2.549_893_906_503_473_6 { // W <= 0.893, X_1 let x = zc.sqrt(); - Ok((-0.999_999_999_999_999_9 + Some((-0.999_999_999_999_999_9 + x * (-2.739_966_866_820_366 + x * (0.026_164_207_726_990_4 + x * (6.370_916_807_894_901 @@ -41,7 +39,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 43.613_924_462_669_37 { // W <= 2.754, X_2 let x = zc.sqrt(); - Ok((-0.999_978_018_005_789_1 + Some((-0.999_978_018_005_789_1 + x * (-0.704_157_515_904_836 + x * (2.123_226_083_280_252_8 + x * (2.389_676_070_293_572 @@ -60,7 +58,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 598.453_533_718_782_8 { // W <= 4.821, X_3 let x = zc.sqrt(); - Ok((-0.989_674_203_372_735 + Some((-0.989_674_203_372_735 + x * (0.595_876_806_063_943_8 + x * (1.422_508_301_815_194_3 + x * (0.448_828_891_683_238_1 @@ -79,7 +77,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 8_049.491_985_075_761_5 { // W <= 7.041, X_4 let x = zc.sqrt(); - Ok((-0.773_164_919_972_062_3 + Some((-0.773_164_919_972_062_3 + x * (1.139_133_350_429_670_3 + x * (0.431_161_172_552_170_74 + x * (0.035_773_078_319_037_505 @@ -98,7 +96,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 111_124.954_121_217_82 { // W <= 9.380, X_5 let x = zc.sqrt(); - Ok((0.120_071_016_715_536_88 + Some((0.120_071_016_715_536_88 + x * (0.833_526_408_299_128_3 + x * (0.070_142_775_916_948_34 + x * (0.001_484_635_798_547_512_4 @@ -117,7 +115,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 1.587_042_981_208_229_7e6 { // W <= 11.809, X_6 let x = zc.sqrt(); - Ok((1.722_110_443_993_771_1 + Some((1.722_110_443_993_771_1 + x * (0.399_195_942_864_842_8 + x * (0.007_988_554_014_068_503 + x * (0.000_042_889_742_253_257_923 @@ -136,7 +134,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 2.341_470_840_187_546e7 { // W <= 14.308, X_7 let x = zc.sqrt(); - Ok((3.752_931_402_343_454_3 + Some((3.752_931_402_343_454_3 + x * (0.154_913_426_903_578_07 + x * (0.000_756_631_406_759_007_9 + x * (1.027_160_923_596_997_8e-6 @@ -155,7 +153,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 3.557_647_430_800_996_4e8 { // W <= 16.865, X_8 let x = zc.sqrt(); - Ok((6.019_654_205_560_656 + Some((6.019_654_205_560_656 + x * (0.053_496_672_841_797_86 + x * (0.000_064_340_849_275_316_5 + x * (2.196_909_010_009_596_7e-8 @@ -174,7 +172,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 5.550_171_629_616_363e9 { // W <= 19.468, X_9 let x = zc.sqrt(); - Ok((8.428_026_850_098_97 + Some((8.428_026_850_098_97 + x * (0.017_155_758_546_279_713 + x * (5.083_662_066_982_932e-6 + x * (4.335_490_369_183_258_4e-10 @@ -193,7 +191,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 8.867_470_483_965_778e10 { // W <= 22.112, X_10 let x = zc.sqrt(); - Ok((10.931_063_230_472_498 + Some((10.931_063_230_472_498 + x * (0.005_222_423_454_024_553_5 + x * (3.799_610_571_181_013e-7 + x * (8.030_579_353_341_036e-12 @@ -212,7 +210,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 1.447_779_186_527_290_3e12 { // W <= 24.791, X_11 let x = zc.sqrt(); - Ok((13.502_943_080_893_871 + Some((13.502_943_080_893_871 + x * (0.001_528_463_650_634_626_6 + x * (2.715_696_735_826_234_5e-8 + x * (1.411_039_405_124_216_2e-13 @@ -231,7 +229,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 2.411_145_863_251_185e13 { // W <= 27.500, X_12 let x = zc.sqrt(); - Ok((16.128_076_167_439_016 + Some((16.128_076_167_439_016 + x * (0.000_433_603_851_764_670_7 + x * (1.869_640_387_182_092e-9 + x * (2.369_179_576_690_148_7e-15 @@ -250,7 +248,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 4.089_703_644_260_084_4e14 { // W <= 30.236, X_13 let x = zc.sqrt(); - Ok((18.796_301_105_534_486 + Some((18.796_301_105_534_486 + x * (0.000_119_894_433_396_464_69 + x * (1.246_337_752_867_686_3e-10 + x * (3.821_945_685_801_037e-17 @@ -269,7 +267,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 7.055_590_147_678_997e15 { // W <= 32.996, X_14 let x = zc.sqrt(); - Ok((21.500_582_830_667_334 + Some((21.500_582_830_667_334 + x * (0.000_032_441_943_237_735_277 + x * (8.076_496_341_683_755e-12 + x * (5.948_844_550_612_289e-19 @@ -288,7 +286,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 1.236_660_755_797_672_8e17 { // W <= 35.779, X_15 let x = zc.sqrt(); - Ok((24.235_812_532_416_976 + Some((24.235_812_532_416_976 + x * (8.616_150_599_577_68e-6 + x * (5.103_343_156_186_827e-13 + x * (8.964_239_366_584_964e-21 @@ -307,7 +305,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 2.199_937_348_793_1e18 { // W <= 38.582, X_16 let x = zc.sqrt(); - Ok((26.998_134_347_987_44 + Some((26.998_134_347_987_44 + x * (2.251_225_776_757_228_4e-6 + x * (3.152_123_075_986_696_7e-14 + x * (1.311_403_571_979_063e-22 @@ -326,7 +324,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 3.968_539_219_834_401_6e19 { // W <= 41.404, X_17 let x = zc.sqrt(); - Ok((29.784_546_702_831_97 + Some((29.784_546_702_831_97 + x * (5.797_176_439_217_133e-7 + x * (1.906_987_279_260_195e-15 + x * (1.866_870_087_085_876_3e-24 @@ -345,7 +343,7 @@ pub fn dw0c(zc: f64) -> Result { } else if zc <= 1.412_707_514_527_465_2e104 { // W <= 234.358, U_18 let y = zc.ln(); - Ok((0.744_134_994_601_267_8 + Some((0.744_134_994_601_267_8 + y * (0.414_032_436_180_059_14 + y * (0.260_125_641_667_734_16 + y * (0.021_450_457_095_960_294 @@ -364,7 +362,7 @@ pub fn dw0c(zc: f64) -> Result { } else { // U_19 let y = zc.ln(); - Ok((-0.615_144_128_127_297_6 + Some((-0.615_144_128_127_297_6 + y * (0.679_793_101_336_309_3 + y * (0.089_685_353_704_585_82 + y * (0.001_564_494_148_398_938 diff --git a/src/accurate/dwm1c.rs b/src/accurate/dwm1c.rs index b344386..5c8dd07 100644 --- a/src/accurate/dwm1c.rs +++ b/src/accurate/dwm1c.rs @@ -1,7 +1,6 @@ //! The original implementation of the secondary branch of the Lambert W function by Toshio Fukushima, accurate to 50 bits, ported to Rust. use super::super::{X0, Z0}; -use crate::{LambertWm1Error, LambertWm1ErrorReason}; /// 50-bit accuracy computation of secondary branch of the Lambert W function, W_-1(z), /// by piecewise minimax rational function approximation @@ -16,15 +15,13 @@ use crate::{LambertWm1Error, LambertWm1ErrorReason}; /// rational function approximation with variable transformation" // Formatting this function takes a lot of time, so I have ran `cargo fmt` on it once, and now no one else has to / Johanna. #[rustfmt::skip] -pub fn dwm1c(z: f64, zc: f64) -> Result { +pub fn dwm1c(z: f64, zc: f64) -> Option { if zc < 0.0 { - Err(LambertWm1Error::new( - LambertWm1ErrorReason::TooSmallArgument, - )) + None } else if z <= -0.3542913309442164 { // W >= -1.3, X_-1 let x = zc.sqrt(); - Ok((-1.0000000000000001110 + Some((-1.0000000000000001110 + x * (4.296_301_617_877_713 + x * (-4.099_140_792_400_746 + x * (-6.844_284_220_083_331 @@ -43,7 +40,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result { } else if z <= -0.188_726_882_822_894_35 { // W >= -2.637, Y_-1 let x = -z / (X0 + (z - Z0).sqrt()); - Ok((-8.225_315_526_444_685 + Some((-8.225_315_526_444_685 + x * (-813.207_067_320_014_9 + x * (-15_270.113_237_678_51 + x * (-79_971.585_089_674_15 @@ -61,7 +58,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result { } else if z <= -0.060_497_597_226_958_34 { // W >= -4.253, Y_-2 let x = -z / (X0 + (z - Z0).sqrt()); - Ok((-9.618_412_744_335_403 + Some((-9.618_412_744_335_403 + x * (-3_557.856_904_301_800_6 + x * (-254_015.593_112_843_8 + x * (-5.392_389_363_067_063_5e6 @@ -80,7 +77,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result { } else if z <= -0.017_105_334_740_676_01 { // W >= -5.832, Y_-3 let x = -z / (X0 + (z - Z0).sqrt()); - Ok((-11.038_489_462_297_466 + Some((-11.038_489_462_297_466 + x * (-15_575.812_882_656_619 + x * (-4.249_294_730_489_777e6 + x * (-3.517_024_593_880_342e8 @@ -99,7 +96,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result { } else if z <= -0.004_595_496_212_794_371 { // W >= -7.382, Y_-4 let x = -z / (X0 + (z - Z0).sqrt()); - Ok((-12.474_405_916_395_746 + Some((-12.474_405_916_395_746 + x * (-68_180.335_575_543_78 + x * (-7.184_659_984_562_01e7 + x * (-2.314_268_822_175_918_2e10 @@ -118,7 +115,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result { } else if z <= -0.001_200_161_067_219_772_4 { // W >= -8.913, Y_-5 let x = -z / (X0 + (z - Z0).sqrt()); - Ok((-13.921_651_376_890_072 + Some((-13.921_651_376_890_072 + x * (-298_789.564_823_880_7 + x * (-1.231_301_993_732_209_2e9 + x * (-1.555_614_908_189_951e12 @@ -137,7 +134,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result { } else if z <= -0.000_307_288_059_321_915 { // W >= -10.433, Y_-6 let x = -z / (X0 + (z - Z0).sqrt()); - Ok((-15.377_894_224_591_557 + Some((-15.377_894_224_591_557 + x * (-1.312_231_200_509_698e6 + x * (-2.140_815_702_211_173_6e10 + x * (-1.071_828_743_155_781_3e14 @@ -156,7 +153,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result { } else if z <= -0.000_077_447_159_838_062_18 { // W >= -11.946, Y_-7 let x = -z / (X0 + (z - Z0).sqrt()); - Ok((-16.841_701_411_264_98 + Some((-16.841_701_411_264_98 + x * (-5.779_082_325_757_714e6 + x * (-3.775_723_079_125_64e11 + x * (-7.571_213_374_258_986e15 @@ -175,7 +172,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result { } else if z <= -4.580_811_969_815_817_5e-17 { // W >= -41.344, V_-8 let u = (-z).ln(); - Ok((-2.083_626_038_401_644 + Some((-2.083_626_038_401_644 + u * (1.612_243_624_227_149_6 + u * (5.446_426_495_963_720_5 + u * (-3.088_633_112_831_716 @@ -194,7 +191,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result { } else if z <= -6.107_367_223_659_479e-79 { // W >= -185.316, V_-9 let u = (-z).ln(); - Ok((0.160_453_837_665_705_42 + Some((0.160_453_837_665_705_42 + u * (2.221_418_252_446_151_4 + u * (-0.941_196_624_920_508_9 + u * (0.091_921_523_818_747_87 @@ -213,7 +210,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result { } else if z < 0.0 { // V_-10 let u = (-z).ln(); - Ok((-1.274_217_970_307_544 + Some((-1.274_217_970_307_544 + u * (1.369_665_880_542_138_4 + u * (-0.125_193_453_875_587_83 + u * (0.002_515_572_246_076_384_3 @@ -230,8 +227,6 @@ pub fn dwm1c(z: f64, zc: f64) -> Result { + u * (4.641_976_809_305_971e-15 + u * (-1.360_871_393_694_260_3e-23))))))))) } else { - Err(LambertWm1Error::new( - LambertWm1ErrorReason::PositiveArgument, - )) + None } } diff --git a/src/accurate/mod.rs b/src/accurate/mod.rs index 8fedaa0..edcb337 100644 --- a/src/accurate/mod.rs +++ b/src/accurate/mod.rs @@ -7,54 +7,41 @@ use dw0c::dw0c; use dwm1c::dwm1c; use super::Z0; -use crate::error::{LambertW0Error, LambertWm1Error}; -/// Computes the principal branch of the Lambert W function, W_0(`z`) to 50 bits of accuracy by piecewise minimax rational function approximation and variable transformations. +/// Computes the principal branch of the Lambert W function, W_0(`z`), to 50 bits of accuracy, if `z` is larger than -1/e. /// -/// Uses the [method of Toshio Fukushima](https://www.researchgate.net/publication/346309410_Precise_and_fast_computation_of_Lambert_W_function_by_piecewise_minimax_rational_function_approximation_with_variable_transformation). -/// -/// # Errors -/// -/// Returns an error if `z` is less than -1/e. +/// Uses the piecewise minimax rational function approximation method of Toshio Fukushima. /// /// # Example /// /// ``` -/// # use lambert_w::LambertW0Error; /// use lambert_w::accurate::lambert_w_0; /// use approx::assert_abs_diff_eq; /// use core::f64::consts::PI; /// -/// let w = lambert_w_0(PI)?; +/// let w = lambert_w_0(PI).unwrap(); /// /// assert_abs_diff_eq!(w, 1.0736581947961492); -/// # Ok::<(), LambertW0Error >(()) /// ``` -pub fn lambert_w_0(z: f64) -> Result { +pub fn lambert_w_0(z: f64) -> Option { dw0c(z - Z0) } -/// Computes the secondary branch of the Lambert W function, W_-1(`z`), to 50 bits of accuracy by piecewise minimax rational function approximation and variable transformations. -/// -/// Uses the [method of Toshio Fukushima](https://www.researchgate.net/publication/346309410_Precise_and_fast_computation_of_Lambert_W_function_by_piecewise_minimax_rational_function_approximation_with_variable_transformation). -/// -/// # Errors +/// Computes the secondary branch of the Lambert W function, W_-1(`z`), to 50 bits of accuracy, if -1/e < `z` <= 0. /// -/// Returns an error if `z` is positive or less than -1/e. +/// Uses the piecewise minimax rational function approximation method of Toshio Fukushima. /// /// # Example /// /// ``` -/// # use lambert_w::LambertWm1Error; /// use lambert_w::accurate::lambert_w_m1; /// use approx::assert_abs_diff_eq; /// use core::f64::consts::PI; /// -/// let w = lambert_w_m1(-1.0/PI)?; +/// let w = lambert_w_m1(-1.0/PI).unwrap(); /// /// assert_abs_diff_eq!(w, -1.6385284199703634); -/// # Ok::<(), LambertWm1Error>(()) /// ``` -pub fn lambert_w_m1(z: f64) -> Result { +pub fn lambert_w_m1(z: f64) -> Option { dwm1c(z, z - Z0) } diff --git a/src/error.rs b/src/error.rs deleted file mode 100644 index 4d6c8d5..0000000 --- a/src/error.rs +++ /dev/null @@ -1,76 +0,0 @@ -use core::fmt; -use std::backtrace::Backtrace; -use std::error::Error; - -/// The error returned by the Lambert W_0 functions when the input is less than -1/e. -#[derive(Debug)] -pub struct LambertW0Error(Backtrace); - -impl LambertW0Error { - pub(crate) fn new() -> Self { - Self(Backtrace::capture()) - } - - /// Returns a [`Backtrace`] to where the error was created. - /// - /// This backtrace was captured with [`Backtrace::capture`], - /// see it for more information about how to make this display information when printed. - pub fn backtrace(&self) -> &Backtrace { - &self.0 - } -} - -impl fmt::Display for LambertW0Error { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "the argument was less than -1/e") - } -} - -impl Error for LambertW0Error {} - -/// The error returned by the Lambert W_-1 functions when the input is positive or less than -1/e. -#[derive(Debug)] -pub struct LambertWm1Error { - backtrace: Backtrace, - reason: LambertWm1ErrorReason, -} - -/// The reason for the error in the Lambert W_-1 functions. -#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] -pub enum LambertWm1ErrorReason { - TooSmallArgument, - PositiveArgument, -} - -impl LambertWm1Error { - /// Returns a [`Backtrace`] to where the error was created. - /// - /// This backtrace was captured with [`Backtrace::capture`], - /// see it for more information about how to make this display information when printed. - pub fn backtrace(&self) -> &Backtrace { - &self.backtrace - } - - /// Returns the reason for the error. - pub fn reason(&self) -> LambertWm1ErrorReason { - self.reason - } - - pub(crate) fn new(reason: LambertWm1ErrorReason) -> Self { - Self { - backtrace: Backtrace::capture(), - reason, - } - } -} - -impl fmt::Display for LambertWm1Error { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - match self.reason { - LambertWm1ErrorReason::TooSmallArgument => write!(f, "the argument was less than -1/e"), - LambertWm1ErrorReason::PositiveArgument => write!(f, "the argument was positive"), - } - } -} - -impl Error for LambertWm1Error {} diff --git a/src/fast/mod.rs b/src/fast/mod.rs index 9e6be93..37a03f9 100644 --- a/src/fast/mod.rs +++ b/src/fast/mod.rs @@ -6,54 +6,40 @@ mod swm1; use sw0::sw0; use swm1::swm1; -use crate::{LambertW0Error, LambertWm1Error}; - -/// Computes the principal branch of the Lambert W function, W_0(`z`), to 24 bits of accuracy by piecewise minimax rational function approximation and variable transformations. -/// -/// Uses the [method of Toshio Fukushima](https://www.researchgate.net/publication/346309410_Precise_and_fast_computation_of_Lambert_W_function_by_piecewise_minimax_rational_function_approximation_with_variable_transformation). -/// -/// # Errors +/// Computes the principal branch of the Lambert W function, W_0(`z`), to 24 bits of accuracy, if `z` > -1/e. /// -/// Returns an error if `z` is less than -1/e. +/// Uses the piecewise minimax rational function approximation method of Toshio Fukushima. /// /// # Example /// /// ``` -/// # use lambert_w::LambertW0Error; /// use lambert_w::fast::lambert_w_0; /// use approx::assert_abs_diff_eq; /// use core::f64::consts::PI; /// -/// let w = lambert_w_0(PI)?; +/// let w = lambert_w_0(PI).unwrap(); /// /// assert_abs_diff_eq!(w, 1.0736581947961492, epsilon = 1e-7); -/// # Ok::<(), LambertW0Error>(()) /// ``` -pub fn lambert_w_0(z: f64) -> Result { +pub fn lambert_w_0(z: f64) -> Option { sw0(z) } -/// Computes the secondary branch of the Lambert W function, W_-1(`z`), to 24 bits of accuracy by piecewise minimax rational function approximation and variable transformations. -/// -/// Uses the [method of Toshio Fukushima](https://www.researchgate.net/publication/346309410_Precise_and_fast_computation_of_Lambert_W_function_by_piecewise_minimax_rational_function_approximation_with_variable_transformation). -/// -/// # Errors +/// Computes the secondary branch of the Lambert W function, W_-1(`z`), to 24 bits of accuracy, if -1/e < `z` <= 0. /// -/// Returns an error if `z` is positive or less than -1/e. +/// Uses the piecewise minimax rational function approximation method of Toshio Fukushima. /// /// # Example /// /// ``` -/// # use lambert_w::LambertWm1Error; /// use lambert_w::fast::lambert_w_m1; /// use approx::assert_abs_diff_eq; /// use core::f64::consts::PI; /// -/// let w = lambert_w_m1(-1.0/PI)?; +/// let w = lambert_w_m1(-1.0/PI).unwrap(); /// /// assert_abs_diff_eq!(w, -1.6385284199703634, epsilon = 1e-7); -/// # Ok::<(), LambertWm1Error>(()) /// ``` -pub fn lambert_w_m1(z: f64) -> Result { +pub fn lambert_w_m1(z: f64) -> Option { swm1(z) } diff --git a/src/fast/sw0.rs b/src/fast/sw0.rs index 550fcfa..bb92305 100644 --- a/src/fast/sw0.rs +++ b/src/fast/sw0.rs @@ -1,7 +1,6 @@ //! The original implementation of the principal branch of the Lambert W function by Toshio Fukushima, accurate to 24 bits, ported to Rust. use super::super::Z0; -use crate::LambertW0Error; /// 24-bit accuracy computation of principal branch of the Lambert W function, W_0(z), /// by piecewise minimax rational function approximation @@ -12,180 +11,218 @@ use crate::LambertW0Error; /// Reference: T. Fukushima (2020) to be submitted /// "Precise and fast computation of Lambert W-functions by piecewise /// rational function approximation with variable transformation" -pub fn sw0(z: f64) -> Result { +pub fn sw0(z: f64) -> Option { if z < Z0 { - Err(LambertW0Error::new()) + None } else if z <= 2.008_217_811_584_472_7 { // W <= 0.854, X_1 let x = (z - Z0).sqrt(); - Ok((-0.999_999_940_395_401_9 - + x * (0.055_730_052_161_777_8 - + x * (2.126_973_249_105_317_3 - + x * (0.813_511_236_783_528_8 + x * 0.016_324_880_146_070_16)))) - / (1. - + x * (2.275_906_559_863_465 - + x * (1.367_597_013_868_904 + x * 0.186_158_234_528_316_23)))) + Some( + (-0.999_999_940_395_401_9 + + x * (0.055_730_052_161_777_8 + + x * (2.126_973_249_105_317_3 + + x * (0.813_511_236_783_528_8 + x * 0.016_324_880_146_070_16)))) + / (1. + + x * (2.275_906_559_863_465 + + x * (1.367_597_013_868_904 + x * 0.186_158_234_528_316_23))), + ) } else if z <= 30.539_142_109_510_895 { // W <= 2.502, X_2 let x = (z - Z0).sqrt(); - Ok((-0.985_519_709_059_991 - + x * (1.077_497_573_381_351_7 - + x * (0.871_751_030_681_775 + x * 0.054_352_728_608_275_766))) - / (1. - + x * (1.186_101_403_701_543_4 - + x * (0.249_962_984_308_281_62 + x * 0.006_881_368_648_675_912)))) + Some( + (-0.985_519_709_059_991 + + x * (1.077_497_573_381_351_7 + + x * (0.871_751_030_681_775 + x * 0.054_352_728_608_275_766))) + / (1. + + x * (1.186_101_403_701_543_4 + + x * (0.249_962_984_308_281_62 + x * 0.006_881_368_648_675_912))), + ) } else if z <= 371.669_843_713_757_76 { // W <= 4.430, X_3 let x = (z - Z0).sqrt(); - Ok((-0.762_397_113_463_688_9 - + x * (1.231_773_161_336_359_6 - + x * (0.243_424_471_130_566_95 + x * 0.004_320_601_393_878_236))) - / (1. - + x * (0.579_386_215_035_869_1 - + x * (0.046_601_427_736_078_775 + x * 0.000_435_128_175_674_741_1)))) + Some( + (-0.762_397_113_463_688_9 + + x * (1.231_773_161_336_359_6 + + x * (0.243_424_471_130_566_95 + x * 0.004_320_601_393_878_236))) + / (1. + + x * (0.579_386_215_035_869_1 + + x * (0.046_601_427_736_078_775 + x * 0.000_435_128_175_674_741_1))), + ) } else if z <= 4_705.918_954_265_969 { // W <= 6.574, X_4 let x = (z - Z0).sqrt(); - Ok((0.085_801_247_434_391_38 - + x * (0.825_397_980_997_483_4 - + x * (0.039_781_960_760_329_076 + x * 0.000_187_855_789_275_838))) - / (1. - + x * (0.213_380_768_170_801_4 - + x * (0.005_462_672_039_792_693_5 + x * 0.000_015_449_534_481_294_754)))) + Some( + (0.085_801_247_434_391_38 + + x * (0.825_397_980_997_483_4 + + x * (0.039_781_960_760_329_076 + x * 0.000_187_855_789_275_838))) + / (1. + + x * (0.213_380_768_170_801_4 + + x * (0.005_462_672_039_792_693_5 + x * 0.000_015_449_534_481_294_754))), + ) } else if z <= 64_640.797_355_310_09 { // W <= 8.892, X_5 let x = (z - Z0).sqrt(); - Ok((1.621_924_538_347_016_9 - + x * (0.388_691_451_325_166_64 - + x * (0.004_575_064_267_850_351_5 + x * 5.538_467_214_864_45e-6))) - / (1. - + x * (0.065_219_460_735_182_41 - + x * (0.000_478_827_607_890_225_1 + x * 3.809_482_814_629_24e-7)))) + Some( + (1.621_924_538_347_016_9 + + x * (0.388_691_451_325_166_64 + + x * (0.004_575_064_267_850_351_5 + x * 5.538_467_214_864_45e-6))) + / (1. + + x * (0.065_219_460_735_182_41 + + x * (0.000_478_827_607_890_225_1 + x * 3.809_482_814_629_24e-7))), + ) } else if z <= 965_649.030_871_163_2 { // W <= 11.351, X_6 let x = (z - Z0).sqrt(); - Ok((3.621_899_608_569_592 - + x * (0.148_846_467_548_801_6 - + x * (0.000_424_696_224_099_984 + x * 1.279_017_971_037_421_7e-7))) - / (1. - + x * (0.017_985_659_319_608_747 - + x * (0.000_035_446_449_757_357_85 + x * 7.506_249_296_303_705e-9)))) + Some( + (3.621_899_608_569_592 + + x * (0.148_846_467_548_801_6 + + x * (0.000_424_696_224_099_984 + x * 1.279_017_971_037_421_7e-7))) + / (1. + + x * (0.017_985_659_319_608_747 + + x * (0.000_035_446_449_757_357_85 + x * 7.506_249_296_303_705e-9))), + ) } else if z <= 1.559_333_422_803_816_6e7 { // W <= 13.928, X_7 let x = (z - Z0).sqrt(); - Ok((5.907_336_973_960_809 - + x * (0.050_053_653_594_737_11 - + x * (0.000_034_072_148_625_204_7 + x * 2.481_206_469_365_548_5e-9))) - / (1. - + x * (0.004_655_899_001_684_321 - + x * (2.344_944_586_080_881e-6 + x * 1.263_142_996_480_846_2e-10)))) + Some( + (5.907_336_973_960_809 + + x * (0.050_053_653_594_737_11 + + x * (0.000_034_072_148_625_204_7 + x * 2.481_206_469_365_548_5e-9))) + / (1. + + x * (0.004_655_899_001_684_321 + + x * (2.344_944_586_080_881e-6 + x * 1.263_142_996_480_846_2e-10))), + ) } else if z <= 2.702_564_027_724_190_4e8 { // W <= 16.605, X_8 let x = (z - Z0).sqrt(); - Ok((8.382_600_584_819_551 - + x * (0.015_360_346_475_232_501 - + x * (2.443_338_439_786_936_6e-6 + x * 4.185_680_326_411_855e-11))) - / (1. - + x * (0.001_150_742_322_378_586_9 - + x * (1.422_142_847_481_351_6e-7 + x * 1.873_917_202_662_012_4e-12)))) + Some( + (8.382_600_584_819_551 + + x * (0.015_360_346_475_232_501 + + x * (2.443_338_439_786_936_6e-6 + x * 4.185_680_326_411_855e-11))) + / (1. + + x * (0.001_150_742_322_378_586_9 + + x * (1.422_142_847_481_351_6e-7 + x * 1.873_917_202_662_012_4e-12))), + ) } else if z <= 4.995_018_739_704_195e9 { // W <= 19.368, X_9 let x = (z - Z0).sqrt(); - Ok((10.996_674_803_992_551 - + x * (0.004_394_213_889_867_383 - + x * (1.596_666_535_484_678e-7 + x * 6.266_538_284_496_873e-13))) - / (1. - + x * (0.000_273_837_576_757_036_47 - + x * (8.015_706_231_969_03e-9 + x * 2.495_698_215_887_173e-14)))) + Some( + (10.996_674_803_992_551 + + x * (0.004_394_213_889_867_383 + + x * (1.596_666_535_484_678e-7 + x * 6.266_538_284_496_873e-13))) + / (1. + + x * (0.000_273_837_576_757_036_47 + + x * (8.015_706_231_969_03e-9 + x * 2.495_698_215_887_173e-14))), + ) } else if z <= 9.791_115_441_672_696e10 { // W <= 22.207, X_10 let x = (z - Z0).sqrt(); - Ok((13.719_833_802_350_86 - + x * (0.001_187_444_380_520_229_2 - + x * (9.630_338_120_016_468e-9 + x * 8.443_452_423_226_163e-15))) - / (1. - + x * (0.000_063_056_372_424_395_35 - + x * (4.235_876_603_109_884_3e-10 + x * 3.020_540_500_543_447_3e-16)))) + Some( + (13.719_833_802_350_86 + + x * (0.001_187_444_380_520_229_2 + + x * (9.630_338_120_016_468e-9 + x * 8.443_452_423_226_163e-15))) + / (1. + + x * (0.000_063_056_372_424_395_35 + + x * (4.235_876_603_109_884_3e-10 + x * 3.020_540_500_543_447_3e-16))), + ) } else if z <= 2.025_975_385_630_21e12 { // W <= 25.114, X_11 let x = (z - Z0).sqrt(); - Ok((16.533_119_481_561_616 - + x * (0.000_305_831_257_519_080_4 - + x * (5.411_294_663_372_01e-10 + x * 1.034_713_033_370_471_1e-16))) - / (1. - + x * (0.000_014_099_161_212_376_34 - + x * (2.112_109_541_235_469_5e-11 + x * 3.352_692_715_745_247e-18)))) + Some( + (16.533_119_481_561_616 + + x * (0.000_305_831_257_519_080_4 + + x * (5.411_294_663_372_01e-10 + x * 1.034_713_033_370_471_1e-16))) + / (1. + + x * (0.000_014_099_161_212_376_34 + + x * (2.112_109_541_235_469_5e-11 + x * 3.352_692_715_745_247e-18))), + ) } else if z <= 4.407_744_425_147_794e13 { // W <= 28.082, X_12 let x = (z - Z0).sqrt(); - Ok((19.423_519_260_478_578 - + x * (0.000_075_559_269_761_977_81 - + x * (2.853_002_312_078_307_5e-11 + x * 1.162_962_709_646_357_9e-18))) - / (1. - + x * (3.069_209_278_972_785_8e-6 - + x * (9.986_661_305_031_147e-13 + x * 3.437_671_711_698_392e-20)))) + Some( + (19.423_519_260_478_578 + + x * (0.000_075_559_269_761_977_81 + + x * (2.853_002_312_078_307_5e-11 + x * 1.162_962_709_646_357_9e-18))) + / (1. + + x * (3.069_209_278_972_785_8e-6 + + x * (9.986_661_305_031_147e-13 + x * 3.437_671_711_698_392e-20))), + ) } else if z <= 1.004_838_215_057_150_5e15 { // W <= 31.106, X_13 let x = (z - Z0).sqrt(); - Ok((22.381_576_050_041_915 - + x * (0.000_017_994_724_029_162_553 - + x * (1.419_487_642_040_223e-12 + x * 1.207_110_515_438_583e-20))) - / (1. - + x * (6.518_396_280_665_677e-7 - + x * (4.495_866_571_281_253_6e-14 + x * 3.275_542_924_502_358e-22)))) + Some( + (22.381_576_050_041_915 + + x * (0.000_017_994_724_029_162_553 + + x * (1.419_487_642_040_223e-12 + x * 1.207_110_515_438_583e-20))) + / (1. + + x * (6.518_396_280_665_677e-7 + + x * (4.495_866_571_281_253_6e-14 + x * 3.275_542_924_502_358e-22))), + ) } else if z <= 2.393_255_260_235_983_6e16 { // W <= 34.182, X_14 let x = (z - Z0).sqrt(); - Ok((25.400_105_417_092_067 - + x * (4.146_737_838_657_924e-6 - + x * (6.696_269_721_968_18e-14 + x * 1.163_790_515_950_647_6e-22))) - / (1. - + x * (1.352_980_135_703_041_8e-7 - + x * (1.933_608_178_532_575e-15 + x * 2.914_939_619_981_625_7e-24)))) + Some( + (25.400_105_417_092_067 + + x * (4.146_737_838_657_924e-6 + + x * (6.696_269_721_968_18e-14 + x * 1.163_790_515_950_647_6e-22))) + / (1. + + x * (1.352_980_135_703_041_8e-7 + + x * (1.933_608_178_532_575e-15 + x * 2.914_939_619_981_625_7e-24))), + ) } else if z <= 5.939_799_659_746_575e17 { // W <= 37.306, X_15 let x = (z - Z0).sqrt(); - Ok((28.473_455_626_379_916 - + x * (9.274_682_469_309_406e-7 - + x * (3.006_899_015_933_680_6e-15 + x * 1.047_355_759_182_202_7e-24))) - / (1. - + x * (2.748_648_970_452_172_8e-8 - + x * (7.967_898_707_103_613e-17 + x * 2.433_166_636_706_153e-26)))) + Some( + (28.473_455_626_379_916 + + x * (9.274_682_469_309_406e-7 + + x * (3.006_899_015_933_680_6e-15 + x * 1.047_355_759_182_202_7e-24))) + / (1. + + x * (2.748_648_970_452_172_8e-8 + + x * (7.967_898_707_103_613e-17 + x * 2.433_166_636_706_153e-26))), + ) } else if z <= 1.532_693_858_990_176_7e19 { // W <= 40.475, X_16 let x = (z - Z0).sqrt(); - Ok((31.597_055_437_846_36 - + x * (2.018_422_527_678_632_4e-7 - + x * (1.289_578_819_651_291e-16 + x * 8.836_117_471_410_11e-27))) - / (1. - + x * (5.472_394_512_609_85e-9 - + x * (3.153_772_917_992_919e-18 + x * 1.912_203_513_257_167e-28)))) + Some( + (31.597_055_437_846_36 + + x * (2.018_422_527_678_632_4e-7 + + x * (1.289_578_819_651_291e-16 + x * 8.836_117_471_410_11e-27))) + / (1. + + x * (5.472_394_512_609_85e-9 + + x * (3.153_772_917_992_919e-18 + x * 1.912_203_513_257_167e-28))), + ) } else if z <= 4.103_565_939_888_539_6e20 { // W <= 43.687, X_17 let x = (z - Z0).sqrt(); - Ok((34.767_124_490_414_52 - + x * (4.283_079_924_069_894_4e-8 - + x * (5.297_588_412_103_653e-18 + x * 7.014_551_539_215_878e-29))) - / (1. - + x * (1.068_930_112_769_633_2e-9 - + x * (1.201_669_906_178_942_4e-19 + x * 1.419_524_481_080_098_4e-30)))) + Some( + (34.767_124_490_414_52 + + x * (4.283_079_924_069_894_4e-8 + + x * (5.297_588_412_103_653e-18 + x * 7.014_551_539_215_878e-29))) + / (1. + + x * (1.068_930_112_769_633_2e-9 + + x * (1.201_669_906_178_942_4e-19 + x * 1.419_524_481_080_098_4e-30))), + ) } else if z <= 2.172_370_661_049_060_6e141 { // W <= 319.673, U_18 let y = z.ln(); - Ok((-0.607_023_733_718_462 - + y * (0.698_287_163_225_269_8 - + y * (0.075_795_135_081_824_75 + y * 0.000_516_692_560_817_372_5))) - / (1. - + y * (0.079_048_429_972_306_02 - + y * (0.000_517_609_908_992_059_8 + y * (-4.243_840_393_198_107e-10))))) + Some( + (-0.607_023_733_718_462 + + y * (0.698_287_163_225_269_8 + + y * (0.075_795_135_081_824_75 + y * 0.000_516_692_560_817_372_5))) + / (1. + + y * (0.079_048_429_972_306_02 + + y * (0.000_517_609_908_992_059_8 + y * (-4.243_840_393_198_107e-10)))), + ) } else { // U_19 let y = z.ln(); - Ok((-3.132_005_602_886_366 - + y * (0.948_894_657_265_326 - + y * (0.008_317_815_296_164_44 + y * 5.558_784_815_783_349e-6))) - / (1. - + y * (0.008_365_681_867_773_006 - + y * (5.559_715_493_597_327e-6 + y * (-3.748_153_583_315_12e-14))))) + Some( + (-3.132_005_602_886_366 + + y * (0.948_894_657_265_326 + + y * (0.008_317_815_296_164_44 + y * 5.558_784_815_783_349e-6))) + / (1. + + y * (0.008_365_681_867_773_006 + + y * (5.559_715_493_597_327e-6 + y * (-3.748_153_583_315_12e-14)))), + ) } } diff --git a/src/fast/swm1.rs b/src/fast/swm1.rs index 4de09a6..ff0a92d 100644 --- a/src/fast/swm1.rs +++ b/src/fast/swm1.rs @@ -1,7 +1,6 @@ //! The original implementation of the secondary branch of the Lambert W function by Toshio Fukushima, accurate to 24 bits, ported to Rust. use super::super::{X0, Z0}; -use crate::{LambertWm1Error, LambertWm1ErrorReason}; /// 24-bit accuracy computation of secondary branch of the Lambert W function, W_-1(z), /// defined as the solution of nonlinear equation, W exp(W) = z, when W < -1 @@ -13,104 +12,120 @@ use crate::{LambertWm1Error, LambertWm1ErrorReason}; /// Reference: T. Fukushima (2020) to be submitted /// "Precise and fast computation of Lambert W-functions by piecewise /// rational function approximation with variable transformation" -pub fn swm1(z: f64) -> Result { +pub fn swm1(z: f64) -> Option { if z < Z0 { - Err(LambertWm1Error::new( - LambertWm1ErrorReason::TooSmallArgument, - )) + None } else if z <= -0.207_293_777_640_384_15 { // W >= -2.483, Y_-1 let x = -z / (X0 + (z - Z0).sqrt()); - Ok((-6.383_722_822_801_905 - + x * (-74.968_653_259_594_05 - + x * (-19.714_821_552_432_483 + x * 70.677_326_667_809_24))) - / (1. - + x * (24.295_836_951_878_69 - + x * (64.112_460_611_386_04 + x * 17.994_497_369_039_312)))) + Some( + (-6.383_722_822_801_905 + + x * (-74.968_653_259_594_05 + + x * (-19.714_821_552_432_483 + x * 70.677_326_667_809_24))) + / (1. + + x * (24.295_836_951_878_69 + + x * (64.112_460_611_386_04 + x * 17.994_497_369_039_312))), + ) } else if z <= -0.071_507_705_083_841_95 { // W >= -4.032, Y_-2 let x = -z / (X0 + (z - Z0).sqrt()); - Ok((-7.723_328_481_229_978 - + x * (-352.484_690_970_429_14 - + x * (-1_242.008_890_368_567_8 + x * 1_171.647_596_062_049_8))) - / (1. - + x * (77.681_242_588_997_4 - + x * (648.564_312_140_752_5 + x * 566.701_549_764_361_6)))) + Some( + (-7.723_328_481_229_978 + + x * (-352.484_690_970_429_14 + + x * (-1_242.008_890_368_567_8 + x * 1_171.647_596_062_049_8))) + / (1. + + x * (77.681_242_588_997_4 + + x * (648.564_312_140_752_5 + x * 566.701_549_764_361_6))), + ) } else if z <= -0.020_704_412_621_717_48 { // W >= -5.600, Y_-3 let x = -z / (X0 + (z - Z0).sqrt()); - Ok((-9.137_773_141_758_155 - + x * (-1_644.724_479_150_889 - + x * (-28_105.096_098_779_683 + x * 3_896.079_810_390_921_5))) - / (1. - + x * (272.375_261_351_239_7 - + x * (7_929.224_261_291_35 + x * 23_980.122_860_821_313)))) + Some( + (-9.137_773_141_758_155 + + x * (-1_644.724_479_150_889 + + x * (-28_105.096_098_779_683 + x * 3_896.079_810_390_921_5))) + / (1. + + x * (272.375_261_351_239_7 + + x * (7_929.224_261_291_35 + x * 23_980.122_860_821_313))), + ) } else if z <= -0.005_480_012_945_209_444 { // W >= -7.178, Y_-4 let x = -z / (X0 + (z - Z0).sqrt()); - Ok((-10.603_388_239_566_373 - + x * (-7_733.348_521_498_648 - + x * (-575_482.407_079_644_3 + x * (-2.154_552_604_188_978e6)))) - / (1. - + x * (1_021.793_856_606_681_7 - + x * (111_300.229_154_865_21 + x * 1.261_425_640_008_844_2e6)))) + Some( + (-10.603_388_239_566_373 + + x * (-7_733.348_521_498_648 + + x * (-575_482.407_079_644_3 + x * (-2.154_552_604_188_978e6)))) + / (1. + + x * (1_021.793_856_606_681_7 + + x * (111_300.229_154_865_21 + x * 1.261_425_640_008_844_2e6))), + ) } else if z <= -0.001_367_466_989_250_804_2 { // W >= -8.766, Y_-5 let x = -z / (X0 + (z - Z0).sqrt()); - Ok((-12.108_699_273_343_438 - + x * (-36_896.535_108_166_376 - + x * (-1.183_112_672_010_605_4e7 + x * (-2.756_583_081_394_092_4e8)))) - / (1. - + x * (4_044.975_306_488_07 - + x * (1.741_827_761_903_000_5e6 + x * 7.843_690_738_080_69e7)))) + Some( + (-12.108_699_273_343_438 + + x * (-36_896.535_108_166_376 + + x * (-1.183_112_672_010_605_4e7 + x * (-2.756_583_081_394_092_4e8)))) + / (1. + + x * (4_044.975_306_488_07 + + x * (1.741_827_761_903_000_5e6 + x * 7.843_690_738_080_69e7))), + ) } else if z <= -0.000_326_142_267_310_725_66 { // W >= -10.367, Y_-6 let x = -z / (X0 + (z - Z0).sqrt()); - Ok((-13.646_761_936_746_191 - + x * (-179_086.115_857_151_48 - + x * (-2.508_463_493_521_464_2e8 + x * (-2.934_370_049_483_371_4e10)))) - / (1. - + x * (16_743.826_607_737_14 - + x * (2.980_965_094_601_174_4e7 + x * 5.573_951_481_695_8e9)))) + Some( + (-13.646_761_936_746_191 + + x * (-179_086.115_857_151_48 + + x * (-2.508_463_493_521_464_2e8 + x * (-2.934_370_049_483_371_4e10)))) + / (1. + + x * (16_743.826_607_737_14 + + x * (2.980_965_094_601_174_4e7 + x * 5.573_951_481_695_8e9))), + ) } else if z <= -0.000_074_906_612_036_101_44 { // W >= -11.983, Y_-7 let x = -z / (X0 + (z - Z0).sqrt()); - Ok((-15.212_958_142_001_646 - + x * (-884_954.687_981_689_6 - + x * (-5.529_815_437_863_348e9 + x * (-3.093_418_743_531_467e12)))) - / (1. - + x * (72_009.255_525_210_12 - + x * (5.505_900_767_187_58e8 + x * 4.432_489_486_700_034e11)))) + Some( + (-15.212_958_142_001_646 + + x * (-884_954.687_981_689_6 + + x * (-5.529_815_437_863_348e9 + x * (-3.093_418_743_531_467e12)))) + / (1. + + x * (72_009.255_525_210_12 + + x * (5.505_900_767_187_58e8 + x * 4.432_489_486_700_034e11))), + ) } else if z <= -1.096_244_452_641_099_5e-19 { // W >= -47.518, V_-8 let u = (-z).ln(); - Ok((-0.032_401_163_177_791_084 - + u * (2.028_194_214_474_250_4 - + u * (-0.527_524_312_425_927_1 + u * 0.017_340_294_772_717_584))) - / (1. - + u * (-0.450_042_744_438_917_4 - + u * (0.017_154_705_753_566_295 + u * (-5.243_819_620_271_836e-7))))) + Some( + (-0.032_401_163_177_791_084 + + u * (2.028_194_214_474_250_4 + + u * (-0.527_524_312_425_927_1 + u * 0.017_340_294_772_717_584))) + / (1. + + u * (-0.450_042_744_438_917_4 + + u * (0.017_154_705_753_566_295 + u * (-5.243_819_620_271_836e-7)))), + ) } else if z <= -2.509_609_929_994_59e-136 { // W >= -317.993, V_-9 let u = (-z).ln(); - Ok((-1.441_124_659_581_209_7 - + u * (1.281_926_963_998_047_7 - + u * (-0.074_979_356_113_812_33 + u * 0.000_476_363_091_620_691_5))) - / (1. - + u * (-0.072_000_873_723_868_65 - + u * (0.000_475_489_329_895_970_3 + u * (-4.171_497_924_754_684e-10))))) + Some( + (-1.441_124_659_581_209_7 + + u * (1.281_926_963_998_047_7 + + u * (-0.074_979_356_113_812_33 + u * 0.000_476_363_091_620_691_5))) + / (1. + + u * (-0.072_000_873_723_868_65 + + u * (0.000_475_489_329_895_970_3 + u * (-4.171_497_924_754_684e-10)))), + ) } else if z < 0.0 { // V_-10 let u = (-z).ln(); - Ok((-3.310_876_091_171_045 - + u * (1.050_067_880_993_517_6 - + u * (-0.008_236_749_582_134_32 + u * 5.528_956_159_491_019e-6))) - / (1. - + u * (-0.008_189_272_743_331_552 - + u * (5.528_007_600_971_195e-6 + u * (-3.922_277_308_457_406_3e-14))))) + Some( + (-3.310_876_091_171_045 + + u * (1.050_067_880_993_517_6 + + u * (-0.008_236_749_582_134_32 + u * 5.528_956_159_491_019e-6))) + / (1. + + u * (-0.008_189_272_743_331_552 + + u * (5.528_007_600_971_195e-6 + u * (-3.922_277_308_457_406_3e-14)))), + ) } else { - Err(LambertWm1Error::new( - LambertWm1ErrorReason::PositiveArgument, - )) + None } } diff --git a/src/lib.rs b/src/lib.rs index 9695ab3..25be34a 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -4,41 +4,36 @@ //! //! This method uses a piecewise minimax rational approximation of the function. //! -//! This crate is a Rust port of the original Fortran 90 code by T. Fukushima. -//! //! ## Examples //! //! Evaluate the principal branch of the Lambert W function to 50 bits of accuracy: #![cfg_attr( - feature = "50", + feature = "50bits", doc = r##" ``` -# use lambert_w::LambertW0Error; use lambert_w::accurate::lambert_w_0; use core::f64::consts::PI; use approx::assert_abs_diff_eq; -let w = lambert_w_0(PI)?; +let w = lambert_w_0(PI).unwrap(); + assert_abs_diff_eq!(w, 1.0736581947961492); -# Ok::<(), LambertW0Error>(()) ``` "## )] //! //! or to only 24 bits of accuracy, but with faster execution time: #![cfg_attr( - feature = "24", + feature = "24bits", doc = r##" ``` -# use lambert_w::LambertW0Error; use lambert_w::fast::lambert_w_0; use core::f64::consts::PI; use approx::assert_abs_diff_eq; -let w = lambert_w_0(PI)?; +let w = lambert_w_0(PI).unwrap(); assert_abs_diff_eq!(w, 1.0736581947961492, epsilon = 1e-7); -# Ok::<(), LambertW0Error>(()) ``` "## )] @@ -52,42 +47,41 @@ assert_abs_diff_eq!(w, 1.0736581947961492, epsilon = 1e-7); //! //! You can disable one of these feature flags to potentially save a little bit of binary size. //! -//! `50` *(enabled by default)*: enables the function versions with 50 bits of accuracy. +//! `50bits` *(enabled by default)*: enables the more accurate function versions with 50 bits of accuracy. //! -//! `24` *(enabled by default)*: enables the function versions with 24 bits of accuracy. +//! `24bits` *(enabled by default)*: enables the faster function versions with 24 bits of accuracy. //! //! It is a compile error to disable both features. #![cfg_attr(docsrs, feature(doc_auto_cfg))] -#[cfg(not(any(feature = "50", feature = "24")))] -compile_error!("one or both of the '24' and '50' features must be enabled"); +#[cfg(not(any(feature = "50bits", feature = "24bits")))] +compile_error!("one or both of the '24bits' and '50bits' features must be enabled"); -#[cfg(feature = "50")] +#[cfg(feature = "50bits")] pub mod accurate; -mod error; -#[cfg(feature = "24")] +#[cfg(feature = "24bits")] pub mod fast; -pub use error::{LambertW0Error, LambertWm1Error, LambertWm1ErrorReason}; - // -1/e const Z0: f64 = -0.367_879_441_171_442_33; // 1/sqrt(e) const X0: f64 = 0.606_530_659_712_633_4; -#[cfg(all(test, any(feature = "24", feature = "50")))] +#[cfg(all(test, any(feature = "24bits", feature = "50bits")))] mod tets { - #[cfg(feature = "50")] + #[cfg(feature = "50bits")] use super::accurate::{lambert_w_0 as lambert_w_0_50, lambert_w_m1 as lambert_w_m1_50}; - #[cfg(feature = "24")] + #[cfg(feature = "24bits")] use super::fast::{lambert_w_0 as lambert_w_0_24, lambert_w_m1 as lambert_w_m1_24}; use approx::assert_abs_diff_eq; + use core::f64::consts::E; - #[cfg(feature = "50")] + #[cfg(feature = "50bits")] #[test] fn test_lambert_w_0_50() { + assert_eq!(lambert_w_0_50(-1.0 / E - f64::EPSILON), None); assert_abs_diff_eq!( lambert_w_0_50(-2.678794411714424e-01).unwrap(), -3.993824525397807e-01 @@ -207,9 +201,10 @@ mod tets { ); } - #[cfg(feature = "24")] + #[cfg(feature = "24bits")] #[test] fn test_lambert_w_0_24() { + assert_eq!(lambert_w_0_24(-1.0 / E - f64::EPSILON), None); assert_abs_diff_eq!( lambert_w_0_24(-2.678794411714424e-01).unwrap(), -3.993824525397807e-01, @@ -342,9 +337,10 @@ mod tets { ); } - #[cfg(feature = "50")] + #[cfg(feature = "50bits")] #[test] fn test_lambert_w_m1_50() { + assert_eq!(lambert_w_m1_50(-1.0 / E - f64::EPSILON), None); assert_abs_diff_eq!( lambert_w_m1_50(-3.578794411714423e-01).unwrap(), -1.253493791367214, @@ -407,11 +403,13 @@ mod tets { lambert_w_m1_50(-1.000000000000008e-145).unwrap(), -3.397029099254290e+02 ); + assert_eq!(lambert_w_m1_50(f64::EPSILON), None); } - #[cfg(feature = "24")] + #[cfg(feature = "24bits")] #[test] fn test_lambert_w_m1_24() { + assert_eq!(lambert_w_m1_24(-1.0 / E - f64::EPSILON), None); assert_abs_diff_eq!( lambert_w_m1_24(-3.578794411714423e-01).unwrap(), -1.253493791367214, @@ -477,5 +475,6 @@ mod tets { -3.397029099254290e+02, epsilon = 1e-4 ); + assert_eq!(lambert_w_m1_24(f64::EPSILON), None); } }