I am trying to implement an equation evaluator, but making a nth root function is surprisingly difficult to do properly..
I'm using the typical "f64.powf(1.0 / f64)" to do the root, but it doesn't give the correct results for negative numbers.
E.g.
(-27f64).powf(1.0 / 3.0) should be -3.0, but it's NaN
(-32f64).powf(1.0 / 5.0) should be -2.0, but it is also NaN
Does anyone know of a performant way to calculate these? (The best I've come up with is taking the root.trunc() % 2.0 to determine whether or not I need to .abs() before doing the .powf and then multiplying by signum, but I hope someone has an easier and more performant way to do this..)
Bonus points, does anyone know why rust gives NaN in these cases?
This is consistent with other languages (C++, JavaScript, C#). The reason is because of floating point numbers aren't real numbers with infinite precision, so we can't approximate non-integer powers of negative numbers. i.e. (-27)^0.33333333 = 1.5 + 2.59807621 i? Increasing the precision does not get us closer to the correct answer, so all non-integer powers of negative numbers are NaN.
Tldr: non-integeral roots of negative numbers are complex (pun intended)
What if you unconditionally do the abs-and-multiply-with-sign? That should yield literally something like 2 bit flips (Godbolt):
pub fn root(x: f64, n: u32) -> f64 {
let exp = 1.0 / n as f64;
if (n & 1) == 0 {
x.powf(exp)
} else {
let absroot = x.abs().powf(exp);
if x < 0.0 {
-absroot
} else {
absroot
}
}
}
Actually, I wouldn't even worry about the performance given the context. The single floating point exponentiation is likely to dominate everything else, if not there's still the division (that you also need unconditionally). So my bet is that shaving off a bunch of negative or is-even checks won't get you very far.
I'm sure you are right, I suppose it will just be expensive in any case, thank you for your response! (And thank you for showing me the compiler explorer! I am going to spend some quality time with it)
That is very interesting that other languages behave the same way, I suppose I shouldn't have expected the computer wouldn't treat the rational number in my code as any more precise than the result it's division will yield.