--- /dev/null
+#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;
+}
--- /dev/null
+#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
--- /dev/null
+#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;
+}
--- /dev/null
+#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
--- /dev/null
+#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;
+}
--- /dev/null
+#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
--- /dev/null
+#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;
+}
--- /dev/null
+#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
projects / ccan / commitdiff