home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / gnu / g__lib / integer.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1993-07-23  |  41.2 KB  |  2,113 lines
  1. /* 
  2. Copyright (C) 1988 Free Software Foundation
  3.     written by Doug Lea (dl@rocky.oswego.edu)
  4.  
  5. This file is part of GNU CC.
  6.  
  7. GNU CC is distributed in the hope that it will be useful,
  8. but WITHOUT ANY WARRANTY.  No author or distributor
  9. accepts responsibility to anyone for the consequences of using it
  10. or for whether it serves any particular purpose or works at all,
  11. unless he says so in writing.  Refer to the GNU CC General Public
  12. License for full details.
  13.  
  14. Everyone is granted permission to copy, modify and redistribute
  15. GNU CC, but only under the conditions described in the
  16. GNU CC General Public License.   A copy of this license is
  17. supposed to have been given to you along with GNU CC so you
  18. can know your rights and responsibilities.  It should be in a
  19. file named COPYING.  Among other things, the copyright notice
  20. and this notice must be preserved on all copies.  
  21. */
  22.  
  23. /*
  24.   Some of the following algorithms are very loosely based on those from 
  25.   MIT C-Scheme bignum.c, which is
  26.       Copyright (c) 1987 Massachusetts Institute of Technology
  27.  
  28.   with other guidance from Knuth, vol. 2
  29.  
  30.   Thanks to the creators of the algorithms.
  31. */
  32.  
  33. #include <Integer.h>
  34. #include <Obstack.h>
  35. #include <std.h>
  36. #include <ctype.h>
  37. #include <math.h>
  38. #include "libconfig.h"
  39.  
  40. /*
  41.   minimum and maximum sizes for an IntRep
  42. */
  43.  
  44. #define MINIntRep_SIZE   16
  45. #define MAXIntRep_SIZE   I_MAXNUM
  46. #define MALLOC_OVERHEAD 4
  47.  
  48.  
  49. // static globals
  50.  
  51. static unsigned short ones[1]   = { 1 }; // to copy a one into a rep
  52.  
  53.  
  54. // utilities to extract and transfer bits 
  55.  
  56. // get low bits
  57.  
  58. inline static unsigned short extract(unsigned long x)
  59. {
  60.   return x & I_MAXNUM;
  61. }
  62.  
  63. // transfer high bits to low
  64.  
  65. inline static unsigned long down(unsigned long x)
  66. {
  67.   return (x >> I_SHIFT) & I_MAXNUM;
  68. }
  69.  
  70. // transfer low bits to high
  71.  
  72. inline static unsigned long up(unsigned long x)
  73. {
  74.   return x << I_SHIFT;
  75. }
  76.  
  77. // compare two equal-length reps
  78.  
  79. inline static int docmp(unsigned short* x, unsigned short* y, int l)
  80. {
  81.   int diff = 0;
  82.   const unsigned short* xs = &(x[l]);
  83.   const unsigned short* ys = &(y[l]);
  84.   while (l-- > 0 && (diff = (*--xs) - (*--ys)) == 0);
  85.   return diff;
  86. }
  87.  
  88. // figure out max length of result of +, -, etc.
  89.  
  90. inline static int calc_len(int len1, int len2, int pad)
  91. {
  92.   return (len1 >= len2)? len1 + pad : len2 + pad;
  93. }
  94.  
  95. // ensure len & sgn are correct
  96.  
  97. inline static void Icheck(IntRep* rep)
  98. {
  99.   int l = rep->len;
  100.   const unsigned short* p = &(rep->s[l]);
  101.   while (l > 0 && *--p == 0) --l;
  102.   if ((rep->len = l) == 0) rep->sgn = I_POSITIVE;
  103. }
  104.  
  105. inline static void Iclear(IntRep* rep)
  106. {
  107.   bzero(rep->s, rep->sz * sizeof(short));
  108. }
  109.  
  110. inline static void Iclear_last(IntRep* rep)
  111. {
  112.   rep->s[rep->len - 1] = 0;
  113. }
  114.  
  115. inline static void Iclear_from(IntRep* rep, int p)
  116. {
  117.   bzero(&(rep->s[p]), (rep->sz - p) * sizeof(short));
  118. }
  119.  
  120. static void nonnil(IntRep* rep)
  121. {
  122.   if (rep == 0) 
  123.     (*lib_error_handler)("Integer", "operation on uninitialized Integer");
  124. }
  125.  
  126. // allocate a new Irep
  127.  
  128. inline static IntRep* Inew(int newlen)
  129. {
  130.   unsigned int siz = sizeof(IntRep) + newlen * sizeof(short) + MALLOC_OVERHEAD;
  131.   unsigned int allocsiz = MINIntRep_SIZE;;
  132.   while (allocsiz < siz) allocsiz <<= 1;
  133.   allocsiz -= MALLOC_OVERHEAD;
  134.   if (allocsiz >= MAXIntRep_SIZE * sizeof(short))
  135.     (*lib_error_handler)("Integer", "Requested length out of range");
  136.     
  137.   IntRep* rep = (IntRep *) new char[allocsiz];
  138.   bzero(rep, allocsiz);
  139.   rep->sz = (allocsiz - sizeof(IntRep) + sizeof(short)) / sizeof(short);
  140.   return rep;
  141. }
  142.  
  143. // allocate: use the bits in src if non-null,
  144.  
  145. IntRep* Ialloc(IntRep* old, const unsigned short* src, int srclen, int newsgn,
  146.               int newlen)
  147. {
  148.   IntRep* rep;
  149.   if (old == 0 || newlen > old->sz)
  150.     rep = Inew(newlen);
  151.   else
  152.     rep = old;
  153.  
  154.   rep->len = newlen;
  155.   rep->sgn = newsgn;
  156.  
  157.   if (srclen != 0) bcopy(src, rep->s, srclen * sizeof(short));
  158.  
  159.   if (old != rep && old != 0) delete old;
  160.   return rep;
  161. }
  162.  
  163. IntRep* Iresize(IntRep* old, int newlen)
  164. {
  165.   IntRep* rep;
  166.   if (old == 0)
  167.   {
  168.     rep = Inew(newlen);
  169.     rep->sgn = I_POSITIVE;
  170.   }
  171.   else if (newlen > old->sz)
  172.   {
  173.     rep = Inew(newlen);
  174.     bcopy(old->s, rep->s, old->len * sizeof(short));
  175.     rep->sgn = old->sgn;
  176.     delete old;
  177.   }
  178.   else
  179.     rep = old;
  180.  
  181.   rep->len = newlen;
  182.  
  183.   return rep;
  184. }
  185.  
  186. // same, for straight copy
  187.  
  188. IntRep* Icopy(IntRep* old, IntRep* src)
  189. {
  190.   IntRep* rep;
  191.   if (old == src) return old; 
  192.   if (src == 0)
  193.   {
  194.     if (old == 0)
  195.       rep = Inew(0);
  196.     rep->len = 0;
  197.     rep->sgn = I_POSITIVE;
  198.   }
  199.   else 
  200.   {
  201.     int newlen = src->len;
  202.     if (old == 0 || newlen > old->sz)
  203.     {
  204.       rep = Inew(newlen);
  205.       if (old != 0) delete old;
  206.     }
  207.     else
  208.       rep = old;
  209.  
  210.     bcopy(src->s, rep->s, newlen * sizeof(short));
  211.     rep->len = newlen;
  212.     rep->sgn = src->sgn;
  213.   }
  214.   return rep;
  215. }
  216.  
  217.  
  218. IntRep* Icopy_long(IntRep* rep, long x)
  219. {
  220.   unsigned short tmp[SHORT_PER_LONG];
  221.   unsigned long u;
  222.   int newsgn;
  223.   if (newsgn = (x >= 0))
  224.     u = x;
  225.   else
  226.     u = -x;
  227.   
  228.   int l = 0;
  229.   while (u != 0)
  230.   {
  231.     tmp[l++] = extract(u);
  232.     u = down(u);
  233.   }
  234.   return Ialloc(rep, tmp, l, newsgn, l);
  235. }
  236.  
  237. inline static IntRep* copy_zero(IntRep* r)
  238. {
  239.   return Ialloc(r, 0, 0, I_POSITIVE, 0);
  240. }
  241.  
  242. inline static IntRep* copy_one(IntRep* r, int sgn)
  243. {
  244.   return Ialloc(r, ones, 1, sgn, 1);
  245. }
  246.  
  247.  
  248. // convert to a legal two's complement long if possible
  249. // if too big, return most negative/positive value
  250.  
  251. long Itolong(IntRep* rep)
  252. { 
  253.   if (rep->len > SHORT_PER_LONG)
  254.     return (rep->sgn == I_POSITIVE) ? MAXLONG : MINLONG;
  255.   else
  256.   {
  257.     unsigned long a = rep->s[SHORT_PER_LONG - 1];
  258.     if (a >= I_MINNUM)
  259.       return (rep->sgn == I_POSITIVE) ? MAXLONG : MINLONG;
  260.     else
  261.     {
  262.       a = up(a) | rep->s[SHORT_PER_LONG - 2];
  263.       if (SHORT_PER_LONG - 3 >= 0)
  264.       {
  265.         for (int i = SHORT_PER_LONG - 3; i >= 0; --i)
  266.           a = up(a) | rep->s[i];
  267.       }
  268.       return (rep->sgn == I_POSITIVE)? a : -((long)a);
  269.     }
  270.   }
  271. }
  272.  
  273.  
  274.  
  275. // Integer support
  276.  
  277.  
  278.  
  279. // test whether op long() will work.
  280. // careful about asymmetry between MINLONG & MAXLONG
  281.  
  282. int Iislong(IntRep* rep)
  283. {
  284.   int l = rep->len;
  285.   if (l < SHORT_PER_LONG)
  286.     return 1;
  287.   else if (l > SHORT_PER_LONG)
  288.     return 0;
  289.   else if (rep->s[SHORT_PER_LONG - 1] < I_MINNUM)
  290.     return 1;
  291.   else if (rep->sgn == I_NEGATIVE && rep->s[SHORT_PER_LONG - 1] == I_MINNUM)
  292.   {
  293.     for (int i = 0; i < SHORT_PER_LONG - 1; ++i)
  294.       if (rep->s[i] != 0)
  295.         return 0;
  296.     return 1;
  297.   }
  298.   else
  299.     return 0;
  300. }
  301.  
  302. // convert to a double 
  303.  
  304. double Itodouble(IntRep* rep)
  305. { 
  306.   double d = 0.0;
  307.   double bound = HUGE / 2.0;
  308.   for (int i = rep->len - 1; i >= 0; --i)
  309.   {
  310.     unsigned short a = I_RADIX >> 1;
  311.     while (a != 0)
  312.     {
  313.       if (d >= bound)
  314.         return (rep->sgn == I_NEGATIVE) ? -HUGE : HUGE;
  315.       d *= 2.0;
  316.       if (rep->s[i] & a)
  317.         d += 1.0;
  318.       a >>= 1;
  319.     }
  320.   }
  321.   if (rep->sgn == I_NEGATIVE)
  322.     return -d;
  323.   else
  324.     return d;
  325. }
  326.  
  327. // see whether op double() will work-
  328. // have to actually try it in order to find out
  329. // since otherwise might trigger fp exception
  330.  
  331. int Iisdouble(IntRep* rep)
  332. {
  333.   double d = 0.0;
  334.   double bound = HUGE / 2.0;
  335.   for (int i = rep->len - 1; i >= 0; --i)
  336.   {
  337.     unsigned short a = I_RADIX >> 1;
  338.     while (a != 0)
  339.     {
  340.       if (d > bound || (d == bound && (i > 0 || (rep->s[i] & a))))
  341.         return 0;
  342.       d *= 2.0;
  343.       if (rep->s[i] & a)
  344.         d += 1.0;
  345.       a >>= 1;
  346.     }
  347.   }
  348.   return 1;
  349. }
  350.  
  351. // real division of num / den
  352.  
  353. double ratio(Integer& num, Integer& den)
  354. {
  355.   Integer q, r;
  356.   divide(num, den, q, r);
  357.   double d1 = double(q);
  358.  
  359.   if (d1 == HUGE || d1 == -HUGE || sign(r) == 0)
  360.     return d1;
  361.   else      // use as much precision as available for fractional part
  362.   {
  363.     double  d2 = 0.0;
  364.     double  d3 = 0.0; 
  365.     double bound = HUGE / 2.0;
  366.     int cont = 1;
  367.     for (int i = den.rep->len - 1; i >= 0 && cont; --i)
  368.     {
  369.       unsigned short a = I_RADIX >> 1;
  370.       while (a != 0)
  371.       {
  372.         if (d2 >= bound)
  373.         {
  374.           cont = 0;
  375.           break;
  376.         }
  377.         d2 *= 2.0;
  378.         if (den.rep->s[i] & a)
  379.           d2 += 1.0;
  380.  
  381.         if (i < r.rep->len)
  382.         {
  383.           d3 *= 2.0;
  384.           if (r.rep->s[i] & a)
  385.             d3 += 1.0;
  386.         }
  387.         a >>= 1;
  388.       }
  389.     }
  390.  
  391.     if (sign(r) < 0)
  392.       d3 = -d3;
  393.     return d1 + d3 / d2;
  394.   }
  395. }
  396.  
  397. // comparison functions
  398.   
  399. int compare(IntRep* x, IntRep* y)
  400. {
  401.   int diff  = x->sgn - y->sgn;
  402.   if (diff == 0)
  403.   {
  404.     diff = x->len - y->len;
  405.     if (diff == 0)
  406.       diff = docmp(x->s, y->s, x->len);
  407.     if (x->sgn == I_NEGATIVE)
  408.       diff = -diff;
  409.   }
  410.   return diff;
  411. }
  412.  
  413. int ucompare(IntRep* x, IntRep* y)
  414. {
  415.   int diff = x->len - y->len;
  416.   if (diff == 0)
  417.   {
  418.     int l = x->len;
  419.     const unsigned short* xs = &(x->s[l]);
  420.     const unsigned short* ys = &(y->s[l]);
  421.     while (l-- > 0 && (diff = (*--xs) - (*--ys)) == 0);
  422.   }
  423.   return diff;
  424. }
  425.  
  426. int compare(IntRep* x, long  y)
  427. {
  428.   int xl = x->len;
  429.   int xsgn = x->sgn;
  430.   if (y == 0)
  431.   {
  432.     if (xl == 0)
  433.       return 0;
  434.     else if (xsgn == I_NEGATIVE)
  435.       return -1;
  436.     else
  437.       return 1;
  438.   }
  439.   else
  440.   {
  441.     int ysgn = y >= 0;
  442.     unsigned long uy = (ysgn)? y : -y;
  443.     int diff = xsgn - ysgn;
  444.     if (diff == 0)
  445.     {
  446.       diff = xl - SHORT_PER_LONG;
  447.       if (diff <= 0)
  448.       {
  449.         unsigned short tmp[SHORT_PER_LONG];
  450.         int yl = 0;
  451.         while (uy != 0)
  452.         {
  453.           tmp[yl++] = extract(uy);
  454.           uy = down(uy);
  455.         }
  456.         diff = xl - yl;
  457.         if (diff == 0)
  458.           diff = docmp(x->s, tmp, xl);
  459.       }
  460.     }
  461.     if (xsgn == I_NEGATIVE)
  462.       diff = -diff;
  463.     return diff;
  464.   }
  465. }
  466.  
  467. int ucompare(IntRep* x, long  y)
  468. {
  469.   int xl = x->len;
  470.   if (y == 0)
  471.     return xl;
  472.   else
  473.   {
  474.     unsigned long uy = (y >= 0)? y : -y;
  475.     int diff = xl - SHORT_PER_LONG;
  476.     if (diff <= 0)
  477.     {
  478.       unsigned short tmp[SHORT_PER_LONG];
  479.       int yl = 0;
  480.       while (uy != 0)
  481.       {
  482.         tmp[yl++] = extract(uy);
  483.         uy = down(uy);
  484.       }
  485.       diff = xl - yl;
  486.       if (diff == 0)
  487.         diff = docmp(x->s, tmp, xl);
  488.     }
  489.     return diff;
  490.   }
  491. }
  492.  
  493.  
  494.  
  495. // arithmetic functions
  496.  
  497. IntRep* add(IntRep* x, int negatex, IntRep* y, int negatey, IntRep* r)
  498. {
  499.   nonnil(x);
  500.   nonnil(y);
  501.  
  502.   int xl = x->len;
  503.   int yl = y->len;
  504.  
  505.   int xsgn = (negatex && xl != 0) ? !x->sgn : x->sgn;
  506.   int ysgn = (negatey && yl != 0) ? !y->sgn : y->sgn;
  507.  
  508.   int xrsame = x == r;
  509.   int yrsame = y == r;
  510.  
  511.   if (yl == 0)
  512.     r = Ialloc(r, x->s, xl, xsgn, xl);
  513.   else if (xl == 0)
  514.     r = Ialloc(r, y->s, yl, ysgn, yl);
  515.   else if (xsgn == ysgn)
  516.   {
  517.     r = Iresize(r, calc_len(xl, yl, 1));
  518.     r->sgn = xsgn;
  519.     Iclear_last(r);
  520.     unsigned short* rs = r->s;
  521.     unsigned short* as;
  522.     unsigned short* bs;
  523.     unsigned short* topa;
  524.     unsigned short* topb;
  525.     if (xl >= yl)
  526.     {
  527.       as =  (xrsame)? r->s : x->s;
  528.       topa = &(as[xl]);
  529.       bs =  (yrsame)? r->s : y->s;
  530.       topb = &(bs[yl]);
  531.     }
  532.     else
  533.     {
  534.       bs =  (xrsame)? r->s : x->s;
  535.       topb = &(bs[xl]);
  536.       as =  (yrsame)? r->s : y->s;
  537.       topa = &(as[yl]);
  538.     }
  539.     unsigned long sum = 0;
  540.     while (bs < topb)
  541.     {
  542.       sum += (unsigned long)(*as++) + (unsigned long)(*bs++);
  543.       *rs++ = extract(sum);
  544.       sum = down(sum);
  545.     }
  546.     while (sum != 0 && as < topa)
  547.     {
  548.       sum += (unsigned long)(*as++);
  549.       *rs++ = extract(sum);
  550.       sum = down(sum);
  551.     }
  552.     if (sum != 0)
  553.       *rs = extract(sum);
  554.     else if (rs != as)
  555.       while (as < topa)
  556.         *rs++ = *as++;
  557.   }
  558.   else
  559.   {
  560.     int comp = ucompare(x, y);
  561.     if (comp == 0)
  562.       r = copy_zero(r);
  563.     else
  564.     {
  565.       r = Iresize(r, calc_len(xl, yl, 0));
  566.       unsigned short* rs = r->s;
  567.       unsigned short* as;
  568.       unsigned short* bs;
  569.       unsigned short* topa;
  570.       unsigned short* topb;
  571.       if (comp > 0)
  572.       {
  573.         as =  (xrsame)? r->s : x->s;
  574.         topa = &(as[xl]);
  575.         bs =  (yrsame)? r->s : y->s;
  576.         topb = &(bs[yl]);
  577.         r->sgn = xsgn;
  578.       }
  579.       else
  580.       {
  581.         bs =  (xrsame)? r->s : x->s;
  582.         topb = &(bs[xl]);
  583.         as =  (yrsame)? r->s : y->s;
  584.         topa = &(as[yl]);
  585.         r->sgn = ysgn;
  586.       }
  587.       unsigned long hi = 1;
  588.       while (bs < topb)
  589.       {
  590.         hi += (unsigned long)(*as++) + I_MAXNUM - (unsigned long)(*bs++);
  591.         *rs++ = extract(hi);
  592.         hi = down(hi);
  593.       }
  594.       while (hi == 0 && as < topa)
  595.       {
  596.         hi = (unsigned long)(*as++) + I_MAXNUM;
  597.         *rs++ = extract(hi);
  598.         hi = down(hi);
  599.       }
  600.       if (rs != as)
  601.         while (as < topa)
  602.           *rs++ = *as++;
  603.     }
  604.   }
  605.   Icheck(r);
  606.   return r;
  607. }
  608.  
  609.  
  610. IntRep* add(IntRep* x, int negatex, long y, IntRep* r)
  611. {
  612.   nonnil(x);
  613.   int xl = x->len;
  614.   int xsgn = (negatex && xl != 0) ? !x->sgn : x->sgn;
  615.   int xrsame = x == r;
  616.  
  617.   int ysgn = (y >= 0);
  618.   unsigned long uy = (ysgn)? y : -y;
  619.  
  620.   if (y == 0)
  621.     r = Ialloc(r, x->s, xl, xsgn, xl);
  622.   else if (xl == 0)
  623.     r = Icopy_long(r, y);
  624.   else if (xsgn == ysgn)
  625.   {
  626.     r = Iresize(r, calc_len(xl, SHORT_PER_LONG, 1));
  627.     r->sgn = xsgn;
  628.     Iclear_last(r);
  629.     unsigned short* rs = r->s;
  630.     unsigned short* as =  (xrsame)? r->s : x->s;
  631.     unsigned short* topa = &(as[xl]);
  632.     unsigned long sum = 0;
  633.     while (as < topa && uy != 0)
  634.     {
  635.       unsigned long u = extract(uy);
  636.       uy = down(uy);
  637.       sum += (unsigned long)(*as++) + u;
  638.       *rs++ = extract(sum);
  639.       sum = down(sum);
  640.     }
  641.     while (sum != 0 && as < topa)
  642.     {
  643.       sum += (unsigned long)(*as++);
  644.       *rs++ = extract(sum);
  645.       sum = down(sum);
  646.     }
  647.     if (sum != 0)
  648.       *rs = extract(sum);
  649.     else if (rs != as)
  650.       while (as < topa)
  651.         *rs++ = *as++;
  652.   }
  653.   else
  654.   {
  655.     unsigned short tmp[SHORT_PER_LONG];
  656.     int yl = 0;
  657.     while (uy != 0)
  658.     {
  659.       tmp[yl++] = extract(uy);
  660.       uy = down(uy);
  661.     }
  662.     int comp = xl - yl;
  663.     if (comp == 0)
  664.       comp = docmp(x->s, tmp, yl);
  665.     if (comp == 0)
  666.       r = copy_zero(r);
  667.     else
  668.     {
  669.       r = Iresize(r, calc_len(xl, yl, 0));
  670.       unsigned short* rs = r->s;
  671.       unsigned short* as;
  672.       unsigned short* bs;
  673.       unsigned short* topa;
  674.       unsigned short* topb;
  675.       if (comp > 0)
  676.       {
  677.         as =  (xrsame)? r->s : x->s;
  678.         topa = &(as[xl]);
  679.         bs =  tmp;
  680.         topb = &(bs[yl]);
  681.         r->sgn = xsgn;
  682.       }
  683.       else
  684.       {
  685.         bs =  (xrsame)? r->s : x->s;
  686.         topb = &(bs[xl]);
  687.         as =  tmp;
  688.         topa = &(as[yl]);
  689.         r->sgn = ysgn;
  690.       }
  691.       unsigned long hi = 1;
  692.       while (bs < topb)
  693.       {
  694.         hi += (unsigned long)(*as++) + I_MAXNUM - (unsigned long)(*bs++);
  695.         *rs++ = extract(hi);
  696.         hi = down(hi);
  697.       }
  698.       while (hi == 0 && as < topa)
  699.       {
  700.         hi = (unsigned long)(*as++) + I_MAXNUM;
  701.         *rs++ = extract(hi);
  702.         hi = down(hi);
  703.       }
  704.       if (rs != as)
  705.         while (as < topa)
  706.           *rs++ = *as++;
  707.     }
  708.   }
  709.   Icheck(r);
  710.   return r;
  711. }
  712.  
  713.  
  714. IntRep* multiply(IntRep* x, IntRep* y, IntRep* r)
  715. {
  716.   nonnil(x);
  717.   nonnil(y);
  718.   int xl = x->len;
  719.   int yl = y->len;
  720.   int rl = xl + yl;
  721.   int xrsame = x == r;
  722.   int yrsame = y == r;
  723.   int xysame = x == y;
  724.   int rsgn = x->sgn == y->sgn;
  725.   
  726.   if (xl == 0 || yl == 0)
  727.     r = copy_zero(r);
  728.   else if (xl == 1 && x->s[0] == 1)
  729.     r = Icopy(r, y);
  730.   else if (yl == 1 && y->s[0] == 1)
  731.     r = Icopy(r, x);
  732.   else if (!xysame)
  733.   {
  734.     r = Iresize(r, rl);
  735.     unsigned short* rs = r->s;
  736.     unsigned short* topr = &(rs[rl]);
  737.  
  738.     // use best inner/outer loop params given constraints
  739.     unsigned short* currentr;
  740.     unsigned short* bota;
  741.     unsigned short* as;
  742.     unsigned short* botb;
  743.     unsigned short* topb;
  744.     if (xrsame)                 
  745.     { 
  746.       Iclear_from(r, xl);
  747.       currentr = &(rs[xl-1]);
  748.       bota = rs;
  749.       as = currentr;
  750.       botb = y->s;
  751.       topb = &(botb[yl]);
  752.     }
  753.     else if (yrsame)
  754.     {
  755.       Iclear_from(r, yl);
  756.       currentr = &(rs[yl-1]);
  757.       bota = rs;
  758.       as = currentr;
  759.       botb = x->s;
  760.       topb = &(botb[xl]);
  761.     }
  762.     else if (xl <= yl)
  763.     {
  764.       Iclear(r);
  765.       currentr = &(rs[xl-1]);
  766.       bota = x->s;
  767.       as = &(bota[xl-1]);
  768.       botb = y->s;
  769.       topb = &(botb[yl]);
  770.     }
  771.     else
  772.     {
  773.       Iclear(r);
  774.       currentr = &(rs[yl-1]);
  775.       bota = y->s;
  776.       as = &(bota[yl-1]);
  777.       botb = x->s;
  778.       topb = &(botb[xl]);
  779.     }
  780.  
  781.     while (as >= bota)
  782.     {
  783.       unsigned long ai = (unsigned long)(*as--);
  784.       unsigned short* rs = currentr--;
  785.       *rs = 0;
  786.       if (ai != 0)
  787.       {
  788.         unsigned long sum = 0;
  789.         unsigned short* bs = botb;
  790.         while (bs < topb)
  791.         {
  792.           sum += ai * (unsigned long)(*bs++) + (unsigned long)(*rs);
  793.           *rs++ = extract(sum);
  794.           sum = down(sum);
  795.         }
  796.         while (sum != 0 && rs < topr)
  797.         {
  798.           sum += (unsigned long)(*rs);
  799.           *rs++ = extract(sum);
  800.           sum = down(sum);
  801.         }
  802.       }
  803.     }
  804.   }
  805.   else                          // x & y same; compute over diagonals
  806.   {
  807.     r = Iresize(r, rl);
  808.     Iclear_last(r);
  809.     unsigned short* botr = r->s;
  810.     unsigned short* topr = &(botr[rl]);
  811.     unsigned short* rs =   &(botr[rl - 2]);
  812.  
  813.     unsigned short* bota = (xrsame)? botr : x->s;
  814.     unsigned short* loa =  &(bota[xl - 1]);
  815.     unsigned short* hia =  loa;
  816.  
  817.     for (; rs >= botr; --rs)
  818.     {
  819.       unsigned short* h = hia;
  820.       unsigned short* l = loa;
  821.       unsigned long prod = (unsigned long)(*h) * (unsigned long)(*l);
  822.       *rs = 0;
  823.  
  824.       for(;;)
  825.       {
  826.         unsigned short* rt = rs;
  827.         unsigned long sum = prod + (unsigned long)(*rt);
  828.         *rt++ = extract(sum);
  829.         sum = down(sum);
  830.         while (sum != 0 && rt < topr)
  831.         {
  832.           sum += (unsigned long)(*rt);
  833.           *rt++ = extract(sum);
  834.           sum = down(sum);
  835.         }
  836.         if (h > l)
  837.         {
  838.           rt = rs;
  839.           sum = prod + (unsigned long)(*rt);
  840.           *rt++ = extract(sum);
  841.           sum = down(sum);
  842.           while (sum != 0 && rt < topr)
  843.           {
  844.             sum += (unsigned long)(*rt);
  845.             *rt++ = extract(sum);
  846.             sum = down(sum);
  847.           }
  848.           if (--h >= ++l)
  849.             prod = (unsigned long)(*h) * (unsigned long)(*l);
  850.           else
  851.             break;
  852.         }
  853.         else
  854.           break;
  855.       }
  856.       if (loa > bota)
  857.         --loa;
  858.       else
  859.         --hia;
  860.     }
  861.   }
  862.   r->sgn = rsgn;
  863.   Icheck(r);
  864.   return r;
  865. }
  866.  
  867.  
  868. IntRep* multiply(IntRep* x, long y, IntRep* r)
  869. {
  870.   nonnil(x);
  871.   int xl = x->len;
  872.     
  873.   if (xl == 0 || y == 0)
  874.     r = copy_zero(r);
  875.   else if (y == 1)
  876.     r = Icopy(r, x);
  877.   else
  878.   {
  879.     int ysgn = y >= 0;
  880.     int rsgn = x->sgn == ysgn;
  881.     unsigned long uy = (ysgn)? y : -y;
  882.     unsigned short tmp[SHORT_PER_LONG];
  883.     int yl = 0;
  884.     while (uy != 0)
  885.     {
  886.       tmp[yl++] = extract(uy);
  887.       uy = down(uy);
  888.     }
  889.  
  890.     int xrsame = x == r;
  891.     int rl = xl + yl;
  892.     r = Iresize(r, rl);
  893.  
  894.     unsigned short* rs = r->s;
  895.     unsigned short* topr = &(rs[rl]);
  896.     unsigned short* currentr;
  897.     unsigned short* bota;
  898.     unsigned short* as;
  899.     unsigned short* botb;
  900.     unsigned short* topb;
  901.  
  902.     if (xrsame)
  903.     { 
  904.       Iclear_from(r, xl);
  905.       currentr = &(rs[xl-1]);
  906.       bota = rs;
  907.       as = currentr;
  908.       botb = tmp;
  909.       topb = &(botb[yl]);
  910.     }
  911.     else if (xl <= yl)
  912.     {
  913.       Iclear(r);
  914.       currentr = &(rs[xl-1]);
  915.       bota = x->s;
  916.       as = &(bota[xl-1]);
  917.       botb = tmp;
  918.       topb = &(botb[yl]);
  919.     }
  920.     else
  921.     {
  922.       Iclear(r);
  923.       currentr = &(rs[yl-1]);
  924.       bota = tmp;
  925.       as = &(bota[yl-1]);
  926.       botb = x->s;
  927.       topb = &(botb[xl]);
  928.     }
  929.  
  930.     while (as >= bota)
  931.     {
  932.       unsigned long ai = (unsigned long)(*as--);
  933.       unsigned short* rs = currentr--;
  934.       *rs = 0;
  935.       if (ai != 0)
  936.       {
  937.         unsigned long sum = 0;
  938.         unsigned short* bs = botb;
  939.         while (bs < topb)
  940.         {
  941.           sum += ai * (unsigned long)(*bs++) + (unsigned long)(*rs);
  942.           *rs++ = extract(sum);
  943.           sum = down(sum);
  944.         }
  945.         while (sum != 0 && rs < topr)
  946.         {
  947.           sum += (unsigned long)(*rs);
  948.           *rs++ = extract(sum);
  949.           sum = down(sum);
  950.         }
  951.       }
  952.     }
  953.     r->sgn = rsgn;
  954.   }
  955.   Icheck(r);
  956.   return r;
  957. }
  958.  
  959.  
  960. // main division routine
  961.  
  962. static void do_divide(unsigned short* rs,
  963.                unsigned short* ys, int yl,
  964.                unsigned short* qs, int ql)
  965. {
  966.   unsigned short* topy = &(ys[yl]);
  967.   unsigned short d1 = ys[yl - 1];
  968.   unsigned short d2 = ys[yl - 2];
  969.  
  970.   int l = ql - 1;
  971.   int i = l + yl;
  972.   
  973.   for (; l >= 0; --l, --i)
  974.   {
  975.     unsigned short qhat;       // guess q
  976.     if (d1 == rs[i])
  977.       qhat = I_MAXNUM;
  978.     else
  979.     {
  980.       unsigned long lr = up((unsigned long)rs[i]) | rs[i-1];
  981.       qhat = lr / d1;
  982.     }
  983.  
  984.     for(;;)     // adjust q, use docmp to avoid overflow problems
  985.     {
  986.       unsigned short ts[3];
  987.       unsigned long prod = (unsigned long)d2 * (unsigned long)qhat;
  988.       ts[0] = extract(prod);
  989.       prod = down(prod) + (unsigned long)d1 * (unsigned long)qhat;
  990.       ts[1] = extract(prod);
  991.       ts[2] = extract(down(prod));
  992.       if (docmp(ts, &(rs[i-2]), 3) > 0)
  993.         --qhat;
  994.       else
  995.         break;
  996.     };
  997.     
  998.     // multiply & subtract
  999.     
  1000.     unsigned short* yt = ys;
  1001.     unsigned short* rt = &(rs[l]);
  1002.     unsigned long prod = 0;
  1003.     unsigned long hi = 1;
  1004.     while (yt < topy)
  1005.     {
  1006.       prod = (unsigned long)qhat * (unsigned long)(*yt++) + down(prod);
  1007.       hi += (unsigned long)(*rt) + I_MAXNUM - (unsigned long)(extract(prod));
  1008.       *rt++ = extract(hi);
  1009.       hi = down(hi);
  1010.     }
  1011.     hi += (unsigned long)(*rt) + I_MAXNUM - (unsigned long)(down(prod));
  1012.     *rt = extract(hi);
  1013.     hi = down(hi);
  1014.     
  1015.     // off-by-one, add back
  1016.     
  1017.     if (hi == 0)
  1018.     {
  1019.       --qhat;
  1020.       yt = ys;
  1021.       rt = &(rs[l]);
  1022.       hi = 0;
  1023.       while (yt < topy)
  1024.       {
  1025.         hi = (unsigned long)(*rt) + (unsigned long)(*yt++) + down(hi);
  1026.         *rt++ = extract(hi);
  1027.       }
  1028.       *rt = 0;
  1029.     }
  1030.     if (qs != 0)
  1031.       qs[l] = qhat;
  1032.   }
  1033. }
  1034.  
  1035. // divide by single digit, return remainder
  1036. // if q != 0, then keep the result in q, else just compute rem
  1037.  
  1038. static int unscale(unsigned short* x, int xl, unsigned short y,
  1039.                    unsigned short* q)
  1040. {
  1041.   if (xl == 0 || y == 1)
  1042.     return 0;
  1043.   else if (q != 0)
  1044.   {
  1045.     unsigned short* botq = q;
  1046.     unsigned short* qs = &(botq[xl - 1]);
  1047.     unsigned short* xs = &(x[xl - 1]);
  1048.     unsigned long rem = 0;
  1049.     while (qs >= botq)
  1050.     {
  1051.       rem = up(rem) | *xs--;
  1052.       unsigned long u = rem / y;
  1053.       *qs-- = extract(u);
  1054.       rem -= u * y;
  1055.     }
  1056.     int r = extract(rem);
  1057.     return r;
  1058.   }
  1059.   else                          // same loop, a bit faster if just need rem
  1060.   {
  1061.     unsigned short* botx = x;
  1062.     unsigned short* xs = &(botx[xl - 1]);
  1063.     unsigned long rem = 0;
  1064.     while (xs >= botx)
  1065.     {
  1066.       rem = up(rem) | *xs--;
  1067.       unsigned long u = rem / y;
  1068.       rem -= u * y;
  1069.     }
  1070.     int r = extract(rem);
  1071.     return r;
  1072.   }
  1073. }
  1074.  
  1075.  
  1076. IntRep* div(IntRep* x, IntRep* y, IntRep* q)
  1077. {
  1078.   nonnil(x);
  1079.   nonnil(y);
  1080.   int xl = x->len;
  1081.   int yl = y->len;
  1082.   if (yl == 0) (*lib_error_handler)("Integer", "attempted division by zero");
  1083.  
  1084.   int comp = ucompare(x, y);
  1085.   int xsgn = x->sgn;
  1086.   int ysgn = y->sgn;
  1087.  
  1088.   int samesign = xsgn == ysgn;
  1089.  
  1090.   if (comp < 0)
  1091.     q = copy_zero(q);
  1092.   else if (comp == 0)
  1093.     q = copy_one(q, samesign);
  1094.   else if (yl == 1)
  1095.   {
  1096.     q = Icopy(q, x);
  1097.     unscale(q->s, q->len, y->s[0], q->s);
  1098.   }
  1099.   else
  1100.   {
  1101.     IntRep* yy = 0;
  1102.     IntRep* r  = 0;
  1103.     unsigned short prescale = (I_RADIX / (1 + y->s[yl - 1]));
  1104.     if (prescale != 1 || y == q)
  1105.     {
  1106.       yy = multiply(y, ((long)prescale & I_MAXNUM), yy);
  1107.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1108.     }
  1109.     else
  1110.     {
  1111.       yy = y;
  1112.       r = Icopy(r, x);
  1113.     }
  1114.  
  1115.     int ql = xl - yl + 1;
  1116.       
  1117.     q = Iresize(q, ql);
  1118.     Iclear(q);
  1119.     do_divide(r->s, yy->s, yl, q->s, ql);
  1120.  
  1121.     if (yy != y) delete yy;
  1122.     delete r;
  1123.   }
  1124.   q->sgn = samesign;
  1125.   Icheck(q);
  1126.   return q;
  1127. }
  1128.  
  1129. IntRep* div(IntRep* x, long y, IntRep* q)
  1130. {
  1131.   nonnil(x);
  1132.   int xl = x->len;
  1133.   if (y == 0) (*lib_error_handler)("Integer", "attempted division by zero");
  1134.  
  1135.   unsigned short ys[SHORT_PER_LONG];
  1136.   unsigned long u;
  1137.   int ysgn = y >= 0;
  1138.   if (ysgn)
  1139.     u = y;
  1140.   else
  1141.     u = -y;
  1142.   int yl = 0;
  1143.   while (u != 0)
  1144.   {
  1145.     ys[yl++] = extract(u);
  1146.     u = down(u);
  1147.   }
  1148.  
  1149.   int comp = xl - yl;
  1150.   if (comp == 0) comp = docmp(x->s, ys, xl);
  1151.  
  1152.   int xsgn = x->sgn;
  1153.   int samesign = xsgn == ysgn;
  1154.  
  1155.   if (comp < 0)
  1156.     q = copy_zero(q);
  1157.   else if (comp == 0)
  1158.   {
  1159.     q = copy_one(q, samesign);
  1160.   }
  1161.   else if (yl == 1)
  1162.   {
  1163.     q = Icopy(q, x);
  1164.     unscale(q->s, q->len, ys[0], q->s);
  1165.   }
  1166.   else
  1167.   {
  1168.     IntRep* r  = 0;
  1169.     unsigned short prescale = (I_RADIX / (1 + ys[yl - 1]));
  1170.     if (prescale != 1)
  1171.     {
  1172.       unsigned long prod = (unsigned long)prescale * (unsigned long)ys[0];
  1173.       ys[0] = extract(prod);
  1174.       prod = down(prod) + (unsigned long)prescale * (unsigned long)ys[1];
  1175.       ys[1] = extract(prod);
  1176.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1177.     }
  1178.     else
  1179.     {
  1180.       r = Icopy(r, x);
  1181.     }
  1182.  
  1183.     int ql = xl - yl + 1;
  1184.       
  1185.     q = Iresize(q, ql);
  1186.     Iclear(q);
  1187.     do_divide(r->s, ys, yl, q->s, ql);
  1188.  
  1189.     delete r;
  1190.   }
  1191.   q->sgn = samesign;
  1192.   Icheck(q);
  1193.   return q;
  1194. }
  1195.  
  1196.  
  1197. void divide(Integer& Ix, long y, Integer& Iq, long& rem)
  1198. {
  1199.   IntRep* x = Ix.rep;
  1200.   nonnil(x);
  1201.   IntRep* q = Iq.rep;
  1202.   int xl = x->len;
  1203.   if (y == 0) (*lib_error_handler)("Integer", "attempted division by zero");
  1204.  
  1205.   unsigned short ys[SHORT_PER_LONG];
  1206.   unsigned long u;
  1207.   int ysgn = y >= 0;
  1208.   if (ysgn)
  1209.     u = y;
  1210.   else
  1211.     u = -y;
  1212.   int yl = 0;
  1213.   while (u != 0)
  1214.   {
  1215.     ys[yl++] = extract(u);
  1216.     u = down(u);
  1217.   }
  1218.  
  1219.   int comp = xl - yl;
  1220.   if (comp == 0) comp = docmp(x->s, ys, xl);
  1221.  
  1222.   int xsgn = x->sgn;
  1223.   int samesign = xsgn == ysgn;
  1224.  
  1225.   if (comp < 0)
  1226.   {
  1227.     rem = Itolong(x);
  1228.     q = copy_zero(q);
  1229.   }
  1230.   else if (comp == 0)
  1231.   {
  1232.     q = copy_one(q, samesign);
  1233.     rem = 0;
  1234.   }
  1235.   else if (yl == 1)
  1236.   {
  1237.     q = Icopy(q, x);
  1238.     rem = unscale(q->s, q->len, ys[0], q->s);
  1239.   }
  1240.   else
  1241.   {
  1242.     IntRep* r  = 0;
  1243.     unsigned short prescale = (I_RADIX / (1 + ys[yl - 1]));
  1244.     if (prescale != 1)
  1245.     {
  1246.       unsigned long prod = (unsigned long)prescale * (unsigned long)ys[0];
  1247.       ys[0] = extract(prod);
  1248.       prod = down(prod) + (unsigned long)prescale * (unsigned long)ys[1];
  1249.       ys[1] = extract(prod);
  1250.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1251.     }
  1252.     else
  1253.     {
  1254.       r = Icopy(r, x);
  1255.     }
  1256.  
  1257.     int ql = xl - yl + 1;
  1258.       
  1259.     q = Iresize(q, ql);
  1260.     Iclear(q);
  1261.     
  1262.     do_divide(r->s, ys, yl, q->s, ql);
  1263.  
  1264.     if (prescale != 1)
  1265.     {
  1266.       Icheck(r);
  1267.       unscale(r->s, r->len, prescale, r->s);
  1268.     }
  1269.     Icheck(r);
  1270.     rem = Itolong(r);
  1271.     delete r;
  1272.   }
  1273.   rem = abs(rem);
  1274.   if (xsgn == I_NEGATIVE) rem = -rem;
  1275.   q->sgn = samesign;
  1276.   Icheck(q);
  1277.   Iq.rep = q;
  1278. }
  1279.  
  1280.  
  1281. void divide(Integer& Ix, Integer& Iy, Integer& Iq, Integer& Ir)
  1282. {
  1283.   IntRep* x = Ix.rep;
  1284.   nonnil(x);
  1285.   IntRep* y = Iy.rep;
  1286.   nonnil(y);
  1287.   IntRep* q = Iq.rep;
  1288.   IntRep* r = Ir.rep;
  1289.  
  1290.   int xl = x->len;
  1291.   int yl = y->len;
  1292.   if (yl == 0)
  1293.     (*lib_error_handler)("Integer", "attempted division by zero");
  1294.  
  1295.   int comp = ucompare(x, y);
  1296.   int xsgn = x->sgn;
  1297.   int ysgn = y->sgn;
  1298.  
  1299.   int samesign = xsgn == ysgn;
  1300.  
  1301.   if (comp < 0)
  1302.   {
  1303.     q = copy_zero(q);
  1304.     r = Icopy(r, x);
  1305.   }
  1306.   else if (comp == 0)
  1307.   {
  1308.     q = copy_one(q, samesign);
  1309.     r = copy_zero(r);
  1310.   }
  1311.   else if (yl == 1)
  1312.   {
  1313.     q = Icopy(q, x);
  1314.     int rem =  unscale(q->s, q->len, y->s[0], q->s);
  1315.     r = Icopy_long(r, rem);
  1316.     if (rem != 0)
  1317.       r->sgn = xsgn;
  1318.   }
  1319.   else
  1320.   {
  1321.     IntRep* yy = 0;
  1322.     unsigned short prescale = (I_RADIX / (1 + y->s[yl - 1]));
  1323.     if (prescale != 1 || y == q || y == r)
  1324.     {
  1325.       yy = multiply(y, ((long)prescale & I_MAXNUM), yy);
  1326.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1327.     }
  1328.     else
  1329.     {
  1330.       yy = y;
  1331.       r = Icopy(r, x);
  1332.     }
  1333.  
  1334.     int ql = xl - yl + 1;
  1335.       
  1336.     q = Iresize(q, ql);
  1337.     Iclear(q);
  1338.     do_divide(r->s, yy->s, yl, q->s, ql);
  1339.  
  1340.     if (yy != y) delete yy;
  1341.     if (prescale != 1)
  1342.     {
  1343.       Icheck(r);
  1344.       unscale(r->s, r->len, prescale, r->s);
  1345.     }
  1346.   }
  1347.   q->sgn = samesign;
  1348.   Icheck(q);
  1349.   Iq.rep = q;
  1350.   Icheck(r);
  1351.   Ir.rep = r;
  1352. }
  1353.  
  1354. IntRep* mod(IntRep* x, IntRep* y, IntRep* r)
  1355. {
  1356.   nonnil(x);
  1357.   nonnil(y);
  1358.   int xl = x->len;
  1359.   int yl = y->len;
  1360.   if (yl == 0) (*lib_error_handler)("Integer", "attempted division by zero");
  1361.  
  1362.   int comp = ucompare(x, y);
  1363.   int xsgn = x->sgn;
  1364.  
  1365.   if (comp < 0)
  1366.     r = Icopy(r, x);
  1367.   else if (comp == 0)
  1368.     r = copy_zero(r);
  1369.   else if (yl == 1)
  1370.   {
  1371.     int rem = unscale(x->s, xl, y->s[0], 0);
  1372.     r = Icopy_long(r, rem);
  1373.     if (rem != 0)
  1374.       r->sgn = xsgn;
  1375.   }
  1376.   else
  1377.   {
  1378.     IntRep* yy = 0;
  1379.     unsigned short prescale = (I_RADIX / (1 + y->s[yl - 1]));
  1380.     if (prescale != 1 || y == r)
  1381.     {
  1382.       yy = multiply(y, ((long)prescale & I_MAXNUM), yy);
  1383.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1384.     }
  1385.     else
  1386.     {
  1387.       yy = y;
  1388.       r = Icopy(r, x);
  1389.     }
  1390.       
  1391.     do_divide(r->s, yy->s, yl, 0, xl - yl + 1);
  1392.  
  1393.     if (yy != y) delete yy;
  1394.  
  1395.     if (prescale != 1)
  1396.     {
  1397.       Icheck(r);
  1398.       unscale(r->s, r->len, prescale, r->s);
  1399.     }
  1400.   }
  1401.   Icheck(r);
  1402.   return r;
  1403. }
  1404.  
  1405. IntRep* mod(IntRep* x, long y, IntRep* r)
  1406. {
  1407.   nonnil(x);
  1408.   int xl = x->len;
  1409.   if (y == 0) (*lib_error_handler)("Integer", "attempted division by zero");
  1410.  
  1411.   unsigned short ys[SHORT_PER_LONG];
  1412.   unsigned long u;
  1413.   int ysgn = y >= 0;
  1414.   if (ysgn)
  1415.     u = y;
  1416.   else
  1417.     u = -y;
  1418.   int yl = 0;
  1419.   while (u != 0)
  1420.   {
  1421.     ys[yl++] = extract(u);
  1422.     u = down(u);
  1423.   }
  1424.  
  1425.   int comp = xl - yl;
  1426.   if (comp == 0) comp = docmp(x->s, ys, xl);
  1427.  
  1428.   int xsgn = x->sgn;
  1429.  
  1430.   if (comp < 0)
  1431.     r = Icopy(r, x);
  1432.   else if (comp == 0)
  1433.     r = copy_zero(r);
  1434.   else if (yl == 1)
  1435.   {
  1436.     int rem = unscale(x->s, xl, ys[0], 0);
  1437.     r = Icopy_long(r, rem);
  1438.     if (rem != 0)
  1439.       r->sgn = xsgn;
  1440.   }
  1441.   else
  1442.   {
  1443.     unsigned short prescale = (I_RADIX / (1 + ys[yl - 1]));
  1444.     if (prescale != 1)
  1445.     {
  1446.       unsigned long prod = (unsigned long)prescale * (unsigned long)ys[0];
  1447.       ys[0] = extract(prod);
  1448.       prod = down(prod) + (unsigned long)prescale * (unsigned long)ys[1];
  1449.       ys[1] = extract(prod);
  1450.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1451.     }
  1452.     else
  1453.     {
  1454.       r = Icopy(r, x);
  1455.     }
  1456.       
  1457.     do_divide(r->s, ys, yl, 0, xl - yl + 1);
  1458.  
  1459.     if (prescale != 1)
  1460.     {
  1461.       Icheck(r);
  1462.       unscale(r->s, r->len, prescale, r->s);
  1463.     }
  1464.   }
  1465.   Icheck(r);
  1466.   return r;
  1467. }
  1468.  
  1469. IntRep* lshift(IntRep* x, long y, IntRep* r)
  1470. {
  1471.   nonnil(x);
  1472.   int xl = x->len;
  1473.   if (xl == 0 || y == 0)
  1474.   {
  1475.     r = Icopy(r, x);
  1476.     return r;
  1477.   }
  1478.  
  1479.   int xrsame = x == r;
  1480.   int rsgn = x->sgn;
  1481.  
  1482.   long ay = (y < 0)? -y : y;
  1483.   int bw = ay / I_SHIFT;
  1484.   int sw = ay % I_SHIFT;
  1485.  
  1486.   if (y > 0)
  1487.   {
  1488.     int rl = bw + xl + 1;
  1489.     r = Iresize(r, rl);
  1490.    
  1491.     unsigned short* botr = r->s;
  1492.     unsigned short* rs = &(botr[rl - 1]);
  1493.     unsigned short* botx = (xrsame)? botr : x->s;
  1494.     unsigned short* xs = &(botx[xl - 1]);
  1495.     unsigned long a = 0;
  1496.     while (xs >= botx)
  1497.     {
  1498.       a = up(a) | ((unsigned long)(*xs--) << sw);
  1499.       *rs-- = extract(down(a));
  1500.     }
  1501.     *rs-- = extract(a);
  1502.     while (rs >= botr)
  1503.       *rs-- = 0;
  1504.   }
  1505.   else
  1506.   {
  1507.     int rl = xl - bw;
  1508.     if (rl < 0)
  1509.       r = copy_zero(r);
  1510.     else
  1511.     {
  1512.       r = Iresize(r, rl);
  1513.       int rw = I_SHIFT - sw;
  1514.       unsigned short* rs = r->s;
  1515.       unsigned short* topr = &(rs[rl]);
  1516.       unsigned short* botx = (xrsame)? rs : x->s;
  1517.       unsigned short* xs =  &(botx[bw]);
  1518.       unsigned short* topx = &(botx[xl]);
  1519.       unsigned long a = (unsigned long)(*xs++) >> sw;
  1520.       while (xs < topx)
  1521.       {
  1522.         a |= (unsigned long)(*xs++) << rw;
  1523.         *rs++ = extract(a);
  1524.         a = down(a);
  1525.       }
  1526.       *rs++ = extract(a);
  1527.       if (xrsame) topr = topx;
  1528.       while (rs < topr)
  1529.         *rs++ = 0;
  1530.     }
  1531.   }
  1532.   r->sgn = rsgn;
  1533.   Icheck(r);
  1534.   return r;
  1535. }
  1536.  
  1537. IntRep* lshift(IntRep* x, IntRep* yy, int negatey, IntRep* r)
  1538. {
  1539.   long y = Itolong(yy);
  1540.   if (negatey)
  1541.     y = -y;
  1542.  
  1543.   return lshift(x, y, r);
  1544. }
  1545.  
  1546. IntRep* bitop(IntRep* x, IntRep* y, IntRep* r, char op)
  1547. {
  1548.   nonnil(x);
  1549.   nonnil(y);
  1550.   int xl = x->len;
  1551.   int yl = y->len;
  1552.   int xrsame = x == r;
  1553.   int yrsame = y == r;
  1554.   r = Iresize(r, calc_len(xl, yl, 0));
  1555.   r->sgn = x->sgn;
  1556.   unsigned short* rs = r->s;
  1557.   unsigned short* topr = &(rs[r->len]);
  1558.   unsigned short* as;
  1559.   unsigned short* bs;
  1560.   unsigned short* topb;
  1561.   if (xl >= yl)
  1562.   {
  1563.     as = (xrsame)? rs : x->s;
  1564.     bs = (yrsame)? rs : y->s;
  1565.     topb = &(bs[yl]);
  1566.   }
  1567.   else
  1568.   {
  1569.     bs = (xrsame)? rs : x->s;
  1570.     topb = &(bs[xl]);
  1571.     as = (yrsame)? rs : y->s;
  1572.   }
  1573.  
  1574.   switch (op)
  1575.   {
  1576.   case '&':
  1577.     while (bs < topb) *rs++ = *as++ & *bs++;
  1578.     while (rs < topr) *rs++ = 0;
  1579.     break;
  1580.   case '|':
  1581.     while (bs < topb) *rs++ = *as++ | *bs++;
  1582.     while (rs < topr) *rs++ = *as++;
  1583.     break;
  1584.   case '^':
  1585.     while (bs < topb) *rs++ = *as++ ^ *bs++;
  1586.     while (rs < topr) *rs++ = *as++;
  1587.     break;
  1588.   }
  1589.   Icheck(r);
  1590.   return r;
  1591. }
  1592.  
  1593. IntRep* bitop(IntRep* x, long y, IntRep* r, char op)
  1594. {
  1595.   nonnil(x);
  1596.   unsigned short tmp[SHORT_PER_LONG];
  1597.   unsigned long u;
  1598.   int newsgn;
  1599.   if (newsgn = (y >= 0))
  1600.     u = y;
  1601.   else
  1602.     u = -y;
  1603.   
  1604.   int l = 0;
  1605.   while (u != 0)
  1606.   {
  1607.     tmp[l++] = extract(u);
  1608.     u = down(u);
  1609.   }
  1610.  
  1611.   int xl = x->len;
  1612.   int yl = l;
  1613.   int xrsame = x == r;
  1614.   r = Iresize(r, calc_len(xl, yl, 0));
  1615.   r->sgn = x->sgn;
  1616.   unsigned short* rs = r->s;
  1617.   unsigned short* topr = &(rs[r->len]);
  1618.   unsigned short* as;
  1619.   unsigned short* bs;
  1620.   unsigned short* topb;
  1621.   if (xl >= yl)
  1622.   {
  1623.     as = (xrsame)? rs : x->s;
  1624.     bs = tmp;
  1625.     topb = &(bs[yl]);
  1626.   }
  1627.   else
  1628.   {
  1629.     bs = (xrsame)? rs : x->s;
  1630.     topb = &(bs[xl]);
  1631.     as = tmp;
  1632.   }
  1633.  
  1634.   switch (op)
  1635.   {
  1636.   case '&':
  1637.     while (bs < topb) *rs++ = *as++ & *bs++;
  1638.     while (rs < topr) *rs++ = 0;
  1639.     break;
  1640.   case '|':
  1641.     while (bs < topb) *rs++ = *as++ | *bs++;
  1642.     while (rs < topr) *rs++ = *as++;
  1643.     break;
  1644.   case '^':
  1645.     while (bs < topb) *rs++ = *as++ ^ *bs++;
  1646.     while (rs < topr) *rs++ = *as++;
  1647.     break;
  1648.   }
  1649.   Icheck(r);
  1650.   return r;
  1651. }
  1652.  
  1653.  
  1654.  
  1655. IntRep*  compl(IntRep* src, IntRep* r)
  1656. {
  1657.   nonnil(src);
  1658.   r = Icopy(r, src);
  1659.   unsigned short* s = r->s;
  1660.   unsigned short* top = &(s[r->len - 1]);
  1661.   while (s < top)
  1662.   {
  1663.     unsigned short cmp = ~(*s);
  1664.     *s++ = cmp;
  1665.   }
  1666.   unsigned short a = *s;
  1667.   unsigned short b = 0;
  1668.   while (a != 0)
  1669.   {
  1670.     b <<= 1;
  1671.     if (!(a & 1)) b |= 1;
  1672.     a >>= 1;
  1673.   }
  1674.   *s = b;
  1675.   Icheck(r);
  1676.   return r;
  1677. }
  1678.  
  1679. void setbit(Integer& x, long b)
  1680. {
  1681.   if (b >= 0)
  1682.   {
  1683.     int bw = (unsigned long)b / I_SHIFT;
  1684.     int sw = (unsigned long)b % I_SHIFT;
  1685.     if (x.rep == 0)
  1686.       x.rep = Ialloc(0, 0, 0, I_POSITIVE, bw + 1);
  1687.     else if (x.rep->len < bw)
  1688.     {
  1689.       int xl = x.rep->len;
  1690.       x.rep = Iresize(x.rep, calc_len(xl, bw+1, 0));
  1691.       Iclear_from(x.rep, xl);
  1692.     }
  1693.     x.rep->s[bw] |= (1 << sw);
  1694.     Icheck(x.rep);
  1695.   }
  1696. }
  1697.  
  1698. void clearbit(Integer& x, long b)
  1699. {
  1700.   if (b >= 0)
  1701.   {
  1702.     int bw = (unsigned long)b / I_SHIFT;
  1703.     int sw = (unsigned long)b % I_SHIFT;
  1704.     if (x.rep == 0)
  1705.       x.rep = Ialloc(0, 0, 0, I_POSITIVE, bw + 1);
  1706.     else if (x.rep->len < bw)
  1707.     {
  1708.       int xl = x.rep->len;
  1709.       x.rep = Iresize(x.rep, calc_len(xl, bw+1, 0));
  1710.       Iclear_from(x.rep, xl);
  1711.     }
  1712.     x.rep->s[bw] &= ~(1 << sw);
  1713.     Icheck(x.rep);
  1714.   }
  1715. }
  1716.  
  1717. int testbit(Integer& x, long b)
  1718. {
  1719.   if (x.rep != 0 && b >= 0)
  1720.   {
  1721.     int bw = (unsigned long)b / I_SHIFT;
  1722.     int sw = (unsigned long)b % I_SHIFT;
  1723.     return (bw < x.rep->len && (x.rep->s[bw] & (1 << sw)) != 0);
  1724.   }
  1725.   else
  1726.     return 0;
  1727. }
  1728.  
  1729.  
  1730. // A  version of knuth's algorithm B / ex. 4.5.3.34
  1731. // A better version that doesn't bother shifting all of `t' forthcoming
  1732.  
  1733. IntRep* gcd(IntRep* x, IntRep* y)
  1734. {
  1735.   nonnil(x);
  1736.   nonnil(y);
  1737.   int ul = x->len;
  1738.   int vl = y->len;
  1739.   
  1740.   if (vl == 0)
  1741.     return Ialloc(0, x->s, ul, I_POSITIVE, ul);
  1742.   else if (ul == 0)
  1743.     return Ialloc(0, y->s, vl, I_POSITIVE, vl);
  1744.  
  1745.   IntRep* u = Ialloc(0, x->s, ul, I_POSITIVE, ul);
  1746.   IntRep* v = Ialloc(0, y->s, vl, I_POSITIVE, vl);
  1747.  
  1748. // find shift so that both not even
  1749.  
  1750.   long k = 0;
  1751.   int l = ul >? vl;
  1752.   int cont = 1;
  1753.   for (int i = 0;  i < l && cont; ++i)
  1754.   {
  1755.     unsigned long a =  (i < ul)? u->s[i] : 0;
  1756.     unsigned long b =  (i < vl)? v->s[i] : 0;
  1757.     for (int j = 0; j < I_SHIFT; ++j)
  1758.     {
  1759.       if ((a | b) & 1)
  1760.       {
  1761.         cont = 0;
  1762.         break;
  1763.       }
  1764.       else
  1765.       {
  1766.         ++k;
  1767.         a >>= 1;
  1768.         b >>= 1;
  1769.       }
  1770.     }
  1771.   }
  1772.  
  1773.   if (k != 0)
  1774.   {
  1775.     u = lshift(u, -k, u);
  1776.     v = lshift(v, -k, v);
  1777.   }
  1778.  
  1779.   IntRep* t;
  1780.   if (u->s[0] & 01)
  1781.     t = Ialloc(0, v->s, v->len, !v->sgn, v->len);
  1782.   else
  1783.     t = Ialloc(0, u->s, u->len, u->sgn, u->len);
  1784.  
  1785.   while (t->len != 0)
  1786.   {
  1787.     long s = 0;                 // shift t until odd
  1788.     cont = 1;
  1789.     int tl = t->len;
  1790.     for (i = 0; i < tl && cont; ++i)
  1791.     {
  1792.       unsigned long a = t->s[i];
  1793.       for (int j = 0; j < I_SHIFT; ++j)
  1794.       {
  1795.         if (a & 1)
  1796.         {
  1797.           cont = 0;
  1798.           break;
  1799.         }
  1800.         else
  1801.         {
  1802.           ++s;
  1803.           a >>= 1;
  1804.         }
  1805.       }
  1806.     }
  1807.  
  1808.     if (s != 0) t = lshift(t, -s, t);
  1809.  
  1810.     if (t->sgn == I_POSITIVE)
  1811.     {
  1812.       u = Icopy(u, t);
  1813.       t = add(t, 0, v, 1, t);
  1814.     }
  1815.     else
  1816.     {
  1817.       v = Ialloc(v, t->s, t->len, !t->sgn, t->len);
  1818.       t = add(t, 0, u, 0, t);
  1819.     }
  1820.   }
  1821.   delete t;
  1822.   delete v;
  1823.   if (k != 0) u = lshift(u, k, u);
  1824.   return u;
  1825. }
  1826.  
  1827. IntTmp lcm(Integer& x, Integer& y)
  1828. {
  1829.   nonnil(x.rep);
  1830.   nonnil(y.rep);
  1831.   Integer g;
  1832.   if (sign(x) == 0 || sign(y) == 0)
  1833.     g = 1;
  1834.   else 
  1835.     g = gcd(x, y);
  1836.   return x / g * y;
  1837. }
  1838.  
  1839.  
  1840. long lg(IntRep* x)
  1841. {
  1842.   nonnil(x);
  1843.   int xl = x->len;
  1844.   if (xl == 0)
  1845.     return 0;
  1846.  
  1847.   long l = (xl - 1) * I_SHIFT - 1;
  1848.   unsigned short a = x->s[xl-1];
  1849.  
  1850.   while (a != 0)
  1851.   {
  1852.     a = a >> 1;
  1853.     ++l;
  1854.   }
  1855.   return l;
  1856. }
  1857.   
  1858. IntRep* power(IntRep* x, long y, IntRep* r)
  1859. {
  1860.   nonnil(x);
  1861.   int sgn;
  1862.   if (x->sgn == I_POSITIVE || (!(y & 1)))
  1863.     sgn = I_POSITIVE;
  1864.   else
  1865.     sgn = I_NEGATIVE;
  1866.  
  1867.   int xl = x->len;
  1868.  
  1869.   if (y == 0 || (xl == 1 && x->s[0] == 1))
  1870.     r = copy_one(r, sgn);
  1871.   else if (xl == 0 || y < 0)
  1872.     r = copy_zero(r);
  1873.   else if (y == 1 || y == -1)
  1874.     r = Icopy(r, x);
  1875.   else
  1876.   {
  1877.     int maxsize = ((lg(x) + 1) * y) / I_SHIFT + 2;     // pre-allocate space
  1878.     IntRep* b = Ialloc(0, x->s, xl, I_POSITIVE, maxsize);
  1879.     b->len = xl;
  1880.     r = Ialloc(r, 0, 0, I_POSITIVE, maxsize);
  1881.     r = copy_one(r, I_POSITIVE);
  1882.     for(;;)
  1883.     {
  1884.       if (y & 1)
  1885.         r = multiply(r, b, r);
  1886.       if ((y >>= 1) == 0)
  1887.         break;
  1888.       else
  1889.         b = multiply(b, b, b);
  1890.     }
  1891.     delete b;
  1892.   }
  1893.   r->sgn = sgn;
  1894.   Icheck(r);
  1895.   return r;
  1896. }
  1897.  
  1898.  
  1899.  
  1900. Integer sqrt(Integer& x)
  1901. {
  1902.   nonnil(x.rep);
  1903.   Integer r = x;
  1904.   int s = sign(x);
  1905.   if (s < 0) x.error("Attempted square root of negative Integer");
  1906.   if (s != 0)
  1907.   {
  1908.     r >>= (lg(x) / 2); // get close
  1909.     Integer q = x / r;
  1910.     while (q < r)
  1911.     {
  1912.       r += q;
  1913.       r >>= 1;
  1914.       q.rep = div(x.rep, r.rep, q.rep);
  1915.     }
  1916.   }
  1917.   return r;
  1918. }
  1919.  
  1920.  
  1921.  
  1922. extern Obstack _libgxx_io_ob;
  1923. extern char* _libgxx_io_oblast;
  1924.  
  1925. IntRep* atoIntRep(const char* s, int base = 10)
  1926. {
  1927.   int sl = strlen(s);
  1928.   IntRep* r = Ialloc(0, 0, 0, 0, sl * (lg(base) + 1) / I_SHIFT + 1);
  1929.   if (s != 0)
  1930.   {
  1931.     char sgn;
  1932.     while (isspace(*s)) ++s;
  1933.     if (*s == '-')
  1934.     {
  1935.       sgn = I_NEGATIVE;
  1936.       s++;
  1937.     }
  1938.     else if (*s == '+')
  1939.     {
  1940.       sgn = I_POSITIVE;
  1941.       s++;
  1942.     }
  1943.     else
  1944.       sgn = I_POSITIVE;
  1945.     for (;;)
  1946.     {
  1947.       long digit;
  1948.       if (*s >= '0' && *s <= '9') digit = *s - '0';
  1949.       else if (*s >= 'a' && *s <= 'z') digit = *s - 'a' + 10;
  1950.       else if (*s >= 'A' && *s <= 'Z') digit = *s - 'A' + 10;
  1951.       else break;
  1952.       if (digit >= base) break;
  1953.       r = multiply(r, base, r);
  1954.       r = add(r, 0, digit, r);
  1955.       r->sgn = sgn;
  1956.       ++s;
  1957.     }
  1958.   }
  1959.   return r;
  1960. }
  1961.  
  1962.  
  1963. char* Itoa(IntRep* x, int base = 10, int width = 0)
  1964. {
  1965.   if (_libgxx_io_oblast)  _libgxx_io_ob.free(_libgxx_io_oblast);
  1966.  
  1967.   int wrksiz = (x->len + 1) * I_SHIFT / lg(base) + 2 + width;
  1968.   _libgxx_io_ob.blank(wrksiz);
  1969.   _libgxx_io_oblast = (char*)(_libgxx_io_ob.finish());
  1970.   char* obbase = _libgxx_io_oblast;
  1971.   char* e = obbase + wrksiz - 1;
  1972.   char* s = e;
  1973.   *--s = 0;
  1974.  
  1975.   if (x->len == 0)
  1976.     *--s = '0';
  1977.   else
  1978.   {
  1979.     IntRep* z = Icopy(0, x);
  1980.  
  1981.     // split division by base into two parts: 
  1982.     // first divide by biggest power of base that fits in an unsigned short,
  1983.     // then use straight signed div/mods from there. 
  1984.  
  1985.     // find power
  1986.     int bpower = 1;
  1987.     unsigned short b = base;
  1988.     unsigned short maxb = I_MAXNUM / base;
  1989.     while (b < maxb)
  1990.     {
  1991.       b *= base;
  1992.       ++bpower;
  1993.     }
  1994.     for(;;)
  1995.     {
  1996.       int rem = unscale(z->s, z->len, b, z->s);
  1997.       Icheck(z);
  1998.       if (z->len == 0)
  1999.       {
  2000.         while (rem != 0)
  2001.         {
  2002.           char ch = rem % base;
  2003.           rem /= base;
  2004.           if (ch >= 10)
  2005.             ch += 'a' - 10;
  2006.           else
  2007.             ch += '0';
  2008.           *--s = ch;
  2009.         }
  2010.         delete z;
  2011.         break;
  2012.       }
  2013.       else
  2014.       {
  2015.         for (int i = 0; i < bpower; ++i)
  2016.         {
  2017.           char ch = rem % base;
  2018.           rem /= base;
  2019.           if (ch >= 10)
  2020.             ch += 'a' - 10;
  2021.           else
  2022.             ch += '0';
  2023.           *--s = ch;
  2024.         }
  2025.       }
  2026.     }
  2027.   }
  2028.  
  2029.   if (x->sgn == I_NEGATIVE) *--s = '-';
  2030.   int w = e - s - 1;
  2031.   while (w++ < width)
  2032.     *--s = ' ';
  2033.   return s;
  2034. }
  2035.  
  2036. char* dec(Integer& x, int width = 0)
  2037. {
  2038.   return Itoa(x, 10, width);
  2039. }
  2040.  
  2041. char* oct(Integer& x, int width = 0)
  2042. {
  2043.   return Itoa(x, 8, width);
  2044. }
  2045.  
  2046. char* hex(Integer& x, int width = 0)
  2047. {
  2048.   return Itoa(x, 16, width);
  2049. }
  2050.  
  2051. istream& operator >> (istream& s, Integer& y)
  2052. {
  2053.   if (!s.readable())
  2054.   {
  2055.     s.clear(_bad);
  2056.     return s;
  2057.   }
  2058.  
  2059.   int got_one = 0;
  2060.   char sgn = 0;
  2061.   char ch;
  2062.   y.rep = copy_zero(y.rep);
  2063.   s >> WS;
  2064.   while (s.good())
  2065.   {
  2066.     s.get(ch);
  2067.     if (ch == '-')
  2068.     {
  2069.       if (sgn == 0)
  2070.         sgn = '-';
  2071.       else
  2072.         break;
  2073.     }
  2074.     else if (ch >= '0' && ch <= '9')
  2075.     {
  2076.       long digit = ch - '0';
  2077.       y *= 10;
  2078.       y += digit;
  2079.       got_one = 1;
  2080.     }
  2081.     else
  2082.       break;
  2083.   }
  2084.   if (s.good())
  2085.     s.unget(ch);
  2086.   if (!got_one)
  2087.     s.clear(_bad);
  2088.  
  2089.   if (sgn == '-')
  2090.     y.negate();
  2091.  
  2092.   return s;
  2093. }
  2094.  
  2095. int Integer::OK()
  2096. {
  2097.   int v = rep != 0;             // have a rep
  2098.   int l = rep->len;
  2099.   int s = rep->sgn;
  2100.   v &= l <= rep->sz;            // length with bounds
  2101.   v &= s == 0 || s == 1;        // legal sign
  2102.   Icheck(rep);                  // and correctly adjusted
  2103.   v &= rep->len == l;
  2104.   v &= rep->sgn == s;
  2105.   if (!v) error("invariant failure");
  2106.   return v;
  2107. }
  2108.  
  2109. void Integer::error(char* msg)
  2110. {
  2111.   (*lib_error_handler)("Integer", msg);
  2112. }
  2113.