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