ccan/junkcode/tterribe@email.unc.edu-nmbrthry/crt.c

#include "crt.h"
#include "egcd.h"

/*Computes the solution to a system of simple linear congruences via the
   Chinese Remainder Theorem.
  This function solves the system of equations
   x = a_i (mod m_i)
  A value x satisfies this equation if there exists an integer k such that
   x = a_i + k*m_i
  Note that under this definition, negative moduli are equivalent to positive
   moduli, and a modulus of 0 demands exact equality.
  x: Returns the solution, if it exists.
     Otherwise, the value is unchanged.
  a: An array of the a_i's.
  m: An array of the m_i's.
     These do not have to be relatively prime.
  n: The number of equations in the system.
  Return: -1 if the system is not consistent, otherwise the modulus by which
           the solution is unique.
          This modulus is the LCM of the m_i's, except in the case where one of
           them is 0, in which case the return value is also 0.*/
int crt(int *_x,const int _a[],const int _m[],int _n){
  int i;
  int m0;
  int a0;
  int x0;
  /*If there are no equations, everything is a solution.*/
  if(_n<1){
    *_x=0;
    return 1;
  }
  /*Start with the first equation.*/
  /*Force all moduli to be positive.*/
  m0=_m[0]<0?-_m[0]:_m[0];
  a0=_a[0];
  /*Add in each additional equation, and use it to derive a new, single
     equation.*/
  for(i=1;i<_n;i++){
    int d;
    int mi;
    int xi;
    /*Force all moduli to be positive.*/
    mi=_m[i]<0?-_m[i]:_m[i];
    /*Compute the inverse of m0 (mod mi) and of mi (mod m0).
      These are only inverses if m0 and mi are relatively prime.*/
    d=egcd(m0,mi,&xi,&x0);
    if(d>1){
      /*The hard case: m0 and mi are not relatively prime.*/
      /*First: check for consistency.*/
      if((a0-_a[i])%d!=0)return -1;
      /*If m0 divides mi, the old equation was completely redundant.*/
      else if(d==m0){
        a0=_a[i];
        m0=mi;
      }
      /*If mi divides m0, the new equation is completely redundant.*/
      else if(d!=mi){
        /*Otherwise the two have a non-trivial combination.
          The system is equivalent to the system
            x == a0 (mod m0/d)
            x == a0 (mod d)
            x == ai (mod mi/d)
            x == ai (mod d)
          Note that a0+c*(ai-a0) == a0 (mod d) == ai (mod d) for any c, since
           d|ai-a0; thus any such c gives solutions that satisfy eqs. 2 and 4.
          Choosing c as a multiple of (m0/d) ensures
            a0+c*(ai-a0) == a0 (mod m0/d), satisfying eq. 1.
          But (m0/d) and (mi/d) are relatively prime, so we can choose
            c = (m0/d)*((m0/d)^-1 mod mi/d),
          Hence c == 1 (mod mi/d), and
            a0+c*(ai-a0) == ai (mod mi/d), satisfying eq. 3.
          The inverse of (m0/d) mod (mi/d) can be computed with the egcd().*/
        m0/=d;
        egcd(m0,mi/d,&xi,&x0);
        a0+=(_a[i]-a0)*m0*xi;
        m0*=mi;
        a0=a0%m0;
      }
    }
    else if(d==1){
      /*m0 and mi are relatively prime, so xi and x0 are the inverses of m0 and
         mi modulo mi and m0, respectively.
        The Chinese Remainder Theorem now states that the solution is given by*/
      a0=a0*mi*x0+_a[i]*m0*xi;
      /* modulo the LCM of m0 and mi.*/
      m0*=mi;
      a0%=m0;
    }
    /*Special case: mi and m0 are both 0.
      Check for consistency.*/
    else if(a0!=_a[i])return -1;
    /*Otherwise, this equation was redundant.*/
  }
  /*If the final modulus was not 0, then constrain the answer to be
     non-negative and less than that modulus.*/
  if(m0!=0){
    x0=a0%m0;
    *_x=x0<0?x0+m0:x0;
  }
  /*Otherwise, there is only one solution.*/
  else *_x=a0;
  return m0;
}