Skip to content

Commit

Permalink
Return Option instead of Result for a cleaner API (#2)
Browse files Browse the repository at this point in the history
* Remove all results

* Test border values

* More descriptive feature flag names

* More description in Cargo.toml

* Add speed-accuracy trade-off info to readme

* Update README and docstring
  • Loading branch information
JSorngard authored Jul 28, 2024
1 parent a605cab commit da48d97
Show file tree
Hide file tree
Showing 12 changed files with 328 additions and 384 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/rust.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion Cargo.lock

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

12 changes: 6 additions & 6 deletions Cargo.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
[package]
name = "lambert_w"
version = "0.1.1"
version = "0.2.0"
edition = "2021"
authors = ["Johanna Sörngård <jsorngard@gmail.com>"]
categories = ["mathematics"]
Expand All @@ -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 <https://stackoverflow.com/a/61417700/>.
[package.metadata.docs.rs]
Expand Down
9 changes: 6 additions & 3 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand All @@ -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);
```
```

## 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`.
44 changes: 21 additions & 23 deletions src/accurate/dw0c.rs
Original file line number Diff line number Diff line change
@@ -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
///
Expand All @@ -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<f64, LambertW0Error> {
pub fn dw0c(zc: f64) -> Option<f64> {
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
Expand All @@ -41,7 +39,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -60,7 +58,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -79,7 +77,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -98,7 +96,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -117,7 +115,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -136,7 +134,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -155,7 +153,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -174,7 +172,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -193,7 +191,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -212,7 +210,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -231,7 +229,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -250,7 +248,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -269,7 +267,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -288,7 +286,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -307,7 +305,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -326,7 +324,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -345,7 +343,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand All @@ -364,7 +362,7 @@ pub fn dw0c(zc: f64) -> Result<f64, LambertW0Error> {
} 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
Expand Down
33 changes: 14 additions & 19 deletions src/accurate/dwm1c.rs
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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<f64, LambertWm1Error> {
pub fn dwm1c(z: f64, zc: f64) -> Option<f64> {
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
Expand All @@ -43,7 +40,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result<f64, LambertWm1Error> {
} 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
Expand All @@ -61,7 +58,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result<f64, LambertWm1Error> {
} 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
Expand All @@ -80,7 +77,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result<f64, LambertWm1Error> {
} 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
Expand All @@ -99,7 +96,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result<f64, LambertWm1Error> {
} 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
Expand All @@ -118,7 +115,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result<f64, LambertWm1Error> {
} 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
Expand All @@ -137,7 +134,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result<f64, LambertWm1Error> {
} 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
Expand All @@ -156,7 +153,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result<f64, LambertWm1Error> {
} 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
Expand All @@ -175,7 +172,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result<f64, LambertWm1Error> {
} 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
Expand All @@ -194,7 +191,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result<f64, LambertWm1Error> {
} 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
Expand All @@ -213,7 +210,7 @@ pub fn dwm1c(z: f64, zc: f64) -> Result<f64, LambertWm1Error> {
} 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
Expand All @@ -230,8 +227,6 @@ pub fn dwm1c(z: f64, zc: f64) -> Result<f64, LambertWm1Error> {
+ u * (4.641_976_809_305_971e-15
+ u * (-1.360_871_393_694_260_3e-23)))))))))
} else {
Err(LambertWm1Error::new(
LambertWm1ErrorReason::PositiveArgument,
))
None
}
}
Loading

0 comments on commit da48d97

Please sign in to comment.