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