#![no_std]
pub mod predicates;
#[cfg(feature = "transpiled")]
mod transpiled;
pub use predicates::{
incircle,
incircle_fast,
insphere,
insphere_fast,
orient2d,
orient2d_fast,
orient3d,
orient3d_fast,
};
#[cfg(test)]
mod tests {
extern crate rand;
use self::rand::{Rng, SeedableRng, StdRng};
use super::*;
use predicates::{
incircle_exact, incircle_slow, insphere_exact, insphere_slow, orient2d_exact,
orient2d_slow, orient3d_exact, orient3d_slow,
};
const SEED: &[usize] = &[1, 2, 3, 4];
const EXP_BOUNDS: [i32; 2] = [-142, 201];
fn tol(rng: &mut StdRng) -> f64 {
let max_exp = (EXP_BOUNDS[0] + 1) as f64;
10.0_f64.powi((max_exp * rng.gen::<f64>()).round() as i32) * (rng.gen::<f64>() - 0.5)
}
macro_rules! quick_perm_impl {
($n:expr, $struct_name:ident) => {
struct $struct_name<T> {
perm: [T; $n],
p: [usize; $n + 1],
index: usize,
}
impl<T> $struct_name<T> {
fn new(v: [T; $n]) -> Self {
let mut p = [0; $n + 1];
for i in 0..$n + 1 {
p[i] = i;
}
$struct_name {
perm: v,
p,
index: 0,
}
}
}
impl<T: Clone> Iterator for $struct_name<T> {
type Item = [T; $n];
fn next(&mut self) -> Option<Self::Item> {
let $struct_name { perm, p, index } = self;
let mut i = *index;
let n = p.len() - 1;
if i == 0 {
*index += 1;
return Some(perm.clone());
} else if i >= n {
return None;
}
p[i] -= 1;
let j = if i % 2 == 0 { 0 } else { p[i] };
perm.swap(j, i);
i = 1;
while p[i] == 0 {
p[i] = i;
i += 1;
}
*index = i;
Some(perm.clone())
}
}
};
}
quick_perm_impl!(3, QuickPerm3);
quick_perm_impl!(4, QuickPerm4);
quick_perm_impl!(5, QuickPerm5);
#[test]
fn orient2d_test() {
let a = [0.0, 1.0];
let b = [2.0, 3.0];
let c = [4.0, 5.0];
assert_eq!(orient2d(a, b, c), 0.0);
}
#[cfg(feature = "transpiled")]
#[test]
fn orient2d_transpiled_regression_test() {
unsafe {
transpiled::exactinit();
}
let mut rng: StdRng = SeedableRng::from_seed(SEED);
let n = 99999;
for _ in 0..n {
let a = [tol(&mut rng), tol(&mut rng)];
let b = [12.0, 12.0];
let c = [24.0, 24.0];
let o2d = orient2d;
let o2d_transpiled = transpiled::orient2d;
for p in QuickPerm3::new([a, b, c]) {
assert_eq!(
o2d(p[0], p[1], p[2]),
o2d_transpiled(p[0], p[1], p[2]),
"{:?}",
a
);
}
}
}
#[test]
fn orient2d_robustness_test() {
let mut rng: StdRng = SeedableRng::from_seed(SEED);
let n = 99999;
for o2d in &[orient2d, orient2d_exact, orient2d_slow] {
for _ in 0..n {
let pa = [tol(&mut rng), tol(&mut rng)];
let pb = [12.0, 12.0];
let pc = [24.0, 24.0];
let main = o2d(pa, pb, pc);
if main == 0.0 {
for [a, b, c] in QuickPerm3::new([pa, pb, pc]) {
assert_eq!(o2d(a, b, c), 0.0);
}
}
let pred = main > 0.0;
for (i, [a, b, c]) in QuickPerm3::new([pa, pb, pc]).enumerate() {
let t = o2d(a, b, c);
assert_eq!(pred, if i % 2 == 0 { t > 0.0 } else { t < 0.0 });
}
}
}
}
#[test]
fn orient2d_regression_test() {
let mut rng: StdRng = SeedableRng::from_seed(SEED);
let n = 99999;
for _ in 0..n {
let pa = [tol(&mut rng), tol(&mut rng)];
let pb = [12.0, 12.0];
let pc = [24.0, 24.0];
let o2d = predicates::orient2d;
let o2de = predicates::orient2d_exact;
let o2ds = predicates::orient2d_slow;
for [a, b, c] in QuickPerm3::new([pa, pb, pc]) {
assert_eq!(
o2d(a, b, c).partial_cmp(&0.0),
o2de(a, b, c).partial_cmp(&0.0),
"{:?}",
pa
);
assert_eq!(
o2d(a, b, c).partial_cmp(&0.0),
o2ds(a, b, c).partial_cmp(&0.0),
"{:?}",
pa
);
}
}
}
#[test]
fn orient2d_fast_test() {
let a = [0.0, 1.0];
let b = [2.0, 3.0];
let c = [4.0, 5.0];
assert_eq!(orient2d_fast(a, b, c), 0.0);
let mut rng: StdRng = SeedableRng::from_seed(SEED);
let tol = 5.0e-10;
for _ in 0..999 {
let a = [0.5 + tol * rng.gen::<f64>(), 0.5 + tol * rng.gen::<f64>()];
let b = [12.0, 12.0];
let c = [24.0, 24.0];
assert_eq!(orient2d_fast(a, b, c) > 0.0, orient2d_fast(b, c, a) > 0.0);
assert_eq!(orient2d_fast(b, c, a) > 0.0, orient2d_fast(c, a, b) > 0.0);
assert_eq!(orient2d_fast(a, b, c) > 0.0, orient2d_fast(b, a, c) < 0.0);
assert_eq!(orient2d_fast(a, b, c) > 0.0, orient2d_fast(a, c, b) < 0.0);
}
}
#[test]
fn orient3d_test() {
let a = [0.0, 1.0, 6.0];
let b = [2.0, 3.0, 4.0];
let c = [4.0, 5.0, 1.0];
let d = [6.0, 2.0, 5.3];
assert_eq!(orient3d(a, b, c, d), 10.0);
}
#[cfg(feature = "transpiled")]
#[test]
fn orient3d_transpiled_regression_test() {
unsafe {
transpiled::exactinit();
}
let mut rng: StdRng = SeedableRng::from_seed(SEED);
let n = 9999;
for _ in 0..n {
let pa = [tol(&mut rng), tol(&mut rng), tol(&mut rng)];
let pb = [12.0, 12.0, 12.0];
let pc = [24.0, 24.0, 24.0];
let pd = [48.0, 48.0, 48.0];
let o3d = predicates::orient3d;
let o3d_transpiled = transpiled::orient3d;
for [a, b, c, d] in QuickPerm4::new([pa, pb, pc, pd]) {
assert_eq!(o3d(a, c, d, b), o3d_transpiled(a, c, d, b), "{:?}", pa);
}
}
}
#[test]
fn orient3d_regression_test() {
let mut rng: StdRng = SeedableRng::from_seed(SEED);
let n = 5000;
for _ in 0..n {
let pa = [tol(&mut rng), tol(&mut rng), tol(&mut rng)];
let pb = [12.0, 12.0, 12.0];
let pc = [24.0, 24.0, 24.0];
let pd = [48.0, 48.0, 48.0];
let o3d = predicates::orient3d;
let o3de = predicates::orient3d_exact;
let o3ds = predicates::orient3d_slow;
for [a, b, c, d] in QuickPerm4::new([pa, pb, pc, pd]) {
assert_eq!(o3d(a, c, d, b), o3de(a, c, d, b), "{:?}", pa);
assert_eq!(o3d(a, c, d, b), o3ds(a, c, d, b), "{:?}", pa);
}
}
}
#[test]
fn orient3d_robustness_test() {
let mut rng: StdRng = SeedableRng::from_seed(SEED);
let n = 999;
for o3d in &[orient3d, orient3d_exact, orient3d_slow] {
for _ in 0..n {
let pa = [tol(&mut rng), tol(&mut rng), tol(&mut rng)];
let pb = [12.0, 12.0, 12.0];
let pc = [24.0, 24.0, 24.0];
let pd = [48.0, 48.0, 48.0];
let main = o3d(pa, pb, pc, pd);
if main == 0.0 {
for [a, b, c, d] in QuickPerm4::new([pa, pb, pc, pd]) {
assert_eq!(o3d(a, b, c, d), 0.0);
}
}
let pred = main > 0.0;
for (i, [a, b, c, d]) in QuickPerm4::new([pa, pb, pc, pd]).enumerate() {
let t = o3d(a, b, c, d);
assert_eq!(
pred,
if i % 2 == 0 { t > 0.0 } else { t < 0.0 },
"{} vs. {} at {:?}",
t,
-main,
pa
);
}
}
}
}
#[test]
fn orient3d_fast_test() {
let a = [0.0, 1.0, 6.0];
let b = [2.0, 3.0, 4.0];
let c = [4.0, 5.0, 1.0];
let d = [6.0, 2.0, 5.3];
assert!(f64::abs(orient3d_fast(a, b, c, d) - 10.0) < 1.0e-4);
let mut rng: StdRng = SeedableRng::from_seed(SEED);
let tol = 5.0;
let a = [
0.5 + tol * rng.gen::<f64>(),
0.5 + tol * rng.gen::<f64>(),
0.5 + tol * rng.gen::<f64>(),
];
let b = [12.0, 12.0, 12.0];
let c = [24.0, 24.0, 24.0];
let d = [48.0, 48.0, 48.0];
assert_eq!(
orient3d_fast(a, b, c, d) > 0.0,
orient3d_fast(b, c, a, d) > 0.0
);
assert_eq!(
orient3d_fast(b, a, c, d) < 0.0,
orient3d_fast(c, b, d, a) < 0.0
);
assert_ne!(
orient3d_fast(b, c, d, a) > 0.0,
orient3d_fast(c, d, a, b) > 0.0
);
assert_ne!(
orient3d_fast(b, d, c, a) < 0.0,
orient3d_fast(c, d, b, a) < 0.0
);
}
#[test]
fn incircle_test() {
let a = [0.0, 1.0];
let b = [1.0, 0.0];
let c = [1.0, 1.0];
let d = [0.0, 0.0];
assert_eq!(incircle(a, b, c, d), 0.0);
let d = [0.1, 0.1];
assert!(incircle(a, b, c, d) > 0.0);
let d = [-0.1, -0.1];
assert!(incircle(a, b, c, d) < 0.0);
}
#[test]
fn incircle_robustness_test() {
let mut rng: StdRng = SeedableRng::from_seed(SEED);
let n = 999;
for ic in &[incircle, incircle_exact, incircle_slow] {
for _ in 0..n {
let pa = [0.0, 1.0];
let pb = [1.0, 0.0];
let pc = [1.0, 1.0];
let pd = [tol(&mut rng), tol(&mut rng)];
let main = ic(pa, pb, pc, pd);
if main == 0.0 {
for [a, b, c, d] in QuickPerm4::new([pa, pb, pc, pd]) {
assert_eq!(ic(a, b, c, d), 0.0);
}
}
let pred = main > 0.0;
for (i, [a, b, c, d]) in QuickPerm4::new([pa, pb, pc, pd]).enumerate() {
let t = ic(a, b, c, d);
assert_eq!(
pred,
if i % 2 == 0 { t > 0.0 } else { t < 0.0 },
"{} vs. {} at {:?}",
t,
-main,
pd
);
}
}
}
}
#[test]
fn insphere_test() {
let a = [0.0, 1.0, 0.0];
let b = [1.0, 0.0, 0.0];
let c = [0.0, 0.0, 1.0];
let d = [1.0, 1.0, 1.0];
let e = [0.0, 0.0, 0.0];
assert_eq!(insphere(a, b, c, d, e), 0.0);
let e = [0.1, 0.1, 0.1];
assert!(insphere(a, b, c, d, e) > 0.0);
let e = [-0.1, -0.1, -0.1];
assert!(insphere(a, b, c, d, e) < 0.0);
}
#[test]
fn insphere_robustness_test() {
let mut rng: StdRng = SeedableRng::from_seed(SEED);
let n = 99;
for is in &[insphere, insphere_exact, insphere_slow] {
for _ in 0..n {
let pa = [0.0, 1.0, 0.0];
let pb = [1.0, 0.0, 0.0];
let pc = [0.0, 0.0, 1.0];
let pd = [1.0, 1.0, 1.0];
let pe = [tol(&mut rng), tol(&mut rng), tol(&mut rng)];
let main = is(pa, pb, pc, pd, pe);
if main == 0.0 {
for [a, b, c, d, e] in QuickPerm5::new([pa, pb, pc, pd, pe]) {
assert_eq!(is(a, b, c, d, e), 0.0);
}
}
let pred = main > 0.0;
for (i, [a, b, c, d, e]) in QuickPerm5::new([pa, pb, pc, pd, pe]).enumerate() {
let t = is(a, b, c, d, e);
assert_eq!(
pred,
if i % 2 == 0 { t > 0.0 } else { t < 0.0 },
"{} vs. {} at {:?}",
t,
-main,
pe
);
}
}
}
}
}