New junkcode from Tim.

junkcode/tterribe@email.unc.edu-nmbrthry/choose.c

+#include "choose.h"
+#include "gcd.h"
+
+/*Computes the number of combinations of _n items, taken _m at a time without
+   overflow.
+  _n: The total number of items.
+  _m: The number taken at a time.
+  Return: The number of combinations of _n items taken _m at a time.*/
+unsigned choose(int _n,int _m){
+  unsigned ret;
+  int      i;
+  ret=1;
+  for(i=1;i<=_m;_n--,i++){
+    int nmi;
+    nmi=_n%i;
+    if(nmi==0)ret*=_n/i;
+    else if(ret%i==0)ret=(ret/i)*_n;
+    else{
+      int d;
+      d=gcd(i,nmi);
+      ret=(ret/(i/d))*(_n/d);
+    }
+  }
+  return ret;
+}

junkcode/tterribe@email.unc.edu-nmbrthry/choose.h

+#if !defined(_choose_H)
+# define _choose_H (1)
+
+/*Computes the number of combinations of _n items, taken _m at a time without
+   overflow.
+  _n: The total number of items.
+  _b: The number taken at a time.
+  Return: The number of combinations of _n items taken _m at a time.*/
+unsigned choose(int _n,int _m);
+
+#endif

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;
+}

junkcode/tterribe@email.unc.edu-nmbrthry/crt.h

+#if !defined(_crt_H)
+# define _crt_H (1)
+
+/*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);
+
+#endif

junkcode/tterribe@email.unc.edu-nmbrthry/egcd.c

+#include "egcd.h"
+
+/*Computes the coefficients of the smallest positive linear combination of two
+   integers _a and _b.
+  These are a solution (u,v) to the equation _a*u+_b*v==gcd(_a,_b).
+  _a: The first integer of which to compute the extended gcd.
+  _b: The second integer of which to compute the extended gcd.
+  _u: Returns the coefficient of _a in the smallest positive linear
+       combination.
+  _v: Returns the coefficient _of b in the smallest positive linear
+       combination.
+  Return: The non-negative gcd of _a and _b.
+          If _a and _b are both 0, then 0 is returned, though in reality the
+           gcd is undefined, as any integer, no matter how large, will divide 0
+           evenly.
+          _a*u+_b*v will not be positive in this case.
+          Note that the solution (u,v) of _a*u+_b*v==gcd(_a,_b) returned is not
+           unique.
+          (u+(_b/gcd(_a,_b))*k,v-(_a/gcd(_a,_b))*k) is also a solution for all
+           k.
+          The coefficients (u,v) might not be positive.*/
+int egcd(int _a,int _b,int *_u,int *_v){
+  int a;
+  int b;
+  int s;
+  int t;
+  int u;
+  int v;
+  /*Make both arguments non-negative.
+    This forces the return value to be non-negative.*/
+  a=_a<0?-_a:_a;
+  b=_b<0?-_b:_b;
+  /*Simply use the extended Euclidean algorithm.*/
+  s=v=0;
+  t=u=1;
+  while(b){
+    int q;
+    int r;
+    int w;
+    q=a/b;
+    r=a%b;
+    a=b;
+    b=r;
+    w=s;
+    s=u-q*s;
+    u=w;
+    w=t;
+    t=v-q*t;
+    v=w;
+  }
+  /*u and v were computed for non-negative a and b.
+    If the arguments passed in were negative, flip the sign.*/
+  *_u=_a<0?-u:u;
+  *_v=_b<0?-v:v;
+  return a;
+}

junkcode/tterribe@email.unc.edu-nmbrthry/egcd.h

+#if !defined(_egcd_H)
+# define _egcd_H (1)
+
+/*Computes the coefficients of the smallest positive linear combination of two
+   integers _a and _b.
+  These are a solution (u,v) to the equation _a*u+_b*v==gcd(_a,_b).
+  _a: The first integer of which to compute the extended gcd.
+  _b: The second integer of which to compute the extended gcd.
+  _u: Returns the coefficient of _a in the smallest positive linear
+       combination.
+  _v: Returns the coefficient _of b in the smallest positive linear
+       combination.
+  Return: The non-negative gcd of _a and _b.
+          If _a and _b are both 0, then 0 is returned, though in reality the
+           gcd is undefined, as any integer, no matter how large, will divide 0
+           evenly.
+          _a*u+_b*v will not be positive in this case.
+          Note that the solution (u,v) of _a*u+_b*v==gcd(_a,_b) returned is not
+           unique.
+          (u+(_b/gcd(_a,_b))*k,v-(_a/gcd(_a,_b))*k) is also a solution for all
+           k.
+          The coefficients (u,v) might not be positive.*/
+int egcd(int _a,int _b,int *_u,int *_v);
+
+#endif

junkcode/tterribe@email.unc.edu-nmbrthry/gcd.c

+#include "gcd.h"
+
+/*Computes the gcd of two integers, _a and _b.
+  _a: The first integer of which to compute the gcd.
+  _b: The second integer of which to compute the gcd.
+  Return: The non-negative gcd of _a and _b.
+          If _a and _b are both 0, then 0 is returned, though in reality the
+           gcd is undefined, as any integer, no matter how large, will divide 0
+           evenly.*/
+int gcd(int _a,int _b){
+  /*Make both arguments non-negative.
+    This forces the return value to be non-negative.*/
+  if(_a<0)_a=-_a;
+  if(_b<0)_b=-_b;
+  /*Simply use the Euclidean algorithm.*/
+  while(_b){
+    int r;
+    r=_a%_b;
+    _a=_b;
+    _b=r;
+  }
+  return _a;
+}

junkcode/tterribe@email.unc.edu-nmbrthry/gcd.h

+#if !defined(_gcd_H)
+# define _gcd_H (1)
+
+/*Computes the gcd of two integers, _a and _b.
+  _a: The first integer of which to compute the gcd.
+  _b: The second integer of which to compute the gcd.
+  Return: The non-negative gcd of _a and _b.
+          If _a and _b are both 0, then 0 is returned, though in reality the
+           gcd is undefined, as any integer, no matter how large, will divide 0
+           evenly.*/
+int gcd(int _a,int _b);
+
+#endif