/* 
Copyright (C) 1988 Free Software Foundation
    written by Doug Lea (dl@rocky.oswego.edu)

This file is part of GNU CC.

GNU CC is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY.  No author or distributor
accepts responsibility to anyone for the consequences of using it
or for whether it serves any particular purpose or works at all,
unless he says so in writing.  Refer to the GNU CC General Public
License for full details.

Everyone is granted permission to copy, modify and redistribute
GNU CC, but only under the conditions described in the
GNU CC General Public License.   A copy of this license is
supposed to have been given to you along with GNU CC so you
can know your rights and responsibilities.  It should be in a
file named COPYING.  Among other things, the copyright notice
and this notice must be preserved on all copies.  
*/

/*
  Some of the following algorithms are very loosely based on those from 
  MIT C-Scheme bignum.c, which is
      Copyright (c) 1987 Massachusetts Institute of Technology

  with other guidance from Knuth, vol. 2

  Thanks to the creators of the algorithms.
*/

#ifdef __GNUG__
#pragma implementation "Integer.h"
#endif
#include <Integer.h>
#include <std.h>
#include <ctype.h>
#include <math.h>
#include <values.h>
#include <_Obstack.h>
#include <AllocRing.h>
#include <new.h>
/*
 Sizes of shifts for multiple-precision arithmetic.
 These should not be changed unless Integer representation
 as unsigned shorts is changed in the implementation files.
*/

#define I_SHIFT         SHORTBITS
#define I_RADIX         ((unsigned long)(1L << I_SHIFT))
#define I_MAXNUM        ((unsigned long)((I_RADIX - 1)))
#define I_MINNUM        ((unsigned long)(I_RADIX >> 1))
#define I_POSITIVE      1
#define I_NEGATIVE      0

/* All routines assume SHORT_PER_LONG > 1 */
#define SHORT_PER_LONG  ((LONGBITS + SHORTBITS - 1) / SHORTBITS)
#define CHAR_PER_LONG   ((LONGBITS + CHARBITS - 1) / CHARBITS)

/*
  minimum and maximum sizes for an IntRep
*/

#define MINIntRep_SIZE   16
#define MAXIntRep_SIZE   I_MAXNUM

#ifndef MALLOC_MIN_OVERHEAD
#define MALLOC_MIN_OVERHEAD 4
#endif


// utilities to extract and transfer bits 

// get low bits

inline static unsigned short extract(unsigned long x)
{
  return x & I_MAXNUM;
}

// transfer high bits to low

inline static unsigned long down(unsigned long x)
{
  return (x >> I_SHIFT) & I_MAXNUM;
}

// transfer low bits to high

inline static unsigned long up(unsigned long x)
{
  return x << I_SHIFT;
}

// compare two equal-length reps

inline static int docmp(const unsigned short* x, const unsigned short* y, int l)
{
  int diff = 0;
  const unsigned short* xs = &(x[l]);
  const unsigned short* ys = &(y[l]);
  while (l-- > 0 && (diff = (*--xs) - (*--ys)) == 0);
  return diff;
}

// figure out max length of result of +, -, etc.

inline static int calc_len(int len1, int len2, int pad)
{
  return (len1 >= len2)? len1 + pad : len2 + pad;
}

// ensure len & sgn are correct

inline static void Icheck(IntRep* rep)
{
  int l = rep->len;
  const unsigned short* p = &(rep->s[l]);
  while (l > 0 && *--p == 0) --l;
  if ((rep->len = l) == 0) rep->sgn = I_POSITIVE;
}


// zero out the end of a rep

inline static void Iclear_from(IntRep* rep, int p)
{
  unsigned short* cp = &(rep->s[p]);
  const unsigned short* cf = &(rep->s[rep->len]);
  while(cp < cf) *cp++ = 0;
}

// copy parts of a rep

static inline void scpy(const unsigned short* src, unsigned short* dest,int nb)
{
  while (--nb >= 0) *dest++ = *src++;
}

// make sure an argument is valid

static inline void nonnil(const IntRep* rep)
{
  if (rep == 0) 
    (*lib_error_handler)("Integer", "operation on uninitialized Integer");
}

// allocate a new Irep. Pad to something close to a power of two.

inline static IntRep* Inew(int newlen)
{
  unsigned int siz = sizeof(IntRep) + newlen * sizeof(short) + 
    MALLOC_MIN_OVERHEAD;
  unsigned int allocsiz = MINIntRep_SIZE;
  while (allocsiz < siz) allocsiz <<= 1;  // find a power of 2
  allocsiz -= MALLOC_MIN_OVERHEAD;
  if (allocsiz >= MAXIntRep_SIZE * sizeof(short))
    (*lib_error_handler)("Integer", "Requested length out of range");
    
  IntRep* rep = (IntRep *) new char[allocsiz];
  rep->sz = (allocsiz - sizeof(IntRep) + sizeof(short)) / sizeof(short);
  return rep;
}

// allocate: use the bits in src if non-null, clear the rest

IntRep* Ialloc(IntRep* old, const unsigned short* src, int srclen, int newsgn,
              int newlen)
{
  IntRep* rep;
  if (old == 0 || newlen > old->sz)
    rep = Inew(newlen);
  else
    rep = old;

  rep->len = newlen;
  rep->sgn = newsgn;

  scpy(src, rep->s, srclen);
  Iclear_from(rep, srclen);

  if (old != rep && old != 0) delete old;
  return rep;
}

// allocate and clear

IntRep* Icalloc(IntRep* old, int newlen)
{
  IntRep* rep;
  if (old == 0 || newlen > old->sz)
  {
    if (old != 0) delete old;
    rep = Inew(newlen);
  }
  else
    rep = old;

  rep->len = newlen;
  rep->sgn = I_POSITIVE;
  Iclear_from(rep, 0);

  return rep;
}

// reallocate

IntRep* Iresize(IntRep* old, int newlen)
{
  IntRep* rep;
  unsigned short oldlen;
  if (old == 0)
  {
    oldlen = 0;
    rep = Inew(newlen);
    rep->sgn = I_POSITIVE;
  }
  else 
  {
    oldlen = old->len;
    if (newlen > old->sz)
    {
      rep = Inew(newlen);
      scpy(old->s, rep->s, oldlen);
      rep->sgn = old->sgn;
      delete old;
    }
    else
      rep = old;
  }

  rep->len = newlen;
  Iclear_from(rep, oldlen);

  return rep;
}


// same, for straight copy

IntRep* Icopy(IntRep* old, const IntRep* src)
{
  if (old == src) return old; 
  IntRep* rep;
  if (src == 0)
  {
    if (old == 0)
      rep = Inew(0);
    else
    {
      rep = old;
      Iclear_from(rep, 0);
    }
    rep->len = 0;
    rep->sgn = I_POSITIVE;
  }
  else 
  {
    int newlen = src->len;
    if (old == 0 || newlen > old->sz)
    {
      if (old != 0) delete old;
      rep = Inew(newlen);
    }
    else
      rep = old;

    rep->len = newlen;
    rep->sgn = src->sgn;

    scpy(src->s, rep->s, newlen);
  }

  return rep;
}

// allocate & copy space for a long

IntRep* Icopy_long(IntRep* old, long x)
{
  unsigned short src[SHORT_PER_LONG];
  unsigned long u;
  int newsgn;
  if (newsgn = (x >= 0))
    u = x;
  else
    u = -x;
  
  unsigned short srclen = 0;
  while (u != 0)
  {
    src[srclen++] = extract(u);
    u = down(u);
  }

  IntRep* rep;
  if (old == 0 || srclen > old->sz)
  {
    if (old != 0) delete old;
    rep = Inew(srclen);
  }
  else
    rep = old;

  rep->len = srclen;
  rep->sgn = newsgn;

  scpy(src, rep->s, srclen);

  return rep;

}

// special case for zero -- it's worth it!

IntRep* Icopy_zero(IntRep* old)
{
  IntRep* rep;
  if (old == 0)
    rep = Inew(0);
  else
    rep = old;

  rep->len = 0;
  rep->sgn = I_POSITIVE;

  return rep;
}

// special case for 1 or -1

IntRep* Icopy_one(IntRep* old, int newsgn)
{
  IntRep* rep;
  if (old == 0 || 1 > old->sz)
  {
    if (old != 0) delete old;
    rep = Inew(1);
  }
  else
    rep = old;

  rep->sgn = newsgn;
  rep->len = 1;
  rep->s[0] = 1;

  return rep;
}

// convert to a legal two's complement long if possible
// if too big, return most negative/positive value

long Itolong(const IntRep* rep)
{ 
  if (rep->len > SHORT_PER_LONG)
    return (rep->sgn == I_POSITIVE) ? MAXLONG : MINLONG;
  else if (rep->len == 0)
    return 0;
  else if (rep->len < SHORT_PER_LONG)
  {
    unsigned long a = rep->s[rep->len-1];
    if (SHORT_PER_LONG > 2) // normally optimized out
    {
      for (int i = rep->len - 2; i >= 0; --i)
        a = up(a) | rep->s[i];
    }
    return (rep->sgn == I_POSITIVE)? a : -((long)a);
  }
  else 
  {
    unsigned long a = rep->s[SHORT_PER_LONG - 1];
    if (a >= I_MINNUM)
      return (rep->sgn == I_POSITIVE) ? MAXLONG : MINLONG;
    else
    {
      a = up(a) | rep->s[SHORT_PER_LONG - 2];
      if (SHORT_PER_LONG > 2)
      {
        for (int i = SHORT_PER_LONG - 3; i >= 0; --i)
          a = up(a) | rep->s[i];
      }
      return (rep->sgn == I_POSITIVE)? a : -((long)a);
    }
  }
}

// test whether op long() will work.
// careful about asymmetry between MINLONG & MAXLONG

int Iislong(const IntRep* rep)
{
  int l = rep->len;
  if (l < SHORT_PER_LONG)
    return 1;
  else if (l > SHORT_PER_LONG)
    return 0;
  else if (rep->s[SHORT_PER_LONG - 1] < I_MINNUM)
    return 1;
  else if (rep->sgn == I_NEGATIVE && rep->s[SHORT_PER_LONG - 1] == I_MINNUM)
  {
    for (int i = 0; i < SHORT_PER_LONG - 1; ++i)
      if (rep->s[i] != 0)
        return 0;
    return 1;
  }
  else
    return 0;
}

// convert to a double 

double Itodouble(const IntRep* rep)
{ 
  double d = 0.0;
  double bound = HUGE / 2.0;
  for (int i = rep->len - 1; i >= 0; --i)
  {
    unsigned short a = I_RADIX >> 1;
    while (a != 0)
    {
      if (d >= bound)
        return (rep->sgn == I_NEGATIVE) ? -HUGE : HUGE;
      d *= 2.0;
      if (rep->s[i] & a)
        d += 1.0;
      a >>= 1;
    }
  }
  if (rep->sgn == I_NEGATIVE)
    return -d;
  else
    return d;
}

// see whether op double() will work-
// have to actually try it in order to find out
// since otherwise might trigger fp exception

int Iisdouble(const IntRep* rep)
{
  double d = 0.0;
  double bound = HUGE / 2.0;
  for (int i = rep->len - 1; i >= 0; --i)
  {
    unsigned short a = I_RADIX >> 1;
    while (a != 0)
    {
      if (d > bound || (d == bound && (i > 0 || (rep->s[i] & a))))
        return 0;
      d *= 2.0;
      if (rep->s[i] & a)
        d += 1.0;
      a >>= 1;
    }
  }
  return 1;
}

// real division of num / den

double ratio(const Integer& num, const Integer& den)
{
  Integer q, r;
  divide(num, den, q, r);
  double d1 = double(q);
 
  if (d1 == HUGE || d1 == -HUGE || sign(r) == 0)
    return d1;
  else      // use as much precision as available for fractional part
  {
    double  d2 = 0.0;
    double  d3 = 0.0; 
    double bound = HUGE / 2.0;
    int cont = 1;
    for (int i = den.rep->len - 1; i >= 0 && cont; --i)
    {
      unsigned short a = I_RADIX >> 1;
      while (a != 0)
      {
        if (d2 >= bound)
        {
          cont = 0;
          break;
        }
        d2 *= 2.0;
        if (den.rep->s[i] & a)
          d2 += 1.0;

        if (i < r.rep->len)
        {
          d3 *= 2.0;
          if (r.rep->s[i] & a)
            d3 += 1.0;
        }
        a >>= 1;
      }
    }

    if (sign(r) < 0)
      d3 = -d3;
    return d1 + d3 / d2;
  }
}

// comparison functions
  
int compare(const IntRep* x, const IntRep* y)
{
  int diff  = x->sgn - y->sgn;
  if (diff == 0)
  {
    diff = x->len - y->len;
    if (diff == 0)
      diff = docmp(x->s, y->s, x->len);
    if (x->sgn == I_NEGATIVE)
      diff = -diff;
  }
  return diff;
}

int ucompare(const IntRep* x, const IntRep* y)
{
  int diff = x->len - y->len;
  if (diff == 0)
  {
    int l = x->len;
    const unsigned short* xs = &(x->s[l]);
    const unsigned short* ys = &(y->s[l]);
    while (l-- > 0 && (diff = (*--xs) - (*--ys)) == 0);
  }
  return diff;
}

int compare(const IntRep* x, long  y)
{
  int xl = x->len;
  int xsgn = x->sgn;
  if (y == 0)
  {
    if (xl == 0)
      return 0;
    else if (xsgn == I_NEGATIVE)
      return -1;
    else
      return 1;
  }
  else
  {
    int ysgn = y >= 0;
    unsigned long uy = (ysgn)? y : -y;
    int diff = xsgn - ysgn;
    if (diff == 0)
    {
      diff = xl - SHORT_PER_LONG;
      if (diff <= 0)
      {
        unsigned short tmp[SHORT_PER_LONG];
        int yl = 0;
        while (uy != 0)
        {
          tmp[yl++] = extract(uy);
          uy = down(uy);
        }
        diff = xl - yl;
        if (diff == 0)
          diff = docmp(x->s, tmp, xl);
      }
    }
    if (xsgn == I_NEGATIVE)
      diff = -diff;
    return diff;
  }
}

int ucompare(const IntRep* x, long  y)
{
  int xl = x->len;
  if (y == 0)
    return xl;
  else
  {
    unsigned long uy = (y >= 0)? y : -y;
    int diff = xl - SHORT_PER_LONG;
    if (diff <= 0)
    {
      unsigned short tmp[SHORT_PER_LONG];
      int yl = 0;
      while (uy != 0)
      {
        tmp[yl++] = extract(uy);
        uy = down(uy);
      }
      diff = xl - yl;
      if (diff == 0)
        diff = docmp(x->s, tmp, xl);
    }
    return diff;
  }
}



