extract complex approximation

Expander · Jan 27, 2024 · 18910aa · 18910aa
1 parent 8e47f05
commit 18910aa
Showing 1 changed file with 54 additions and 45 deletions.
diff --git a/src/li2.rs b/src/li2.rs
@@ -2,11 +2,11 @@ use num::complex::Complex;
 use crate::cln::CLn;
 
 trait Li2Approx<T> {
-    /// rational function approximation of Re[Li2(x)] for x in [0, 1/2]
     fn approx(&self) -> T;
 }
 
 impl Li2Approx<f32> for f32 {
+    /// rational function approximation of Re[Li2(x)] for x in [0, 1/2]
     fn approx(&self) -> f32 {
         let cp = [ 1.00000020_f32, -0.780790946_f32, 0.0648256871_f32 ];
         let cq = [ 1.00000000_f32, -1.03077545_f32, 0.211216710_f32 ];
@@ -20,6 +20,7 @@ impl Li2Approx<f32> for f32 {
 }
 
 impl Li2Approx<f64> for f64 {
+    /// rational function approximation of Re[Li2(x)] for x in [0, 1/2]
     fn approx(&self) -> f64 {
         let cp = [
             0.9999999999999999502e+0,
@@ -51,6 +52,56 @@ impl Li2Approx<f64> for f64 {
     }
 }
 
+impl Li2Approx<Complex<f32>> for Complex<f32> {
+    fn approx(&self) -> Complex<f32> {
+        // bf[1..N-1] are the even Bernoulli numbers / (2 n + 1)!
+        // generated by: Table[BernoulliB[2 n]/(2 n + 1)!, {n, 1, 19}]
+        let bf = [
+            - 1.0_f32/4.0_f32,
+              1.0_f32/36.0_f32,
+            - 1.0_f32/3600.0_f32,
+              1.0_f32/211680.0_f32,
+        ];
+        let x = *self;
+        let x2 = x*x;
+
+        x + x2*(bf[0] + x*(bf[1] + x2*(bf[2] + x2*bf[3])))
+    }
+}
+
+impl Li2Approx<Complex<f64>> for Complex<f64> {
+    fn approx(&self) -> Complex<f64> {
+        // bf[1..N-1] are the even Bernoulli numbers / (2 n + 1)!
+        // generated by: Table[BernoulliB[2 n]/(2 n + 1)!, {n, 1, 19}]
+        let bf = [
+            - 1.0_f64/4.0_f64,
+              1.0_f64/36.0_f64,
+            - 1.0_f64/3600.0_f64,
+              1.0_f64/211680.0_f64,
+            - 1.0_f64/10886400.0_f64,
+              1.0_f64/526901760.0_f64,
+            - 4.0647616451442255e-11_f64,
+              8.9216910204564526e-13_f64,
+            - 1.9939295860721076e-14_f64,
+              4.5189800296199182e-16_f64,
+        ];
+        let x = *self;
+        let x2 = x*x;
+        let x4 = x2*x2;
+
+        x + x2*(bf[0] +
+                x*(bf[1] +
+                   x2*(
+                       bf[2] +
+                       x2*bf[3] +
+                       x4*(bf[4] + x2*bf[5]) +
+                       x4*x4*(bf[6] + x2*bf[7] + x4*(bf[8] + x2*bf[9]))
+                   )
+                )
+        )
+    }
+}
+
 /// Provides the 2nd order polylogarithm (dilogarithm) function
 /// `li2()` of a number of type `T`.
 pub trait Li2<T> {
@@ -158,16 +209,6 @@ impl Li2<Complex<f32>> for Complex<f32> {
     /// ```
     fn li2(&self) -> Complex<f32> {
         let pi = std::f32::consts::PI;
-
-        // bf[1..N-1] are the even Bernoulli numbers / (2 n + 1)!
-        // generated by: Table[BernoulliB[2 n]/(2 n + 1)!, {n, 1, 19}]
-        let bf = [
-            - 1.0_f32/4.0_f32,
-              1.0_f32/36.0_f32,
-            - 1.0_f32/3600.0_f32,
-              1.0_f32/211680.0_f32,
-        ];
-
         let rz = self.re;
         let iz = self.im;
 
@@ -200,10 +241,7 @@ impl Li2<Complex<f32>> for Complex<f32> {
                     }
                 };
 
-                let u2 = u*u;
-                let sum = u + u2*(bf[0] + u*(bf[1] + u2*(bf[2] + u2*bf[3])));
-
-                sgn*sum + rest
+                sgn*u.approx() + rest
             }
         }
     }
@@ -225,22 +263,6 @@ impl Li2<Complex<f64>> for Complex<f64> {
     /// ```
     fn li2(&self) -> Complex<f64> {
         let pi = std::f64::consts::PI;
-
-        // bf[1..N-1] are the even Bernoulli numbers / (2 n + 1)!
-        // generated by: Table[BernoulliB[2 n]/(2 n + 1)!, {n, 1, 19}]
-        let bf = [
-            - 1.0_f64/4.0_f64,
-              1.0_f64/36.0_f64,
-            - 1.0_f64/3600.0_f64,
-              1.0_f64/211680.0_f64,
-            - 1.0_f64/10886400.0_f64,
-              1.0_f64/526901760.0_f64,
-            - 4.0647616451442255e-11_f64,
-              8.9216910204564526e-13_f64,
-            - 1.9939295860721076e-14_f64,
-              4.5189800296199182e-16_f64,
-        ];
-
         let rz = self.re;
         let iz = self.im;
 
@@ -273,20 +295,7 @@ impl Li2<Complex<f64>> for Complex<f64> {
                     }
                 };
 
-                let u2 = u*u;
-                let u4 = u2*u2;
-                let sum =
-                    u +
-                    u2*(bf[0] +
-                    u *(bf[1] +
-                    u2*(
-                        bf[2] +
-                        u2*bf[3] +
-                        u4*(bf[4] + u2*bf[5]) +
-                        u4*u4*(bf[6] + u2*bf[7] + u4*(bf[8] + u2*bf[9]))
-                    )));
-
-                sgn*sum + rest
+                sgn*u.approx() + rest
             }
         }
     }