Address georust#6: polyval.

Slightly different from my original proposal.
stonylohr · Jan 30, 2021 · 6fe6f60 · 6fe6f60
1 parent d4162ad
commit 6fe6f60
Show file tree

Hide file tree

Showing 2 changed files with 28 additions and 29 deletions.
diff --git a/src/geodesic.rs b/src/geodesic.rs
@@ -122,7 +122,7 @@ impl Geodesic {
         let mut k = 0;
         for j in (0..GEODESIC_ORDER).rev() {
             let m = j.min(GEODESIC_ORDER as i64 - j - 1);
-            _A3x[k as usize] = geomath::polyval(m, &COEFF_A3, o as usize, _n)
+            _A3x[k as usize] = geomath::polyval(m as isize, &COEFF_A3[o as usize..], _n)
                 / COEFF_A3[(o + m + 1) as usize] as f64;
             k += 1;
             o += m + 2;
@@ -135,7 +135,7 @@ impl Geodesic {
         for l in 1..GEODESIC_ORDER {
             for j in (l..GEODESIC_ORDER).rev() {
                 let m = j.min(GEODESIC_ORDER as i64 - j - 1);
-                _C3x[k as usize] = geomath::polyval(m, &COEFF_C3, o as usize, _n)
+                _C3x[k as usize] = geomath::polyval(m as isize, &COEFF_C3[o as usize..], _n)
                     / COEFF_C3[(o + m + 1) as usize] as f64;
                 k += 1;
                 o += m + 2;
@@ -149,7 +149,7 @@ impl Geodesic {
         for l in 0..GEODESIC_ORDER {
             for j in (l..GEODESIC_ORDER).rev() {
                 let m = GEODESIC_ORDER as i64 - j - 1;
-                _C4x[k as usize] = geomath::polyval(m, &COEFF_C4, o as usize, _n)
+                _C4x[k as usize] = geomath::polyval(m as isize, &COEFF_C4[o as usize..], _n)
                     / COEFF_C4[(o + m + 1) as usize] as f64;
                 k += 1;
                 o += m + 2;
@@ -186,27 +186,27 @@ impl Geodesic {
     }
 
     pub fn _A3f(&self, eps: f64) -> f64 {
-        geomath::polyval(self.GEODESIC_ORDER - 1, &self._A3x, 0, eps)
+        geomath::polyval(self.GEODESIC_ORDER as isize - 1, &self._A3x, eps)
     }
 
     pub fn _C3f(&self, eps: f64, c: &mut [f64]) {
         let mut mult = 1.0;
-        let mut o = 0.0;
-        for l in 1..self.GEODESIC_ORDER {
-            let m = self.GEODESIC_ORDER - l - 1;
+        let mut o = 0;
+        for l in 1..self.GEODESIC_ORDER as usize {
+            let m = self.GEODESIC_ORDER as usize - l - 1;
             mult *= eps;
-            c[l as usize] = mult * geomath::polyval(m, &self._C3x, o as usize, eps);
-            o += m as f64 + 1.0;
+            c[l] = mult * geomath::polyval(m as isize, &self._C3x[o..], eps);
+            o += m + 1;
         }
     }
 
     pub fn _C4f(&self, eps: f64, c: &mut [f64]) {
         let mut mult = 1.0;
-        let mut o = 0.0;
-        for l in 0..self.GEODESIC_ORDER {
-            let m = self.GEODESIC_ORDER - l - 1;
-            c[l as usize] = mult * geomath::polyval(m, &self._C4x, o as usize, eps);
-            o += m as f64 + 1.0;
+        let mut o = 0;
+        for l in 0..self.GEODESIC_ORDER as usize {
+            let m = self.GEODESIC_ORDER as usize - l - 1;
+            c[l] = mult * geomath::polyval(m as isize, &self._C4x[o..], eps);
+            o += m + 1;
             mult *= eps;
         }
     }

diff --git a/src/geomath.rs b/src/geomath.rs
@@ -51,17 +51,16 @@ pub fn sum(u: f64, v: f64) -> (f64, f64) {
 }
 
 // Evaluate a polynomial
-pub fn polyval(n: i64, p: &[f64], s: usize, x: f64) -> f64 {
-    let mut s = s;
-    let mut n = n;
-    let mut y = if n < 0 { 0.0 } else { p[s] };
-    assert!((n as usize) < (std::usize::MAX - s));
-    while n > 0 {
-        n -= 1;
-        s += 1;
-        y = y * x + p[s];
+pub fn polyval(n: isize, p: &[f64], x: f64) -> f64 {
+    if n < 0 {
+        0.0
+    } else {
+        let mut y = p[0];
+        for val in &p[1..=n as usize] {
+            y = y * x + val;
+        }
+        y
     }
-    y
 }
 
 // Round an angle so taht small values underflow to 0
@@ -300,7 +299,7 @@ pub fn astroid(x: f64, y: f64) -> f64 {
 pub fn _A1m1f(eps: f64, geodesic_order: i64) -> f64 {
     const COEFF: [f64; 5] = [1.0, 4.0, 64.0, 0.0, 256.0];
     let m: i64 = geodesic_order / 2;
-    let t = polyval(m, &COEFF, 0, sq(eps)) / COEFF[(m + 1) as usize] as f64;
+    let t = polyval(m as isize, &COEFF, sq(eps)) / COEFF[(m + 1) as usize] as f64;
     (t + eps) / (1.0 - eps)
 }
 
@@ -315,7 +314,7 @@ pub fn _C1f(eps: f64, c: &mut [f64], geodesic_order: i64) {
     for l in 1..=geodesic_order {
         let m = ((geodesic_order - l) / 2) as i64;
         c[l as usize] =
-            d * polyval(m, &COEFF, o as usize, eps2) / COEFF[(o + m + 1) as usize] as f64;
+            d * polyval(m as isize, &COEFF[o as usize..], eps2) / COEFF[(o + m + 1) as usize] as f64;
         o += m + 2;
         d *= eps;
     }
@@ -332,7 +331,7 @@ pub fn _C1pf(eps: f64, c: &mut [f64], geodesic_order: i64) {
     for l in 1..=geodesic_order {
         let m = (geodesic_order - l) / 2;
         c[l as usize] =
-            d * polyval(m as i64, &COEFF, o as usize, eps2) / COEFF[(o + m + 1) as usize] as f64;
+            d * polyval(m as isize, &COEFF[o as usize..], eps2) / COEFF[(o + m + 1) as usize] as f64;
         o += m + 2;
         d *= eps;
     }
@@ -341,7 +340,7 @@ pub fn _C1pf(eps: f64, c: &mut [f64], geodesic_order: i64) {
 pub fn _A2m1f(eps: f64, geodesic_order: i64) -> f64 {
     const COEFF: [f64; 5] = [-11.0, -28.0, -192.0, 0.0, 256.0];
     let m: i64 = geodesic_order / 2;
-    let t = polyval(m, &COEFF, 0, sq(eps)) / COEFF[(m + 1) as usize] as f64;
+    let t = polyval(m as isize, &COEFF, sq(eps)) / COEFF[(m + 1) as usize] as f64;
     (t - eps) / (1.0 + eps)
 }
 
@@ -356,7 +355,7 @@ pub fn _C2f(eps: f64, c: &mut [f64], geodesic_order: i64) {
     for l in 1..=geodesic_order {
         let m = (geodesic_order - l) / 2;
         c[l as usize] =
-            d * polyval(m as i64, &COEFF, o as usize, eps2) / COEFF[(o + m + 1) as usize] as f64;
+            d * polyval(m as isize, &COEFF[o as usize..], eps2) / COEFF[(o + m + 1) as usize] as f64;
         o += m + 2;
         d *= eps;
     }