// arithmetic functions

IntRep* add(const IntRep* x, int negatex, 
            const IntRep* y, int negatey, IntRep* r)
{
  nonnil(x);
  nonnil(y);

  int xl = x->len;
  int yl = y->len;

  int xsgn = (negatex && xl != 0) ? !x->sgn : x->sgn;
  int ysgn = (negatey && yl != 0) ? !y->sgn : y->sgn;

  int xrsame = x == r;
  int yrsame = y == r;

  if (yl == 0)
    r = Ialloc(r, x->s, xl, xsgn, xl);
  else if (xl == 0)
    r = Ialloc(r, y->s, yl, ysgn, yl);
  else if (xsgn == ysgn)
  {
    if (xrsame || yrsame)
      r = Iresize(r, calc_len(xl, yl, 1));
    else
      r = Icalloc(r, calc_len(xl, yl, 1));
    r->sgn = xsgn;
    unsigned short* rs = r->s;
    const unsigned short* as;
    const unsigned short* bs;
    const unsigned short* topa;
    const unsigned short* topb;
    if (xl >= yl)
    {
      as =  (xrsame)? r->s : x->s;
      topa = &(as[xl]);
      bs =  (yrsame)? r->s : y->s;
      topb = &(bs[yl]);
    }
    else
    {
      bs =  (xrsame)? r->s : x->s;
      topb = &(bs[xl]);
      as =  (yrsame)? r->s : y->s;
      topa = &(as[yl]);
    }
    unsigned long sum = 0;
    while (bs < topb)
    {
      sum += (unsigned long)(*as++) + (unsigned long)(*bs++);
      *rs++ = extract(sum);
      sum = down(sum);
    }
    while (sum != 0 && as < topa)
    {
      sum += (unsigned long)(*as++);
      *rs++ = extract(sum);
      sum = down(sum);
    }
    if (sum != 0)
      *rs = extract(sum);
    else if (rs != as)
      while (as < topa)
        *rs++ = *as++;
  }
  else
  {
    int comp = ucompare(x, y);
    if (comp == 0)
      r = Icopy_zero(r);
    else
    {
      if (xrsame || yrsame)
        r = Iresize(r, calc_len(xl, yl, 0));
      else
        r = Icalloc(r, calc_len(xl, yl, 0));
      unsigned short* rs = r->s;
      const unsigned short* as;
      const unsigned short* bs;
      const unsigned short* topa;
      const unsigned short* topb;
      if (comp > 0)
      {
        as =  (xrsame)? r->s : x->s;
        topa = &(as[xl]);
        bs =  (yrsame)? r->s : y->s;
        topb = &(bs[yl]);
        r->sgn = xsgn;
      }
      else
      {
        bs =  (xrsame)? r->s : x->s;
        topb = &(bs[xl]);
        as =  (yrsame)? r->s : y->s;
        topa = &(as[yl]);
        r->sgn = ysgn;
      }
      unsigned long hi = 1;
      while (bs < topb)
      {
        hi += (unsigned long)(*as++) + I_MAXNUM - (unsigned long)(*bs++);
        *rs++ = extract(hi);
        hi = down(hi);
      }
      while (hi == 0 && as < topa)
      {
        hi = (unsigned long)(*as++) + I_MAXNUM;
        *rs++ = extract(hi);
        hi = down(hi);
      }
      if (rs != as)
        while (as < topa)
          *rs++ = *as++;
    }
  }
  Icheck(r);
  return r;
}


IntRep* add(const IntRep* x, int negatex, long y, IntRep* r)
{
  nonnil(x);
  int xl = x->len;
  int xsgn = (negatex && xl != 0) ? !x->sgn : x->sgn;
  int xrsame = x == r;

  int ysgn = (y >= 0);
  unsigned long uy = (ysgn)? y : -y;

  if (y == 0)
    r = Ialloc(r, x->s, xl, xsgn, xl);
  else if (xl == 0)
    r = Icopy_long(r, y);
  else if (xsgn == ysgn)
  {
    if (xrsame)
      r = Iresize(r, calc_len(xl, SHORT_PER_LONG, 1));
    else
      r = Icalloc(r, calc_len(xl, SHORT_PER_LONG, 1));
    r->sgn = xsgn;
    unsigned short* rs = r->s;
    const unsigned short* as =  (xrsame)? r->s : x->s;
    const unsigned short* topa = &(as[xl]);
    unsigned long sum = 0;
    while (as < topa && uy != 0)
    {
      unsigned long u = extract(uy);
      uy = down(uy);
      sum += (unsigned long)(*as++) + u;
      *rs++ = extract(sum);
      sum = down(sum);
    }
    while (sum != 0 && as < topa)
    {
      sum += (unsigned long)(*as++);
      *rs++ = extract(sum);
      sum = down(sum);
    }
    if (sum != 0)
      *rs = extract(sum);
    else if (rs != as)
      while (as < topa)
        *rs++ = *as++;
  }
  else
  {
    unsigned short tmp[SHORT_PER_LONG];
    int yl = 0;
    while (uy != 0)
    {
      tmp[yl++] = extract(uy);
      uy = down(uy);
    }
    int comp = xl - yl;
    if (comp == 0)
      comp = docmp(x->s, tmp, yl);
    if (comp == 0)
      r = Icopy_zero(r);
    else
    {
      if (xrsame)
        r = Iresize(r, calc_len(xl, yl, 0));
      else
        r = Icalloc(r, calc_len(xl, yl, 0));
      unsigned short* rs = r->s;
      const unsigned short* as;
      const unsigned short* bs;
      const unsigned short* topa;
      const unsigned short* topb;
      if (comp > 0)
      {
        as =  (xrsame)? r->s : x->s;
        topa = &(as[xl]);
        bs =  tmp;
        topb = &(bs[yl]);
        r->sgn = xsgn;
      }
      else
      {
        bs =  (xrsame)? r->s : x->s;
        topb = &(bs[xl]);
        as =  tmp;
        topa = &(as[yl]);
        r->sgn = ysgn;
      }
      unsigned long hi = 1;
      while (bs < topb)
      {
        hi += (unsigned long)(*as++) + I_MAXNUM - (unsigned long)(*bs++);
        *rs++ = extract(hi);
        hi = down(hi);
      }
      while (hi == 0 && as < topa)
      {
        hi = (unsigned long)(*as++) + I_MAXNUM;
        *rs++ = extract(hi);
        hi = down(hi);
      }
      if (rs != as)
        while (as < topa)
          *rs++ = *as++;
    }
  }
  Icheck(r);
  return r;
}


