Floating point precision for log10

[Jegge]

Hi there,

while experimenting with Uiua, I found that e.g. ₙ10 1000 gives the incorrect result of 2.9999999999999996 instead of 3, caused by a rounding error:

ₙ10 1 = 0
ₙ10 10 = 1 
ₙ10 100 = 2
ₙ10 1000 = 2.9999999999999996
ₙ10 10000 = 4
ₙ10 100000 = 5
ₙ10 1000000 = 5.999999999999999
ₙ10 10000000 = 7
ₙ10 100000000 = 8
ₙ10 1000000000 = 8.999999999999998

The following test case should uncover the problem:

∵(⍤"ₙ10 failed!" =ₙ10 ⁿ∶10 .) ⇡11

This error can be reproduced in the pad as well as in a locally installed interpreter.

It seems that there is an Issue with Rust's underlying f64/f32 log and that there exists an integer variant that circumvents/fixes the problem. However,

uiua/src/algorithm/pervade.rs

Line 884 in f19ec8b

b.log(a),

seems to define the binary operation based on the floating point type f64. I suppose a possible solution would be a special implementation that explicitly uses ilog instead of log in the case that the supplied value is integral? Unfortunately I am a complete noob concerning Rust, so I can not provide you with a pull request. But I hope my analysis helps 😁

Jegge

[SybelBlue]

Hi Jegge!

tl;dr: Short of implementing a new log function, implementing symbolic algebra optimizations, or switching to a new floating/fixed point type, the root problem cannot be fixed. A practical solution might be to round displayed outputs to an agreeable decimal place, and have slightly more flexible floating point equality checking behind the scenes.

I've done some investigating, and it looks like the default number type is f64. This means that a drop-in replacement won't work because ilog only works on integer types. I did attempt your proposed solution though, I've copied my work down below.

As referenced in the rust version of this issue, the inaccuracy boils down to the reality that IEEE doesn't have a standard for this because there are many trade-offs when designing intrinsic math functions. Rust's implementation used b.ln()/a.ln(), which ends in an inherently inaccurate floating point division.

I think a practical solution to this would be rounding floats in the output to a reasonable decimal place (and adding a feature/flag to disable this behavior). Python took this approach, and standardized rounding displayed output as of v3.1 according to the docs. They also provide a math.isclose to avoid inequalities behind the scenes. Rust's own assert_float_eq module is equipped with many near-equality checkers.

A Partial Solution based on Your Suggestion

I expanded the macro you referenced and added ilog for the u8 byte type, but this is useless unless byte arithmetic is enabled for log which I don't yet know how to do. What's more, a single byte wouldn't even be able to do log10(1000).

My implementation of the log module is below:

pub mod log {
    use super::*;
    #[cfg(feature = "bytes")]
    pub fn byte_byte(a: u8, b: u8) -> f64 {
        b.ilog(a) as f64
    }
    #[cfg(feature = "bytes")]
    pub fn byte_num(a: u8, b: f64) -> f64 {
        let a = f64::from(a);
        b.log(a)
    }
    #[cfg(feature = "bytes")]
    pub fn num_byte(a: f64, b: u8) -> f64 {
        let b = f64::from(b);
        b.log(a)
    }
    pub fn num_num(a: f64, b: f64) -> f64 {
        b.log(a)
    }
    #[cfg(feature = "complex")]
    pub fn com_x(a: Complex, b: impl Into<Complex>) -> Complex {
        let b = b.into();
        b.log(a)
    }
    #[cfg(feature = "complex")]
    pub fn x_com(a: impl Into<Complex>, b: Complex) -> Complex {
        let a = a.into();
        b.log(a)
    }
    pub fn error<T: Display>(a: T, b: T, env: &Uiua) -> UiuaError {
        env.error(format!("Cannot get the log base {0} of {1}", b, a))
    }
}

[Jegge]

Hi SybeBlue,

thank you for your reply. Is #[cfg(feature = "bytes")] a compile-time option? To be honest, I do not completely understand your proposed solution. But in the meantime I've also given this a thought and I whipped up two variants of my initial solution. I've tried to format them rust-like, but since I've never coded in this language, don't expect them to compile ^^ :

The first variant checks for an integral a and casts to i64 if necessary. For the sake of simplicity I've used == for the f64, I am aware that maybe something as float_eq might be needed.

pub fn num_num(a: f64, b: f64) -> f64 {
  if (a.fract() > 0) {
    a.log(b)
  } else if (b == 2) {
   (a as i64).ilog2() as f64
  } else if (b == 10) {
   (a as i64).ilog10() as f64
  } else {
   (a as i64).ilog(b as i64) as f64
  }
}

The second variant should be slightly faster than the first. According to the docs f64::log2 and f64::log10 improve the precision for these special bases and should be used instead of log:

The result might not be correctly rounded owing to implementation details; self.log2() can produce more accurate results for base 2, and self.log10() can produce more accurate results for base 10.

This variant might be a reasonable middle ground:

pub fn num_num(a: f64, b: f64) -> f64 {
  if (b == 2) {
   a.log2()
  } else if (b == 10) {
   a.log10()
  } else {
   a.log(b)
  }
}

[kaikalii]

One thing I would be worried about with solutions that do checks like this is speed.
Adding branches to a pervasive operation will slow it down, since it has to check on every element.
How much slower it would get I'm not sure.

[Jegge]

Well I think that the speed penalty should not be an issue, since it's a simple test, especially in variant two. If I were to choose, I would prefer correctness over speed.

[kaikalii]

I've noticed another problem with some of the implementations above.
ilog rounds, which isn't appropriate.