Bugfix: Imprecision of complex li2 for small arguments

src/cln.rs

@@ -10,6 +10,7 @@ use num::{Float, FromPrimitive};
 /// the imaginary part is within [-pi, pi].
 pub trait CLn<T> {
     fn cln(&self) -> T;
+    fn cln_1p(&self) -> T;
 }
 
 impl<T: Float + FromPrimitive> CLn<Complex<T>> for Complex<T> {
@@ -22,6 +23,24 @@ impl<T: Float + FromPrimitive> CLn<Complex<T>> for Complex<T> {
             self.ln()
         }
     }
+
+    fn cln_1p(&self) -> Complex<T> {
+        let z = Complex::new(
+            self.re,
+            // convert -0.0 to 0.0
+            if self.im == T::zero() { T::zero() } else { self.im },
+        );
+        let u = Complex::new(T::one() + z.re, z.im);
+        if z.im == T::zero() && z.re > -T::one() {
+            Complex::new(z.re.ln_1p(), T::zero())
+        } else if u == Complex::new(T::one(), T::zero()) {
+            z
+        } else if u.re <= T::zero() {
+            u.ln()
+        } else {
+            u.ln()*(z/Complex::new(u.re - T::one(), u.im))
+        }
+    }
 }
 
 #[test]
@@ -34,3 +53,18 @@ fn test_cln() {
     assert!((-x64).cln() == Complex::new(0.0_f64,  std::f64::consts::PI));
     assert!((-x64).ln()  == Complex::new(0.0_f64, -std::f64::consts::PI));
 }
+
+#[test]
+fn test_cln_1p() {
+    let x32 = Complex::new(0.0_f32, 0.0_f32);
+    let x64 = Complex::new(0.0_f64, 0.0_f64);
+
+    assert!(x32.cln_1p().norm() == 0.0_f32);
+    assert!(x64.cln_1p().norm() == 0.0_f64);
+
+    let o32 = Complex::new(2.0_f32, 0.0_f32);
+    let o64 = Complex::new(2.0_f64, 0.0_f64);
+
+    assert!(((-o32).cln_1p() - Complex::new(0.0_f32,  std::f32::consts::PI)).norm() < 4.0_f32*std::f32::EPSILON);
+    assert!((-o64).cln_1p() == Complex::new(0.0_f64,  std::f64::consts::PI));
+}

src/li1.rs

@@ -41,6 +41,6 @@ impl Li1<Complex<f64>> for Complex<f64> {
     /// assert!((Complex::new(1.0_f64, 1.0_f64).li1() - Complex::new(0.0_f64, 1.5707963267948966_f64)).norm() < std::f64::EPSILON);
     /// ```
     fn li1(&self) -> Complex<f64> {
-        -(1.0 - self).cln()
+        -(-self).cln_1p()
     }
 }

src/li2.rs

@@ -229,17 +229,17 @@ impl Li2<Complex<f32>> for Complex<f32> {
             } else if rz <= 0.5_f32 {
                 if nz > 1.0_f32 {
                     let l = (-self).cln();
-                    -(-(1.0_f32 - 1.0_f32/self).cln()).approx() - 0.5_f32*l*l - pi*pi/6.0_f32
+                    -(-(-1.0_f32/self).cln_1p()).approx() - 0.5_f32*l*l - pi*pi/6.0_f32
                 } else { // nz <= 1
-                    (-(1.0_f32 - self).cln()).approx()
+                    (-(-self).cln_1p()).approx()
                 }
             } else { // rz > 0.5
                 if nz <= 2.0_f32*rz {
                     let l = -(self).cln();
-                    -l.approx() + l*(1.0_f32 - self).cln() + pi*pi/6.0_f32
+                    -l.approx() + l*(-self).cln_1p() + pi*pi/6.0_f32
                 } else { // nz > 2*rz
                     let l = (-self).cln();
-                    -(-(1.0_f32 - 1.0_f32/self).cln()).approx() - 0.5_f32*l*l - pi*pi/6.0_f32
+                    -(-(-1.0_f32/self).cln_1p()).approx() - 0.5_f32*l*l - pi*pi/6.0_f32
                 }
             }
         }
@@ -279,17 +279,17 @@ impl Li2<Complex<f64>> for Complex<f64> {
             } else if rz <= 0.5_f64 {
                 if nz > 1.0_f64 {
                     let l = (-self).cln();
-                    -(-(1.0_f64 - 1.0_f64/self).cln()).approx() - 0.5_f64*l*l - pi*pi/6.0_f64
+                    -(-(-1.0_f64/self).cln_1p()).approx() - 0.5_f64*l*l - pi*pi/6.0_f64
                 } else { // nz <= 1
-                    (-(1.0_f64 - self).cln()).approx()
+                    (-(-self).cln_1p()).approx()
                 }
             } else { // rz > 0.5
                 if nz <= 2.0_f64*rz {
                     let l = -(self).cln();
-                    -l.approx() + l*(1.0_f64 - self).cln() + pi*pi/6.0_f64
+                    -l.approx() + l*(-self).cln_1p() + pi*pi/6.0_f64
                 } else { // nz > 2*rz
                     let l = (-self).cln();
-                    -(-(1.0_f64 - 1.0_f64/self).cln()).approx() - 0.5_f64*l*l - pi*pi/6.0_f64
+                    -(-(-1.0_f64/self).cln_1p()).approx() - 0.5_f64*l*l - pi*pi/6.0_f64
                 }
             }
         }

src/li3.rs

