From f283059c8d8f8ff4267e4622ab3ff40fed39d0d7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= Date: Tue, 30 Jul 2024 12:07:16 +0200 Subject: [PATCH 1/9] Add Lambert W function --- Cargo.toml | 1 + src/prelude/mod.rs | 7 +++++-- src/prelude/simpler.rs | 30 ++++++++++++++++++++++++++++ src/special/function.rs | 43 +++++++++++++++++++++++++++++++++++++++++ 4 files changed, 79 insertions(+), 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 3ba8b886..4e8de363 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -39,6 +39,7 @@ serde = { version = "1.0", features = ["derive"], optional = true } json = { version = "0.12", optional = true } arrow2 = { version = "0.18", features = ["io_parquet", "io_parquet_compression"], optional = true } num-complex = { version = "0.4", optional = true } +lambert_w = {version = "0.3.0", default-features = false, features = ["24bits", "50bits"]} [package.metadata.docs.rs] rustdoc-args = [ "--html-in-header", "katex-header.html", "--cfg", "docsrs"] diff --git a/src/prelude/mod.rs b/src/prelude/mod.rs index 7836b46d..5d37c9a0 100644 --- a/src/prelude/mod.rs +++ b/src/prelude/mod.rs @@ -193,7 +193,10 @@ pub use crate::util::{api::*, low_level::*, non_macro::*, print::*, useful::*, w pub use crate::statistics::{dist::*, ops::*, rand::*, stat::*}; #[allow(unused_imports)] -pub use crate::special::function::*; +pub use crate::special::function::{ + beta, erf, erfc, gamma, gaussian, inc_beta, inc_gamma, inv_erf, inv_erfc, inv_inc_gamma, + inv_inv_beta, ln_gamma, phi, poch, +}; #[allow(unused_imports)] pub use crate::numerical::{ @@ -206,7 +209,7 @@ pub use crate::numerical::{ utils::*, }; -pub use simpler::{eigen, integrate, chebyshev_polynomial, cubic_hermite_spline}; +pub use simpler::{eigen, integrate, chebyshev_polynomial, cubic_hermite_spline, lambert_w0, lambert_w_m1}; #[allow(unused_imports)] pub use crate::statistics::stat::Metric::*; diff --git a/src/prelude/simpler.rs b/src/prelude/simpler.rs index 2b1851f9..435cfd5d 100644 --- a/src/prelude/simpler.rs +++ b/src/prelude/simpler.rs @@ -152,6 +152,36 @@ pub fn cubic_hermite_spline(node_x: &[f64], node_y: &[f64]) -> anyhow::Result f64 { + lambert_w0_flex(z, LambertWAccuracyMode::Precise) +} + +/// The secondary branch of the Lambert W function, W_-1(`z`). +/// +/// Returns [`NAN`](f64::NAN) if the given input is positive or smaller than -1/e (≈ -0.36787944117144233). +/// +/// Wrapper of `lambert_w_m1` function of `lambert_w` crate. +/// +/// # Reference +/// +/// [Toshio Fukushima, Precise and fast computation of Lambert W function by piecewise minimax rational function approximation with variable transformation](https://www.researchgate.net/publication/346309410_Precise_and_fast_computation_of_Lambert_W_function_by_piecewise_minimax_rational_function_approximation_with_variable_transformation) +pub fn lambert_w_m1(z: f64) -> f64 { + lambert_w_m1_flex(z, LambertWAccuracyMode::Precise) +} + /// Simple handle parquet #[cfg(feature="parquet")] pub trait SimpleParquet: Sized { diff --git a/src/special/function.rs b/src/special/function.rs index 1fe0787b..01ef92e0 100644 --- a/src/special/function.rs +++ b/src/special/function.rs @@ -126,3 +126,46 @@ pub fn phi(x: f64) -> f64 { // special_fun::unsafe_cephes_double::hyp2f1(a, b, c, x) // } // } + +/// The principal branch of the Lambert W function, W_0(`z`). +/// +/// Returns [`NAN`](f64::NAN) if the given input is smaller than -1/e (≈ -0.36787944117144233). +/// +/// Wrapper of the `lambert_w_0` and `sp_lambert_w_0` functions of the `lambert_w` crate. +/// +/// # Reference +/// +/// [Toshio Fukushima, Precise and fast computation of Lambert W function by piecewise minimax rational function approximation with variable transformation](https://www.researchgate.net/publication/346309410_Precise_and_fast_computation_of_Lambert_W_function_by_piecewise_minimax_rational_function_approximation_with_variable_transformation) +pub fn lambert_w0(z: f64, mode: LambertWAccuracyMode) -> f64 { + match mode { + LambertWAccuracyMode::Precise => lambert_w::lambert_w_0(z), + LambertWAccuracyMode::Simple => lambert_w::sp_lambert_w_0(z), + } + .unwrap_or(f64::NAN) +} + +/// The secondary branch of the Lambert W function, W_-1(`z`). +/// +/// Returns [`NAN`](f64::NAN) if the given input is positive or smaller than -1/e (≈ -0.36787944117144233). +/// +/// Wrapper of the `lambert_w_m1` and `sp_lambert_w_m1` functions of the `lambert_w` crate. +/// +/// # Reference +/// +/// [Toshio Fukushima, Precise and fast computation of Lambert W function by piecewise minimax rational function approximation with variable transformation](https://www.researchgate.net/publication/346309410_Precise_and_fast_computation_of_Lambert_W_function_by_piecewise_minimax_rational_function_approximation_with_variable_transformation) +pub fn lambert_wm1(z: f64, mode: LambertWAccuracyMode) -> f64 { + match mode { + LambertWAccuracyMode::Precise => lambert_w::lambert_w_m1(z), + LambertWAccuracyMode::Simple => lambert_w::sp_lambert_w_m1(z), + } + .unwrap_or(f64::NAN) +} + +/// Decides the accuracy mode of the Lambert W functions. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum LambertWAccuracyMode { + /// Faster, 24 bits of accuracy. + Simple, + /// Slower, 50 bits of accuracy. + Precise, +} \ No newline at end of file From 75f4604e57fa8203bb23afd54e13767c8130f320 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= Date: Tue, 30 Jul 2024 12:09:14 +0200 Subject: [PATCH 2/9] Unify formatting of the lambert_w import with all others --- Cargo.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index 4e8de363..c31c407b 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -39,7 +39,7 @@ serde = { version = "1.0", features = ["derive"], optional = true } json = { version = "0.12", optional = true } arrow2 = { version = "0.18", features = ["io_parquet", "io_parquet_compression"], optional = true } num-complex = { version = "0.4", optional = true } -lambert_w = {version = "0.3.0", default-features = false, features = ["24bits", "50bits"]} +lambert_w = { version = "0.3.0", default-features = false, features = ["24bits", "50bits"] } [package.metadata.docs.rs] rustdoc-args = [ "--html-in-header", "katex-header.html", "--cfg", "docsrs"] From 15796f6f5b89451c2ff989709bd011335074c834 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= Date: Tue, 30 Jul 2024 12:15:20 +0200 Subject: [PATCH 3/9] Add documentation information about the accuracy mode --- src/special/function.rs | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/special/function.rs b/src/special/function.rs index 01ef92e0..cebbfd89 100644 --- a/src/special/function.rs +++ b/src/special/function.rs @@ -131,6 +131,9 @@ pub fn phi(x: f64) -> f64 { /// /// Returns [`NAN`](f64::NAN) if the given input is smaller than -1/e (≈ -0.36787944117144233). /// +/// Use [`LambertWAccuracyMode::Precise`] for 50 bits of accuracy and the [`Simple`](LambertWAccuracyMode::Simple) mode +/// for only 24 bits, but with faster execution time. +/// /// Wrapper of the `lambert_w_0` and `sp_lambert_w_0` functions of the `lambert_w` crate. /// /// # Reference @@ -148,8 +151,11 @@ pub fn lambert_w0(z: f64, mode: LambertWAccuracyMode) -> f64 { /// /// Returns [`NAN`](f64::NAN) if the given input is positive or smaller than -1/e (≈ -0.36787944117144233). /// +/// Use [`LambertWAccuracyMode::Precise`] for 50 bits of accuracy and the [`Simple`](LambertWAccuracyMode::Simple) mode +/// for only 24 bits, but with faster execution time. +/// /// Wrapper of the `lambert_w_m1` and `sp_lambert_w_m1` functions of the `lambert_w` crate. -/// +/// /// # Reference /// /// [Toshio Fukushima, Precise and fast computation of Lambert W function by piecewise minimax rational function approximation with variable transformation](https://www.researchgate.net/publication/346309410_Precise_and_fast_computation_of_Lambert_W_function_by_piecewise_minimax_rational_function_approximation_with_variable_transformation) From 17c14dc780255690ec923225459ecb79bf5d5567 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= Date: Tue, 30 Jul 2024 12:15:59 +0200 Subject: [PATCH 4/9] Add accuracy information to the simple version in prelude --- src/prelude/simpler.rs | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/prelude/simpler.rs b/src/prelude/simpler.rs index 435cfd5d..d8cb0cb6 100644 --- a/src/prelude/simpler.rs +++ b/src/prelude/simpler.rs @@ -160,6 +160,8 @@ use crate::special::function::{ /// /// Returns [`NAN`](f64::NAN) if the given input is smaller than -1/e (≈ -0.36787944117144233). /// +/// Accurate to 50 bits. +/// /// Wrapper of `lambert_w_0` function of `lambert_w` crate. /// /// # Reference @@ -173,6 +175,8 @@ pub fn lambert_w0(z: f64) -> f64 { /// /// Returns [`NAN`](f64::NAN) if the given input is positive or smaller than -1/e (≈ -0.36787944117144233). /// +/// Accurate to 50 bits. +/// /// Wrapper of `lambert_w_m1` function of `lambert_w` crate. /// /// # Reference From aba9cfe038900eb1bd8be10adb6df0f1f4f1b2cc Mon Sep 17 00:00:00 2001 From: Tae-Geun Kim Date: Tue, 30 Jul 2024 20:22:41 +0900 Subject: [PATCH 5/9] Fix typo of comments [`LambertWAccuracyMode::Precise`] -> [`Precise](LambertWAccuracyMode::Precise) --- src/special/function.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/special/function.rs b/src/special/function.rs index cebbfd89..1850aaff 100644 --- a/src/special/function.rs +++ b/src/special/function.rs @@ -131,7 +131,7 @@ pub fn phi(x: f64) -> f64 { /// /// Returns [`NAN`](f64::NAN) if the given input is smaller than -1/e (≈ -0.36787944117144233). /// -/// Use [`LambertWAccuracyMode::Precise`] for 50 bits of accuracy and the [`Simple`](LambertWAccuracyMode::Simple) mode +/// Use [`Precise`](LambertWAccuracyMode::Precise) for 50 bits of accuracy and the [`Simple`](LambertWAccuracyMode::Simple) mode /// for only 24 bits, but with faster execution time. /// /// Wrapper of the `lambert_w_0` and `sp_lambert_w_0` functions of the `lambert_w` crate. @@ -151,7 +151,7 @@ pub fn lambert_w0(z: f64, mode: LambertWAccuracyMode) -> f64 { /// /// Returns [`NAN`](f64::NAN) if the given input is positive or smaller than -1/e (≈ -0.36787944117144233). /// -/// Use [`LambertWAccuracyMode::Precise`] for 50 bits of accuracy and the [`Simple`](LambertWAccuracyMode::Simple) mode +/// Use [`Precise`](LambertWAccuracyMode::Precise) for 50 bits of accuracy and the [`Simple`](LambertWAccuracyMode::Simple) mode /// for only 24 bits, but with faster execution time. /// /// Wrapper of the `lambert_w_m1` and `sp_lambert_w_m1` functions of the `lambert_w` crate. @@ -174,4 +174,4 @@ pub enum LambertWAccuracyMode { Simple, /// Slower, 50 bits of accuracy. Precise, -} \ No newline at end of file +} From aafaa3ecd5ac13f48e28fdc4d7ada1820d8e5ace Mon Sep 17 00:00:00 2001 From: Axect Date: Tue, 30 Jul 2024 20:34:14 +0900 Subject: [PATCH 6/9] FIX: Fix typo for inverse incomplete beta function - `inv_inv_beta` -> `inv_inc_beta` --- src/special/function.rs | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/src/special/function.rs b/src/special/function.rs index 1850aaff..e81c86ac 100644 --- a/src/special/function.rs +++ b/src/special/function.rs @@ -107,7 +107,7 @@ pub fn inc_beta(a: f64, b: f64, x: f64) -> f64 { /// Inverse regularized incomplete beta function /// /// Wrapper of `invbetai` function of `puruspe` crate -pub fn inv_inv_beta(p: f64, a: f64, b: f64) -> f64 { +pub fn inv_inc_beta(p: f64, a: f64, b: f64) -> f64 { puruspe::invbetai(p, a, b) } @@ -118,15 +118,6 @@ pub fn phi(x: f64) -> f64 { 0.5 * (1f64 + erf(x / 2f64.sqrt())) } -// /// Hypergeometric function 2F1 -// /// -// /// Wrapper of `hyp2f1` function of `special-fun` crate -// pub fn hyp2f1(a: f64, b: f64, c: f64, x: f64) -> f64 { -// unsafe { -// special_fun::unsafe_cephes_double::hyp2f1(a, b, c, x) -// } -// } - /// The principal branch of the Lambert W function, W_0(`z`). /// /// Returns [`NAN`](f64::NAN) if the given input is smaller than -1/e (≈ -0.36787944117144233). From 6c85064b317cd2490ebf4145aa19f9783066682f Mon Sep 17 00:00:00 2001 From: Axect Date: Tue, 30 Jul 2024 20:41:45 +0900 Subject: [PATCH 7/9] FIX: Fix typo for inv_inc_beta in prelude --- src/prelude/mod.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/prelude/mod.rs b/src/prelude/mod.rs index 5d37c9a0..ce5dd035 100644 --- a/src/prelude/mod.rs +++ b/src/prelude/mod.rs @@ -195,7 +195,7 @@ pub use crate::statistics::{dist::*, ops::*, rand::*, stat::*}; #[allow(unused_imports)] pub use crate::special::function::{ beta, erf, erfc, gamma, gaussian, inc_beta, inc_gamma, inv_erf, inv_erfc, inv_inc_gamma, - inv_inv_beta, ln_gamma, phi, poch, + inv_inc_beta, ln_gamma, phi, poch, }; #[allow(unused_imports)] From 82d0a52f47930d217c7dd8ef4c535ba6f0de7150 Mon Sep 17 00:00:00 2001 From: Axect Date: Tue, 30 Jul 2024 20:42:06 +0900 Subject: [PATCH 8/9] IMPL: Implement very simple test for lambert_w --- tests/special.rs | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 tests/special.rs diff --git a/tests/special.rs b/tests/special.rs new file mode 100644 index 00000000..3954b4bf --- /dev/null +++ b/tests/special.rs @@ -0,0 +1,7 @@ +use peroxide::fuga::{*, LambertWAccuracyMode::*}; + +#[test] +fn lambert_w_test() { + assert_eq!(lambert_w0(1.0, Precise), 0.567143290409784); + assert!(nearly_eq(lambert_w0(1.0, Simple), 0.567143290409784)); +} From 9247f21eca2e8d84f32fea530a627604b2fe111a Mon Sep 17 00:00:00 2001 From: Axect Date: Tue, 30 Jul 2024 20:54:02 +0900 Subject: [PATCH 9/9] RLSE: Ver 0.37.8 - Integrate with Lambert W crate (#63) --- Cargo.toml | 2 +- RELEASES.md | 23 +++++++++++++++++++++++ 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index c31c407b..26bc9d99 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "peroxide" -version = "0.37.7" +version = "0.37.8" authors = ["axect "] edition = "2018" description = "Rust comprehensive scientific computation library contains linear algebra, numerical analysis, statistics and machine learning tools with farmiliar syntax" diff --git a/RELEASES.md b/RELEASES.md index 4ef2435d..8bc5a2dd 100644 --- a/RELEASES.md +++ b/RELEASES.md @@ -1,3 +1,26 @@ +# Release 0.37.8 (2024-07-30) + +- Integrate with [lambert_w](https://crates.io/crates/lambert_w) crate ([#63](https://github.com/Axect/Peroxide/pull/63)) (Thanks to [@JSorngard](https://github.com/JSorngard)) + - Write flexible wrapper for [lambert_w](https://crates.io/crates/lambert_w) + ```rust + pub enum LambertWAccuracyMode { + Simple, // Faster, 24 bits of accuracy + Precise, // Slower, 50 bits of accuracy + } + + pub fn lambert_w0(z: f64, mode: LambertWAccuracyMode) -> f64; + pub fn lambert_wm1(z: f64, mode: LambertWAccuracyMode) -> f64; + ``` + + - Write default Lambert W function for `prelude` (Precise as default) + ```rust + use peroxide::prelude::*; + + fn main() { + lambert_w0(1.0).print(); // Same as fuga::lambert_w0(1.0, LambertWAccuracyMode::Simple) + } + ``` + # Release 0.37.7 (2024-07-05) - Bump `pyo3` dependency to `0.22`