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.

About these ads

Pages: 1 2

8 Responses to “Modular Arithmetic”

  1. [...] Praxis – Modular Arithmetic By Remco Niemeijer Yesterday’s Programming Praxis problem is about making a convenient way to do modular arithmetic. While the [...]

  2. Remco Niemeijer said

    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]
    
  3. [...] 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 [...]

  4. programmingpraxis said

    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)))))))

  5. Lii said

    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

  6. Matthew said

    I am wondering, how do you deal with integer overflow when doing modular arithmetic? Especially when the size of modulus is close to the size of a c integer.

    For example, for addition in c I am trying the following with unsigned integers

    unsigned int res, a, b, p;

    res = a+b;
    if(res < a)
    while (res < (res = res + (0-p) % p));

    but I still need to figure out if it is correct.
    I just want a result that's congruent to the correct result – it doesn't have to be the smallest remainder mod p.

  7. Matthew said

    To answer my own question, here is what I came up with:
    uint res = a+b;
    if(res < a){
    while((p2 = (p < p) p = p2;
    if (p > res) res = res – p; else res = (res – p) – p;
    }

    But a better strategy may be to make sure the modulus p is always one bit less than the range of the integer, use signed ints and then reduce mod p before each addition or subtraction and not have any extra code.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 609 other followers

%d bloggers like this: