#[derive(Copy, Clone, Debug)]
pub struct Ball {
pub m: f64,
pub r: f64,
}
#[inline]
fn two_sum(a: f64, b: f64) -> (f64, f64) {
let s = a + b;
let bb = s - a;
let err = (a - (s - bb)) + (b - bb);
(s, err)
}
#[inline]
fn two_prod(a: f64, b: f64) -> (f64, f64) {
let p = a * b;
let err = f64::mul_add(a, b, -p);
(p, err)
}
impl Ball {
#[inline]
pub fn from_f64(x: f64) -> Self {
Ball { m: x, r: 0.0 }
}
#[inline]
pub fn unknown() -> Self {
Ball {
m: 0.0,
r: f64::INFINITY,
}
}
#[inline]
fn widen_ulp(x: f64) -> f64 {
let up = f64::from_bits(x.to_bits() + 1);
(up - x).abs()
}
#[inline]
pub fn add(self, o: Self) -> Self {
let (s, e) = two_sum(self.m, o.m);
Ball {
m: s,
r: self.r + o.r + e.abs() + Self::widen_ulp(s),
}
}
#[inline]
pub fn sub(self, o: Self) -> Self {
self.add(Ball { m: -o.m, r: o.r })
}
#[inline]
pub fn neg(self) -> Self {
Ball {
m: -self.m,
r: self.r,
}
}
#[inline]
pub fn mul(self, o: Self) -> Self {
let (p, e) = two_prod(self.m, o.m);
Ball {
m: p,
r: self.m.abs() * o.r
+ o.m.abs() * self.r
+ self.r * o.r
+ e.abs()
+ Self::widen_ulp(p),
}
}
#[inline]
fn recip(self) -> Self {
if self.m.abs() <= self.r {
return Ball::unknown();
}
let denom = self.m.abs() * (self.m.abs() - self.r);
let m = 1.0 / self.m;
let r = self.r / denom; Ball { m, r }
}
#[inline]
pub fn div(self, o: Self) -> Self {
self.mul(o.recip())
}
#[inline]
pub fn sign_if_certain(self) -> Option<i8> {
if self.r.is_infinite() {
return None;
}
if self.m > self.r {
Some(1)
} else if self.m < -self.r {
Some(-1)
} else {
None
}
}
}