ChineseRemainderTheorem.frink

Download or view ChineseRemainderTheorem.frink in plain text format


/** This contains a function to solve a system of modular integer congruences
    using the Chinese Remainder Theorem.

    See:
    https://brilliant.org/wiki/chinese-remainder-theorem/

    Given pairwise coprime positive integer moduli m1, m2, m3, ... mk
    and arbitrary integer remainders r1, r2, r3, ... rk,
    this solves the system of simultaneous congruences

    x ≡ r1 (mod m1)
    x ≡ r2 (mod m2)
    x ≡ r3 (mod m3)
    ...
    x ≡ rk (mod mk)

    And returns the unique solution (x mod N) where N is the product of all the
    m terms.

    arguments:
       [r, m]  where r and m are arrays of the remainder terms r and the
               modulus terms m respectively.  These must be of the same length.
    returns
       x, the unique solution mod N where N is the product of all the M terms.

    Example:
       ChineseRemainder[[1,4,6], [3,5,7]] == 34

     The "comet" puzzle on the Brilliant site above can be solved by:
        ChineseRemainder[[2017,2014,2008],[3,8,13],2017] == 2086

    Note:
       This does not verify that the moduli are co-prime.  But maybe that can
       be done? See
        https://en.wikipedia.org/wiki/Chinese_remainder_theorem#Generalization_to_non-coprime_moduli

       This does not currently seem to work for negative numbers.  Fix that.
*/

ChineseRemainder[r, m] :=
{
   return ChineseRemainder[r, m, 0]
}

/** Returns a solution using the Chinese Remainder Theorem as above, but making
    the solution x be the smallest solution >= d. */

ChineseRemainder[r, m, d] :=
{
   if length[r] != length[m]
   {
      println["ChineseRemainder:  r and m must be arrays of the same length."]
      return undef
   }
   
   N = product[m]

   y = new array
   z = new array
   x = 0
   for i = rangeOf[m]
   {
      y@i = N / m@i
      z@i = modInverse[y@i, m@i]
      if z@i == undef
      {
         println["ChineseRemainder:  modInverse returned undef for modInverse[" + y@i + ", " + m@i + "]"]
         return undef
      }
      
      x = x + r@i y@i z@i
   }

   xp = x mod N
   f = d div N
   r = f * N + xp
   if r < d
      r = r + N

   return r
}


Download or view ChineseRemainderTheorem.frink in plain text format


This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 19355 days, 23 hours, 43 minutes ago.