home *** CD-ROM | disk | FTP | other *** search
/ Skunkware 5 / Skunkware 5.iso / lib / calc / chrem.cal < prev    next >
Encoding:
Text File  |  1995-07-17  |  4.2 KB  |  182 lines
  1. /*
  2.  * chrem - Chinese remainder theorem/problem solver
  3.  *
  4.  * When possible, chrem finds solutions for x of a set of congruences
  5.  * of the form:
  6.  *
  7.  *            x = r1 (mod m1)
  8.  *            x = r2 (mod m2)
  9.  *               ...
  10.  *
  11.  * where the residues r1, r2, ... and the moduli m1, m2, ... are
  12.  * given integers.  The Chinese remainder theorem states that if
  13.  * m1, m2, ... are relatively prime in pairs, the above congruences
  14.  * have a unique solution modulo  m1 * m2 * ...   If m1, m2, ...
  15.  * are not relatively prime in pairs, it is possible that no solution
  16.  * exists.  If solutions exist, the general solution is expressible as:
  17.  *
  18.  *                   x = r (mod m)
  19.  *
  20.  * where m = lcm(m1,m2,...), and if m > 0, 0 <= r < m.  This
  21.  * solution may be interpreted as:
  22.  *
  23.  *                  x = r + k * m            [[NOTE 1]]
  24.  *
  25.  * where k is an arbitrary integer.
  26.  *
  27.  ***
  28.  *
  29.  * usage:
  30.  *
  31.  *    chrem(r1,m1 [,r2,m2, ...])
  32.  *
  33.  *         r1, r2, ...      remainder integers or null values
  34.  *         m1, m2, ...      moduli integers
  35.  *
  36.  *    chrem(r_list, [m_list])
  37.  *
  38.  *        r_list    list (r1,r2, ...)
  39.  *        m_list    list (m1,m2, ...)
  40.  *
  41.  *         If m_list is omitted, then 'defaultmlist' is used.
  42.  *        This default list is a global value that may be changed
  43.  *        by the user.  Initially it is the first 8 primes.
  44.  *
  45.  * If a remainder is null(), then the corresponding congruence is 
  46.  * ignored.  This is useful when working with a fixed list of moduli.  
  47.  * 
  48.  * If there are more remainders than moduli, then the later moduli are 
  49.  * ignored.
  50.  *
  51.  * The moduli may be any integers, not necessarily relatively prime in 
  52.  * pairs (as required for the Chinese remainder theorem).  Any moduli 
  53.  * may be zero;  x = r (mod 0) has the meaning of x = r.
  54.  *
  55.  * returns:
  56.  *
  57.  *    If args were integer pairs:
  58.  *
  59.  *          r        ('r' is defined above, see [[NOTE 1]])
  60.  *
  61.  *    If 1 or 2 list args were given:
  62.  *
  63.  *      (r, m)    ('r' and 'm' are defined above, see [[NOTE 1]])
  64.  *
  65.  * NOTE: In all cases, null() is returned if there is no solution.
  66.  *
  67.  ***
  68.  *
  69.  * This function may be used to solve the following historical problems:
  70.  *
  71.  *   Sun-Tsu, 1st century A.D.  
  72.  *
  73.  *    To find a number for which the reminders after division by 3, 5, 7 
  74.  *    are 2, 3, 2, respectively:
  75.  *
  76.  *        chrem(2,3,3,5,2,7) ---> 23
  77.  *
  78.  *   Fibonacci, 13th century A.D.
  79.  *
  80.  *    To find a number divisible by 7 which leaves remainder 1 when 
  81.  *    divided by 2, 3, 4, 5, or 6:
  82.  *
  83.  *
  84.  *        chrem(list(0,1,1,1,1,1),list(7,2,3,4,5,6)) ---> (301,420)
  85.  *
  86.  *    i.e., any value that is 301 mod 420.
  87.  *
  88.  * Written by: Ernest W Bowen <ernie@neumann.une.edu.au>
  89.  * Interface by: Landon Curt Noll <chongo@toad.com>
  90.  */
  91.  
  92. static defaultmlist = list(2,3,5,7,11,13,17,19); /* The first eight primes */
  93.  
  94. define chrem()
  95. {
  96.     local argc;        /* number of args given */
  97.     local rlist;    /* reminder list - ri */
  98.     local mlist;    /* modulus list - mi */
  99.     local list_args;      /* true => args given are lists, not r1,m1, ... */
  100.     local m,z,r,y,d,t,x,u,i;
  101.  
  102.     /* 
  103.      * parse args 
  104.      */
  105.     argc = param(0);
  106.     if (argc == 0) {
  107.     quit "usage: chrem(r1,m1 [,r2,m2 ...]) or chrem(r_list, m_list)";
  108.     }
  109.     list_args = islist(param(1));
  110.     if (list_args) {
  111.     rlist = param(1);
  112.     mlist = (argc == 1) ? defaultmlist : param(2);
  113.     if (size(rlist) > size(mlist)) {
  114.         quit "too many residues";
  115.     }
  116.     } else {
  117.     if (argc % 2 == 1) {
  118.         quit "odd number integers given";
  119.     }
  120.     rlist = list();
  121.     mlist = list();
  122.     for (i=1; i <= argc; i+=2) {
  123.         push(rlist, param(i));
  124.             push(mlist, param(i+1));
  125.     }
  126.     }
  127.  
  128.     /* 
  129.      * solve the problem found in rlist & mlist 
  130.      */
  131.     m = 1; 
  132.     z = 0;
  133.     while (size(rlist)) { 
  134.     r=pop(rlist); 
  135.     y=abs(pop(mlist));
  136.     if (r==null()) 
  137.         continue;
  138.     if (m) {
  139.         if (y) {
  140.         d = t = z - r;
  141.         m = lcm(x=m, y);
  142.         while (d % y) { 
  143.             u = x; 
  144.             x %= y; 
  145.             swap(x,y);
  146.             if (y==0) 
  147.             return;
  148.             z += (t *= -u/y);
  149.         }
  150.         } else { 
  151.         if ((r % m) != (z % m)) 
  152.             return;
  153.         else {
  154.             m = 0; 
  155.             z = r;
  156.         }
  157.         }
  158.     } else if (((y) && (r % y != z % y)) || (r != z)) 
  159.         return;
  160.     }
  161.     if (m) { 
  162.     z %= m; 
  163.     if (z < 0) 
  164.         z += m;
  165.     }
  166.  
  167.     /* 
  168.      * return information as required 
  169.      */
  170.     if (list_args) {
  171.     return list(z,m);
  172.     } else {
  173.     return z;
  174.     }
  175. }
  176.  
  177. global lib_debug;
  178. if (lib_debug >= 0) {
  179.     print "chrem(r1,m1 [,r2,m2 ...]) defined";
  180.     print "chrem(rlist [,mlist]) defined";
  181. }
  182.