Skip to content

Commit

Permalink
Cap the precision estimate to the max precision
Browse files Browse the repository at this point in the history 
- Print one unstable digit by default
  • Loading branch information
@benruijl
benruijl committed Jun 10, 2024
1 parent a9bf19c commit fbe5c9a
Showing 1 changed file with 67 additions and 1 deletion.
68 changes: 67 additions & 1 deletion src/domains/float.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ pub trait NumericalFloatLike:

/// Get the number of precise binary digits.
fn get_precision(&self) -> u32;
fn get_epsilon(&self) -> f64;

/// Sample a point on the interval [0, 1].
fn sample_unit<R: Rng + ?Sized>(&self, rng: &mut R) -> Self;
Expand Down Expand Up @@ -150,10 +151,16 @@ impl NumericalFloatLike for f64 {
a as f64
}

#[inline(always)]
fn get_precision(&self) -> u32 {
53
}

#[inline(always)]
fn get_epsilon(&self) -> f64 {
f64::EPSILON / 2.
}

fn sample_unit<R: Rng + ?Sized>(&self, rng: &mut R) -> Self {
rng.gen()
}
Expand Down Expand Up @@ -668,6 +675,11 @@ impl NumericalFloatLike for Float {
self.prec()
}

#[inline(always)]
fn get_epsilon(&self) -> f64 {
2.0f64.powi(-(self.prec() as i32))
}

