2^65 modulo 101 incorrect answer

Question

This code checks that the value a maps uniquely for the values 1 to 100 using the formula (a^x) % 101

local function f(a)
    found = {}
    bijective = true

    for x = 1, 100 do
        value = (a^x) % 101
        if found[value] then
            bijective = false
            break
        else
            found[value] = x
        end
    end

    return bijective
end

However does not produce the expected result. it maps 2^65 % 101 to 56, which matches the value produced by 2^12 % 101 and I get a false result, however the correct value for 2^65 % 101 is 57 and 2 actually should produce all unique values resulting in a true result.

The error described above is specifically on Lua 5.1, is this just a quirk of Lua's number typing? Is there a way to make this function work correctly in 5.1?

Answer 1

The error described above is specifically on Lua 5.1, is this just a quirk of Lua's number typing? Is there a way to make this function work correctly in 5.1?

First of all, this is not an issue with Lua's number typing since 2^65, being a (rather small) power of two, can be represented exactly by the double precision since it uses an exponent-mantissa representation. The mantissa can simply be set to all zeroes (leading one is implicit) and the exponent must be set to 65 (+ offset).

I tried this on different Lua versions and PUC Lua 5.1 & 5.2 as well as LuaJIT have the issue; Lua 5.3 (and presumably later versions as well) are fine. Interestingly, using math.fmod(2^65, 101) returns the correct result on the older Lua versions but 2^65 % 101 does not (it returns 0 instead).

This surprised me so I dug in the Lua 5.1 sources. This is the implementation of math.fmod:

#include <math.h>

...

static int math_fmod (lua_State *L) {
  lua_pushnumber(L, fmod(luaL_checknumber(L, 1), luaL_checknumber(L, 2)));
  return 1;
}

this also is the only place where fmod from math.h appears to be used. The % operator on the other hand is implemented as documented in the reference manual:

#define luai_nummod(a,b)    ((a) - floor((a)/(b))*(b)) 

in src/luaconf.h. You could trivially redefine it as fmod(a,b) to fix your issue. In fact Lua 5.4 does something similar and even provides an elaborate explanation in its sources!

/*
** modulo: defined as 'a - floor(a/b)*b'; the direct computation
** using this definition has several problems with rounding errors,
** so it is better to use 'fmod'. 'fmod' gives the result of
** 'a - trunc(a/b)*b', and therefore must be corrected when
** 'trunc(a/b) ~= floor(a/b)'. That happens when the division has a
** non-integer negative result: non-integer result is equivalent to
** a non-zero remainder 'm'; negative result is equivalent to 'a' and
** 'b' with different signs, or 'm' and 'b' with different signs
** (as the result 'm' of 'fmod' has the same sign of 'a').
*/
#if !defined(luai_nummod)
#define luai_nummod(L,a,b,m)  \
  { (void)L; (m) = l_mathop(fmod)(a,b); \
    if (((m) > 0) ? (b) < 0 : ((m) < 0 && (b) > 0)) (m) += (b); }
#endif

---

Is there a way to make this function work correctly in 5.1?

Yes: The easy way is to use fmod. This may work for these particular numbers since they still fit in doubles due to the base being 2 and the exponent being moderately small, but it won't work in the general case. The better approach is to leverage modular arithmetics to keep your intermediate results small, never storing numbers significantly larger than 101^2 since (a * b) % c == (a % c) * (b % c).

local function f(a)
    found = {}
    bijective = true

    local value = 1
    for _ = 1, 100 do
        value = (value * a) % 101 -- a^x % 101
        if found[value] then
            bijective = false
            break
        else
            found[value] = x
        end
    end

    return bijective
end