ChineseRemainderTheorem.frink

View or download 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
//   println["N is $N"]
   f = d div N
//   println["f is $f"]
   r = f * N + xp
   if r < d
      r = r + N

   return r
}


View or download 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 19100 days, 8 hours, 58 minutes ago.