home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / hugs101.zip / hugs101sc.zip / hugsdist / src / bignums.c next >
C/C++ Source or Header  |  1995-03-02  |  18KB  |  690 lines
  1. /* --------------------------------------------------------------------------
  2.  * bignums.c:   Copyright (c) Mark P Jones 1991-1994.   All rights reserved.
  3.  *              See NOTICE for details and conditions of use etc...
  4.  *              Hugs version 1.0 August 1994, derived from Gofer 2.30a
  5.  *
  6.  * Functions for manipulating Haskell Integers (bignums).
  7.  * ------------------------------------------------------------------------*/
  8.  
  9. /*#define DEBUG_BIGNUMS*/
  10.  
  11. Bignum bigInt     Args((Int));
  12. Bignum bigDouble  Args((double));
  13. Cell   bigToInt   Args((Bignum));
  14. Cell   bigToFloat Args((Bignum));
  15. Bignum bigStr      Args((String));
  16. Bignum bigNeg     Args((Bignum));
  17. Bool   bigEven    Args((Bignum));
  18. Int    bigCmp     Args((Bignum,Bignum));
  19. Bignum bigAdd     Args((Bignum,Bignum));
  20. Bignum bigSub     Args((Bignum,Bignum));
  21. Bignum bigMul     Args((Bignum,Bignum));
  22. Bignum bigQrm     Args((Bignum,Bignum));
  23.  
  24. #ifdef DEBUG_BIGNUMS
  25. static Void local bigDump(List ds, Int n) {
  26.     while (nonNull(ds) && 0<n--) {
  27.     printf(":%04d",digitOf(hd(ds)));
  28.     ds = tl(ds);
  29.     }
  30. }
  31. #endif
  32.  
  33. /* --------------------------------------------------------------------------
  34.  * Local function prototypes:
  35.  * ------------------------------------------------------------------------*/
  36.  
  37. static Int    local digitsCmp Args((List,List));
  38. static Bignum local digitsAdd Args((Cell,List,List));
  39. static Bignum local digitsSub Args((List,List));
  40. static List   local digitsQrm Args((List,List));
  41.  
  42. /*---------------------------------------------------------------------------
  43.  * Simple bignum primitives:
  44.  *-------------------------------------------------------------------------*/
  45.  
  46. Bignum bigInt(n)            /* convert Int to bignum       */
  47. Int n; {
  48.     if (n==0)
  49.     return ZERONUM;
  50.     else {
  51.     unsigned long no;
  52.     Cell bn, nx;
  53.     if (n<0) {
  54.         no = (unsigned long)(-n);
  55.         bn = pair(NEGNUM,NIL);
  56.     }
  57.     else {
  58.         no = (unsigned long)(n);
  59.         bn = pair(POSNUM,NIL);
  60.     }
  61.     for (nx=bn; no>0; no/=BIGBASE, nx=tl(nx)) 
  62.         tl(nx) = pair(mkDigit(no%BIGBASE),NIL);
  63.     return bn;
  64.     }
  65. }
  66.  
  67. Bignum bigDouble(d)            /* convert double to bignum       */
  68. double d; {
  69.     if (d==0)
  70.     return ZERONUM;
  71.     else {
  72.     Cell bn, nx;
  73.     if (d<0) {
  74.         d  = (-d);
  75.         bn = pair(NEGNUM,NIL);
  76.     }
  77.     else
  78.         bn = pair(POSNUM,NIL);
  79.         for (d=floor(d), nx=bn; d>0; nx=tl(nx)) {
  80.         double n = fmod(d,(double)(BIGBASE));
  81.         tl(nx)   = pair(mkDigit(((Int)(n))),NIL);
  82.         d        = (d - n) / BIGBASE;
  83.     }
  84.     return bn;
  85.     }
  86. }
  87.  
  88. Cell bigToInt(n)            /* convert bignum to Int       */
  89. Bignum n; {
  90.     if (n==ZERONUM)
  91.     return mkInt(0);
  92.     else {
  93.     Int  m  = 0;
  94.     List ds = snd(n);
  95.     List rs = NIL;
  96.     for (; nonNull(ds); ds=tl(ds))
  97.         rs = cons(hd(ds),rs);
  98.     for (; nonNull(rs); rs=tl(rs)) {
  99.         Int d = digitOf(hd(rs));
  100.         if (m > ((MAXPOSINT - d)/BIGBASE))
  101.         return NIL;
  102.         else
  103.         m = m * BIGBASE + d;
  104.     }
  105.     return mkInt(fst(n)==POSNUM ? m : (-m));
  106.     }
  107. }
  108.  
  109. Cell bigToFloat(n)            /* convert bignum to float       */
  110. Bignum n; {
  111.     if (n==ZERONUM)
  112.     return mkFloat((Float)0);
  113.     else {
  114.     Float m  = 0;
  115.     List  ds = snd(n);
  116.     List  rs = NIL;
  117.     for (; nonNull(ds); ds=tl(ds))
  118.         rs = cons(hd(ds),rs);
  119.     for (; nonNull(rs); rs=tl(rs))
  120.         m = m * BIGBASE + digitOf(hd(rs));
  121.     return mkFloat(fst(n)==POSNUM ? m : (-m));
  122.     }
  123. }
  124.  
  125. Bignum bigStr(s)            /* convert String to bignum       */
  126. String s; {
  127.     List   ds = NIL;
  128.     String t  = s;
  129.     Int    i;
  130.  
  131.     if (*t == '-')            /* look for leading minus       */
  132.     t++;                /* and ignore for time being       */
  133.     if (i=strlen(t)%BIGEXP) {
  134.     Int d = 0;
  135.     while (0<i--)
  136.         d = d*10 + (*t++ - '0');
  137.     if (d)
  138.         ds = cons(mkDigit(d),NIL);
  139.     }
  140.     while (*t) {            /* now scan rest of string       */
  141.     Int d = 0;
  142.     for (i=BIGEXP; i>0; i--)
  143.         d = d*10 + (*t++ - '0');
  144.     if (nonNull(ds) || d)
  145.         ds = cons(mkDigit(d),ds);
  146.     }
  147.  
  148.     return isNull(ds) ? ZERONUM : pair((*s=='-' ? NEGNUM : POSNUM), ds);
  149. }
  150.  
  151. Cell bigOut(a,s,b)            /* bignum output, prepend digits to*/
  152. Bignum a;                /* stream s               */
  153. Cell   s;
  154. Bool   b; {                /* TRUE => wrap neg int in parens  */
  155.     if (a==ZERONUM)
  156.     return ap(consChar('0'),s);
  157.     else {
  158.     Bool isNeg = fst(a)==NEGNUM;    /* keep record of sign           */
  159.     List ds    = snd(a);        /* list of digits           */
  160.  
  161.     if (b && isNeg)            /* print closing paren           */
  162.         s = ap(consChar(')'),s);
  163.     
  164.     for (;;) {
  165.         Int d = digitOf(hd(ds));    /* get next digit           */
  166.         ds    = tl(ds);        /* move to next digit           */
  167.         if (nonNull(ds)) {        /* more digits to come after this  */
  168.         Int i = BIGEXP;
  169.         for (; i>0; i--, d/=10)
  170.             s = ap(consChar('0'+(d%10)),s);
  171.         }
  172.         else {            /* leading (non-zero) digits       */
  173.         for (; d; d/=10)
  174.             s = ap(consChar('0'+(d%10)),s);
  175.         break;
  176.         }
  177.     }
  178.  
  179.     if (isNeg)            /* print minus sign           */
  180.         s = ap(consChar('-'),s);
  181.     if (b && isNeg)            /* print open paren           */
  182.         s = ap(consChar('('),s);
  183.     return s;
  184.     }
  185. }
  186.  
  187. Bignum bigNeg(a)            /* unary negation           */
  188. Bignum a; {
  189.     if (a==ZERONUM)
  190.     return ZERONUM;
  191.     else
  192.     return pair(((fst(a)==POSNUM) ? NEGNUM : POSNUM), snd(a));
  193. }
  194.  
  195. Bool bigEven(a)                /* test for even number           */
  196. Bignum a; {                /* ASSUMES BIGBASE is EVEN!       */
  197.     return a==ZERONUM || (digitOf(hd(snd(a))) % 2 == 0) ;
  198. }
  199.  
  200. /*---------------------------------------------------------------------------
  201.  * Bignum comparison routines:
  202.  *-------------------------------------------------------------------------*/
  203.  
  204. Int bigCmp(a,b)                /* Compare bignums returning:       */
  205. Bignum a, b; {                /* -1 if a<b,  +1 if a>b,  0 o/w   */
  206.     if (a==ZERONUM)
  207.     return (b==ZERONUM) ? 0 : ((fst(b)==POSNUM) ? (-1) : 1);
  208.     else if (fst(a)==NEGNUM)
  209.     if (b==ZERONUM || fst(b)==POSNUM)
  210.         return (-1);
  211.     else
  212.         return digitsCmp(snd(b),snd(a));
  213.     else
  214.     if (b==ZERONUM || fst(b)==NEGNUM)
  215.         return 1;
  216.     else
  217.         return digitsCmp(snd(a),snd(b));
  218. }
  219.  
  220. static Int local digitsCmp(xs,ys)    /* Compare positive digit streams  */
  221. List xs, ys; {                /* -1 if xs<ys, +1 if xs>ys, 0 if= */
  222.     Int s = 0;
  223.     for (; nonNull(xs) && nonNull(ys); xs=tl(xs), ys=tl(ys)) {
  224.     Int x = hd(xs);
  225.     Int y = hd(ys);
  226.     if (x<y)
  227.         s = (-1);
  228.     else if (x>y)
  229.         s = 1;
  230.     }
  231.     return (nonNull(xs) ? 1 : (nonNull(ys) ? (-1) : s));
  232. }
  233.  
  234. /*---------------------------------------------------------------------------
  235.  * Addition and subtraction:
  236.  *-------------------------------------------------------------------------*/
  237.  
  238. static Bignum local digitsAdd(sign,xs,ys)/* Addition of digit streams       */
  239. Cell sign;
  240. List xs, ys; {
  241.     Cell bn = pair(sign,NIL);
  242.     Cell nx = bn;
  243.     Int  c  = 0;
  244.     for (;;) {
  245.     if (nonNull(xs)) {            /* Add any digits to carry */
  246.         if (nonNull(ys)) {
  247.         c += digitOf(hd(xs)) + digitOf(hd(ys));
  248.         xs = tl(xs);
  249.         ys = tl(ys);
  250.         }
  251.         else if (c==0) {            /* look for short cut when */
  252.         tl(nx) = xs;            /* a stream ends and there */
  253.         break;                /* is no outstanding carry */
  254.         }
  255.         else {
  256.         c += digitOf(hd(xs));
  257.         xs = tl(xs);
  258.         }
  259.     }
  260.     else if (c==0) {
  261.         tl(nx) = ys;
  262.         break;
  263.     }
  264.     else if (nonNull(ys)) {
  265.         c += digitOf(hd(ys));
  266.         ys = tl(ys);
  267.     }
  268.  
  269.     if (c>=BIGBASE) {            /* Calculate output digit  */
  270.         nx = tl(nx) = cons(mkDigit(c-BIGBASE),NIL);
  271.         c  = 1;
  272.         }
  273.     else {                    /* Carry will always be >0 */
  274.         nx = tl(nx) = cons(mkDigit(c),NIL);    /* at this point       */
  275.         c  = 0;
  276.     }
  277.     }
  278.     return bn;
  279. }
  280.  
  281. static Bignum local digitsSub(xs,ys)    /* Subtraction of digit streams    */
  282. List xs, ys; {
  283.     Cell bn;
  284.     Cell nx;
  285.     Int  b  = 0;
  286.     Int  lz = 0;
  287.  
  288.     switch (digitsCmp(xs,ys)) {
  289.      case (-1) : nx = xs;        /* if xs<ys, return -(ys-xs)       */
  290.             xs = ys;
  291.             ys = nx;
  292.             bn = pair(NEGNUM,NIL);
  293.             break;
  294.  
  295.     case   0  : return ZERONUM;    /* if xs=ys, return 0           */
  296.  
  297.     case   1  : bn = pair(POSNUM,NIL);
  298.             break;        /* if xs>ys, return +(xs-ys)       */
  299.     }
  300.  
  301.     nx = bn;                /* Now we can assume that xs>ys    */
  302.     for (;;) {                /* Scan each digit           */
  303.     Int y = b;
  304.     if (nonNull(ys)) {
  305.         y += digitOf(hd(xs)) - digitOf(hd(ys));
  306.         xs = tl(xs);
  307.         ys = tl(ys);
  308.     }
  309.     else if (y==0) {
  310.         if (nonNull(xs))
  311.         for (; lz>0; lz--)
  312.             nx = tl(nx) = cons(mkDigit(0),NIL);
  313.         tl(nx) = xs;
  314.         break;
  315.     }
  316.     else {
  317.         y += digitOf(hd(xs));    /* xs>ys, so we can't run out of   */
  318.         xs = tl(xs);        /* digits of xs while y!=0       */
  319.     }
  320.  
  321.     if (y<0) {            /* Calculate output digit       */
  322.         y += BIGBASE;
  323.         b  = (-1);
  324.         }
  325.     else
  326.         b  = 0;
  327.  
  328.     if (y==0)            /* Don't insert leading zeros       */
  329.         lz++;
  330.     else {
  331.         for (; lz>0; lz--)
  332.         nx = tl(nx) = cons(mkDigit(0),NIL);
  333.         nx = tl(nx) = cons(mkDigit(y),NIL);
  334.     }
  335.     }
  336.     return bn;
  337. }
  338.  
  339. Bignum bigAdd(a,b)            /* Bignum addition           */
  340. Bignum a, b; {
  341.     if (a==ZERONUM)
  342.     return b;
  343.     else if (b==ZERONUM)
  344.     return a;
  345.     else if (fst(a)==POSNUM)
  346.     if (fst(b)==POSNUM)
  347.         return digitsAdd(POSNUM,snd(a),snd(b));
  348.     else
  349.         return digitsSub(snd(a),snd(b));
  350.     else /* fst(a)==NEGNUM */
  351.     if (fst(b)==NEGNUM)
  352.         return digitsAdd(NEGNUM,snd(a),snd(b));
  353.     else
  354.         return digitsSub(snd(b),snd(a));
  355. }
  356.  
  357. Bignum bigSub(a,b)            /* Bignum subtraction           */
  358. Bignum a, b; {
  359.     if (a==ZERONUM)
  360.     return bigNeg(b);
  361.     else if (b==ZERONUM)
  362.     return a;
  363.     else if (fst(a)==POSNUM)
  364.     if (fst(b)==NEGNUM)
  365.         return digitsAdd(POSNUM,snd(a),snd(b));
  366.     else
  367.         return digitsSub(snd(a),snd(b));
  368.     else /* fst(a)==NEGNUM */
  369.     if (fst(b)==POSNUM)
  370.         return digitsAdd(NEGNUM,snd(a),snd(b));
  371.     else
  372.         return digitsSub(snd(b),snd(a));
  373. }
  374.  
  375. /*---------------------------------------------------------------------------
  376.  * Multiplication:
  377.  *-------------------------------------------------------------------------*/
  378.  
  379. Bignum bigMul(a,b)        /* Calculate product of bignums a and b       */
  380. Bignum a, b; {
  381.     if (a==ZERONUM || b==ZERONUM)    /* if either operand is zero, then */
  382.     return ZERONUM;            /* so is the result ...           */
  383.     else {                /* otherwise, use rule of signs:   */
  384.     Bignum bn = ap((hd(a)==hd(b) ? POSNUM : NEGNUM), NIL);
  385.     List   nx = bn;
  386.     for (; nonNull(b=tl(b)); nx=tl(nx)) {    /* loop through digits of b*/
  387.         List as = nx;        /* At each stage of the loop, add  */
  388.         List xs = tl(a);        /* y * xs to the value in result,  */
  389.         Int  y  = digitOf(hd(b));    /* using c as carry           */
  390.         Int  c  = 0;
  391.         for (; nonNull(xs); xs=tl(xs)) {    /* loop through digits of a*/
  392.         c += digitOf(hd(xs)) * y;
  393.         if (nonNull(tl(as))) {
  394.             as = tl(as);
  395.             c += digitOf(hd(as));
  396.         }
  397.         else
  398.             as = tl(as) = cons(NIL,NIL);
  399.         hd(as) = mkDigit(c % BIGBASE);
  400.         c     /= BIGBASE;
  401.         }
  402.         if (c>0)            /* add carry digit, if required       */
  403.         tl(as) = cons(mkDigit(c),NIL);
  404.     }
  405.     return bn;
  406.     }
  407. }
  408.  
  409. /*---------------------------------------------------------------------------
  410.  * Division:
  411.  *-------------------------------------------------------------------------*/
  412.  
  413. static List local bigRem = NIL;        /* used to return remainder       */
  414.  
  415. static List local digitsQrm(us,vs)    /* digits quotient and remainder   */
  416. List us, vs; {
  417.     if (isNull(tl(vs))) {        /* single digit divisor           */
  418.     Int  v   = digitOf(hd(vs));
  419.     Int  r   = 0;
  420.     List us1 = NIL;            /* first, copy and reverse us       */
  421.     for (; nonNull(us); us=tl(us))
  422.         us1 = cons(hd(us),us1);
  423.     while (nonNull(us1)) {        /* then do division, MSD first       */
  424.         Cell tmp = tl(us1);
  425.         Int  u   = r * BIGBASE + digitOf(hd(us1));
  426.         r        = u % v;
  427.         u         = u / v;
  428.         if (nonNull(us) || u) {    /* output quotient digit       */
  429.         hd(us1) = mkDigit(u);
  430.         tl(us1) = us;
  431.         us     = us1;
  432.         }
  433.         us1         = tmp;
  434.     }
  435.     bigRem = r ? singleton(mkDigit(r)) : NIL;
  436.     return us;
  437.     }
  438.     else {                /* at least two digits in divisor  */
  439.  
  440.     /* The division algorithm used here is, inevitably, based on the
  441.      * description in Knuth's volume 2 on Seminumerical algorithms,
  442.      * and is probably at least as comprehensible as the MIX
  443.      * implementation given there.  It took me a long time to
  444.      * figure that out, so let's just hope that I got it right!
  445.      */
  446.  
  447.     List us1 = NIL;
  448.     List vs1 = NIL;
  449.     List ds  = us;
  450.     Int  v1  = 0, v2  = 0;
  451.     Int  uj  = 0, uj1 = 0, uj2 = 0;
  452.      Int  n   = 0;
  453.     List qs  = NIL;
  454.     Int  sc;
  455.  
  456.     while (nonNull(us) && nonNull(vs)) {
  457.         v2  = v1;
  458.         v1  = digitOf(hd(vs));
  459.         vs1 = cons(hd(vs),vs1);
  460.         vs  = tl(vs);
  461.  
  462.         uj2 = uj1;
  463.         uj1 = digitOf(hd(us));
  464.         us1 = cons(hd(us),us1);
  465.         us  = tl(us);
  466.  
  467.         n++;
  468.     }
  469.  
  470.     if (nonNull(vs)) {        /* if us is shorter than vs, then  */
  471.         bigRem = ds;        /* remainder is us, quotient zero  */
  472.         return NIL;
  473.     }
  474.     vs = rev(vs1);
  475.  
  476.     /* Now we have:
  477.      *   n      = number of digits in vs which is at least two (we
  478.      *          also know that us has at least n digits),
  479.      *   v1, v2 = most significant digits of vs
  480.      *   vs     = digits of vs with least significant digit first
  481.      */
  482.  
  483. #ifdef DEBUG_BIGNUMS
  484. printf("initial vs (n=%d, v1=%d, v2=%d): ",n,v1,v2);
  485. bigDump(vs,n);
  486. putchar('\n');
  487. #endif
  488.     while (nonNull(us)) {
  489.         uj2 = uj1;
  490.         uj1 = digitOf(hd(us));
  491.         us1 = cons(hd(us),us1);
  492.         us  = tl(us);
  493.     }
  494.  
  495.     us = cons(mkDigit(uj=0),NIL);
  496.     ds = us1;
  497.     for (vs1=tl(vs); nonNull(vs1); vs1=tl(vs1)) {
  498.         us1    = tl(ds);
  499.         tl(ds) = us;
  500.         us     = ds;
  501.         ds     = us1;
  502.     }
  503.  
  504.     /* And, at this point, we have:
  505.      *   us = first (n-1) significant digits of original numerator,
  506.      *      with least significant digit first, and a zero at the
  507.      *      end (MSD) of the list, so that length us == n.
  508.      *   ds = remaining digits of the numerator, most significant
  509.      *      digit first.
  510.      *   uj, uj1, uj2
  511.      *      = first three significant digits of us.  (At this point,
  512.      *      uj is actually zero.)
  513.      */
  514.  
  515. #ifdef DEBUG_BIGNUMS
  516. printf("initial us (uj=%d, uj1=%d, uj2=%d): ",uj,uj1,uj2);
  517. bigDump(us,n);
  518. putchar('\n');
  519. printf("initial ds: ");
  520. bigDump(ds,1000);
  521. putchar('\n');
  522. #endif
  523.         sc = BIGBASE / (v1+1);
  524. #ifdef DEBUG_BIGNUMS
  525. printf("scaling factor %d\n",sc);
  526. #endif
  527.     if (sc!=1) {            /* Scale numerator and denominator */
  528.         Int c = 0;
  529.         v1 = v2 = 0;
  530.         for (vs1=vs; nonNull(vs1); vs1=tl(vs1)) {
  531.         v2 = v1;
  532.         v1 = sc * digitOf(hd(vs1)) + c;
  533.         c  = v1 / BIGBASE;
  534.         hd(vs1) = mkDigit(v1%=BIGBASE);
  535.         }                /* no carry here, guaranteed       */
  536.  
  537.         c = uj = uj1 = uj2 = 0;
  538.         for (us1=ds=rev(ds); nonNull(us1); us1=tl(us1)) {
  539.         uj2 = uj1;
  540.         uj1 = uj;
  541.         uj  = sc * digitOf(hd(us1)) + c;
  542.         c   = uj / BIGBASE;
  543.         hd(us1) = mkDigit(uj%=BIGBASE);
  544.         }
  545.         for (ds=rev(ds), us1=us; nonNull(us1); us1=tl(us1)) {
  546.         uj2 = uj1;
  547.         uj1 = uj;
  548.         uj  = sc * digitOf(hd(us1)) + c;
  549.         c   = uj / BIGBASE;
  550.         hd(us1) = mkDigit(uj%=BIGBASE);
  551.         }                /* no carry here, guaranteed       */
  552.     }
  553.  
  554. #ifdef DEBUG_BIGNUMS
  555. printf("scaled vs (n=%d, v1=%d, v2=%d): ",n,v1,v2);
  556. bigDump(vs,n);
  557. putchar('\n');
  558. printf("scaled us (uj=%d, uj1=%d, uj2=%d): ",uj,uj1,uj2);
  559. bigDump(us,n);
  560. putchar('\n');
  561. printf("scaled ds: ");
  562. bigDump(ds,1000);
  563. putchar('\n');
  564. #endif
  565.     /* Most of the conditions above are still valid, except that both
  566.      * the numerator and denominator have been multiplied by the scaling
  567.      * factor sc, and the values of the various digit positions have been
  568.      * updated accordingly.
  569.      *
  570.      * Now we can start the division algorithm proper:
  571.      */
  572.  
  573.     while (nonNull(ds)) {
  574.         Int c;            /* Guess a value for quotient digit*/
  575.         Int qhat = (uj==v1) ? (BIGBASE-1) : (uj*BIGBASE+uj1)/v1;
  576.         while (v2*qhat > (uj*BIGBASE+uj1-qhat*v1)*BIGBASE+uj2)
  577.         qhat--;            /* and then try to improve it       */
  578.  
  579.         us1    = tl(ds);        /* take digit off off front of ds  */
  580.         tl(ds) = us;        /* and add to front of us       */
  581.         us     = ds;
  582.         ds       = us1;
  583.  
  584. #ifdef DEBUG_BIGNUMS
  585. printf("To divide us (uj=%d, uj1=%d, uj2=%d): ",uj,uj1,uj2);
  586. bigDump(us,n+1);
  587. printf(" by vs (v1=%d, v2=%d): ",v1,v2);
  588. bigDump(vs,n);
  589. printf(", guess qhat=%d\n",qhat);
  590. #endif
  591.         uj  = uj1 = uj2 = c = 0;    /* us := us - qhat * vs           */
  592.         us1 = us;
  593.         vs1 = vs;
  594.         do {
  595.         uj2 = uj1;
  596.         uj1 = uj;
  597.         uj  = digitOf(hd(us1)) - qhat*digitOf(hd(vs1)) - c;
  598.         if (uj>=0)
  599.             c = 0;
  600.         else {
  601.             c   = (BIGBASE - 1 - uj) / BIGBASE;
  602.             uj += c*BIGBASE;
  603.         }
  604.         hd(us1) = mkDigit(uj);
  605.         us1     = tl(us1);
  606.         vs1     = tl(vs1);
  607.         } while (nonNull(vs1));
  608.  
  609.         if (digitOf(hd(us1))<c) {    /* maybe we overestimated qhat?       */
  610. #ifdef DEBUG_BIGNUMS
  611. printf("overestimated qhat by one!\n");
  612. #endif
  613.         qhat--;            /* so reduce guess           */
  614.         uj  = uj1 = uj2 = c = 0;/* and set us := us + vs       */
  615.         us1 = us;        /* (we can't have overestimated by */
  616.         vs1 = vs;        /* more than one thanks to the       */
  617.         do {            /* initial scaling of us, vs by sc)*/
  618.             uj2 = uj1;
  619.             uj1 = uj;
  620.             uj  = digitOf(hd(us1)) + digitOf(hd(vs1)) + c;
  621.             if (uj>=BIGBASE)
  622.             c = 1, uj -= BIGBASE;
  623.             else
  624.             c = 0;
  625.             hd(us1) = mkDigit(uj);
  626.             us1     = tl(us1);
  627.             vs1     = tl(vs1);
  628.         } while (nonNull(vs1));
  629.         }
  630. #ifdef DEBUG_BIGNUMS
  631. printf("There remains (uj=%d, uj1=%d, uj2=%d): ",uj,uj1,uj2);
  632. bigDump(us,n);
  633. putchar('\n');
  634. #endif
  635.  
  636.         if (nonNull(qs) || qhat)    /* output quotient digit, without  */
  637.         qs = cons(mkDigit(qhat),qs);    /* leading zeros       */
  638.     }
  639.  
  640. #ifdef DEBUG_BIGNUMS
  641. printf("done quotient\n");
  642. #endif
  643.     /* Now we have the quotient digits (if any) with least significant
  644.      * digit first in qs, and sc times the remainder is the first n
  645.      * digits of us.  All that remains is to adjust the remainder:
  646.      */
  647.  
  648.     us1 = rev(take(n,us));
  649.     us  = NIL;
  650.     uj  = 0;            /* reuse variable uj as a carry       */
  651.     while (nonNull(us1)) {
  652.         Int y = uj * BIGBASE + digitOf(hd(us1));
  653.         uj    = y % sc;
  654.         y    /= sc;
  655.         if (nonNull(us) || y) {
  656.         vs1     = tl(us1);
  657.         tl(us1) = us;
  658.         hd(us1) = mkDigit(y);
  659.         us      = us1;
  660.         us1    = vs1;
  661.         }
  662.         else
  663.         us1 = tl(us1);
  664.     }
  665.     bigRem = us;
  666.     return qs;
  667.     }
  668. }
  669.  
  670. Cell bigQrm(a,b)            /* bignum quotient and remainder   */
  671. Bignum a, b; {
  672.     if (a==ZERONUM)
  673.     return ap(ap(mkTuple(2),ZERONUM),ZERONUM);
  674.     else if (b!=ZERONUM) {
  675.     /* The sign of the remainder is positive if numerator and denominator
  676.      * have the same sign, negative if the signs differ.  The sign of the
  677.      * quotient is always the same as the sign of the numerator.
  678.      */
  679.     Cell qs = digitsQrm(snd(a),snd(b));
  680.     Cell rs = bigRem;
  681.     qs      = isNull(qs) ? ZERONUM
  682.                  : pair((fst(a)==fst(b) ? POSNUM : NEGNUM), qs);
  683.     rs      = isNull(rs) ? ZERONUM : pair(hd(a),rs);
  684.     return ap(ap(mkTuple(2),qs),rs);
  685.     }
  686.     return NIL;                /* indicates division by zero       */
  687. }
  688.  
  689. /*-------------------------------------------------------------------------*/
  690.