Chinese Remainders

August 27, 2010

The solution is easily computed by the Chinese remainder theorem, which states that given a set of positive, pairwise coprime moduli mi, 1 < i < r with product p = m1 · … · mr and corresponding residues ai mod mi, then the r relations nai mod mi have a unique solution computed as the sum over r of ai · vi · pi modulo p, where pi = p / mi and the vi are the inverses vi · pi ≡ 1 (mod mi).

Here is our version of the calculation:

(define (chinese-remainder as ms)
  (let ((p (apply * ms)))
    (define (f a m)
      (let* ((v (/ p m)) (b (inverse v m)))
        (* a b v)))
    (modulo (apply + (map f as ms)) p)))

The Chinese general had a thousand troups:

> (chinese-remainder '(10 4 2) '(11 12 13))
1000

We used the modular inversion function from the exercise on modular arithmetic. You can run the program at http://programmingpraxis.codepad.org/5LUl9JTa.

In testing this program, there were some calculations that failed with a divide-by-zero error. The problem arose in a modular inverse function that was stolen from a leading textbook. The book clearly states that the modulus need not be prime, as indeed it will not be in some uses of the Chinese remainder theorem, but in fact that is incorrect, and in that version of the modular inverse function the modulus must indeed by prime. The textbook authors were aware of the error, and noted the error in the errata on their website. Unfortunately, that algorithm was used in several exercises on Programming Praxis, including all recent uses of the prime? function, all of which will have to be modified to use the correct modular inverse function. This is, of course, a strong reminder that borrowing code or algorithms, even from a trusted source, makes the code or algorithm your own, and you, not the original author, bear responsibility for any errors that appear.

About these ads

Pages: 1 2