fn sample_unit<R: Rng + ?Sized>(&self, rng: &mut R) -> Self {
let f: f64 = rng.gen();
Float::with_val(self.prec(), f)
Expand Down Expand Up @@ -997,6 +1009,13 @@ impl<T: NumericalFloatLike> ErrorPropagatingFloat<T> {
pub fn get_precision(&self) -> f64 {
-self.prec.log10()
}

/// Truncate the precision to the maximal number of stable decimal digits
/// of the underlying float.
pub fn truncate(mut self) -> Self {
self.prec = self.prec.max(self.value.get_epsilon());
self
}
}

impl<T: NumericalFloatLike> fmt::Display for ErrorPropagatingFloat<T> {
Expand Down Expand Up @@ -1108,6 +1127,10 @@ impl<T: NumericalFloatComparison> NumericalFloatLike for ErrorPropagatingFloat<T
self.value.get_precision()
}

fn get_epsilon(&self) -> f64 {
2.0f64.powi(-(self.get_precision() as i32))
}

fn sample_unit<R: Rng + ?Sized>(&self, rng: &mut R) -> Self {
ErrorPropagatingFloat {
value: self.value.sample_unit(rng),
Expand Down Expand Up @@ -1154,6 +1177,7 @@ impl<T: Real + NumericalFloatComparison> Real for ErrorPropagatingFloat<T> {
value: self.value.sqrt(),
prec: self.prec / 2.,
}
.truncate()
}

fn log(&self) -> Self {
Expand All @@ -1162,27 +1186,31 @@ impl<T: Real + NumericalFloatComparison> Real for ErrorPropagatingFloat<T> {
prec: self.prec / r.clone().to_f64().abs(),
value: r,
}
.truncate()
}

fn exp(&self) -> Self {
ErrorPropagatingFloat {
value: self.value.exp(),
prec: self.value.clone().to_f64().abs() * self.prec,
}
.truncate()
}

fn sin(&self) -> Self {
ErrorPropagatingFloat {
prec: self.prec * self.value.clone().to_f64().abs() / self.value.tan().to_f64().abs(),
value: self.value.sin(),
}
.truncate()
}

fn cos(&self) -> Self {
ErrorPropagatingFloat {
prec: self.prec * self.value.clone().to_f64().abs() * self.value.tan().to_f64().abs(),
value: self.value.cos(),
}
.truncate()
}

fn tan(&self) -> Self {
Expand All @@ -1192,6 +1220,7 @@ impl<T: Real + NumericalFloatComparison> Real for ErrorPropagatingFloat<T> {
prec: self.prec * self.value.clone().to_f64().abs() * (tt.inv() + tt),
value: t,
}
.truncate()
}

fn asin(&self) -> Self {
Expand All @@ -1202,6 +1231,7 @@ impl<T: Real + NumericalFloatComparison> Real for ErrorPropagatingFloat<T> {
prec: self.prec * v.abs() / tt,
value: t,
}
.truncate()
}

fn acos(&self) -> Self {
Expand All @@ -1212,6 +1242,7 @@ impl<T: Real + NumericalFloatComparison> Real for ErrorPropagatingFloat<T> {
prec: self.prec * v.abs() / tt,
value: t,
}
.truncate()
}

fn atan2(&self, x: &Self) -> Self {
Expand All @@ -1224,20 +1255,23 @@ impl<T: Real + NumericalFloatComparison> Real for ErrorPropagatingFloat<T> {
prec: r.prec * r2 / tt,
value: t,
}
.truncate()
}

fn sinh(&self) -> Self {
ErrorPropagatingFloat {
prec: self.prec * self.value.clone().to_f64().abs() / self.value.tanh().to_f64().abs(),
value: self.value.sinh(),
}
.truncate()
}

fn cosh(&self) -> Self {
ErrorPropagatingFloat {
prec: self.prec * self.value.clone().to_f64().abs() * self.value.tanh().to_f64().abs(),
value: self.value.cosh(),
}
.truncate()
}

fn tanh(&self) -> Self {
Expand All @@ -1247,6 +1281,7 @@ impl<T: Real + NumericalFloatComparison> Real for ErrorPropagatingFloat<T> {
prec: self.prec * self.value.clone().to_f64().abs() * (tt.inv() - tt),
value: t,
}
.truncate()
}

fn asinh(&self) -> Self {
Expand All @@ -1257,6 +1292,7 @@ impl<T: Real + NumericalFloatComparison> Real for ErrorPropagatingFloat<T> {
prec: self.prec * v.abs() / tt,
value: t,
}
.truncate()
}

fn acosh(&self) -> Self {
Expand All @@ -1267,6 +1303,7 @@ impl<T: Real + NumericalFloatComparison> Real for ErrorPropagatingFloat<T> {
prec: self.prec * v.abs() / tt,
value: t,
}
.truncate()
}

fn atanh(&self) -> Self {
Expand All @@ -1277,6 +1314,7 @@ impl<T: Real + NumericalFloatComparison> Real for ErrorPropagatingFloat<T> {
prec: self.prec * v.abs() / tt,
value: t,
}
.truncate()
}

fn powf(&self, e: Self) -> Self {
Expand All @@ -1285,6 +1323,7 @@ impl<T: Real + NumericalFloatComparison> Real for ErrorPropagatingFloat<T> {
value: self.value.powf(e.value.clone()),
prec: (self.prec + e.prec * v.ln().abs()) * e.value.clone().to_f64().abs(),
}
.truncate()
}
}

Expand Down Expand Up @@ -1343,10 +1382,16 @@ macro_rules! simd_impl {
Self::from(a as f64)
}

#[inline(always)]
fn get_precision(&self) -> u32 {
53
}

#[inline(always)]
fn get_epsilon(&self) -> f64 {
f64::EPSILON / 2.
}

fn sample_unit<R: Rng + ?Sized>(&self, rng: &mut R) -> Self {
Self::from(rng.gen::<f64>())
}
Expand Down Expand Up @@ -1500,10 +1545,16 @@ impl NumericalFloatLike for Rational {
Rational::Natural(a, 1)
}

#[inline(always)]
fn get_precision(&self) -> u32 {
u32::MAX
}

#[inline(always)]
fn get_epsilon(&self) -> f64 {
0.
}

fn sample_unit<R: Rng + ?Sized>(&self, rng: &mut R) -> Self {
let rng1 = rng.gen::<i64>();
let rng2 = rng.gen::<i64>();
Expand Down Expand Up @@ -1892,10 +1943,16 @@ impl<T: Real> NumericalFloatLike for Complex<T> {
}
}

#[inline(always)]
fn get_precision(&self) -> u32 {
self.re.get_precision().min(self.im.get_precision())
}

#[inline(always)]
fn get_epsilon(&self) -> f64 {
(2.0f64).powi(-(self.get_precision() as i32))
}

fn sample_unit<R: Rng + ?Sized>(&self, rng: &mut R) -> Self {
Complex {
re: self.re.sample_unit(rng),
Expand Down Expand Up @@ -2072,7 +2129,16 @@ mod test {
+ b.acosh() / f.atanh()
+ b.powf(a);
assert_eq!(r.value, 17293.219725825093);
assert_eq!(r.get_precision(), 14.836811363436391);
// error is 14.836811363436391 when the f64 could have theoretically grown in between
assert_eq!(r.get_precision(), 14.836795991431746);
}

#[test]
fn error_truncation() {
let a = ErrorPropagatingFloat::new(0.0000000123456789, 9.)
.exp()
.log();
assert_eq!(a.get_precision(), 8.046104745509947);
}

#[test]
Expand Down

0 comments on commit fbe5c9a

Please sign in to comment.