incorrect implementation of prime field $\mathbb{Z}/p$ when p is large

[lampretl]

In source code

Ripserer.jl/src/base/primefield.jl

Line 75 in 3f2d63f

@eval (Base.$op)(i::Mod{M}, j::Mod{M}) where {M} = Mod{M}($op(Int(i), Int(j)))

we see that Int overflow can happen before the modulo operation % is applied, which leads to incorrect results:

using Ripserer, Primes
p=9223372036854775783
print(isprime(p))
x, y = Mod{p}(p-1), Mod{p}(p-2)
(x+y).value, p-3  # we see these are different, even though they should be equal: (p-1)+(p-2) ≡ 2p-3 ≡ p-3 (mod p)

The current implementation of + and * operations gives a correct result when 1<p<sqrt(typemax(Int)), as this inequality makes sure that no overflow happens. Perhaps you could add this assertion somewhere in your code.

Otherwise, a better implementation of integers modulo m can be found in this issue of mine.

Does the original Ripser in C++ have such an implementation of integers modulo m? Would you be so kind and point me to the line in code where I could find this?

[mtsch]

Nice catch thanks! I always assumed people will use integers modulo a small number (maybe up to 17 or so), so integer overflow never even crossed my mind.
I wasn't aware of the Mods package. I'll check to see how fast it is for my specific usecase. It would be good to depend on a high quality implementation rather than rolling my own.

In the meantime, ripserer should work fine with coefficients from that package. You can run it that way by providing the field keyword to it.

ripserer(...; field=Mods.Mod{M})

[mtsch]

In ripser, the inverses get precomputed and stored in a vector, so it only really works for very small p.

I think these should be the relevant lines
https://github.com/Ripser/ripser/blob/01add51ff64aaf40889483260cc5c3b7d0f2a1e7/ripser.cpp#L127
https://github.com/Ripser/ripser/blob/01add51ff64aaf40889483260cc5c3b7d0f2a1e7/ripser.cpp#L762

[lampretl]

Regarding the Mods package, its performance is slow. I benchmarked my implementation with several others in the issue, and my simple code proved to be much faster than the others. Feel free to use my implementation.