IntRep* multiply(const IntRep* x, const IntRep* y, IntRep* r)
{
  nonnil(x);
  nonnil(y);
  int xl = x->len;
  int yl = y->len;
  int rl = xl + yl;
  int rsgn = x->sgn == y->sgn;
  int xrsame = x == r;
  int yrsame = y == r;
  int xysame = x == y;
  
  if (xl == 0 || yl == 0)
    r = Icopy_zero(r);
  else if (xl == 1 && x->s[0] == 1)
    r = Icopy(r, y);
  else if (yl == 1 && y->s[0] == 1)
    r = Icopy(r, x);
  else if (!(xysame && xrsame))
  {
    if (xrsame || yrsame)
      r = Iresize(r, rl);
    else
      r = Icalloc(r, rl);
    unsigned short* rs = r->s;
    unsigned short* topr = &(rs[rl]);

    // use best inner/outer loop params given constraints
    unsigned short* currentr;
    const unsigned short* bota;
    const unsigned short* as;
    const unsigned short* botb;
    const unsigned short* topb;
    if (xrsame)                 
    { 
      currentr = &(rs[xl-1]);
      bota = rs;
      as = currentr;
      botb = y->s;
      topb = &(botb[yl]);
    }
    else if (yrsame)
    {
      currentr = &(rs[yl-1]);
      bota = rs;
      as = currentr;
      botb = x->s;
      topb = &(botb[xl]);
    }
    else if (xl <= yl)
    {
      currentr = &(rs[xl-1]);
      bota = x->s;
      as = &(bota[xl-1]);
      botb = y->s;
      topb = &(botb[yl]);
    }
    else
    {
      currentr = &(rs[yl-1]);
      bota = y->s;
      as = &(bota[yl-1]);
      botb = x->s;
      topb = &(botb[xl]);
    }

    while (as >= bota)
    {
      unsigned long ai = (unsigned long)(*as--);
      unsigned short* rs = currentr--;
      *rs = 0;
      if (ai != 0)
      {
        unsigned long sum = 0;
        const unsigned short* bs = botb;
        while (bs < topb)
        {
          sum += ai * (unsigned long)(*bs++) + (unsigned long)(*rs);
          *rs++ = extract(sum);
          sum = down(sum);
        }
        while (sum != 0 && rs < topr)
        {
          sum += (unsigned long)(*rs);
          *rs++ = extract(sum);
          sum = down(sum);
        }
      }
    }
  }
  else                          // x, y, and r same; compute over diagonals
  {
    r = Iresize(r, rl);
    unsigned short* botr = r->s;
    unsigned short* topr = &(botr[rl]);
    unsigned short* rs =   &(botr[rl - 2]);

    const unsigned short* bota = (xrsame)? botr : x->s;
    const unsigned short* loa =  &(bota[xl - 1]);
    const unsigned short* hia =  loa;

    for (; rs >= botr; --rs)
    {
      const unsigned short* h = hia;
      const unsigned short* l = loa;
      unsigned long prod = (unsigned long)(*h) * (unsigned long)(*l);
      *rs = 0;

      for(;;)
      {
        unsigned short* rt = rs;
        unsigned long sum = prod + (unsigned long)(*rt);
        *rt++ = extract(sum);
        sum = down(sum);
        while (sum != 0 && rt < topr)
        {
          sum += (unsigned long)(*rt);
          *rt++ = extract(sum);
          sum = down(sum);
        }
        if (h > l)
        {
          rt = rs;
          sum = prod + (unsigned long)(*rt);
          *rt++ = extract(sum);
          sum = down(sum);
          while (sum != 0 && rt < topr)
          {
            sum += (unsigned long)(*rt);
            *rt++ = extract(sum);
            sum = down(sum);
          }
          if (--h >= ++l)
            prod = (unsigned long)(*h) * (unsigned long)(*l);
          else
            break;
        }
        else
          break;
      }
      if (loa > bota)
        --loa;
      else
        --hia;
    }
  }
  r->sgn = rsgn;
  Icheck(r);
  return r;
}


IntRep* multiply(const IntRep* x, long y, IntRep* r)
{
  nonnil(x);
  int xl = x->len;
    
  if (xl == 0 || y == 0)
    r = Icopy_zero(r);
  else if (y == 1)
    r = Icopy(r, x);
  else
  {
    int ysgn = y >= 0;
    int rsgn = x->sgn == ysgn;
    unsigned long uy = (ysgn)? y : -y;
    unsigned short tmp[SHORT_PER_LONG];
    int yl = 0;
    while (uy != 0)
    {
      tmp[yl++] = extract(uy);
      uy = down(uy);
    }

    int rl = xl + yl;
    int xrsame = x == r;
    if (xrsame)
      r = Iresize(r, rl);
    else
      r = Icalloc(r, rl);

    unsigned short* rs = r->s;
    unsigned short* topr = &(rs[rl]);
    unsigned short* currentr;
    const unsigned short* bota;
    const unsigned short* as;
    const unsigned short* botb;
    const unsigned short* topb;

    if (xrsame)
    { 
      currentr = &(rs[xl-1]);
      bota = rs;
      as = currentr;
      botb = tmp;
      topb = &(botb[yl]);
    }
    else if (xl <= yl)
    {
      currentr = &(rs[xl-1]);
      bota = x->s;
      as = &(bota[xl-1]);
      botb = tmp;
      topb = &(botb[yl]);
    }
    else
    {
      currentr = &(rs[yl-1]);
      bota = tmp;
      as = &(bota[yl-1]);
      botb = x->s;
      topb = &(botb[xl]);
    }

    while (as >= bota)
    {
      unsigned long ai = (unsigned long)(*as--);
      unsigned short* rs = currentr--;
      *rs = 0;
      if (ai != 0)
      {
        unsigned long sum = 0;
        const unsigned short* bs = botb;
        while (bs < topb)
        {
          sum += ai * (unsigned long)(*bs++) + (unsigned long)(*rs);
          *rs++ = extract(sum);
          sum = down(sum);
        }
        while (sum != 0 && rs < topr)
        {
          sum += (unsigned long)(*rs);
          *rs++ = extract(sum);
          sum = down(sum);
        }
      }
    }
    r->sgn = rsgn;
  }
  Icheck(r);
  return r;
}


// main division routine

static void do_divide(unsigned short* rs,
                      const unsigned short* ys, int yl,
                      unsigned short* qs, int ql)
{
  const unsigned short* topy = &(ys[yl]);
  unsigned short d1 = ys[yl - 1];
  unsigned short d2 = ys[yl - 2];
 
  int l = ql - 1;
  int i = l + yl;
  
  for (; l >= 0; --l, --i)
  {
    unsigned short qhat;       // guess q
    if (d1 == rs[i])
      qhat = I_MAXNUM;
    else
    {
      unsigned long lr = up((unsigned long)rs[i]) | rs[i-1];
      qhat = lr / d1;
    }

    for(;;)     // adjust q, use docmp to avoid overflow problems
    {
      unsigned short ts[3];
      unsigned long prod = (unsigned long)d2 * (unsigned long)qhat;
      ts[0] = extract(prod);
      prod = down(prod) + (unsigned long)d1 * (unsigned long)qhat;
      ts[1] = extract(prod);
      ts[2] = extract(down(prod));
      if (docmp(ts, &(rs[i-2]), 3) > 0)
        --qhat;
      else
        break;
    };
    
    // multiply & subtract
    
    const unsigned short* yt = ys;
    unsigned short* rt = &(rs[l]);
    unsigned long prod = 0;
    unsigned long hi = 1;
    while (yt < topy)
    {
      prod = (unsigned long)qhat * (unsigned long)(*yt++) + down(prod);
      hi += (unsigned long)(*rt) + I_MAXNUM - (unsigned long)(extract(prod));
      *rt++ = extract(hi);
      hi = down(hi);
    }
    hi += (unsigned long)(*rt) + I_MAXNUM - (unsigned long)(down(prod));
    *rt = extract(hi);
    hi = down(hi);
    
    // off-by-one, add back
    
    if (hi == 0)
    {
      --qhat;
      yt = ys;
      rt = &(rs[l]);
      hi = 0;
      while (yt < topy)
      {
        hi = (unsigned long)(*rt) + (unsigned long)(*yt++) + down(hi);
        *rt++ = extract(hi);
      }
      *rt = 0;
    }
    if (qs != 0)
      qs[l] = qhat;
  }
}

