Merge #39

39: Address #6: polyval r=michaelkirk a=stonylohr

Slightly different from my original proposal. There are a few pieces worth highlighting.

First, I think the most idiomatic Rust approach would actually be to drop the `n` argument as well as `s`, and let the size of the `p` argument imply the order of the polynomial. The reason I didn't take this further step is that keeping `n` requires fewer changes, and leaves the code a closer parallel to its C++ and Java counterparts. I figure we can always choose to reconsider later, perhaps if and when we focus on performance. The performance implications are likely minor, and I'm not sure which approach they favor.

Second, while keeping `n`, it's somewhat tempting to change it to `usize`, since its support for the negative special case doesn't appear to be used in any of Karney's code that I've examined or logged. It's certainly not used in any of our current geodesic functionality. It may be used in geoid, magnetic model, or projection special cases that aren't hit by any of the scenarios I've run in my instrumented C++. My current guess is that it's not actually used and was only added for completeness. Anyway, I left it signed for now, but changed it from i64 to isize since we ultimately need to cast it to usize for use as an array index, and any value too large for i32 would strongly hint that we're already off the rails. The largest `n` I've seen in instrumented scenarios is 29.

Third, I went with a slightly different polyval implementation than my original proposal. I felt that my new version is slightly more clear that it never tries to cast a negative `n` to usize (though both should actually be safe), and that it might give the optimizer more of a hint that we don't care about the `i` value except for accessing an element of `p` (though I suspect it would recognize that anyway).

Fourth, I tried to keep changes to polyval's callers minimal in most cases, but there were two I couldn't resist: Both `_C3f` and `_C4f` were using floating point values to represent array indices. I switched these to usize, also tweaking the types of some variables they interact with.

Finally, these changes should be behavior-neutral. I think they'll help me with arcdirect troubleshooting because polyval is called extensively in some sections of code used to calculate Geodesic member variables that end up with suspect values, and these changes will make it easier for me to compare the C++ and Rust code for these sections. However, they don't do anything to resolve the arcdirect problems themselves.

I'd be happy to revise my approach if you'd prefer something different for any of these various items.

- [ ] I agree to follow the project's [code of conduct](https://github.com/georust/geo/blob/master/CODE_OF_CONDUCT.md).
- [ ] I added an entry to `CHANGES.md` if knowledge of this change could be valuable to users.
---



Co-authored-by: Stony <stony.lohr@gmail.com>

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;
         }
     }

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;
     }