@@ -151,11 +151,11 @@ impl Li3<Complex<f64>> for Complex<f64> {
                 u8*(cs[2] + u2*cs[3] + u4*(cs[4] + u2*cs[5])) +
                 u8*u8*cs[6]
             } else if nz <= 1.0 {
-                cli3_unit_circle(-(1.0 - self).cln())
+                cli3_unit_circle(-(-self).cln_1p())
             } else { // nz > 1
                 let arg = if pz > 0.0 { pz - pi } else { pz + pi };
                 let lmz = Complex::new(lnz, arg); // (-self).cln()
-                cli3_unit_circle(-(1.0 - 1.0/self).cln()) - lmz*(lmz*lmz/6.0 + z2)
+                cli3_unit_circle(-(-1.0/self).cln_1p()) - lmz*(1.0/6.0*lmz*lmz + z2)
             }
         }
     }

src/li4.rs

@@ -192,13 +192,13 @@ impl Li4<Complex<f64>> for Complex<f64> {
                    u8*u8*cs[6]
                 )
             } else if nz <= 1.0 {
-                cli4_unit_circle(-(1.0 - self).cln())
+                cli4_unit_circle(-(-self).cln_1p())
             } else { // nz > 1.0
                 let pi4  = pi2*pi2;
                 let arg = if pz > 0.0 { pz - pi } else { pz + pi };
                 let lmz = Complex::new(lnz, arg); // (-self).cln()
                 let lmz2 = lmz*lmz;
-                -cli4_unit_circle(-(1.0 - 1.0/self).cln()) + 1.0/360.0*(-7.0*pi4 + lmz2*(-30.0*pi2 - 15.0*lmz2))
+                -cli4_unit_circle(-(-1.0/self).cln_1p()) + 1.0/360.0*(-7.0*pi4 + lmz2*(-30.0*pi2 - 15.0*lmz2))
             }
         }
     }

src/li5.rs

@@ -59,13 +59,13 @@ impl Li5<Complex<f64>> for Complex<f64> {
                 u2 * (cs[4] +
                 u2 * (cs[5]))))))))
             } else if nz <= 1.0 {
-                cli5_unit_circle(-(1.0 - self).cln())
+                cli5_unit_circle(-(-self).cln_1p())
             } else { // nz > 1.0
                 let pi4  = pi2*pi2;
                 let arg = if pz > 0.0 { pz - pi } else { pz + pi };
                 let lmz = Complex::new(lnz, arg); // (-self).cln()
                 let lmz2 = lmz*lmz;
-                cli5_unit_circle(-(1.0 - 1.0/self).cln()) - 1.0/360.0*lmz*(7.0*pi4 + lmz2*(10.0*pi2 + 3.0*lmz2))
+                cli5_unit_circle(-(-1.0/self).cln_1p()) - 1.0/360.0*lmz*(7.0*pi4 + lmz2*(10.0*pi2 + 3.0*lmz2))
             }
         }
     }

src/li6.rs

@@ -60,14 +60,14 @@ impl Li6<Complex<f64>> for Complex<f64> {
                 u2 * (cs[3] +
                 u2 * (cs[4]))))))))
             } else if nz <= 1.0 {
-                cli6_unit_circle(-(1.0 - self).cln())
+                cli6_unit_circle(-(-self).cln_1p())
             } else { // nz > 1.0
                 let pi4 = pi2*pi2;
                 let pi6 = pi2*pi4;
                 let arg = if pz > 0.0 { pz - pi } else { pz + pi };
                 let lmz = Complex::new(lnz, arg); // (-self).cln()
                 let lmz2 = lmz*lmz;
-                -cli6_unit_circle(-(1.0 - 1.0/self).cln()) - 31.0*pi6/15120.0 + lmz2*(-7.0/720.0*pi4 + lmz2*(-1.0/144.0*pi2 - 1.0/720.0*lmz2))
+                -cli6_unit_circle(-(-1.0/self).cln_1p()) - 31.0*pi6/15120.0 + lmz2*(-7.0/720.0*pi4 + lmz2*(-1.0/144.0*pi2 - 1.0/720.0*lmz2))
             }
         }
     }

tests/common.rs

@@ -6,6 +6,25 @@ use std::path::PathBuf;
 use std::str::FromStr;
 
 
+#[macro_export]
+macro_rules! max_rel_diff {
+    ($a:expr, $b:expr) => {
+        ($a - $b).abs()/$a.abs().max($b.abs())
+    }
+}
+
+
+#[macro_export]
+macro_rules! assert_close_rel {
+    ($a:expr, $b:expr, $eps:expr) => {
+        if max_rel_diff!($a, $b) >= $eps {
+            println!("Numbers differ relatively by more than {:e}: {:e} != {:e}", $eps, $a, $b);
+            assert!(false);
+        }
+    }
+}
+
+
 #[macro_export]
 macro_rules! assert_eq_float {
     ($a:expr, $b:expr, $eps:expr) => {

tests/li2.rs

@@ -85,6 +85,14 @@ fn special_values() {
         assert_eq!(li2.im, li2_expected.im);
     }
 
+    {
+        // test value close to zero
+        let z = Complex::new(4.831285545908206e-6_f64, 0.004396919500211628_f64);
+        let expected = Complex::new(-1.94166578202937687444628936853e-9_f64, 0.00439692067657240512726530759719387623_f64);
+        assert_close_rel!(z.li2().re, expected.re, 1e-12);
+        assert_close_rel!(z.li2().im, expected.im, 1e-15);
+    }
+
     // test value that causes overflow if squared
     assert!(!Complex::new(1e300_f64, 1.0_f64).li2().is_infinite());
     assert!(!Complex::new(1.0_f64, 1e300_f64).li2().is_infinite());