0

I've been working on a deterministic maths library for LuaJIT, and after research I am aware that multiplication, division, addition, subtraction, and square root are deterministic (as long as the same rounding mode is active) due to the IEEE 754 standard. However I also have a hunch that x^y where y is an integer may also be deterministic, as there is no error-- at least, no error when x is an integer-- when I try it, as if an algorithm like exponentiation by squaring was active. It would make sense if that were the case on all systems. I'm wondering if anyone in the know would have any idea if that's true.

By deterministic I mean giving the exact result on all systems.

Tachytaenius
  • 147
  • 8
  • `pow(x,y)` is not required by IEEE 754, some languages may not even provide `pow(double,int)` but always use something equivalent to `pow(x,y) := exp2(y*log2(x))`. – chtz Oct 10 '22 at 20:29
  • 2
    @chtz From first-hand experience: Not all platforms implement `pow(double,int)` via `exp{2}, log{2}` for all values of the `int` argument. – njuffa Oct 10 '22 at 22:09
  • 5
    @Tachytaenius The answer to the question in the subject line is "no". There are no guarantees that exponentiation is implemented in any particular way, even when restricting to integer exponents. – njuffa Oct 10 '22 at 22:12
  • 1
    In particular: if you evaluate `pow(2, 3)` on two different platforms, you are quite likely to discover that one gives you 7.9999999, and one gives you 8.0. – Steve Summit Oct 11 '22 at 12:15
  • Thanks for your answers, everyone. I will avoid using the `pow` and `^` I am provided for deterministic programs, then. – Tachytaenius Oct 11 '22 at 12:22
  • @SteveSummit any suggestions as to which system(s) would actually do that? I'm two bits out after a naive `exp(log(2)*3)` (with binary64s) but wonder what standard libraries are that bad – Sam Mason Oct 13 '22 at 09:54
  • @SamMason Sorry, no suggestions. Back when I was learning C, all `pow()` implementations were "that bad". So the C Standard has never required anything better. I'm annoyed by this; I wish the Standard *did* require something better. These days, good systems *are* better, but I don't think all systems. I just tried MacOS and Linux (Debian Buster), and they both gave 8.0000000000, so perhaps my "if you evaluate on two different platforms" was a notch too pessimistic. – Steve Summit Oct 14 '22 at 01:52
  • @SteveSummit I tried recent versions of glibc and musl before asking and both seem to have good implementations as well. debian presumably uses glibc by default. the C standard remains an annoying lowest common denominator, would be nice if there were some way of saying I want a conventional modern processor (e.g. 8bit bytes, ieee754 floats) as most code doesn't care about breaking on obscure systems anyway (e.g. your interactive GTK+ GUI isn't going to be run on a DSP) – Sam Mason Oct 14 '22 at 13:18

1 Answers1

1

Is raising a float to an integer power guaranteed to be the same on all ordinary systems?

No.

Consider a function z = pow(x, y) where x, y have integer values (x is odd, y > 0). When z is close, yet less than DBL_MAX, the math correct answer is perhaps a 1000 binary digits long. If code use exponentiation by squaring, there are many roundings that occur in the sub-steps. If code used extended math, the final result would sometimes certainly be better using the same algorithm, yet different than a limited double math approach. Of course other algorithms can generate different results.

Languages rarely specified the exact computational approach for transcendental functions.

It is interesting OP is seeking the same result and not the best result. To find the best result can involve some heavy computation and then risks performance complaints. It is a trade-off that as advances occur, tend to go for the better answer - without too much excessive computation time.

Do you want consistency, or the best answer? Careful what you wish for - and this integer example.

chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256
  • 1
    (+1) A nitpicking quibble: The arguments to `pow()` are rational numbers, and if we restrict to the common case where there is no complex result, exponentiating a rational number by another rational number delivers an *algebraic*, not a *transcendental* result. Basically n-th root of x to the m. Not that I know any standard floating-point library that I have encountered that computes it that way :-) – njuffa Nov 02 '22 at 06:24
  • @njuffa Fair point. – chux - Reinstate Monica Nov 02 '22 at 13:46