BSD 4_4_Lite1 development

[unix-history] / usr / src / contrib / calc-2.9.3t6 / lib / chrem.cal

/*
* 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>
*/

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";
}