// divide by single digit, return remainder
// if q != 0, then keep the result in q, else just compute rem

static int unscale(const unsigned short* x, int xl, unsigned short y,
                   unsigned short* q)
{
  if (xl == 0 || y == 1)
    return 0;
  else if (q != 0)
  {
    unsigned short* botq = q;
    unsigned short* qs = &(botq[xl - 1]);
    const unsigned short* xs = &(x[xl - 1]);
    unsigned long rem = 0;
    while (qs >= botq)
    {
      rem = up(rem) | *xs--;
      unsigned long u = rem / y;
      *qs-- = extract(u);
      rem -= u * y;
    }
    int r = extract(rem);
    return r;
  }
  else                          // same loop, a bit faster if just need rem
  {
    const unsigned short* botx = x;
    const unsigned short* xs = &(botx[xl - 1]);
    unsigned long rem = 0;
    while (xs >= botx)
    {
      rem = up(rem) | *xs--;
      unsigned long u = rem / y;
      rem -= u * y;
    }
    int r = extract(rem);
    return r;
  }
}


IntRep* div(const IntRep* x, const IntRep* y, IntRep* q)
{
  nonnil(x);
  nonnil(y);
  int xl = x->len;
  int yl = y->len;
  if (yl == 0) (*lib_error_handler)("Integer", "attempted division by zero");

  int comp = ucompare(x, y);
  int xsgn = x->sgn;
  int ysgn = y->sgn;

  int samesign = xsgn == ysgn;

  if (comp < 0)
    q = Icopy_zero(q);
  else if (comp == 0)
    q = Icopy_one(q, samesign);
  else if (yl == 1)
  {
    q = Icopy(q, x);
    unscale(q->s, q->len, y->s[0], q->s);
  }
  else
  {
    IntRep* yy = 0;
    IntRep* r  = 0;
    unsigned short prescale = (I_RADIX / (1 + y->s[yl - 1]));
    if (prescale != 1 || y == q)
    {
      yy = multiply(y, ((long)prescale & I_MAXNUM), yy);
      r = multiply(x, ((long)prescale & I_MAXNUM), r);
    }
    else
    {
      yy = (IntRep*)y;
      r = Icalloc(r, xl + 1);
      scpy(x->s, r->s, xl);
    }

    int ql = xl - yl + 1;
      
    q = Icalloc(q, ql);
    do_divide(r->s, yy->s, yl, q->s, ql);

    if (yy != y) delete yy;
    delete r;
  }
  q->sgn = samesign;
  Icheck(q);
  return q;
}

IntRep* div(const IntRep* x, long y, IntRep* q)
{
  nonnil(x);
  int xl = x->len;
  if (y == 0) (*lib_error_handler)("Integer", "attempted division by zero");

  unsigned short ys[SHORT_PER_LONG];
  unsigned long u;
  int ysgn = y >= 0;
  if (ysgn)
    u = y;
  else
    u = -y;
  int yl = 0;
  while (u != 0)
  {
    ys[yl++] = extract(u);
    u = down(u);
  }

  int comp = xl - yl;
  if (comp == 0) comp = docmp(x->s, ys, xl);

  int xsgn = x->sgn;
  int samesign = xsgn == ysgn;

  if (comp < 0)
    q = Icopy_zero(q);
  else if (comp == 0)
  {
    q = Icopy_one(q, samesign);
  }
  else if (yl == 1)
  {
    q = Icopy(q, x);
    unscale(q->s, q->len, ys[0], q->s);
  }
  else
  {
    IntRep* r  = 0;
    unsigned short prescale = (I_RADIX / (1 + ys[yl - 1]));
    if (prescale != 1)
    {
      unsigned long prod = (unsigned long)prescale * (unsigned long)ys[0];
      ys[0] = extract(prod);
      prod = down(prod) + (unsigned long)prescale * (unsigned long)ys[1];
      ys[1] = extract(prod);
      r = multiply(x, ((long)prescale & I_MAXNUM), r);
    }
    else
    {
      r = Icalloc(r, xl + 1);
      scpy(x->s, r->s, xl);
    }

    int ql = xl - yl + 1;
      
    q = Icalloc(q, ql);
    do_divide(r->s, ys, yl, q->s, ql);

    delete r;
  }
  q->sgn = samesign;
  Icheck(q);
  return q;
}


void divide(const Integer& Ix, long y, Integer& Iq, long& rem)
{
  const IntRep* x = Ix.rep;
  nonnil(x);
  IntRep* q = Iq.rep;
  int xl = x->len;
  if (y == 0) (*lib_error_handler)("Integer", "attempted division by zero");

  unsigned short ys[SHORT_PER_LONG];
  unsigned long u;
  int ysgn = y >= 0;
  if (ysgn)
    u = y;
  else
    u = -y;
  int yl = 0;
  while (u != 0)
  {
    ys[yl++] = extract(u);
    u = down(u);
  }

  int comp = xl - yl;
  if (comp == 0) comp = docmp(x->s, ys, xl);

  int xsgn = x->sgn;
  int samesign = xsgn == ysgn;

  if (comp < 0)
  {
    rem = Itolong(x);
    q = Icopy_zero(q);
  }
  else if (comp == 0)
  {
    q = Icopy_one(q, samesign);
    rem = 0;
  }
  else if (yl == 1)
  {
    q = Icopy(q, x);
    rem = unscale(q->s, q->len, ys[0], q->s);
  }
  else
  {
    IntRep* r  = 0;
    unsigned short prescale = (I_RADIX / (1 + ys[yl - 1]));
    if (prescale != 1)
    {
      unsigned long prod = (unsigned long)prescale * (unsigned long)ys[0];
      ys[0] = extract(prod);
      prod = down(prod) + (unsigned long)prescale * (unsigned long)ys[1];
      ys[1] = extract(prod);
      r = multiply(x, ((long)prescale & I_MAXNUM), r);
    }
    else
    {
      r = Icalloc(r, xl + 1);
      scpy(x->s, r->s, xl);
    }

    int ql = xl - yl + 1;
      
    q = Icalloc(q, ql);
    
    do_divide(r->s, ys, yl, q->s, ql);

    if (prescale != 1)
    {
      Icheck(r);
      unscale(r->s, r->len, prescale, r->s);
    }
    Icheck(r);
    rem = Itolong(r);
    delete r;
  }
  rem = abs(rem);
  if (xsgn == I_NEGATIVE) rem = -rem;
  q->sgn = samesign;
  Icheck(q);
  Iq.rep = q;
}


