Skip to content

Commit

Permalink
Merge branch 'mkirk/usize'
Browse files Browse the repository at this point in the history
  • Loading branch information
@michaelkirk
michaelkirk committed Jan 25, 2024
2 parents c4de165 + a9b8f92 commit 9f2b3b9
Show file tree
Hide file tree
Showing 4 changed files with 69 additions and 74 deletions.
5 changes: 4 additions & 1 deletion .github/workflows/test.yml
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,8 @@
on: push
name: Run tests
on:
- push
- pull_request
- merge_group
jobs:
linting:
name: formatting and clippy
Expand Down
68 changes: 29 additions & 39 deletions src/geodesic.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,13 +25,13 @@ pub struct Geodesic {
pub _b: f64,
pub _c2: f64,
_etol2: f64,
_A3x: [f64; GEODESIC_ORDER as usize],
_C3x: [f64; _nC3x_ as usize],
_C4x: [f64; _nC4x_ as usize],
_A3x: [f64; GEODESIC_ORDER],
_C3x: [f64; _nC3x_],
_C4x: [f64; _nC4x_],

pub GEODESIC_ORDER: i64,
_nC3x_: i64,
_nC4x_: i64,
pub GEODESIC_ORDER: usize,
_nC3x_: usize,
_nC4x_: usize,
maxit1_: u64,
maxit2_: u64,

Expand Down Expand Up @@ -80,11 +80,11 @@ const COEFF_C4: [f64; 77] = [
-3328.0, 1144.0, 315315.0, -128.0, 135135.0, -2560.0, 832.0, 405405.0, 128.0, 99099.0,
];

pub const GEODESIC_ORDER: i64 = 6;
pub const GEODESIC_ORDER: usize = 6;
#[allow(non_upper_case_globals)]
const _nC3x_: i64 = 15;
const _nC3x_: usize = 15;
#[allow(non_upper_case_globals)]
const _nC4x_: i64 = 21;
const _nC4x_: usize = 21;

impl Geodesic {
pub fn new(a: f64, f: f64) -> Self {
Expand Down Expand Up @@ -113,42 +113,37 @@ impl Geodesic {
/ 2.0;
let _etol2 = 0.1 * _tol2_ / (f.abs().max(0.001) * (1.0 - f / 2.0).min(1.0) / 2.0).sqrt();

let mut _A3x: [f64; GEODESIC_ORDER as usize] = [0.0; GEODESIC_ORDER as usize];
let mut _C3x: [f64; _nC3x_ as usize] = [0.0; _nC3x_ as usize];
let mut _C4x: [f64; _nC4x_ as usize] = [0.0; _nC4x_ as usize];
let mut _A3x: [f64; GEODESIC_ORDER] = [0.0; GEODESIC_ORDER];
let mut _C3x: [f64; _nC3x_] = [0.0; _nC3x_];
let mut _C4x: [f64; _nC4x_] = [0.0; _nC4x_];

// Call a3coeff
let mut o: i64 = 0;
let mut o: usize = 0;
for (k, j) in (0..GEODESIC_ORDER).rev().enumerate() {
let m = j.min(GEODESIC_ORDER - j - 1);
_A3x[k] = geomath::polyval(m as isize, &COEFF_A3[o as usize..], _n)
/ COEFF_A3[(o + m + 1) as usize];
_A3x[k] = geomath::polyval(m, &COEFF_A3[o..], _n) / COEFF_A3[o + m + 1];
o += m + 2;
}

// c3coeff
let mut o: i64 = 0;
let mut o = 0;
let mut k = 0;

for l in 1..GEODESIC_ORDER {
for j in (l..GEODESIC_ORDER).rev() {
let m = j.min(GEODESIC_ORDER - j - 1);
_C3x[k as usize] = geomath::polyval(m as isize, &COEFF_C3[o as usize..], _n)
/ COEFF_C3[(o + m + 1) as usize];
_C3x[k] = geomath::polyval(m, &COEFF_C3[o..], _n) / COEFF_C3[o + m + 1];
k += 1;
o += m + 2;
}
}

// c4coeff
let mut o: i64 = 0;
let mut o = 0;
let mut k = 0;

for l in 0..GEODESIC_ORDER {
for j in (l..GEODESIC_ORDER).rev() {
let m = GEODESIC_ORDER - j - 1;
_C4x[k as usize] = geomath::polyval(m as isize, &COEFF_C4[o as usize..], _n)
/ COEFF_C4[(o + m + 1) as usize];
_C4x[k] = geomath::polyval(m, &COEFF_C4[o..], _n) / COEFF_C4[o + m + 1];
k += 1;
o += m + 2;
}
Expand Down Expand Up @@ -184,31 +179,26 @@ impl Geodesic {
}

pub fn _A3f(&self, eps: f64) -> f64 {
geomath::polyval(self.GEODESIC_ORDER as isize - 1, &self._A3x, eps)
geomath::polyval(self.GEODESIC_ORDER - 1, &self._A3x, eps)
}

pub fn _C3f(&self, eps: f64, c: &mut [f64]) {
let mut mult = 1.0;
let mut o = 0;
for (l, c_item) in c
.iter_mut()
.enumerate()
.take(self.GEODESIC_ORDER as usize)
.skip(1)
{
let m = self.GEODESIC_ORDER as usize - l - 1;
for (l, c_item) in c.iter_mut().enumerate().take(self.GEODESIC_ORDER).skip(1) {
let m = self.GEODESIC_ORDER - l - 1;
mult *= eps;
*c_item = mult * geomath::polyval(m as isize, &self._C3x[o..], eps);
*c_item = mult * geomath::polyval(m, &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;
for (l, c_item) in c.iter_mut().enumerate().take(self.GEODESIC_ORDER as usize) {
let m = self.GEODESIC_ORDER as usize - l - 1;
*c_item = mult * geomath::polyval(m as isize, &self._C4x[o..], eps);
for (l, c_item) in c.iter_mut().enumerate().take(self.GEODESIC_ORDER) {
let m = self.GEODESIC_ORDER - l - 1;
*c_item = mult * geomath::polyval(m, &self._C4x[o..], eps);
o += m + 1;
mult *= eps;
}
Expand Down Expand Up @@ -265,7 +255,7 @@ impl Geodesic {
}
} else if outmask & (caps::REDUCEDLENGTH | caps::GEODESICSCALE) != 0 {
for l in 1..=self.GEODESIC_ORDER {
C2a[l as usize] = A1 * C1a[l as usize] - A2 * C2a[l as usize];
C2a[l] = A1 * C1a[l] - A2 * C2a[l];
}
J12 = m0x * sig12
+ (geomath::sin_cos_series(true, ssig2, csig2, C2a)
Expand Down Expand Up @@ -613,10 +603,10 @@ impl Geodesic {
let dn1 = (1.0 + self._ep2 * geomath::sq(sbet1)).sqrt();
let dn2 = (1.0 + self._ep2 * geomath::sq(sbet2)).sqrt();

const CARR_SIZE: usize = GEODESIC_ORDER as usize + 1;
const CARR_SIZE: usize = GEODESIC_ORDER + 1;
let mut C1a: [f64; CARR_SIZE] = [0.0; CARR_SIZE];
let mut C2a: [f64; CARR_SIZE] = [0.0; CARR_SIZE];
let mut C3a: [f64; GEODESIC_ORDER as usize] = [0.0; GEODESIC_ORDER as usize];
let mut C3a: [f64; GEODESIC_ORDER] = [0.0; GEODESIC_ORDER];

let mut meridian = lat1 == -90.0 || slam12 == 0.0;
let mut calp1 = 0.0;
Expand Down Expand Up @@ -833,7 +823,7 @@ impl Geodesic {
let A4 = geomath::sq(self.a) * calp0 * salp0 * self._e2;
geomath::norm(&mut ssig1, &mut csig1);
geomath::norm(&mut ssig2, &mut csig2);
let mut C4a: [f64; GEODESIC_ORDER as usize] = [0.0; GEODESIC_ORDER as usize];
let mut C4a: [f64; GEODESIC_ORDER] = [0.0; GEODESIC_ORDER];
self._C4f(eps, &mut C4a);
let B41 = geomath::sin_cos_series(false, ssig1, csig1, &C4a);
let B42 = geomath::sin_cos_series(false, ssig2, csig2, &C4a);
Expand Down
20 changes: 10 additions & 10 deletions src/geodesic_line.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,11 +16,11 @@ pub struct GeodesicLine {
_B21: f64,
_B31: f64,
_B41: f64,
_C1a: [f64; GEODESIC_ORDER as usize + 1],
_C1pa: [f64; GEODESIC_ORDER as usize + 1],
_C2a: [f64; GEODESIC_ORDER as usize + 1],
_C3a: [f64; GEODESIC_ORDER as usize],
_C4a: [f64; GEODESIC_ORDER as usize],
_C1a: [f64; GEODESIC_ORDER + 1],
_C1pa: [f64; GEODESIC_ORDER + 1],
_C2a: [f64; GEODESIC_ORDER + 1],
_C3a: [f64; GEODESIC_ORDER],
_C4a: [f64; GEODESIC_ORDER],
_b: f64,
_c2: f64,
_calp0: f64,
Expand Down Expand Up @@ -107,7 +107,7 @@ impl GeodesicLine {
let eps = _k2 / (2.0 * (1.0 + (1.0 + _k2).sqrt()) + _k2);

let mut _A1m1 = 0.0;
let mut _C1a: [f64; GEODESIC_ORDER as usize + 1] = [0.0; GEODESIC_ORDER as usize + 1];
let mut _C1a: [f64; GEODESIC_ORDER + 1] = [0.0; GEODESIC_ORDER + 1];
let mut _B11 = 0.0;
let mut _stau1 = 0.0;
let mut _ctau1 = 0.0;
Expand All @@ -121,21 +121,21 @@ impl GeodesicLine {
_ctau1 = _csig1 * c - _ssig1 * s;
}

let mut _C1pa: [f64; GEODESIC_ORDER as usize + 1] = [0.0; GEODESIC_ORDER as usize + 1];
let mut _C1pa: [f64; GEODESIC_ORDER + 1] = [0.0; GEODESIC_ORDER + 1];
if caps & caps::CAP_C1p != 0 {
geomath::_C1pf(eps, &mut _C1pa, geod.GEODESIC_ORDER);
}

let mut _A2m1 = 0.0;
let mut _C2a: [f64; GEODESIC_ORDER as usize + 1] = [0.0; GEODESIC_ORDER as usize + 1];
let mut _C2a: [f64; GEODESIC_ORDER + 1] = [0.0; GEODESIC_ORDER + 1];
let mut _B21 = 0.0;
if caps & caps::CAP_C2 != 0 {
_A2m1 = geomath::_A2m1f(eps, geod.GEODESIC_ORDER);
geomath::_C2f(eps, &mut _C2a, geod.GEODESIC_ORDER);
_B21 = geomath::sin_cos_series(true, _ssig1, _csig1, &_C2a);
}

let mut _C3a: [f64; GEODESIC_ORDER as usize] = [0.0; GEODESIC_ORDER as usize];
let mut _C3a: [f64; GEODESIC_ORDER] = [0.0; GEODESIC_ORDER];
let mut _A3c = 0.0;
let mut _B31 = 0.0;
if caps & caps::CAP_C3 != 0 {
Expand All @@ -144,7 +144,7 @@ impl GeodesicLine {
_B31 = geomath::sin_cos_series(true, _ssig1, _csig1, &_C3a);
}

let mut _C4a: [f64; GEODESIC_ORDER as usize] = [0.0; GEODESIC_ORDER as usize];
let mut _C4a: [f64; GEODESIC_ORDER] = [0.0; GEODESIC_ORDER];
let mut _A4 = 0.0;
let mut _B41 = 0.0;
if caps & caps::CAP_C4 != 0 {
Expand Down
50 changes: 26 additions & 24 deletions src/geomath.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,16 +52,12 @@ pub fn sum(u: f64, v: f64) -> (f64, f64) {
}

// Evaluate a polynomial
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
pub fn polyval(n: usize, p: &[f64], x: f64) -> f64 {
let mut y = p[0];
for val in &p[1..=n] {
y = y * x + val;
}
y
}

// Round an angle so taht small values underflow to 0
Expand Down Expand Up @@ -291,66 +287,72 @@ pub fn astroid(x: f64, y: f64) -> f64 {
}
}

pub fn _A1m1f(eps: f64, geodesic_order: i64) -> f64 {
pub fn _A1m1f(eps: f64, geodesic_order: usize) -> 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 as isize, &COEFF, sq(eps)) / COEFF[(m + 1) as usize];
let m = geodesic_order / 2;
let t = polyval(m, &COEFF, sq(eps)) / COEFF[m + 1];
(t + eps) / (1.0 - eps)
}

pub fn _C1f(eps: f64, c: &mut [f64], geodesic_order: i64) {
pub fn _C1f(eps: f64, c: &mut [f64], geodesic_order: usize) {
const COEFF: [f64; 18] = [
-1.0, 6.0, -16.0, 32.0, -9.0, 64.0, -128.0, 2048.0, 9.0, -16.0, 768.0, 3.0, -5.0, 512.0,
-7.0, 1280.0, -7.0, 2048.0,
];
let eps2 = sq(eps);
let mut d = eps;
let mut o = 0;
// Clippy wants us to turn this into `c.iter_mut().enumerate().take(geodesic_order + 1).skip(1)`
// but benching (rust-1.75) shows that it would be slower.
#[allow(clippy::needless_range_loop)]
for l in 1..=geodesic_order {
let m = (geodesic_order - l) / 2;
c[l as usize] =
d * polyval(m as isize, &COEFF[o as usize..], eps2) / COEFF[(o + m + 1) as usize];
c[l] = d * polyval(m, &COEFF[o..], eps2) / COEFF[o + m + 1];
o += m + 2;
d *= eps;
}
}

pub fn _C1pf(eps: f64, c: &mut [f64], geodesic_order: i64) {
pub fn _C1pf(eps: f64, c: &mut [f64], geodesic_order: usize) {
const COEFF: [f64; 18] = [
205.0, -432.0, 768.0, 1536.0, 4005.0, -4736.0, 3840.0, 12288.0, -225.0, 116.0, 384.0,
-7173.0, 2695.0, 7680.0, 3467.0, 7680.0, 38081.0, 61440.0,
];
let eps2 = sq(eps);
let mut d = eps;
let mut o = 0;
// Clippy wants us to turn this into `c.iter_mut().enumerate().take(geodesic_order + 1).skip(1)`
// but benching (rust-1.75) shows that it would be slower.
#[allow(clippy::needless_range_loop)]
for l in 1..=geodesic_order {
let m = (geodesic_order - l) / 2;
c[l as usize] =
d * polyval(m as isize, &COEFF[o as usize..], eps2) / COEFF[(o + m + 1) as usize];
c[l] = d * polyval(m, &COEFF[o..], eps2) / COEFF[o + m + 1];
o += m + 2;
d *= eps;
}
}

pub fn _A2m1f(eps: f64, geodesic_order: i64) -> f64 {
pub fn _A2m1f(eps: f64, geodesic_order: usize) -> 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 as isize, &COEFF, sq(eps)) / COEFF[(m + 1) as usize];
let m = geodesic_order / 2;
let t = polyval(m, &COEFF, sq(eps)) / COEFF[m + 1];
(t - eps) / (1.0 + eps)
}

pub fn _C2f(eps: f64, c: &mut [f64], geodesic_order: i64) {
pub fn _C2f(eps: f64, c: &mut [f64], geodesic_order: usize) {
const COEFF: [f64; 18] = [
1.0, 2.0, 16.0, 32.0, 35.0, 64.0, 384.0, 2048.0, 15.0, 80.0, 768.0, 7.0, 35.0, 512.0, 63.0,
1280.0, 77.0, 2048.0,
];
let eps2 = sq(eps);
let mut d = eps;
let mut o = 0;
// Clippy wants us to turn this into `c.iter_mut().enumerate().take(geodesic_order + 1).skip(1)`
// but benching (rust-1.75) shows that it would be slower.
#[allow(clippy::needless_range_loop)]
for l in 1..=geodesic_order {
let m = (geodesic_order - l) / 2;
c[l as usize] =
d * polyval(m as isize, &COEFF[o as usize..], eps2) / COEFF[(o + m + 1) as usize];
c[l] = d * polyval(m, &COEFF[o..], eps2) / COEFF[o + m + 1];
o += m + 2;
d *= eps;
}
Expand Down

0 comments on commit 9f2b3b9

Please sign in to comment.