chrem.cal

早期freebsd实现

/*
 * 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 <[e-mail]>
 * Interface by: Landon Curt Noll <[e-mail]>
 */

static 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 (lib_debug >= 0) {
    print "chrem(r1,m1 [,r2,m2 ...]) defined";
    print "chrem(rlist [,mlist]) defined";
}