Different results with floor(a, b) and a // b

Question

--this function isn't relevant, but I included it for completeness
function gcd(a, b)
    local temp
    while b > 0 do
        temp = a % b
        a = b
        b = temp
    end
    return a
end

function extendedgcd(a, b)
    if b == 0 then
        return 1, 0, a
    end
    local x, y, d = extendedgcd(b, a % b)
    --this assertion fails for a = 568784642688024576, b = 5
    --left side gives 113756928537604915 (correct), while right side gives 113756928537604912 (wrong)
    assert(a // b == math.floor(a / b))
    --so here, I can't use math.floor
    return y, x - a // b * y, d
end

function modularinverse(e, mod)
    if gcd(e, mod) > 1 then
        return nil
    end
    local x, _, _ = extendedgcd(e, mod)
    if x < 0 then
        x = x % mod
    end
    if x * e % mod ~= 1 then
        error(string.format("Modular inverse (%d) doesn't produce 1 (e = %d, mod = %d)", x, e, mod))
    end
    return x
end

modularinverse(5, 568784642688024576)

For a = 5 and b = 568784642688024576, the assertion in extendedgcd fails. I'm not an expert at FP precision, but there's a difference of 3 between the two so I don't think there's a rounding/precision error here. But I might be wrong.

Normally I'd just use // instead, but I can't because the target platform isn't running Lua 5.3, which is when the operator was added.

What am I missing? How do I make it work with floor, or is there another way?

Also of note: The same issue happened when I rewrote it in Python to fiddle (with math.floor and //).

Answer 1

A more precise answer:

As paper and pencil will show: 568784642688024576 / 5 = 113756928537604915.02

The most accurate binary representation of that quotient as a double precision number is:

0100001101111001010000100101010100101110010000111001001100110011

Which, in decimal, is: 1.13756928537604912E17 (notice the ...912 ending)

Now, if you decrease the binary representation by one, like this:

0100001101111001010000100101010100101110010000111001001100110010

Then it equals: 1.13756928537604896E17 (...896 at the end!)

If you were to increase the original binary number by one, like this:

0100001101111001010000100101010100101110010000111001001100110100

Then that equals: 1.13756928537604928E17 (...928 at the end!)

So, there are exact double precision binary representations for these numbers:

113756928537604896

113756928537604912 (closest to actual quotient)

113756928537604928

The above can be verified using online converters, such as here.

So the lesson is:

Integer division will give the correct answer (i.e. the integer portion of the quotient). Flooring floating-point division relies on binary representation of numbers, which is not always precise.

Outside of the scope of this answer, but required reading to really understand any of this:

How do the above binary numbers represent their decimal equivalents?

• Short answer: IEEE-754. I wish you luck and lot’s of coffee if you plan on studying that standards document.

What is the algorithmic difference between / and //?

• In other words, why does // get us the answer we want above?