treat corner cases (inf, nan etc.) as in std::ln

Expander · Nov 25, 2024 · 9999dea · 9999dea
1 parent 37e111b
commit 9999dea
Showing 1 changed file with 26 additions and 7 deletions.
diff --git a/src/cln.rs b/src/cln.rs
@@ -28,7 +28,7 @@ impl<T: Float + Pi<T>> CLn<Complex<T>> for Complex<T> {
     fn cln(&self) -> Complex<T> {
         if self.im == T::zero() && self.re > T::zero() {
             Complex::new(self.re.ln(), T::zero())
-        } else if self.im == T::zero() {
+        } else if self.im == T::zero() && self.re < T::zero() {
             Complex::new((-self.re).ln(), T::pi())
         } else {
             self.ln()
@@ -38,11 +38,30 @@ impl<T: Float + Pi<T>> CLn<Complex<T>> for Complex<T> {
 
 #[test]
 fn test_cln() {
-    let x32 = Complex::new(1.0_f32, 0.0_f32);
-    let x64 = Complex::new(1.0_f64, 0.0_f64);
+    // test positive zero
+    let cnegone32 = Complex::new(-1.0_f32, 0.0_f32);
+    let cnegone64 = Complex::new(-1.0_f64, 0.0_f64);
 
-    assert!(((-x32).cln() - Complex::new(0.0_f32,  std::f32::consts::PI)).norm() < 4.0_f32*std::f32::EPSILON);
-    assert!(((-x32).ln()  - Complex::new(0.0_f32, -std::f32::consts::PI)).norm() < 4.0_f32*std::f32::EPSILON);
-    assert!((-x64).cln() == Complex::new(0.0_f64,  std::f64::consts::PI));
-    assert!((-x64).ln()  == Complex::new(0.0_f64, -std::f64::consts::PI));
+    assert!((cnegone32.cln() - Complex::new(0.0_f32, std::f32::consts::PI)).norm() < 4.0_f32*std::f32::EPSILON);
+    assert!((cnegone32.ln()  - Complex::new(0.0_f32, std::f32::consts::PI)).norm() < 4.0_f32*std::f32::EPSILON);
+    assert!(cnegone64.cln() == Complex::new(0.0_f64, std::f64::consts::PI));
+    assert!(cnegone64.ln()  == Complex::new(0.0_f64, std::f64::consts::PI));
+
+    // test negative zero
+    let cone32 = Complex::new(1.0_f32, 0.0_f32);
+    let cone64 = Complex::new(1.0_f64, 0.0_f64);
+
+    assert!(((-cone32).cln() - Complex::new(0.0_f32,  std::f32::consts::PI)).norm() < 4.0_f32*std::f32::EPSILON);
+    assert!(((-cone32).ln()  - Complex::new(0.0_f32, -std::f32::consts::PI)).norm() < 4.0_f32*std::f32::EPSILON);
+    assert!((-cone64).cln() == Complex::new(0.0_f64,  std::f64::consts::PI));
+    assert!((-cone64).ln()  == Complex::new(0.0_f64, -std::f64::consts::PI));
+
+    // test zero input
+    let rzero32 = 0.0_f32;
+    let rzero64 = 0.0_f64;
+
+    assert!(Complex::new( rzero32, rzero32).cln() == Complex::new(std::f32::NEG_INFINITY, 0.0));
+    assert!(Complex::new(-rzero32, rzero32).cln() == Complex::new(std::f32::NEG_INFINITY, std::f32::consts::PI));
+    assert!(Complex::new( rzero64, rzero64).cln() == Complex::new(std::f64::NEG_INFINITY, 0.0));
+    assert!(Complex::new(-rzero64, rzero64).cln() == Complex::new(std::f64::NEG_INFINITY, std::f64::consts::PI));
 }