Modular Arithmetic
July 7, 2009
We begin with the oldest algorithm still in use today, Euclid’s algorithm for calculating the greatest common divisor, which pre-dates even the grade-school algorithms for basic arithmetic; Knuth calls this the “grand-daddy” of all algorithms. The normal formulation of Euclid’s algorithm uses division via the modulo operator, but Euclid used repeated subtraction:
(define (grand-daddy m n)
(cond ((< m n) (grand-daddy m (- n m)))
((< n m) (grand-daddy (- m n) n))
(else m)))
Étienne Bézout was an eighteenth-century French mathematician known primarily for his contributions to number theory. Bézout's identity is a linear Diophantine equation that states if x and y are non-zero integers with greatest common divisor g, there exist integers a and b such that ax + by = g. An extended version of Euclid’s algorithm can calculate the Bézout coefficients a and b, and at the same time calculate the greatest common divisor g:
(define (euclid x y)
(let loop ((a 1) (b 0) (g x) (u 0) (v 1) (w y))
(if (zero? w) (values a b g)
(let ((q (quotient g w)))
(loop u v w (- a (* q u)) (- b (* q v)) (- g (* q w)))))))
The calculation of the modular inverse uses the extended Euclid’s algorithm:
(define (inverse x m)
(if (not (= (gcd x m) 1))
(error 'inverse "divisor must be coprime to modulus")
(call-with-values
(lambda () (euclid x m))
(lambda (a b g) (modulo a m)))))
We previously used modular exponentiation in the Primality Checking and Monte Carlo Factorization exercises, and it is already part of the Standard Prelude; we repeat the definition here:
(define (expm b e m)
(define (m* x y) (modulo (* x y) m))
(cond ((zero? e) 1)
((even? e) (expm (m* b b) (/ e 2) m))
(else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))
The jacobi symbol is a simple calculation (simpler than calculating the modular square root) that is 0 if a is not coprime to n, +1 if a is a quadratic residue of n, and -1 if a is not a quadratic residue of n:
(define (jacobi a n)
(if (not (and (integer? a) (integer? n) (positive? n) (odd? n)))
(error 'jacobi "modulus must be positive odd integer")
(let jacobi ((a a) (n n))
(cond ((= a 0) 0)
((= a 1) 1)
((= a 2) (case (modulo n 8) ((1 7) 1) ((3 5) -1)))
((even? a) (* (jacobi 2 n) (jacobi (quotient a 2) n)))
((< n a) (jacobi (modulo a n) n))
((and (= (modulo a 4) 3) (= (modulo n 4) 3)) (- (jacobi n a)))
(else (jacobi n a))))))
The modular square root is defined only if the modulus is an odd prime and the target number has a quadratic residue; where it does exist, its two conjugates can be computed as follows:
(define (mod-sqrt a p)
(define (both n) (list n (- p n)))
(cond ((not (and (odd? p) (prime? p)))
(error 'mod-sqrt "modulus must be an odd prime"))
((not (= (jacobi a p) 1))
(error 'mod-sqrt "must be a quadratic residual"))
(else (let ((a (modulo a p)))
(case (modulo p 8)
((3 7) (both (expm a (/ (+ p 1) 4) p)))
((5) (let* ((x (expm a (/ (+ p 3) 8) p))
(c (expm x 2 p)))
(if (= a c) (both x)
(both (modulo (* x (expm 2 (/ (- p 1) 4) p)) p)))))
(else (let* ((d (let loop ((d 2))
(if (= (jacobi d p) -1) d
(loop (+ d 1)))))
(s (let loop ((p (- p 1)) (s 0))
(if (odd? p) s
(loop (quotient p 2) (+ s 1)))))
(t (quotient (- p 1) (expt 2 s)))
(big-a (expm a t p))
(big-d (expm d t p))
(m (let loop ((i 0) (m 0))
(cond ((= i s) m)
((= (- p 1)
(expm (* big-a (expm big-d m p))
(expt 2 (- s 1 i)) p))
(loop (+ i 1) (+ m (expt 2 i))))
(else (loop (+ i 1) m))))))
(both (modulo (* (expm a (/ (+ t 1) 2) p)
(expm big-d (/ m 2) p)) p)))))))))
It is convenient to have an environment in which all calculations are performed with respect to a given modulus. The (with-modulus modulus expr ...) macro puts together all the various modular operators in a single environment, providing == for determining congruence, +, -, * and / with their usual meanings, ^ for modular exponentiation, and sqrt for the modular square root. Just like the regular Scheme operators, +, -, * and / can take any number of operands, and special meanings are give to unary – (the conjugate operator) and unary / (the inverse):
(define-syntax (with-modulus stx)
(syntax-case stx ()
((with-modulus e expr ...)
(with-syntax ((modulus (datum->syntax-object (syntax with-modulus) 'modulus))
(== (datum->syntax-object (syntax with-modulus) '== ))
(+ (datum->syntax-object (syntax with-modulus) '+ ))
(- (datum->syntax-object (syntax with-modulus) '- ))
(* (datum->syntax-object (syntax with-modulus) '* ))
(/ (datum->syntax-object (syntax with-modulus) '/ ))
(^ (datum->syntax-object (syntax with-modulus) '^ ))
(sqrt (datum->syntax-object (syntax with-modulus) 'sqrt )))
(syntax (letrec ((fold (lambda (op base xs)
(if (null? xs) base
(fold op (op base (car xs)) (cdr xs))))))
(let* ((modulus e)
(mod (lambda (x)
(if (not (integer? x))
(error 'with-modulus "all arguments must be integers")
(modulo x modulus))))
(== (lambda (x y) (= (mod x) (mod y))))
(+ (lambda xs (fold (lambda (x y) (mod (+ x (mod y)))) 0 xs)))
(- (lambda (x . xs)
(if (null? xs)
(mod (- 0 x))
(fold (lambda (x y) (mod (- x (mod y)))) x xs))))
(* (lambda xs (fold (lambda (x y) (mod (* x (mod y)))) 1 xs)))
(/ (lambda (x . xs)
(if (null? xs)
(inverse x e)
(fold (lambda (x y) (* x (inverse y e))) x xs))))
(^ (lambda (base exp) (expm base exp e)))
(sqrt (lambda (x) (mod-sqrt x e))))
expr ...)))))))
Here are some examples:
(with-modulus 12
(display (== 17 5)) (newline) ; #t
(display (+ 8 9)) (newline) ; 5
(display (- 4 9)) (newline) ; 7
(display (* 3 7)) (newline) ; 9
(display (/ 9 7)) (newline)) ; 3
(with-modulus 13
(display (^ 6 2)) (newline) ; 10
(display (^ 7 2)) (newline) ; 10
(display (sqrt 10)) (newline)) ; ±6
The specific algorithms given above derive from the book Prime Numbers: A Computational Perspective by Richard Crandall and Carl Pomerance, but there are many other sources. The square-root function uses prime? and its companion check? from the Primality Checking exercise, and rand and randint from the Standard Prelude. You can run this program at http://programmingpraxis.codepad.org/1QA39Fbf.
Pages: 1 2
[...] Praxis – Modular Arithmetic By Remco Niemeijer Yesterday’s Programming Praxis problem is about making a convenient way to do modular arithmetic. While the [...]
My Haskell solution (see http://bonsaicode.wordpress.com/2009/07/08/programming-praxis-modular-arithmetic/ for a version with comments):
import Data.Bits import Data.List euclid :: Integral a => a -> a -> a euclid x y = f 1 0 x 0 1 y where f a b g u v w | w == 0 = mod a y | otherwise = f u v w (a-q*u) (b-q*v) (g-q*w) where q = div g w inv :: Integral a => a -> a -> a inv x m | gcd x m == 1 = euclid x m | otherwise = error "divisor must be coprime to modulus" (<=>) :: Integral a => a -> a -> a -> Bool a <=> b = \m -> mod a m == mod b m (<+>), (<->), (<*>), (</>) :: Integral a => a -> a -> a -> a a <+> b = \m -> mod (a + mod b m) m a <-> b = a <+> (-b) a <*> b = \m -> mod (a * mod b m) m a </> b = \m -> (a <*> inv b m) m expm :: Integer -> Integer -> Integer -> Integer expm b e m = foldl' (\r (b', _) -> mod (r * b') m) 1 . filter (flip testBit 0 . snd) . zip (iterate (flip mod m . (^ 2)) b) $ takeWhile (> 0) $ iterate (`shiftR` 1) e (<^>) :: Integer -> Integer -> Integer -> Integer (<^>) = expm sqrtm :: Integral a => a -> a -> [a] sqrtm x m = case [n | n <- [1..m], (n <*> n) m == mod x m] of (_:_:_:_) -> [] s -> s modulo :: a -> (a -> b) -> b modulo = flip ($) main :: IO () main = do mapM_ (modulo 12) [ print . (17 <=> 5), print . (8 <+> 9), print . (4 <-> 9), print . (3 <*> 7), print . (9 </> 7)] mapM_ (modulo 13) [ print . (6 <^> 2), print . (7 <^> 2), print . sqrtm 10][...] OLE 2 or ActiveX, we are referring to a very specific type of object, sometimes called a … Modular Arithmetic « Programming Praxis – programmingpraxis.com 12/31/1969 Programming Praxis. A collection of etudes, updated weekly, for [...]
Here is an improved version of modular inverse, based on Algorithm 9.4.4 from the book Prime Numbers: A Computational Perspective by Richard Crandall and Carl Pomerance:
(define (inverse x m)(if (not (= (gcd x m) 1))
(error 'inverse "divisor must be coprime to modulus")
(let loop ((z (modulo x m)) (a 1))
(if (= z 1) a
(let ((q (- (quotient m z))))
(loop (+ m (* q z)) (modulo (* q a) m)))))))
Here is a variant of the above code where the operation take and return functions. The intention is to get rid of the ugliness of having to pass the modulo value in every function invocation.
The disadvantage is that you cant use ints directly in the expressions, you have to wrap them in the ugly i functions. Is there a solution to this?
test_func = ((i 5) *% (i 6) +% (i 7)) 13
(=%) :: TModNum -> TModNum -> (Integer -> Bool)
a =% b = \m -> (a m) == (b m)
type TModNum = (Integer -> Integer)
(+%), (-%), (*%), (/%) :: TModNum -> TModNum -> TModNum
a +% b = \m -> mod ((a m) + (b m)) m
a -% b = \m -> mod ((a m) - (b m)) m
a *% b = \m -> mod ((a m) * (b m)) m
a /% b = \m -> mod ((a m) * (inv (b m) m)) m
i :: Integer -> (Integer -> Integer)
i n = \m -> mod n m