void divide(const Integer& Ix, const Integer& Iy, Integer& Iq, Integer& Ir)
{
  const IntRep* x = Ix.rep;
  nonnil(x);
  const IntRep* y = Iy.rep;
  nonnil(y);
  IntRep* q = Iq.rep;
  IntRep* r = Ir.rep;

  int xl = x->len;
  int yl = y->len;
  if (yl == 0)
    (*lib_error_handler)("Integer", "attempted division by zero");

  int comp = ucompare(x, y);
  int xsgn = x->sgn;
  int ysgn = y->sgn;

  int samesign = xsgn == ysgn;

  if (comp < 0)
  {
    q = Icopy_zero(q);
    r = Icopy(r, x);
  }
  else if (comp == 0)
  {
    q = Icopy_one(q, samesign);
    r = Icopy_zero(r);
  }
  else if (yl == 1)
  {
    q = Icopy(q, x);
    int rem =  unscale(q->s, q->len, y->s[0], q->s);
    r = Icopy_long(r, rem);
    if (rem != 0)
      r->sgn = xsgn;
  }
  else
  {
    IntRep* yy = 0;
    unsigned short prescale = (I_RADIX / (1 + y->s[yl - 1]));
    if (prescale != 1 || y == q || y == r)
    {
      yy = multiply(y, ((long)prescale & I_MAXNUM), yy);
      r = multiply(x, ((long)prescale & I_MAXNUM), r);
    }
    else
    {
      yy = (IntRep*)y;
      r = Icalloc(r, xl + 1);
      scpy(x->s, r->s, xl);
    }

    int ql = xl - yl + 1;
      
    q = Icalloc(q, ql);
    do_divide(r->s, yy->s, yl, q->s, ql);

    if (yy != y) delete yy;
    if (prescale != 1)
    {
      Icheck(r);
      unscale(r->s, r->len, prescale, r->s);
    }
  }
  q->sgn = samesign;
  Icheck(q);
  Iq.rep = q;
  Icheck(r);
  Ir.rep = r;
}

IntRep* mod(const IntRep* x, const IntRep* y, IntRep* r)
{
  nonnil(x);
  nonnil(y);
  int xl = x->len;
  int yl = y->len;
  if (yl == 0) (*lib_error_handler)("Integer", "attempted division by zero");

  int comp = ucompare(x, y);
  int xsgn = x->sgn;

  if (comp < 0)
    r = Icopy(r, x);
  else if (comp == 0)
    r = Icopy_zero(r);
  else if (yl == 1)
  {
    int rem = unscale(x->s, xl, y->s[0], 0);
    r = Icopy_long(r, rem);
    if (rem != 0)
      r->sgn = xsgn;
  }
  else
  {
    IntRep* yy = 0;
    unsigned short prescale = (I_RADIX / (1 + y->s[yl - 1]));
    if (prescale != 1 || y == r)
    {
      yy = multiply(y, ((long)prescale & I_MAXNUM), yy);
      r = multiply(x, ((long)prescale & I_MAXNUM), r);
    }
    else
    {
      yy = (IntRep*)y;
      r = Icalloc(r, xl + 1);
      scpy(x->s, r->s, xl);
    }
      
    do_divide(r->s, yy->s, yl, 0, xl - yl + 1);

    if (yy != y) delete yy;

    if (prescale != 1)
    {
      Icheck(r);
      unscale(r->s, r->len, prescale, r->s);
    }
  }
  Icheck(r);
  return r;
}

IntRep* mod(const IntRep* x, long y, IntRep* r)
{
  nonnil(x);
  int xl = x->len;
  if (y == 0) (*lib_error_handler)("Integer", "attempted division by zero");

  unsigned short ys[SHORT_PER_LONG];
  unsigned long u;
  int ysgn = y >= 0;
  if (ysgn)
    u = y;
  else
    u = -y;
  int yl = 0;
  while (u != 0)
  {
    ys[yl++] = extract(u);
    u = down(u);
  }

  int comp = xl - yl;
  if (comp == 0) comp = docmp(x->s, ys, xl);

  int xsgn = x->sgn;

  if (comp < 0)
    r = Icopy(r, x);
  else if (comp == 0)
    r = Icopy_zero(r);
  else if (yl == 1)
  {
    int rem = unscale(x->s, xl, ys[0], 0);
    r = Icopy_long(r, rem);
    if (rem != 0)
      r->sgn = xsgn;
  }
  else
  {
    unsigned short prescale = (I_RADIX / (1 + ys[yl - 1]));
    if (prescale != 1)
    {
      unsigned long prod = (unsigned long)prescale * (unsigned long)ys[0];
      ys[0] = extract(prod);
      prod = down(prod) + (unsigned long)prescale * (unsigned long)ys[1];
      ys[1] = extract(prod);
      r = multiply(x, ((long)prescale & I_MAXNUM), r);
    }
    else
    {
      r = Icalloc(r, xl + 1);
      scpy(x->s, r->s, xl);
    }
      
    do_divide(r->s, ys, yl, 0, xl - yl + 1);

    if (prescale != 1)
    {
      Icheck(r);
      unscale(r->s, r->len, prescale, r->s);
    }
  }
  Icheck(r);
  return r;
}

IntRep* lshift(const IntRep* x, long y, IntRep* r)
{
  nonnil(x);
  int xl = x->len;
  if (xl == 0 || y == 0)
  {
    r = Icopy(r, x);
    return r;
  }

  int xrsame = x == r;
  int rsgn = x->sgn;

  long ay = (y < 0)? -y : y;
  int bw = ay / I_SHIFT;
  int sw = ay % I_SHIFT;

  if (y > 0)
  {
    int rl = bw + xl + 1;
    if (xrsame)
      r = Iresize(r, rl);
    else
      r = Icalloc(r, rl);

    unsigned short* botr = r->s;
    unsigned short* rs = &(botr[rl - 1]);
    const unsigned short* botx = (xrsame)? botr : x->s;
    const unsigned short* xs = &(botx[xl - 1]);
    unsigned long a = 0;
    while (xs >= botx)
    {
      a = up(a) | ((unsigned long)(*xs--) << sw);
      *rs-- = extract(down(a));
    }
    *rs-- = extract(a);
    while (rs >= botr)
      *rs-- = 0;
  }
  else
  {
    int rl = xl - bw;
    if (rl < 0)
      r = Icopy_zero(r);
    else
    {
      if (xrsame)
        r = Iresize(r, rl);
      else
        r = Icalloc(r, rl);
      int rw = I_SHIFT - sw;
      unsigned short* rs = r->s;
      unsigned short* topr = &(rs[rl]);
      const unsigned short* botx = (xrsame)? rs : x->s;
      const unsigned short* xs =  &(botx[bw]);
      const unsigned short* topx = &(botx[xl]);
      unsigned long a = (unsigned long)(*xs++) >> sw;
      while (xs < topx)
      {
        a |= (unsigned long)(*xs++) << rw;
        *rs++ = extract(a);
        a = down(a);
      }
      *rs++ = extract(a);
      if (xrsame) topr = (unsigned short*)topx;
      while (rs < topr)
        *rs++ = 0;
    }
  }
  r->sgn = rsgn;
  Icheck(r);
  return r;
}

IntRep* lshift(const IntRep* x, const IntRep* yy, int negatey, IntRep* r)
{
  long y = Itolong(yy);
  if (negatey)
    y = -y;

  return lshift(x, y, r);
}

