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