« What would Dijkstra do? Proving the associativity of min » Data structure challenge: finding the rightmost empty slot

Competitive Programming in Haskell: modular arithmetic, part 2

Posted on March 3, 2020
Tagged , , ,

In my last post I wrote about modular exponentiation and egcd. In this post, I consider the problem of solving modular equivalences, building on code from the previous post.

Solving linear congruences

A linear congruence is a modular equivalence of the form

\(ax \equiv b \pmod m\).

Let’s write a function to solve such equivalences for \(x\). We want a pair of integers \(y\) and \(k\) such that \(x\) is a solution to \(ax \equiv b \pmod m\) if and only if \(x \equiv y \pmod k\). This isn’t hard to write in the end, but takes a little bit of thought to do it properly.

First of all, if \(a\) and \(m\) are relatively prime (that is, \(\gcd(a,m) = 1\)) then we know from the last post that \(a\) has an inverse modulo \(m\); multiplying both sides by \(a^{-1}\) yields the solution \(x \equiv a^{-1} b \pmod m\).

OK, but what if \(\gcd(a,m) > 1\)? In this case there might not even be any solutions. For example, \(2x \equiv 3 \pmod 4\) has no solutions: any even number will be equivalent to \(0\) or \(2\) modulo \(4\), so there is no value of \(x\) such that double it will be equivalent to \(3\). On the other hand, \(2x \equiv 2 \pmod 4\) is OK: this will be true for any odd value of \(x\), that is, \(x \equiv 1 \pmod 2\). In fact, it is easy to see that any common divisor of \(a\) and \(m\) must also divide \(b\) in order to have any solutions. In case the GCD of \(a\) and \(m\) does divide \(b\), we can simply divide through by the GCD (including dividing the modulus \(m\)!) and then solve the resulting equivalence.

-- solveMod a b m solves ax = b (mod m), returning a pair (y,k) (with
-- 0 <= y < k) such that x is a solution iff x = y (mod k).
solveMod :: Integer -> Integer -> Integer -> Maybe (Integer, Integer)
solveMod a b m
  | g == 1         = Just ((b * inverse m a) `mod` m, m)
  | b `mod` g == 0 = solveMod (a `div` g) (b `div` g) (m `div` g)
  | otherwise      = Nothing
  where
    g = gcd a m

Solving systems of congruences with CRT

In its most basic form, the Chinese remainder theorem (CRT) says that if we have a system of two modular equations

\(\begin{array}{rcl}x &\equiv& a \pmod m \\ x &\equiv& b \pmod n\end{array}\)

then as long as \(m\) and \(n\) are relatively prime, there is a unique solution for \(x\) modulo the product \(mn\); that is, the system of two equations is equivalent to a single equation of the form

\(x \equiv c \pmod {mn}.\)

We first compute the Bézout coefficients \(u\) and \(v\) such that \(mu + nv = 1\) using egcd, and then compute the solution as \(c = anv + bmu\). Indeed,

\(c = anv + bmu = a(1 - mu) + bmu = a - amu + bmu = a + (b-a)mu\)

and hence \(c \equiv a \pmod m\); similarly \(c \equiv b \pmod n\).

However, this is not quite general enough: we want to still be able to say something useful even if \(\gcd(m,n) &gt; 1\). I won’t go through the whole proof, but it turns out that there is a solution if and only if \(a \equiv b \pmod {\gcd(m,n)}\), and we can just divide everything through by \(g = \gcd(m,n)\), as we did for solving linear congruences. Here’s the code:

-- gcrt2 (a,n) (b,m) solves the pair of modular equations
--
--   x = a (mod n)
--   x = b (mod m)
--
-- It returns a pair (c, k) such that all solutions for x satisfy x =
-- c (mod k), that is, solutions are of the form x = kt + c for
-- integer t.
gcrt2 :: (Integer, Integer) -> (Integer, Integer) -> Maybe (Integer, Integer)
gcrt2 (a,n) (b,m)
  | a `mod` g == b `mod` g = Just (((a*v*m + b*u*n) `div` g) `mod` k, k)
  | otherwise              = Nothing
  where
    (g,u,v) = egcd n m
    k = (m*n) `div` g

From here we can bootstrap ourselves into solving systems of more than two equations, by iteratively combining two equations into one.

-- gcrt solves a system of modular equations.  Each equation x = a
-- (mod n) is given as a pair (a,n).  Returns a pair (z, k) such that
-- solutions for x satisfy x = z (mod k), that is, solutions are of
-- the form x = kt + z for integer t.
gcrt :: [(Integer, Integer)] -> Maybe (Integer, Integer)
gcrt []         = Nothing
gcrt [e]        = Just e
gcrt (e1:e2:es) = gcrt2 e1 e2 >>= \e -> gcrt (e:es)

Practice problems

And here are a bunch of problems for you to practice!