use super::{Point2, Point3, Vector3};
use crate::features::FiniteVector3;
use crate::scalar::FiniteReal;
const EPS_UNIT_RESULT: f64 = 8.0 * f64::EPSILON;
#[test]
fn a_normal_power_of_two_is_the_power_powi_forms() {
for exponent in -1100..=1100 {
let formed = super::normal_power_of_two(exponent).map(f64::to_bits);
if (-1022..=1023).contains(&exponent) {
assert_eq!(
formed,
Some(2.0_f64.powi(exponent).to_bits()),
"2^{exponent}"
);
} else {
assert_eq!(formed, None, "2^{exponent}");
}
}
for exponent in [i32::MIN, i32::MAX] {
assert_eq!(super::normal_power_of_two(exponent), None);
}
}
fn power_of_two_product_by_powi(mut value: f64, mut exponent: i32) -> f64 {
while exponent > 1023 {
value *= 2.0_f64.powi(1023);
if !value.is_finite() {
return value;
}
exponent -= 1023;
}
while exponent < -1022 {
value *= 2.0_f64.powi(-969);
if value == 0.0 {
break;
}
exponent += 969;
}
value * 2.0_f64.powi(exponent)
}
fn assert_power_of_two_product_matches(value: f64, exponent: i32) {
assert_eq!(
super::power_of_two_product(value, exponent).to_bits(),
power_of_two_product_by_powi(value, exponent).to_bits(),
"{value:e} * 2^{exponent}"
);
}
#[test]
fn a_power_of_two_product_is_the_product_powi_forms_at_every_exponent() {
let values = [
0.0,
f64::MIN_POSITIVE,
f64::MIN_POSITIVE / 3.0,
f64::from_bits(1),
0.75,
1.0,
f64::MAX,
f64::INFINITY,
f64::NAN,
];
for value in values.into_iter().flat_map(|value| [value, -value]) {
for exponent in -2200..=2200 {
assert_power_of_two_product_matches(value, exponent);
}
}
}
proptest::proptest! {
#[test]
fn a_power_of_two_product_is_the_product_powi_forms(
bits in proptest::prop_oneof![
proptest::num::u64::ANY,
0..1_u64 << 53,
(1_u64 << 63)..(1_u64 << 63) | (1_u64 << 53),
],
exponent in -2200..=2200_i32,
) {
assert_power_of_two_product_matches(f64::from_bits(bits), exponent);
}
}
#[test]
fn numerical_audit_translation_retains_finite_cancellation() {
assert_eq!(
Point3::new(-1e308, 0.0, 0.0).translated(Vector3::new(1e308, 0.0, 0.0), 2.0),
Point3::new(1e308, 0.0, 0.0)
);
}
#[test]
fn numerical_audit_dot_retains_finite_cancellation() {
assert_eq!(
Vector3::new(1e308, 1e308, 0.0).dot(Vector3::new(2.0, -2.0, 0.0)),
0.0
);
}
#[test]
fn numerical_audit_cross_of_parallel_large_vectors_is_zero() {
let vector = Vector3::new(1e200, 1e200, 1e200);
assert_eq!(vector.cross(vector), Vector3::new(0.0, 0.0, 0.0));
}
#[test]
fn numerical_audit_unit_preserves_representable_subnormal_component() {
let tiny = f64::from_bits(1);
let vector = Vector3::new(1.0, tiny, 0.0);
assert_eq!(
FiniteVector3::new(vector)
.and_then(FiniteVector3::unit_nonzero)
.unwrap()
.y
.to_bits(),
tiny.to_bits()
);
assert_eq!(vector.unit().unwrap().y.to_bits(), tiny.to_bits());
}
#[test]
fn numerical_audit_squared_distance_accumulates_subnormal_squares() {
let component = 2.0_f64.powi(-538);
let point = Point3::new(component, component, component);
assert_eq!(
point.distance_squared(Point3::new(0.0, 0.0, 0.0)).to_bits(),
1
);
}
#[test]
fn point_is_finite_rejects_a_nonfinite_coordinate_in_any_slot() {
assert!(Point3::new(1.0, -2.0, 3.0).is_finite());
for value in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
assert!(!Point3::new(value, 1.0, 1.0).is_finite());
assert!(!Point3::new(1.0, value, 1.0).is_finite());
assert!(!Point3::new(1.0, 1.0, value).is_finite());
}
}
#[test]
fn vector_is_finite_rejects_a_nonfinite_component_in_any_slot() {
assert!(Vector3::new(1.0, -2.0, 3.0).is_finite());
for value in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
assert!(!Vector3::new(value, 1.0, 1.0).is_finite());
assert!(!Vector3::new(1.0, value, 1.0).is_finite());
assert!(!Vector3::new(1.0, 1.0, value).is_finite());
}
}
#[test]
fn parameter_point_is_finite_rejects_a_nonfinite_coordinate_in_either_slot() {
assert!(Point2::new(1.0, -2.0).is_finite());
for value in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
assert!(!Point2::new(value, 1.0).is_finite());
assert!(!Point2::new(1.0, value).is_finite());
}
}
#[test]
fn unit_vector_preserves_finite_direction_without_squared_length_overflow() {
for scale in [1.0, 2.0_f64.powi(600), f64::MAX] {
for (input, expected) in [
(Vector3::new(scale, 0.0, 0.0), Vector3::new(1.0, 0.0, 0.0)),
(Vector3::new(0.0, -scale, 0.0), Vector3::new(0.0, -1.0, 0.0)),
(Vector3::new(0.0, 0.0, scale), Vector3::new(0.0, 0.0, 1.0)),
(
Vector3::new(scale, -scale, scale),
Vector3::new(1.0, -1.0, 1.0).scale(1.0 / 3.0_f64.sqrt()),
),
] {
let actual = input.unit().expect("finite nondegenerate direction");
assert!((actual.x - expected.x).abs() <= EPS_UNIT_RESULT);
assert!((actual.y - expected.y).abs() <= EPS_UNIT_RESULT);
assert!((actual.z - expected.z).abs() <= EPS_UNIT_RESULT);
assert!((actual.norm() - 1.0).abs() <= EPS_UNIT_RESULT);
}
}
}
#[test]
fn unit_vector_preserves_the_epsilon_degeneracy_boundary() {
for x in [0.0, f64::from_bits(1), f64::MIN_POSITIVE, f64::EPSILON] {
assert!(Vector3::new(x, 0.0, 0.0).unit().is_none());
assert!(Vector3::new(0.0, -x, 0.0).unit().is_none());
}
let above = f64::from_bits(f64::EPSILON.to_bits() + 1);
assert_eq!(
Vector3::new(above, 0.0, 0.0).unit(),
Some(Vector3::new(1.0, 0.0, 0.0))
);
let diagonal = Vector3::new(f64::EPSILON, f64::EPSILON, 0.0)
.unit()
.unwrap();
assert!((diagonal.x - 1.0 / 2.0_f64.sqrt()).abs() <= EPS_UNIT_RESULT);
assert!((diagonal.y - diagonal.x).abs() <= EPS_UNIT_RESULT);
assert_eq!(diagonal.z, 0.0);
}
#[test]
fn unit_vector_refuses_each_nonfinite_component() {
for value in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
for input in [
Vector3::new(value, 1.0, 1.0),
Vector3::new(1.0, value, 1.0),
Vector3::new(1.0, 1.0, value),
] {
assert!(input.unit().is_none());
}
}
}
#[test]
fn numerical_audit_lengths_and_nonzero_directions_keep_the_finite_range() {
for magnitude in [f64::from_bits(1), 1.0e-200, 1.0, 1.0e200, f64::MAX] {
let vector = Vector3::new(magnitude, 0.0, 0.0);
assert_eq!(vector.norm(), magnitude);
assert_eq!(
Point3::new(magnitude, 0.0, 0.0).distance(Point3::new(0.0, 0.0, 0.0)),
magnitude
);
assert_eq!(
FiniteVector3::new(vector).and_then(FiniteVector3::unit_nonzero),
Some(Vector3::new(1.0, 0.0, 0.0))
);
}
let diagonal = FiniteVector3::new(Vector3::new(f64::MAX, f64::MAX, 0.0))
.and_then(FiniteVector3::unit_nonzero)
.unwrap();
assert!((diagonal.norm() - 1.0).abs() <= EPS_UNIT_RESULT);
assert!(FiniteVector3::ZERO.unit_nonzero().is_none());
assert!(FiniteVector3::new(Vector3::new(f64::INFINITY, 0.0, 0.0))
.and_then(FiniteVector3::unit_nonzero)
.is_none());
}
#[test]
fn numerical_audit_parameter_reflection_preserves_shifted_endpoints() {
let finite = |value: f64| FiniteReal::new(value).unwrap();
for [start, end] in [
[1.0e16, 1.0e16 + 2.0],
[1.0e308, 1.1e308],
[-f64::MAX, f64::MAX],
] {
let [start, end] = [finite(start), finite(end)];
assert_eq!(super::reflect_parameter(start, start, end), Some(end));
assert_eq!(super::reflect_parameter(end, start, end), Some(start));
}
assert_eq!(
super::reflect_parameter(finite(0.0), finite(-f64::MAX), finite(f64::MAX))
.map(FiniteReal::get),
Some(0.0)
);
assert!(super::reflect_parameter(finite(-f64::MAX), finite(0.0), finite(f64::MAX)).is_none());
}
#[test]
fn numerical_audit_product_quotient_preserves_representable_results() {
let finite = |value: f64| FiniteReal::new(value).unwrap();
assert_eq!(
super::multiply_divide(finite(0.0), finite(1.0e200), finite(1.0e-200)).map(FiniteReal::get),
Some(0.0)
);
for value in [f64::from_bits(1), 1.0e-200, 1.0, 1.0e200, f64::MAX] {
let result = super::multiply_divide(finite(value), finite(value), finite(value))
.unwrap()
.get();
assert!((result / value - 1.0).abs() <= 4.0 * f64::EPSILON);
}
assert!(super::multiply_divide(FiniteReal::ONE, FiniteReal::ONE, FiniteReal::ZERO).is_none());
}
#[test]
fn numerical_ranges_products_keep_intermediate_exponents() {
use super::{product_quotient, scaled_sinh_cosh};
assert_eq!(
product_quotient([1e200, 1e200], [4.0, 1e200]).map(FiniteReal::get),
Some(2.5e199)
);
assert_eq!(
product_quotient([2.0, 1e308, 0.5], []).map(FiniteReal::get),
Some(1e308)
);
assert_eq!(product_quotient([1.0], [0.0]), None);
assert_eq!(product_quotient([f64::NAN], [1.0]), None);
assert_eq!(product_quotient([f64::MAX, 2.0], []), None);
for parameter in [-720.0_f64, 720.0] {
let expected = (parameter.abs() + 1e-10_f64.ln()).exp() * 0.5;
let (sinh, cosh) = scaled_sinh_cosh(
FiniteReal::new(1e-10).unwrap(),
FiniteReal::new(parameter).unwrap(),
)
.unwrap();
let (sinh, cosh) = (sinh.get(), cosh.get());
assert!((cosh / expected - 1.0).abs() < 1024.0 * f64::EPSILON);
assert_eq!(sinh, parameter.signum() * cosh);
}
assert_eq!(
scaled_sinh_cosh(FiniteReal::ZERO, FiniteReal::new(2000.0).unwrap())
.ok()
.map(|(sinh, cosh)| (sinh.get(), cosh.get())),
Some((0.0, 0.0))
);
assert_eq!(
scaled_sinh_cosh(FiniteReal::ONE, FiniteReal::new(2000.0).unwrap()).ok(),
None
);
assert_eq!(
scaled_sinh_cosh(FiniteReal::new(2.0).unwrap(), FiniteReal::ZERO)
.ok()
.map(|(sinh, cosh)| (sinh.get(), cosh.get())),
Some((0.0, 2.0))
);
}
#[test]
fn numerical_ranges_interpolation_and_wrapping_avoid_endpoint_subtraction_overflow() {
use super::{interpolate, wrap_parameter};
assert_eq!(
interpolate(-1e308, 1e308, 0.5).map(FiniteReal::get),
Some(0.0)
);
assert_eq!(
interpolate(-1e308, 1e308, 0.0).map(FiniteReal::get),
Some(-1e308)
);
assert_eq!(
interpolate(-1e308, 1e308, 1.0).map(FiniteReal::get),
Some(1e308)
);
assert_eq!(interpolate(0.0, 1.0, f64::NAN), None);
assert_eq!(
wrap_parameter(1e308, -1e308, 0.0).map(FiniteReal::get),
Some(-1e308)
);
assert_eq!(
wrap_parameter(1.5e308, -1e308, 1e308).map(FiniteReal::get),
Some(-5e307)
);
assert_eq!(
wrap_parameter(1.0, -1.0, 1.0).map(FiniteReal::get),
Some(-1.0)
);
assert_eq!(
wrap_parameter(0.5, -1.0, 1.0).map(FiniteReal::get),
Some(0.5)
);
assert_eq!(wrap_parameter(1.0, 0.0, 0.0), None);
}
#[test]
fn a_power_of_two_bound_normalises_without_moving_a_significand() {
use super::{power_of_two_bound, scale_power_of_two};
for value in [
f64::from_bits(1),
-1.0e-200,
0.5,
1.0,
-6.0,
1.0e200,
f64::MAX,
] {
let exponent = power_of_two_bound(value).unwrap();
let normalised = scale_power_of_two(value, -exponent).unwrap().get();
assert!((0.5..1.0).contains(&normalised.abs()));
assert_eq!(
scale_power_of_two(normalised, exponent).map(FiniteReal::get),
Some(value)
);
}
assert_eq!(power_of_two_bound(1.0), Some(1));
assert_eq!(power_of_two_bound(-6.0), Some(3));
for value in [0.0, -0.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
assert_eq!(power_of_two_bound(value), None);
}
}
#[test]
fn audit_regression_power_scaling_preserves_extreme_finite_results() {
use super::scale_power_of_two;
let minimum = f64::from_bits(1);
let maximum_power = 2.0_f64.powi(1023);
assert_eq!(
scale_power_of_two(minimum, 2097).map(FiniteReal::get),
Some(maximum_power)
);
assert_eq!(
scale_power_of_two(maximum_power, -2097).map(FiniteReal::get),
Some(minimum)
);
assert_eq!(
scale_power_of_two(1.5, -1074).map(FiniteReal::get),
Some(2.0 * minimum)
);
assert_eq!(scale_power_of_two(1.0, i32::MAX), None);
assert_eq!(
scale_power_of_two(1.0, i32::MIN).map(FiniteReal::get),
Some(0.0)
);
assert_eq!(
scale_power_of_two(-0.0, i32::MAX).unwrap().get().to_bits(),
(-0.0_f64).to_bits()
);
assert_eq!(scale_power_of_two(f64::NAN, 0), None);
}
#[test]
fn numerical_0922b_parameter_fraction_preserves_finite_charts() {
for (parameter, start, end, expected) in [
(0.0, -1e308, 1e308, 0.5),
(-1e308, -1e308, 1e308, 0.0),
(1e308, -1e308, 1e308, 1.0),
(-2.0, 0.0, 1.0, -2.0),
(3.0, 0.0, 1.0, 3.0),
] {
assert_eq!(
super::parameter_fraction(parameter, start, end).map(FiniteReal::get),
Some(expected)
);
}
}
#[test]
fn numerical_audit_domains_do_not_admit_disjoint_parameter_ranges() {
for domain in [[0.0_f64, 1.], [0., 1e-16], [-1e308, 1e308]] {
for value in [domain[0], domain[0].midpoint(domain[1]), domain[1]] {
assert!(super::parameter_in_domain(
value,
domain,
64. * f64::EPSILON
));
}
}
for d in [1., 1e-16] {
assert!(!super::parameter_in_domain(
2. * d,
[0., d],
64. * f64::EPSILON
));
assert!(!super::parameter_in_domain(-d, [0., d], 64. * f64::EPSILON));
assert!(super::parameter_in_domain(
d.next_up(),
[0., d],
64. * f64::EPSILON
));
}
assert!(!super::parameter_in_domain(f64::NAN, [0., 1.], 0.));
}
#[test]
fn scaled_hyperbolic_products_that_overflow_carry_their_plain_products() {
use super::scaled_sinh_cosh;
// cosh(1) is finite; its product with MAX is not.
let (sinh, cosh) = scaled_sinh_cosh(FiniteReal::new(f64::MAX).unwrap(), FiniteReal::ONE)
.expect_err("the products overflow");
assert_eq!((sinh, cosh), (f64::INFINITY, f64::INFINITY));
// Past the cosh range the products are exponentiated in thirds.
assert_eq!(
scaled_sinh_cosh(FiniteReal::ONE.negated(), FiniteReal::new(-2000.0).unwrap()),
Err((f64::INFINITY, f64::NEG_INFINITY))
);
}