14 Responses to “Chinese Remainders”

  1. kbob said

    I think there is an infinite number of solutions of the form 1000 + k * (11 * 12 * 13), for k in 0, 1, …

  2. Graham said

    Cheating, with Sage:

    sage: crt(10, 4, 11, 12)
    76

    Sage’s crt(a, b, m, n) solves the problem x = a mod m and x = b mod n

  3. programmingpraxis said

    Kbob: Of course. But this is a blog about programming, not mathematics, so we seldom worry about such things.

    Graham: I think something is a little bit wrong there.

  4. Bryce said

    I will admit that I do not understand the mathematics to how you get a single answer to this. But hopefully that is something I can figure out over the weekend.

    In the meantime I can up with this function in Haskell:

    module Main where

    import Data.List

    main :: IO ()
    main = print $ x `intersect` (y `intersect` z)
    where x = [x | x <- [1..50000], mod x 11 == 10]
    y = [y | y <- [1..50000], mod y 12 == 4]
    z = [z | z <- [1..50000], mod z 13 == 12]

    Of which 1000 is the first answer, followed by 28 more.(50000 was just an arbitrarily picked upper bound.)

  5. Ramchip said

    The obvious solution:
    head $ [ a | a <- [1..], a `mod` 11 == 10, a `mod` 12 == 4, a `mod` 13 == 12 ]

  6. Graham said

    Ooops! Thanks for pointing out my mistake. Serves me right for not reading clearly.
    Again in Sage:
    sage: crt([10,4,12], [11,12,13])
    1000

    In Python, an answer similar to Ramchip’s (except we need to provide an upper bound for the list, because Python doesn’t have Haskell’s nice lazy evaluation; note that 1716 = 11*12*13):

    [x for x in xrange(1716) if x % 11 == 10 and x % 12 == 4 and x % 13 == 12][0]
    

    And for good measure, the same brute force logic in Fortran:

    program ppraxis_crt
    implicit none
    integer :: i
    do i = 1, 1716
        if (mod(i,13) == 12 .and. mod(i,12) == 4 .and. mod(i,11) == 10) then
            write(*,*) i
            stop
        end if
    end do
    stop
    end program ppraxis_crt
    

    Mathematically speaking, there are infinitely many solutions. Technically, the Chinese Remainder Theorem finds us the congruence class modulo the product of the (three in this case) moduli which satisfies all three congruences, given certain restrictions. For those interested in the inner workings, I’ve found that Wikipedia is surprisingly strong when it comes to mathematics of a reasonably high level.

  7. kbob said

    Here’s a Python brute force solution with lazy evaluation.

    count() generates 0, 1, …

    (for n in …) generates values of n that meet the remainder criteria.

    The .next() method runs the generator until it returns its next (first) value.

    from itertools import count
    (n for n in count() if all(n % m == r for m, r in ((11,10), (12,4), (13,12)))).next()
    
  8. kbob said

    I still say the problem is not quite specified. Nothing in the problem statement lets us determine whether the general had 1000 soldiers, 2716 soldiers, 4432 soldiers…

  9. Graham said

    kbob, nice job with the lazy evaluation in Python! As for your concerns, you’re right; the CRT gives us a congruence class for a solution which consists of infinitely many integers. Thus, any integer equivalent to 1000 mod 1716 could be an answer (including negative numbers, but that wouldn’t make sense in the context of the question).

  10. Sam W said

    I got 408 using brute force. 408%11 = 10, 408%12=4, 408%13=0

    int i;
    for(i = 21; ; ++ i)
        if(!((i%11) -10) && !((i%12)-4) && !(i%13))
            break;
    fprintf(stdout, "The emperor had %d troops.\n", i);
    
  11. programmingpraxis said

    Sam W: The remainder for the 13-rank case is 12, not 0.

  12. wm said

    This is quite old post, but I think it’s worth to note. Since result is always in form 13*k + 12, we can init k=21 and then increment by 13. First solution is found after 76 iterations, not 1000. :)

  13. David said

    Factor command line:

    car of infinite lazy list.

    
    ( scratchpad ) 1 lfrom [ [ 11 mod 10 = ] [ 12 mod 4 = ] [ 13 mod 12 = ] tri and and ] lfilter
    
    --- Data stack:
    T{ memoized-cons f ~lazy-filter~ ~array~ ~array~ ~array~ }
    ( scratchpad ) car .
    1000
    
  14. David said

    Redid it in FORTH using the Chinese Remainder Theorem. Includes the infrastructure for calculating the modular inverse. Didn’t know about the CRT before this…

    : egcd ( a b -- a b )
        dup 0= IF
            2drop 1 0
        ELSE
            dup -rot /mod               \ -- r=a%b q=a/b b
            -rot recurse                \ -- q (s,t) = egcd(b, r)
            >r swap r@ * - r> swap      \ -- t (s - q*t)
        THEN ;
    
    : egcd>gcd  ( a b x y -- n )  \ calculate gcd from egcd
        rot * -rot * + ;
    
    : mod-inv  ( a m -- a' )     \ modular inverse with coprime check
        2dup egcd over >r egcd>gcd r> swap 1 <> -24 and throw ;
    
    : array-product ( adr count -- n )
        1 -rot  cells bounds DO  i @ *  cell +LOOP ;
    
    : crt-from-array  ( adr1 adr2 count -- n )
        2dup array-product   locals| M count m[] a[] |
        0  \ result
        count 0 DO
            m[] i cells + @
            dup M swap /
            dup rot mod-inv *
            a[] i cells + @ * +
        LOOP  M mod ;
    
    create crt-residues[]  10 cells allot
    create crt-moduli[]    10 cells allot
    
    : crt ( .... n -- n )  \ takes pairs of "n (mod m)" from stack.
        10 min  locals| n |
        n 0 DO
            i cells crt-moduli[] + !
            i cells crt-residues[] + !
        LOOP
        crt-residues[] crt-moduli[] n crt-from-array ;
    
    : troops ( -- )
        ." The emperor had " 10 11  4 12  12 13  3 crt . ." troops" cr ;
    
    Gforth 0.7.2, Copyright (C) 1995-2008 Free Software Foundation, Inc.
    Gforth comes with ABSOLUTELY NO WARRANTY; for details type `license'
    Type `bye' to exit
    include crt.fs  ok
    troops The emperor had 1000 troops
     ok
    

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 574 other followers

%d bloggers like this: