From 6fe6f60f26c8ffb440534fdb1eec1f346afdfaac Mon Sep 17 00:00:00 2001 From: Stony Date: Sat, 30 Jan 2021 14:54:02 -0800 Subject: [PATCH] Address #6: polyval. Slightly different from my original proposal. --- src/geodesic.rs | 28 ++++++++++++++-------------- src/geomath.rs | 29 ++++++++++++++--------------- 2 files changed, 28 insertions(+), 29 deletions(-) diff --git a/src/geodesic.rs b/src/geodesic.rs index a4307ce..34229a8 100644 --- 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 index 56e7771..d4a4b5a 100644 --- 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; }