git.ozlabs.org Git - ccan/blob - junkcode/tterribe@email.unc.edu-nmbrthry/crt.c

   1 #include "crt.h"
   2 #include "egcd.h"
   3
   4 /*Computes the solution to a system of simple linear congruences via the
   5    Chinese Remainder Theorem.
   6   This function solves the system of equations
   7    x = a_i (mod m_i)
   8   A value x satisfies this equation if there exists an integer k such that
   9    x = a_i + k*m_i
  10   Note that under this definition, negative moduli are equivalent to positive
  11    moduli, and a modulus of 0 demands exact equality.
  12   x: Returns the solution, if it exists.
  13      Otherwise, the value is unchanged.
  14   a: An array of the a_i's.
  15   m: An array of the m_i's.
  16      These do not have to be relatively prime.
  17   n: The number of equations in the system.
  18   Return: -1 if the system is not consistent, otherwise the modulus by which
  19            the solution is unique.
  20           This modulus is the LCM of the m_i's, except in the case where one of
  21            them is 0, in which case the return value is also 0.*/
  22 int crt(int *_x,const int _a[],const int _m[],int _n){
  23   int i;
  24   int m0;
  25   int a0;
  26   int x0;
  27   /*If there are no equations, everything is a solution.*/
  28   if(_n<1){
  29     *_x=0;
  30     return 1;
  31   }
  32   /*Start with the first equation.*/
  33   /*Force all moduli to be positive.*/
  34   m0=_m[0]<0?-_m[0]:_m[0];
  35   a0=_a[0];
  36   /*Add in each additional equation, and use it to derive a new, single
  37      equation.*/
  38   for(i=1;i<_n;i++){
  39     int d;
  40     int mi;
  41     int xi;
  42     /*Force all moduli to be positive.*/
  43     mi=_m[i]<0?-_m[i]:_m[i];
  44     /*Compute the inverse of m0 (mod mi) and of mi (mod m0).
  45       These are only inverses if m0 and mi are relatively prime.*/
  46     d=egcd(m0,mi,&xi,&x0);
  47     if(d>1){
  48       /*The hard case: m0 and mi are not relatively prime.*/
  49       /*First: check for consistency.*/
  50       if((a0-_a[i])%d!=0)return -1;
  51       /*If m0 divides mi, the old equation was completely redundant.*/
  52       else if(d==m0){
  53         a0=_a[i];
  54         m0=mi;
  55       }
  56       /*If mi divides m0, the new equation is completely redundant.*/
  57       else if(d!=mi){
  58         /*Otherwise the two have a non-trivial combination.
  59           The system is equivalent to the system
  60             x == a0 (mod m0/d)
  61             x == a0 (mod d)
  62             x == ai (mod mi/d)
  63             x == ai (mod d)
  64           Note that a0+c*(ai-a0) == a0 (mod d) == ai (mod d) for any c, since
  65            d|ai-a0; thus any such c gives solutions that satisfy eqs. 2 and 4.
  66           Choosing c as a multiple of (m0/d) ensures
  67             a0+c*(ai-a0) == a0 (mod m0/d), satisfying eq. 1.
  68           But (m0/d) and (mi/d) are relatively prime, so we can choose
  69             c = (m0/d)*((m0/d)^-1 mod mi/d),
  70           Hence c == 1 (mod mi/d), and
  71             a0+c*(ai-a0) == ai (mod mi/d), satisfying eq. 3.
  72           The inverse of (m0/d) mod (mi/d) can be computed with the egcd().*/
  73         m0/=d;
  74         egcd(m0,mi/d,&xi,&x0);
  75         a0+=(_a[i]-a0)*m0*xi;
  76         m0*=mi;
  77         a0=a0%m0;
  78       }
  79     }
  80     else if(d==1){
  81       /*m0 and mi are relatively prime, so xi and x0 are the inverses of m0 and
  82          mi modulo mi and m0, respectively.
  83         The Chinese Remainder Theorem now states that the solution is given by*/
  84       a0=a0*mi*x0+_a[i]*m0*xi;
  85       /* modulo the LCM of m0 and mi.*/
  86       m0*=mi;
  87       a0%=m0;
  88     }
  89     /*Special case: mi and m0 are both 0.
  90       Check for consistency.*/
  91     else if(a0!=_a[i])return -1;
  92     /*Otherwise, this equation was redundant.*/
  93   }
  94   /*If the final modulus was not 0, then constrain the answer to be
  95      non-negative and less than that modulus.*/
  96   if(m0!=0){
  97     x0=a0%m0;
  98     *_x=x0<0?x0+m0:x0;
  99   }
 100   /*Otherwise, there is only one solution.*/
 101   else *_x=a0;
 102   return m0;
 103 }