/* * chrem - Chinese remainder theorem/problem solver * * When possible, chrem finds solutions for x of a set of congruences * of the form: * * x = r1 (mod m1) * x = r2 (mod m2) * ... * * where the residues r1, r2, ... and the moduli m1, m2, ... are * given integers. The Chinese remainder theorem states that if * m1, m2, ... are relatively prime in pairs, the above congruences * have a unique solution modulo m1 * m2 * ... If m1, m2, ... * are not relatively prime in pairs, it is possible that no solution * exists. If solutions exist, the general solution is expressible as: * * x = r (mod m) * * where m = lcm(m1,m2,...), and if m > 0, 0 <= r < m. This * solution may be interpreted as: * * x = r + k * m [[NOTE 1]] * * where k is an arbitrary integer. * *** * * usage: * * chrem(r1,m1 [,r2,m2, ...]) * * r1, r2, ... remainder integers or null values * m1, m2, ... moduli integers * * chrem(r_list, [m_list]) * * r_list list (r1,r2, ...) * m_list list (m1,m2, ...) * * If m_list is omitted, then 'defaultmlist' is used. * This default list is a global value that may be changed * by the user. Initially it is the first 8 primes. * * If a remainder is null(), then the corresponding congruence is * ignored. This is useful when working with a fixed list of moduli. * * If there are more remainders than moduli, then the later moduli are * ignored. * * The moduli may be any integers, not necessarily relatively prime in * pairs (as required for the Chinese remainder theorem). Any moduli * may be zero; x = r (mod 0) has the meaning of x = r. * * returns: * * If args were integer pairs: * * r ('r' is defined above, see [[NOTE 1]]) * * If 1 or 2 list args were given: * * (r, m) ('r' and 'm' are defined above, see [[NOTE 1]]) * * NOTE: In all cases, null() is returned if there is no solution. * *** * * This function may be used to solve the following historical problems: * * Sun-Tsu, 1st century A.D. * * To find a number for which the reminders after division by 3, 5, 7 * are 2, 3, 2, respectively: * * chrem(2,3,3,5,2,7) ---> 23 * * Fibonacci, 13th century A.D. * * To find a number divisible by 7 which leaves remainder 1 when * divided by 2, 3, 4, 5, or 6: * * * chrem(list(0,1,1,1,1,1),list(7,2,3,4,5,6)) ---> (301,420) * * i.e., any value that is 301 mod 420. * * Written by: Ernest W Bowen <ernie@neumann.une.edu.au> * Interface by: Landon Curt Noll <chongo@toad.com> */ defaultmlist = list(2,3,5,7,11,13,17,19); /* The first eight primes */ define chrem() { local argc; /* number of args given */ local rlist; /* reminder list - ri */ local mlist; /* modulus list - mi */ local list_args; /* true => args given are lists, not r1,m1, ... */ local m,z,r,y,d,t,x,u,i; /* * parse args */ argc = param(0); if (argc == 0) { quit "usage: chrem(r1,m1 [,r2,m2 ...]) or chrem(r_list, m_list)"; } list_args = islist(param(1)); if (list_args) { rlist = param(1); mlist = (argc == 1) ? defaultmlist : param(2); if (size(rlist) > size(mlist)) { quit "too many residues"; } } else { if (argc % 2 == 1) { quit "odd number integers given"; } rlist = list(); mlist = list(); for (i=1; i <= argc; i+=2) { push(rlist, param(i)); push(mlist, param(i+1)); } } /* * solve the problem found in rlist & mlist */ m = 1; z = 0; while (size(rlist)) { r=pop(rlist); y=abs(pop(mlist)); if (r==null()) continue; if (m) { if (y) { d = t = z - r; m = lcm(x=m, y); while (d % y) { u = x; x %= y; swap(x,y); if (y==0) return; z += (t *= -u/y); } } else { if ((r % m) != (z % m)) return; else { m = 0; z = r; } } } else if (((y) && (r % y != z % y)) || (r != z)) return; } if (m) { z %= m; if (z < 0) z += m; } /* * return information as required */ if (list_args) { return list(z,m); } else { return z; } } global lib_debug; if (!isnum(lib_debug)||lib_debug>0) print "chrem(r1,m1 [,r2,m2 ...]) defined"; if (!isnum(lib_debug)||lib_debug>0) print "chrem(rlist [,mlist]) defined";