From e27b98757ca0232fc0d2b35cb8dcba792e9387a3 Mon Sep 17 00:00:00 2001 From: Soumya Date: Thu, 19 Sep 2024 17:18:42 +0200 Subject: [PATCH 1/4] Add complex ln_1p() - Added a complex equivalent of the real function ln_1p() --- src/complex_float.rs | 6 ++++++ src/lib.rs | 6 ++++++ 2 files changed, 12 insertions(+) diff --git a/src/complex_float.rs b/src/complex_float.rs index 873fe73..d1c793f 100644 --- a/src/complex_float.rs +++ b/src/complex_float.rs @@ -144,6 +144,10 @@ pub trait ComplexFloat: Num + NumCast + Copy + Neg + private::Sea /// /// Formula: `a+bi -> a-bi` fn conj(self) -> Self; + + /// Returns ln(1+n) (natural logarithm) more accurately than if the operations + /// were performed separately + fn ln_1p(self) -> Self; } macro_rules! forward { @@ -235,6 +239,7 @@ where Float::acosh(self) -> Self; Float::atanh(self) -> Self; Float::abs(self) -> Self; + Float::ln_1p(self) -> Self; } } @@ -306,6 +311,7 @@ impl ComplexFloat for Complex { Complex::asinh(self) -> Self; Complex::acosh(self) -> Self; Complex::atanh(self) -> Self; + Complex::ln_1p(self) -> Self; } forward_ref! { diff --git a/src/lib.rs b/src/lib.rs index 661b67b..6225cd6 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -624,6 +624,12 @@ impl Complex { pub fn fdiv(self, other: Complex) -> Complex { self * other.finv() } + + #[inline] + pub fn ln_1p(self) -> Self { + let complex_num = Self::one() + self; + complex_num.ln() + } } #[cfg(any(feature = "std", feature = "libm"))] From 7d3927cde475a7065f14701d60c33c320fb48719 Mon Sep 17 00:00:00 2001 From: Soumya Date: Thu, 19 Sep 2024 18:03:51 +0200 Subject: [PATCH 2/4] Added doc for ln_1p - Added doc with formula and test for ln_1p() --- src/complex_float.rs | 8 ++++++-- src/lib.rs | 16 ++++++++++++++++ 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/src/complex_float.rs b/src/complex_float.rs index d1c793f..6dc00dc 100644 --- a/src/complex_float.rs +++ b/src/complex_float.rs @@ -145,8 +145,12 @@ pub trait ComplexFloat: Num + NumCast + Copy + Neg + private::Sea /// Formula: `a+bi -> a-bi` fn conj(self) -> Self; - /// Returns ln(1+n) (natural logarithm) more accurately than if the operations - /// were performed separately + /// Returns ln(1+n) (natural logarithm) more accurately + /// than if the operations were performed separately + /// + /// Formula: ln(1+z) + /// + /// where z = a+bi fn ln_1p(self) -> Self; } diff --git a/src/lib.rs b/src/lib.rs index 6225cd6..0cf634d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -625,6 +625,22 @@ impl Complex { self * other.finv() } + /// The ln_1p() function computes the natural logarithm + /// of (1 + z), which is particularly useful for + /// small values of `z` to avoid loss of precision + /// + /// # Examples + /// + /// ```rust + /// use num_complex::Complex64; + /// use num_complex::ComplexFloat; + /// + /// let a = Complex64::new(2.0e-16, 3.0e-16); + /// + /// let approx_val = a.ln_1p(); + /// let expected_val = Complex64::new(2.24e-16, 2.99e-16); + /// assert!((approx_val - expected_val).norm() < 1e-17); + /// ``` #[inline] pub fn ln_1p(self) -> Self { let complex_num = Self::one() + self; From 29937ed02e3a35d50b538e22c77284c0508da366 Mon Sep 17 00:00:00 2001 From: Soumya Date: Thu, 19 Sep 2024 18:38:03 +0200 Subject: [PATCH 3/4] Add complex exp_m1() - Added a complex equivalent of the real function exp_m1() - Added doc with formula and test for exp_m1 --- src/complex_float.rs | 10 ++++++++++ src/lib.rs | 28 ++++++++++++++++++++++++---- 2 files changed, 34 insertions(+), 4 deletions(-) diff --git a/src/complex_float.rs b/src/complex_float.rs index 6dc00dc..5b1c143 100644 --- a/src/complex_float.rs +++ b/src/complex_float.rs @@ -152,6 +152,14 @@ pub trait ComplexFloat: Num + NumCast + Copy + Neg + private::Sea /// /// where z = a+bi fn ln_1p(self) -> Self; + + /// Returns e^(self) - 1 in a way that is accurate + /// even if the number is close to zero + /// + /// Formaula: e^(z) - 1 + /// + /// where z = a+bi + fn exp_m1(self) -> Self; } macro_rules! forward { @@ -244,6 +252,7 @@ where Float::atanh(self) -> Self; Float::abs(self) -> Self; Float::ln_1p(self) -> Self; + Float::exp_m1(self) -> Self; } } @@ -316,6 +325,7 @@ impl ComplexFloat for Complex { Complex::acosh(self) -> Self; Complex::atanh(self) -> Self; Complex::ln_1p(self) -> Self; + Complex::exp_m1(self) -> Self; } forward_ref! { diff --git a/src/lib.rs b/src/lib.rs index 0cf634d..4bdfd8e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -638,13 +638,33 @@ impl Complex { /// let a = Complex64::new(2.0e-16, 3.0e-16); /// /// let approx_val = a.ln_1p(); - /// let expected_val = Complex64::new(2.24e-16, 2.99e-16); - /// assert!((approx_val - expected_val).norm() < 1e-17); + /// let expected_val = Complex64::new(2.2204e-16, 2.999e-16); + /// assert!((approx_val - expected_val).norm() < 1e-18); /// ``` #[inline] pub fn ln_1p(self) -> Self { - let complex_num = Self::one() + self; - complex_num.ln() + (Self::one() + self).ln() + } + + /// The exp_m1() function computes the exponential + /// of z minus one (`e^(z) - 1`), which is particularly + /// useful for small values of `z` to avoid loss of precision + /// + /// # Examples + /// + /// ```rust + /// use num_complex::Complex64; + /// use num_complex::ComplexFloat; + /// + /// let a = Complex64::new(2.0e-16, 3.0e-16); + /// + /// let approx_val = a.exp_m1(); + /// let expected_val = Complex64::new(2.2204e-16, 3.0e-16); + /// assert!((approx_val - expected_val).norm() < 1e-18); + /// ``` + #[inline] + pub fn exp_m1(self) -> Self { + self.exp() - Self::one() } } From 9bd0dd43248707ec22d1e42a6f3c0b766cf19157 Mon Sep 17 00:00:00 2001 From: Soumya Date: Thu, 19 Sep 2024 18:52:12 +0200 Subject: [PATCH 4/4] Added ln_1p and exp_m1 - Added tests for them --- src/lib.rs | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/lib.rs b/src/lib.rs index 4bdfd8e..3e3732c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -2474,6 +2474,22 @@ pub(crate) mod test { )); } } + + #[test] + fn test_ln_1p() { + for &c in all_consts.iter() { + // ln_1p(z) = ln(1+z) + assert!(close(c.ln_1p(), (1.0 + c).ln())); + } + } + + #[test] + fn test_exp_m1() { + for &c in all_consts.iter() { + // exp_m1(z) = exp(z) - 1 + assert!(close(c.exp_m1(), (c).exp() - 1.0)); + } + } } // Test both a + b and a += b