IntRep* bitop(const IntRep* x, const IntRep* y, IntRep* r, char op)
{
  nonnil(x);
  nonnil(y);
  int xl = x->len;
  int yl = y->len;
  int xsgn = x->sgn;
  int xrsame = x == r;
  int yrsame = y == r;
  if (xrsame || yrsame)
    r = Iresize(r, calc_len(xl, yl, 0));
  else
    r = Icalloc(r, calc_len(xl, yl, 0));
  r->sgn = xsgn;
  unsigned short* rs = r->s;
  unsigned short* topr = &(rs[r->len]);
  const unsigned short* as;
  const unsigned short* bs;
  const unsigned short* topb;
  if (xl >= yl)
  {
    as = (xrsame)? rs : x->s;
    bs = (yrsame)? rs : y->s;
    topb = &(bs[yl]);
  }
  else
  {
    bs = (xrsame)? rs : x->s;
    topb = &(bs[xl]);
    as = (yrsame)? rs : y->s;
  }

  switch (op)
  {
  case '&':
    while (bs < topb) *rs++ = *as++ & *bs++;
    while (rs < topr) *rs++ = 0;
    break;
  case '|':
    while (bs < topb) *rs++ = *as++ | *bs++;
    while (rs < topr) *rs++ = *as++;
    break;
  case '^':
    while (bs < topb) *rs++ = *as++ ^ *bs++;
    while (rs < topr) *rs++ = *as++;
    break;
  }
  Icheck(r);
  return r;
}

IntRep* bitop(const IntRep* x, long y, IntRep* r, char op)
{
  nonnil(x);
  unsigned short tmp[SHORT_PER_LONG];
  unsigned long u;
  int newsgn;
  if (newsgn = (y >= 0))
    u = y;
  else
    u = -y;
  
  int l = 0;
  while (u != 0)
  {
    tmp[l++] = extract(u);
    u = down(u);
  }

  int xl = x->len;
  int yl = l;
  int xsgn = x->sgn;
  int xrsame = x == r;
  if (xrsame)
    r = Iresize(r, calc_len(xl, yl, 0));
  else
    r = Icalloc(r, calc_len(xl, yl, 0));
  r->sgn = xsgn;
  unsigned short* rs = r->s;
  unsigned short* topr = &(rs[r->len]);
  const unsigned short* as;
  const unsigned short* bs;
  const unsigned short* topb;
  if (xl >= yl)
  {
    as = (xrsame)? rs : x->s;
    bs = tmp;
    topb = &(bs[yl]);
  }
  else
  {
    bs = (xrsame)? rs : x->s;
    topb = &(bs[xl]);
    as = tmp;
  }

  switch (op)
  {
  case '&':
    while (bs < topb) *rs++ = *as++ & *bs++;
    while (rs < topr) *rs++ = 0;
    break;
  case '|':
    while (bs < topb) *rs++ = *as++ | *bs++;
    while (rs < topr) *rs++ = *as++;
    break;
  case '^':
    while (bs < topb) *rs++ = *as++ ^ *bs++;
    while (rs < topr) *rs++ = *as++;
    break;
  }
  Icheck(r);
  return r;
}



IntRep*  compl(const IntRep* src, IntRep* r)
{
  nonnil(src);
  r = Icopy(r, src);
  unsigned short* s = r->s;
  unsigned short* top = &(s[r->len - 1]);
  while (s < top)
  {
    unsigned short cmp = ~(*s);
    *s++ = cmp;
  }
  unsigned short a = *s;
  unsigned short b = 0;
  while (a != 0)
  {
    b <<= 1;
    if (!(a & 1)) b |= 1;
    a >>= 1;
  }
  *s = b;
  Icheck(r);
  return r;
}

void setbit(Integer& x, long b)
{
  if (b >= 0)
  {
    int bw = (unsigned long)b / I_SHIFT;
    int sw = (unsigned long)b % I_SHIFT;
    if (x.rep == 0)
      x.rep = Icalloc(0, bw + 1);
    else if (x.rep->len < bw)
    {
      int xl = x.rep->len;
      x.rep = Iresize(x.rep, calc_len(xl, bw+1, 0));
    }
    x.rep->s[bw] |= (1 << sw);
    Icheck(x.rep);
  }
}

void clearbit(Integer& x, long b)
{
  if (b >= 0)
  {
    int bw = (unsigned long)b / I_SHIFT;
    int sw = (unsigned long)b % I_SHIFT;
    if (x.rep == 0)
      x.rep = Icalloc(0, bw + 1);
    else if (x.rep->len < bw)
    {
      int xl = x.rep->len;
      x.rep = Iresize(x.rep, calc_len(xl, bw+1, 0));
    }
    x.rep->s[bw] &= ~(1 << sw);
    Icheck(x.rep);
  }
}

int testbit(const Integer& x, long b)
{
  if (x.rep != 0 && b >= 0)
  {
    int bw = (unsigned long)b / I_SHIFT;
    int sw = (unsigned long)b % I_SHIFT;
    return (bw < x.rep->len && (x.rep->s[bw] & (1 << sw)) != 0);
  }
  else
    return 0;
}

// A  version of knuth's algorithm B / ex. 4.5.3.34
// A better version that doesn't bother shifting all of `t' forthcoming

IntRep* gcd(const IntRep* x, const IntRep* y)
{
  nonnil(x);
  nonnil(y);
  int ul = x->len;
  int vl = y->len;
  
  if (vl == 0)
    return Ialloc(0, x->s, ul, I_POSITIVE, ul);
  else if (ul == 0)
    return Ialloc(0, y->s, vl, I_POSITIVE, vl);

  IntRep* u = Ialloc(0, x->s, ul, I_POSITIVE, ul);
  IntRep* v = Ialloc(0, y->s, vl, I_POSITIVE, vl);

// find shift so that both not even

  long k = 0;
  int l = (ul <= vl)? ul : vl;
  int cont = 1;
  for (int i = 0;  i < l && cont; ++i)
  {
    unsigned long a =  (i < ul)? u->s[i] : 0;
    unsigned long b =  (i < vl)? v->s[i] : 0;
    for (int j = 0; j < I_SHIFT; ++j)
    {
      if ((a | b) & 1)
      {
        cont = 0;
        break;
      }
      else
      {
        ++k;
        a >>= 1;
        b >>= 1;
      }
    }
  }

  if (k != 0)
  {
    u = lshift(u, -k, u);
    v = lshift(v, -k, v);
  }

  IntRep* t;
  if (u->s[0] & 01)
    t = Ialloc(0, v->s, v->len, !v->sgn, v->len);
  else
    t = Ialloc(0, u->s, u->len, u->sgn, u->len);

  while (t->len != 0)
  {
    long s = 0;                 // shift t until odd
    cont = 1;
    int tl = t->len;
    for (i = 0; i < tl && cont; ++i)
    {
      unsigned long a = t->s[i];
      for (int j = 0; j < I_SHIFT; ++j)
      {
        if (a & 1)
        {
          cont = 0;
          break;
        }
        else
        {
          ++s;
          a >>= 1;
        }
      }
    }

    if (s != 0) t = lshift(t, -s, t);

    if (t->sgn == I_POSITIVE)
    {
      u = Icopy(u, t);
      t = add(t, 0, v, 1, t);
    }
    else
    {
      v = Ialloc(v, t->s, t->len, !t->sgn, t->len);
      t = add(t, 0, u, 0, t);
    }
  }
  delete t;
  delete v;
  if (k != 0) u = lshift(u, k, u);
  return u;
}



long lg(const IntRep* x)
{
  nonnil(x);
  int xl = x->len;
  if (xl == 0)
    return 0;

  long l = (xl - 1) * I_SHIFT - 1;
  unsigned short a = x->s[xl-1];

  while (a != 0)
  {
    a = a >> 1;
    ++l;
  }
  return l;
}
  
