1/*
2 * chrem - Chinese remainder theorem/problem solver
3 *
4 * Copyright (C) 1999,2021  Ernest Bowen and Landon Curt Noll
5 *
6 * Primary author:  Ernest Bowen
7 *
8 * Calc is open software; you can redistribute it and/or modify it under
9 * the terms of the version 2.1 of the GNU Lesser General Public License
10 * as published by the Free Software Foundation.
11 *
12 * Calc is distributed in the hope that it will be useful, but WITHOUT
13 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
15 * Public License for more details.
16 *
17 * A copy of version 2.1 of the GNU Lesser General Public License is
18 * distributed with calc under the filename COPYING-LGPL.  You should have
19 * received a copy with calc; if not, write to Free Software Foundation, Inc.
20 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
21 *
22 * Under source code control:	1992/09/26 01:00:47
23 * File existed as early as:	1992
24 *
25 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
26 */
27
28/*
29 * When possible, chrem finds solutions for x of a set of congruence
30 * of the form:
31 *
32 *			x = r1 (mod m1)
33 *			x = r2 (mod m2)
34 *			   ...
35 *
36 * where the residues r1, r2, ... and the moduli m1, m2, ... are
37 * given integers.  The Chinese remainder theorem states that if
38 * m1, m2, ... are relatively prime in pairs, the above congruence
39 * have a unique solution modulo  m1 * m2 * ...	  If m1, m2, ...
40 * are not relatively prime in pairs, it is possible that no solution
41 * exists.  If solutions exist, the general solution is expressible as:
42 *
43 *		     x = r (mod m)
44 *
45 * where m = lcm(m1,m2,...), and if m > 0, 0 <= r < m.	This
46 * solution may be interpreted as:
47 *
48 *		    x = r + k * m			[[NOTE 1]]
49 *
50 * where k is an arbitrary integer.
51 *
52 ***
53 *
54 * usage:
55 *
56 *	chrem(r1,m1 [,r2,m2, ...])
57 *
58 *	    r1, r2, ...		remainder integers or null values
59 *	    m1, m2, ...		moduli integers
60 *
61 *	chrem(r_list, [m_list])
62 *
63 *	    r_list	list (r1,r2, ...)
64 *	    m_list	list (m1,m2, ...)
65 *
66 *	    If m_list is omitted, then 'defaultmlist' is used.
67 *	    This default list is a global value that may be changed
68 *	    by the user.  Initially it is the first 8 primes.
69 *
70 * If a remainder is null(), then the corresponding congruence is
71 * ignored.  This is useful when working with a fixed list of moduli.
72 *
73 * If there are more remainders than moduli, then the later moduli are
74 * ignored.
75 *
76 * The moduli may be any integers, not necessarily relatively prime in
77 * pairs (as required for the Chinese remainder theorem).  Any moduli
78 * may be zero;	 x = r (mod 0) has the meaning of x = r.
79 *
80 * returns:
81 *
82 *    If args were integer pairs:
83 *
84 *	  r		('r' is defined above, see [[NOTE 1]])
85 *
86 *    If 1 or 2 list args were given:
87 *
88 *	  (r, m)	('r' and 'm' are defined above, see [[NOTE 1]])
89 *
90 * NOTE: In all cases, null() is returned if there is no solution.
91 *
92 ***
93 *
94 * This function may be used to solve the following historical problems:
95 *
96 *   Sun-Tsu, 1st century A.D.
97 *
98 *	To find a number for which the reminders after division by 3, 5, 7
99 *	are 2, 3, 2, respectively:
100 *
101 *	    chrem(2,3,3,5,2,7) ---> 23
102 *
103 *   Fibonacci, 13th century A.D.
104 *
105 *	To find a number divisible by 7 which leaves remainder 1 when
106 *	divided by 2, 3, 4, 5, or 6:
107 *
108 *
109 *	    chrem(list(0,1,1,1,1,1),list(7,2,3,4,5,6)) ---> (301,420)
110 *
111 *	i.e., any value that is 301 mod 420.
112 */
113
114
115static defaultmlist = list(2,3,5,7,11,13,17,19); /* The first eight primes */
116
117define chrem()
118{
119    local argc;		/* number of args given */
120    local rlist;	/* reminder list - ri */
121    local mlist;	/* modulus list - mi */
122    local list_args;	/* true => args given are lists, not r1,m1, ... */
123    local m,z,r,y,d,t,x,u,i;
124
125    /*
126     * parse args
127     */
128    argc = param(0);
129    if (argc == 0) {
130	quit "usage: chrem(r1,m1 [,r2,m2 ...]) or chrem(r_list, m_list)";
131    }
132    list_args = islist(param(1));
133    if (list_args) {
134	rlist = param(1);
135	mlist = (argc == 1) ? defaultmlist : param(2);
136	if (size(rlist) > size(mlist)) {
137	    quit "too many residues";
138	}
139    } else {
140	if (argc % 2 == 1) {
141	    quit "odd number integers given";
142	}
143	rlist = list();
144	mlist = list();
145	for (i=1; i <= argc; i+=2) {
146	    push(rlist, param(i));
147	    push(mlist, param(i+1));
148	}
149    }
150
151    /*
152     * solve the problem found in rlist & mlist
153     */
154    m = 1;
155    z = 0;
156    while (size(rlist)) {
157	r=pop(rlist);
158	y=abs(pop(mlist));
159	if (r==null())
160	    continue;
161	if (m) {
162	    if (y) {
163		d = t = z - r;
164		m = lcm(x=m, y);
165		while (d % y) {
166		    u = x;
167		    x %= y;
168		    swap(x,y);
169		    if (y==0)
170			return;
171		    z += (t *= -u/y);
172		}
173	    } else {
174		if ((r % m) != (z % m))
175		    return;
176		else {
177		    m = 0;
178		    z = r;
179		}
180	    }
181	} else if (((y) && (r % y != z % y)) || (r != z))
182	    return;
183    }
184    if (m) {
185	z %= m;
186	if (z < 0)
187	    z += m;
188    }
189
190    /*
191     * return information as required
192     */
193    if (list_args) {
194	return list(z,m);
195    } else {
196	return z;
197    }
198}
199
200if (config("resource_debug") & 3) {
201    print "chrem(r1,m1 [,r2,m2 ...]) defined";
202    print "chrem(rlist [,mlist]) defined";
203}
204