Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Return Option instead of Result for a cleaner API #2

Merged
merged 10 commits into from
Jul 28, 2024
Merged

Return Option instead of Result for a cleaner API #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
ba4c55c
Remove all results
JSorngard Jul 28, 2024
4e05d93
We don't have to be that explicit
JSorngard Jul 28, 2024
19d5afa
Less references
JSorngard Jul 28, 2024
7a2b2f7
Test border values
JSorngard Jul 28, 2024
142dd50
Use the right function
JSorngard Jul 28, 2024
e2160b4
More descriptive feature flag names
JSorngard Jul 28, 2024
b21df2e
More description in Cargo.toml
JSorngard Jul 28, 2024
d1d8788
Use correct feature names in CI
JSorngard Jul 28, 2024
6597076
Add speed-accuracy trade-off info to readme
JSorngard Jul 28, 2024
c1e434c
Update README and docstring
JSorngard Jul 28, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .github/workflows/rust.yml
View file
Open in desktop
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
View file
Open in desktop

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
View file
Open in desktop
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 <[email protected]>"]
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
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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