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
.
Let’s write a function to solve such equivalences for . We want a pair of integers
and
such that
is a solution to
if and only if
. This isn’t hard to write in the end, but takes a little bit of thought to do it properly.
First of all, if and
are relatively prime (that is,
) then we know from the last post that
has an inverse modulo
; multiplying both sides by
yields the solution
.
OK, but what if ? In this case there might not even be any solutions. For example,
has no solutions: any even number will be equivalent to
or
modulo
, so there is no value of
such that double it will be equivalent to
. On the other hand,
is OK: this will be true for any odd value of
, that is,
. In fact, it is easy to see that any common divisor of
and
must also divide
in order to have any solutions. In case the GCD of
and
does divide
, we can simply divide through by the GCD (including dividing the modulus
!) 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
then as long as and
are relatively prime, there is a unique solution for
modulo the product
; that is, the system of two equations is equivalent to a single equation of the form
We first compute the Bézout coefficients and
such that
using
egcd
, and then compute the solution as . Indeed,
and hence ; similarly
.
However, this is not quite general enough: we want to still be able to say something useful even if . I won’t go through the whole proof, but it turns out that there is a solution if and only if
, and we can just divide everything through by
, 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!
Pingback: Resumen de lecturas compartidas del 1 al 7 de marzo de 2020 | Vestigium
Pingback: Competitive programming in Haskell: summer series | blog :: Brent -> [String]
Pingback: Competitive programming in Haskell: permutations | blog :: Brent -> [String]
Pingback: Resumen de lecturas compartidas durante marzo de 2020 | Vestigium