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