avoid local variables

Expander · Jan 28, 2024 · 867ac16 · 867ac16
1 parent 5569f1b
commit 867ac16
Showing 1 changed file with 32 additions and 28 deletions.
diff --git a/src/li5.rs b/src/li5.rs
@@ -22,18 +22,6 @@ impl Li5<Complex<f64>> for Complex<f64> {
         let pi  = std::f64::consts::PI;
         let pi2 = pi*pi;
         let z5  = 1.0369277551433699; // zeta(5)
-        let bf  = [
-            1.0                   , -15.0/32.0             ,
-            1.3953189300411523e-01, -2.8633777006172840e-02,
-            4.0317412551440329e-03, -3.3985018004115226e-04,
-            4.5445184621617666e-06,  2.3916808048569012e-06,
-           -1.2762692600122747e-07, -3.1628984306505932e-08,
-            3.2848118445335192e-09,  4.7613713995660579e-10,
-           -8.0846898171909830e-11, -7.2387648587737207e-12,
-            1.9439760115173968e-12,  1.0256978405977236e-13,
-           -4.6180551009884830e-14, -1.1535857196470580e-15,
-            1.0903545401333394e-15
-        ];
 
         if self.im == 0.0 && self.re == 0.0 {
             Complex::new(0.0, 0.0)
@@ -71,28 +59,44 @@ impl Li5<Complex<f64>> for Complex<f64> {
                 u2 * (cs[4] +
                 u2 * (cs[5]))))))))
             } else {
-                let (u, rest) = if nz <= 1.0 {
-                    (-(1.0 - self).cln(), Complex::new(0.0, 0.0))
+                if nz <= 1.0 {
+                    cli5_unit_circle(-(1.0 - self).cln())
                 } 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;
-                    (-(1.0 - 1.0/self).cln(), -1.0/360.0*lmz*(7.0*pi4 + lmz2*(10.0*pi2 + 3.0*lmz2)))
-                };
-
-                let u2 = u*u;
-                let u4 = u2*u2;
-                let u8 = u4*u4;
-
-                rest +
-                u*bf[0] +
-                u2*(bf[1] + u*bf[2]) +
-                u4*(bf[3] + u*bf[4] + u2*(bf[5] + u*bf[6])) +
-                u8*(bf[7] + u*bf[8] + u2*(bf[9] + u*bf[10]) +
-                    u4*(bf[11] + u*bf[12] + u2*(bf[13] + u*bf[14]))) +
-                u8*u8*(bf[15] + u*bf[16] + u2*(bf[17] + u*bf[18]))
+                    cli5_unit_circle(-(1.0 - 1.0/self).cln()) - 1.0/360.0*lmz*(7.0*pi4 + lmz2*(10.0*pi2 + 3.0*lmz2))
+                }
             }
         }
     }
 }
+
+/// series approximation of Li5(z) for |z| <= 1
+/// in terms of x = -ln(1 - z)
+fn cli5_unit_circle(x: Complex<f64>) -> Complex<f64> {
+    let bf  = [
+        1.0                   , -15.0/32.0             ,
+        1.3953189300411523e-01, -2.8633777006172840e-02,
+        4.0317412551440329e-03, -3.3985018004115226e-04,
+        4.5445184621617666e-06,  2.3916808048569012e-06,
+       -1.2762692600122747e-07, -3.1628984306505932e-08,
+        3.2848118445335192e-09,  4.7613713995660579e-10,
+       -8.0846898171909830e-11, -7.2387648587737207e-12,
+        1.9439760115173968e-12,  1.0256978405977236e-13,
+       -4.6180551009884830e-14, -1.1535857196470580e-15,
+        1.0903545401333394e-15
+    ];
+
+    let x2 = x*x;
+    let x4 = x2*x2;
+    let x8 = x4*x4;
+
+    x*bf[0] +
+    x2*(bf[1] + x*bf[2]) +
+    x4*(bf[3] + x*bf[4] + x2*(bf[5] + x*bf[6])) +
+    x8*(bf[7] + x*bf[8] + x2*(bf[9] + x*bf[10]) +
+        x4*(bf[11] + x*bf[12] + x2*(bf[13] + x*bf[14]))) +
+    x8*x8*(bf[15] + x*bf[16] + x2*(bf[17] + x*bf[18]))
+}