IntRep* power(const IntRep* x, long y, IntRep* r)
{
  nonnil(x);
  int sgn;
  if (x->sgn == I_POSITIVE || (!(y & 1)))
    sgn = I_POSITIVE;
  else
    sgn = I_NEGATIVE;

  int xl = x->len;

  if (y == 0 || (xl == 1 && x->s[0] == 1))
    r = Icopy_one(r, sgn);
  else if (xl == 0 || y < 0)
    r = Icopy_zero(r);
  else if (y == 1 || y == -1)
    r = Icopy(r, x);
  else
  {
    int maxsize = ((lg(x) + 1) * y) / I_SHIFT + 2;     // pre-allocate space
    IntRep* b = Ialloc(0, x->s, xl, I_POSITIVE, maxsize);
    b->len = xl;
    r = Icalloc(r, maxsize);
    r = Icopy_one(r, I_POSITIVE);
    for(;;)
    {
      if (y & 1)
        r = multiply(r, b, r);
      if ((y >>= 1) == 0)
        break;
      else
        b = multiply(b, b, b);
    }
    delete b;
  }
  r->sgn = sgn;
  Icheck(r);
  return r;
}

IntRep* abs(const IntRep* src, IntRep* dest)
{
  nonnil(src);
  if (src != dest)
    dest = Icopy(dest, src);
  dest->sgn = I_POSITIVE;
  return dest;
}

IntRep* negate(const IntRep* src, IntRep* dest)
{
  nonnil(src);
  if (src != dest)
    dest = Icopy(dest, src);
  if (dest->len != 0) 
    dest->sgn = !dest->sgn;
  return dest;
}

#if defined(__GNUG__) && !defined(NO_NRV)

Integer sqrt(const Integer& x) return r(x)
{
  int s = sign(x);
  if (s < 0) x.error("Attempted square root of negative Integer");
  if (s != 0)
  {
    r >>= (lg(x) / 2); // get close
    Integer q;
    div(x, r, q);
    while (q < r)
    {
      r += q;
      r >>= 1;
      div(x, r, q);
    }
  }
  return;
}

Integer lcm(const Integer& x, const Integer& y) return r
{
  if (!x.initialized() || !y.initialized())
    x.error("operation on uninitialized Integer");
  Integer g;
  if (sign(x) == 0 || sign(y) == 0)
    g = 1;
  else 
    g = gcd(x, y);
  div(x, g, r);
  mul(r, y, r);
}

#else 
Integer sqrt(const Integer& x) 
{
  Integer r(x);
  int s = sign(x);
  if (s < 0) x.error("Attempted square root of negative Integer");
  if (s != 0)
  {
    r >>= (lg(x) / 2); // get close
    Integer q;
    div(x, r, q);
    while (q < r)
    {
      r += q;
      r >>= 1;
      div(x, r, q);
    }
  }
  return r;
}

Integer lcm(const Integer& x, const Integer& y) 
{
  Integer r;
  if (!x.initialized() || !y.initialized())
    x.error("operation on uninitialized Integer");
  Integer g;
  if (sign(x) == 0 || sign(y) == 0)
    g = 1;
  else 
    g = gcd(x, y);
  div(x, g, r);
  mul(r, y, r);
  return r;
}

#endif



IntRep* atoIntRep(const char* s, int base)
{
  int sl = strlen(s);
  IntRep* r = Icalloc(0, sl * (lg(base) + 1) / I_SHIFT + 1);
  if (s != 0)
  {
    char sgn;
    while (isspace(*s)) ++s;
    if (*s == '-')
    {
      sgn = I_NEGATIVE;
      s++;
    }
    else if (*s == '+')
    {
      sgn = I_POSITIVE;
      s++;
    }
    else
      sgn = I_POSITIVE;
    for (;;)
    {
      long digit;
      if (*s >= '0' && *s <= '9') digit = *s - '0';
      else if (*s >= 'a' && *s <= 'z') digit = *s - 'a' + 10;
      else if (*s >= 'A' && *s <= 'Z') digit = *s - 'A' + 10;
      else break;
      if (digit >= base) break;
      r = multiply(r, base, r);
      r = add(r, 0, digit, r);
      ++s;
    }
    r->sgn = sgn;
  }
  return r;
}



extern AllocRing _libgxx_fmtq;

char* Itoa(const IntRep* x, int base, int width)
{
  int wrksiz = (x->len + 1) * I_SHIFT / lg(base) + 2 + width;
  char* fmtbase = (char *) _libgxx_fmtq.alloc(wrksiz);
  char* e = fmtbase + wrksiz - 1;
  char* s = e;
  *--s = 0;

  if (x->len == 0)
    *--s = '0';
  else
  {
    IntRep* z = Icopy(0, x);

    // split division by base into two parts: 
    // first divide by biggest power of base that fits in an unsigned short,
    // then use straight signed div/mods from there. 

    // find power
    int bpower = 1;
    unsigned short b = base;
    unsigned short maxb = I_MAXNUM / base;
    while (b < maxb)
    {
      b *= base;
      ++bpower;
    }
    for(;;)
    {
      int rem = unscale(z->s, z->len, b, z->s);
      Icheck(z);
      if (z->len == 0)
      {
        while (rem != 0)
        {
          char ch = rem % base;
          rem /= base;
          if (ch >= 10)
            ch += 'a' - 10;
          else
            ch += '0';
          *--s = ch;
        }
        delete z;
        break;
      }
      else
      {
        for (int i = 0; i < bpower; ++i)
        {
          char ch = rem % base;
          rem /= base;
          if (ch >= 10)
            ch += 'a' - 10;
          else
            ch += '0';
          *--s = ch;
        }
      }
    }
  }

  if (x->sgn == I_NEGATIVE) *--s = '-';
  int w = e - s - 1;
  while (w++ < width)
    *--s = ' ';
  return s;
}

char* dec(const Integer& x, int width)
{
  return Itoa(x, 10, width);
}

char* oct(const Integer& x, int width)
{
  return Itoa(x, 8, width);
}

char* hex(const Integer& x, int width)
{
  return Itoa(x, 16, width);
}

istream& operator >> (istream& s, Integer& y)
{
//#ifndef VMS		// is this really needed?  I think not.
  if (!s.readable())
  {
    s.set(_bad);
    return s;
  }
//#endif

  int got_one = 0;
  char sgn = 0;
  char ch;
  y.rep = Icopy_zero(y.rep);
  s >> WS;
  if (!s.good())
  {
    s.set(_bad);
    return s;
  }
  while (s.get(ch))
  {
    if (ch == '-')
    {
      if (sgn == 0)
        sgn = '-';
      else
        break;
    }
    else if (ch >= '0' && ch <= '9')
    {
      long digit = ch - '0';
      y *= 10;
      y += digit;
      got_one = 1;
    }
    else
      break;
  }
  if (s.good())
    s.unget(ch);
  if (!got_one)
    s.set(_bad);

  if (sgn == '-')
    y.negate();

  return s;
}

int Integer::OK() const
{
  int v = rep != 0;             // have a rep
  int l = rep->len;
  int s = rep->sgn;
  v &= l <= rep->sz;            // length with bounds
  v &= s == 0 || s == 1;        // legal sign
  Icheck(rep);                  // and correctly adjusted
  v &= rep->len == l;
  v &= rep->sgn == s;
  if (!v) error("invariant failure");
  return v;
}

volatile void Integer::error(const char* msg) const
{
  (*lib_error_handler)("Integer", msg);
}