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

This file is part of the GNU C++ Library.  This library is free
software; you can redistribute it and/or modify it under the terms of
the GNU Library General Public License as published by the Free
Software Foundation; either version 2 of the License, or (at your
option) any later version.  This library is distributed in the hope
that it will be useful, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE.  See the GNU Library General Public License for more details.
You should have received a copy of the GNU Library General Public
License along with this library; if not, write to the Free Software
Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/

/*
  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
#endif
#ifdef DEBUG
#include <assert.h>
#endif
#include <math.h>
#include "Integer.h"
// #include <std.h>
#include <ctype.h>
#include <limits.h>
#include <string.h>
#include <iostream.h>
// Chen Li:
// #include <Obstack.h>
// #include <AllocRing.h>
#include <new.h>
#include <assert.h>
// #include <builtin.h>
#include "Integer.hP"

#undef OK
#define M_fft 2013265921
#define W_fft 440564289

IntRep _ZeroRep = {1, 0, 1, {0}};
IntRep _OneRep = {1, 0, 1, {1}};
IntRep _MinusOneRep = {1, 0, 0, {1}};

void my_error_handler(const char *className, const char *mesg) {
  cerr << className << " error: " << mesg << endl;
  abort();
}


// utilities to extract and transfer bits 

// get low bits

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

// transfer high bits to low

inline static unsigned long down(unsigned long x)
{
  // Chen Li: why does it use '&' operator in the expression below???
  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

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) 
    my_error_handler("Integer", "operation on uninitialized Integer");
}

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

int num_IntRep_new = 0;

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))
    my_error_handler("Integer", "Requested length out of range");
    
  num_IntRep_new++ ;
  IntRep* rep = new (operator new (allocsiz)) IntRep;
  rep->sz = ((long)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 && !STATIC_IntRep(old)) delete old;
  return rep;
}

// allocate and clear

IntRep* Icalloc(IntRep* old, int newlen)
{
  IntRep* rep;
  if (old == 0 || newlen > old->sz)
  {
    if (old != 0 && !STATIC_IntRep(old)) 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;
      if (!STATIC_IntRep(old)) 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 && !STATIC_IntRep(old)) 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)
{
  int newsgn = (x >= 0);
  IntRep* rep = Icopy_ulong(old, newsgn ? x : -x);
  rep->sgn = newsgn;
  return rep;
}

IntRep* Icopy_ulong(IntRep* old, unsigned long x)
{
  unsigned short src[SHORT_PER_LONG];
  
  unsigned short srclen = 0;
  while (x != 0)
  {
    src[srclen++] = extract(x);
    x = down(x);
  }

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

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

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

  return rep;
}

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

IntRep* Icopy_zero(IntRep* old)
{
  if (old == 0 || STATIC_IntRep(old))
    return &_ZeroRep;

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

  return old;
}

// special case for 1 or -1

IntRep* Icopy_one(IntRep* old, int newsgn)
{
  if (old == 0 || 1 > old->sz)
  {
    if (old != 0 && !STATIC_IntRep(old)) delete old;
    return newsgn==I_NEGATIVE ? &_MinusOneRep : &_OneRep;
  }

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

  return old;
}

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

long Itolong(const IntRep* rep)
{ 
  if ((unsigned)(rep->len) > (unsigned)(SHORT_PER_LONG))
    return (rep->sgn == I_POSITIVE) ? LONG_MAX : LONG_MIN;
  else if (rep->len == 0)
    return 0;
  else if ((unsigned)(rep->len) < (unsigned)(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)? static_cast<long>(a) : 
      -(static_cast<long>(a));
  }
  else 
  {
    unsigned long a = rep->s[SHORT_PER_LONG - 1];
    if (a >= I_MINNUM)
      return (rep->sgn == I_POSITIVE) ? LONG_MAX : LONG_MIN;
    else
    {
      a = up(a) | rep->s[SHORT_PER_LONG - 2];
      if (SHORT_PER_LONG > 2)
      {
        for (int i = static_cast<int>(SHORT_PER_LONG) - 3; i >= 0; --i)
          a = up(a) | rep->s[i];
      }
      return (rep->sgn == I_POSITIVE)? static_cast<long>(a)
	: -(static_cast<long>(a));
    }
  }
}

// test whether op long() will work.
// careful about asymmetry between LONG_MIN & LONG_MAX

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

// 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;
}

// There is an upper limit on bit size for a value in 
// the dft calc with the number transformed into a L-bit 
// vector.  The NumBits is passed to the function to determine L.
int getLfromNumBits(int NumBits)
{
  // The array values are L*K. 
  // For 32-bit numbers, the array runs from L=2 to L=13 (neither 1 or 14 or useful).
  
  // For each L, K is the max power of 2 such that K <= M/2^(2L + 1) )
  // Another constraint is that 2K must be <= 2^27 because, briefly put, M=15*2^27+1 
 
  // Generalizing: 
  // M = 1 + m * 2^x (M is a prime less than 2^N with N being the # of bits) 
  // w = (N-1)^m (mod M) 
  //    Explanation: w^i (mod M) = (N-1)^(m*i) !==1 for i= 1, 2, .. 2^x -1, 
  //    but w^(2^x) (mod M) = (N-1)^(M-1) (mod M) = 1 by Fermat's Little Theorem,
  //    so w is a 2^x-th primitive root of unity.
    
  static const int LBounds[] = 
  { 2 * (1 << 25),
    3 * (1 << 23),
    4 * (1 << 21),
    5 * (1 << 19),
    6 * (1 << 17),
    7 * (1 << 15),
    8 * (1 << 13),
    9 * (1 << 11),
   10 * (1 <<  9),
   11 * (1 <<  7),
   12 * (1 <<  5),
   13 * (1 <<  3)
  };
   
  static const int c = sizeof(LBounds)/sizeof(int);
  //cout << "Here is c: " << c << endl;

  for (int L=c; L> 0; L--) 
     if (NumBits <= LBounds[L-1])
     {  
        //co << "NumBits = " << NumBits << " can be handled by L = " << L+2 << "bits" << endl;
        return L+1;
     }
  // cout << "Returning L = 0 ?????? " << endl; 
  return 0;
}

// this is the recursive part of the dft algorithm
// which splits the values to process into even and odd
// by multiplying a factor and passing two "remainders" which
// when you get down to the final iteration represent a unique element
// in the A array 
int DFTrecur(const int *A     , int Xj          , int vFactor, int vRemainder, 
             int XoddFactor   , int K      , 
             int zeroHalfFlag , int *dft)
{
  //cout << "In DFT Recur with " << vFactor *2 << " vs. " << K << endl;  

  if (vFactor *2 > K)
  {
    // then set this element and return
    dft[vRemainder] = ( A[vRemainder] * Xj * XoddFactor) % M_fft; 
    return 1;
  }
 
  int x2j = (Xj * Xj) % M_fft;

  // recursive case
  // Pe(X^2j) = 

  DFTrecur(A, x2j,     vFactor*2, vRemainder*2, XoddFactor, 
           K, zeroHalfFlag, dft);  

  //Optimization:  This prevents the program from
  // processing the zero-padded top half K sized vector
  // from being processed 
  // Also when the vector is summed, the 2nd half will be skipped
  if (!zeroHalfFlag)
  {
     // X^j*Po(X^2j) = 
     DFTrecur(A, x2j,     vFactor*2, vRemainder*2+1, (XoddFactor * Xj)% M_fft, 
	      K, 0 /* blank out zero flag */, dft);
  }
}

// this calculates the polynomials for DFT by repeatedly calling DFTrecur
// and summing the elements for A, to get a unique element
int DFT(const int *A, int K, int *scratch, int inverse, int zeroHalfFlag, int *dft)
{
  for (int WpowImodM = 1, i = 0; 
       i < K; 
       i++, WpowImodM = (WpowImodM * WpowImodM) % M_fft)
  {    
    //should I blank scratch subseq. times
    //cout << "**DFT Loop: " << i << " of " << K << endl;
    // the sum of the result from one DFT calculation = Ai
    DFTrecur(A,        // L bit vector
            WpowImodM, // w^i
            1,         // array index (multiple of two)
            0,         // remainder
            1,         // cumul. "odd" factor
            K,
            zeroHalfFlag, // true if 1st 1/2 of vector zero-padded
            scratch);

    //cout << "**DFT Loop after DFTrecur: " << i << " of " << K << endl;
    
    int *scr_ptr = scratch;
    
    int j, dftsum, upper_bound = (zeroHalfFlag ? K/2 : K); 
    
    //cout << "**DFT Loop Upperbound: " << upper_bound << endl;
    for (dftsum = 0, j = 0; j < upper_bound; scr_ptr++, j++) 
    {   
      dftsum = (dftsum + *scr_ptr) % M_fft;
      // iterate thru vector keep mod running total
      //cout << "**DFT Loop after DFTrecur: " << j << " of " << upper_bound << endl;
    }
    if (inverse)      
      dft[(K-i)% K] = dftsum / K; 
    else
      dft[i] = dftsum;
  }

  return 1;
}

// called by "convertToLBitVector" this finds the value of the bits
// between start and end
int getValueNextLBits(const IntRep *x, int StartBit, int EndBit)
{  
   static const int iS = sizeof(x->s[0])*8; //assumes a short int
   int StartDiv = StartBit / iS;
   int StartMod = StartBit % iS;
   int EndDiv = EndBit / iS;
   int EndMod = EndBit % iS;
   
   //cout << endl << "StartDiv, StartMod, EndDiv, EndMod, iS, xStart, xEnd: " << endl;
   //cout << StartDiv << ", " << StartMod << ", " << EndDiv << ", " << EndMod;
   //cout << ", " << iS << ", " << x->s[StartDiv] << ", " << x->s[EndDiv] <<endl;

   // if div is equal, easier case
   if (StartDiv == EndDiv)
   {  // knock it right, StartMod bits then, get the next End-Start bits via 
      // mod. 
      return ( (x->s[StartDiv] >> StartMod) % (1 << (EndMod-StartMod)) ); 
   } 
   // otherwise break into two pieces
   else
   {  // similar to above, except split between two elements in the array
      return ( (x->s[StartDiv] >> StartMod) 
               +
               ( (x->s[EndDiv] % (1 << EndMod))
                        << (iS - StartMod) ) 
             ); 
   }
}

// converts an IntRep struct into an array of numbers with
// values of at most 2^L-1 so it can be used for the FFT calculation 
int convertToLBitVector(const IntRep *x, int L, int K, int *A)
{
   int NumBits = x->len * sizeof(x->s[0]) * 8;
   int StartBit, EndBit;
   cout << "--->NumBits, L, K: " << NumBits << ", " << L << ", " << K << endl; 
   for (int I = 0, StartBit = 0, EndBit = L; 
        StartBit < NumBits; 
        I++, StartBit += L, EndBit += L)
  {
     if (EndBit > NumBits)
        EndBit = NumBits;
     A[I] = getValueNextLBits(x,StartBit,EndBit);
     //cout << "The next L Bits are " << A[I] << " from " << StartBit << " to " << EndBit-1 << endl;
  } 
}

// TO DO: make the calculation mod M
int accumWithCarry(int N, int *Z)
{
  static const int pow2 = 1 << (8 * sizeof(int)-1); // 2^(8*sizeof(int)-1)

  // test if a carry is needed (sum/2 > 2^(8*sizeof(int)-1))
  if ( (N/2 + Z[0]/2) > pow2)
  {
    Z[0] += (N/2 + Z[0]/2 - pow2)*2 + (N % 2) + (Z[0] % 2);  
    accumWithCarry(1, Z+1); // just in case
  }
  else
  {
    Z[0] += N;
  }
  return 1;
}

int accumWithCarryShort(int N, unsigned short *Z)
{
  static const int pow2 = 1 << (8*sizeof(*Z)-1); // 2^(sizeof(short)-1)
 
  cout << "Here we are N = " << N << ", Z = " << *Z << endl;
  // test if a carry is needed (sum/2 > 2^(8*sizeof(short)-1))
  if ( (N/2 + Z[0]/2) > pow2)
  {
    Z[0] += (N/2 + Z[0]/2 - pow2)*2 + (N % 2) + (Z[0] % 2);
    accumWithCarryShort(1, Z+1); // just in case
  }
  else
  {
    Z[0] += N;
  }
  return 1;
}

// this function does component-wise multiplication of
// two vectors of I-bit integers mod M.
int multiply(const int *X, const int *Y, int K, int *v)
{
  static const int MMod = 1 + (~0) - M_fft*(int)(2 * (double)( (~0+1)/2/M_fft ) ); // M mod (2^sizeof(int));  
  static const int Base = 1 << (sizeof(int)/2); // Base=2^(sizeof(int)/2) 
  int temp, xT, yT, xD, xM, yD, yM, *z=v, facI, facIHalf;
  const int *x = X, *y = Y;
  int fac, iModCarry, iMod;

  for (int j = 0; j < K; j++)
  {
    cout << "Before j : " << j << ", X : " << X[j] << ", Y : " << Y[j] << endl;  
  }
 
  //blank out z before starting
  for (int b = 0; b < K; b++, z++)
    *z = 0;
  z = v;

  // what if I multiply this and it exceeds the size
  for (int i=0; i < K; i++, x++, y++, z++)
  {
    //do x * y as (xD*Base + xM)*(yD*Base + yM)
    // xD*yD*Base^2 + (xD*yM+xM*yD)*Base + xM*yM 
    xT = x[0] % M_fft; // take mod just to be safe
    yT = y[0] % M_fft;
 
    xD = xT / Base;
    xM = xT % Base;
    yD = yT / Base;
    yM = yT % Base;

    facI = xD * yD % M_fft; // to be mult by Base^2
    facIHalf = xD*yM + xM*yD % M_fft; // to be mult by Base
    fac = xM * yM % M_fft;
    
    int temp = ((facIHalf % Base) * Base) % M_fft;
    *z = fac + temp; 
    temp = (facIHalf/Base * MMod ) % M_fft;
    temp = (temp + ((facI*MMod) % M_fft) ) % M_fft; 
    *z = (*z + temp) % M_fft; 
     
    /* Old Method
    accumWithCarry( fac,                       z);   
    accumWithCarry( (facIHalf % Base) * Base,  z);
    accumWithCarry( facIHalf / Base,         z+1);
    accumWithCarry( facI,                    z+1);*/
  }
  
  for (int j = 0; j < K; j++)
  {
    cout << "j : " << j << ", X : " << X[j] << ", Y : " << Y[j] << ", V : " << v[j] << endl;  
  }
}

// this will convert the array of I-bit integers back
// to the IntRep* form.
int convertIntVecToIntRep(const int *x, int L, int K, IntRep *r)
{
  // add the numbers, but shift each by L-bits
  static const int size = sizeof(short)*8;
  int offset, remain, xval;
  unsigned short *R = r->s;
  const int *X = x;
  int rc = 0;

  cout << "The length of r is " << r->len << endl;
   
  for (int i = 0, rBit = 0; i < K; i++, rBit +=L, X++)
  {
     xval = *X; 
     cout << "Converting at " << i << " of " << K << ", with rBit % size: " << rBit % size << endl;
     if (!(rBit % size) )
     {      
       accumWithCarryShort(xval, R); // No splitting
     }
     else
     {
        offset = size - 1 - (rBit % size);
        remain = size - offset;
        cout << "Here are offset = " << offset << ", remain = " << remain << endl; 
        accumWithCarryShort( (1 << offset) * (xval % (1 << remain)), R);
        accumWithCarryShort(                 (xval / (1 << remain)), R+1); 
     }
     cout << "R = " << *R << ", at i =  " << i << ", rc = " << rc << endl;

     if ( ((rBit + L) % size) < L)
     {
        R++; // increment the representaton's pointer
        rc++;
     }
  }
  cout << endl << endl << "Final Answer:" << endl;
  for (int i = 0; i < r->len; i++)
  {
    cout << "Block " << i << ": " << r->s[i] << endl;
  }
}

// this multiplies two large numbers via the FFT method.
// It is called by the multiply method if a certain size threshold
// is exceeded.
IntRep *multiplyFFT(const IntRep *x, const IntRep* y, IntRep *r)
{
  cout << "Calling multiplyFFT" << endl;
  // the numbers have been validated by "multiply" already
  int L, K, i, NumBits, Ai;
  int xysame = x == y;
  int *A, *B, *temp, *dftA, *dftB, *iDft, *prodDftADftB; 
  static const int size = 8 * sizeof(x->s[0]);

  //determine L and K
  if (xysame)
     NumBits = size * x->len;
  else 
  {
     if (x->len > y->len)
       NumBits = size * x->len;
     else
       NumBits = size * y->len;
  }

  int rl = x->len + y->len;
  if ((x->len == r->len) || (y->len == r->len))
      r = Iresize(r, rl);
    else
      r = Icalloc(r, rl);

  L = getLfromNumBits(NumBits);
  K = int(NumBits/L + .5);
  cout << "L = " << L << ", K = " << K << endl;

  temp = new int[2*K+1]; // reusable array  
  A = new int[2*K+1];

  for (Ai = 0; Ai < 2*K +1; Ai++)
  {
     A[Ai] = 0;
     temp[Ai] = 0;
  }

  convertToLBitVector(x,L,2*K,A); //convert to L-bit vector of length 2K 

  dftA = new int[2*K+1];
  for (Ai = 0; Ai < 2*K +1; Ai++)
  {
     dftA[Ai] = 0;
  }

  DFT(A,2*K,temp,0,1,dftA); // normal DFT with "zero" flag

  cout << "Done with DFT" << endl;

  if (xysame)
     // don't need to do DFT calculation twice
     dftB = dftA;
  else  
  {
     B = new int[2*K+1];
     convertToLBitVector(y,L,2*K,B); //convert to L-bit vector of length 2K
     dftB = new int[2*K+1]; 
     for (Ai = 0; Ai < 2*K +1; Ai++)
     {
        dftB[Ai] = 0;
     }
     DFT(B,2*K,temp,0,1,dftB); // normal DFT with "zero" flag passed
  }
 
  //calculate the component-wise product of the two dfts
  prodDftADftB = new int[2*K+1];
  for (Ai = 0; Ai < 2*K +1; Ai++)
  {
     prodDftADftB[Ai] = 0;
  }

  multiply(dftA, dftB, 2*K, prodDftADftB);
  cout << "*** Done multiplying the two" << endl;   
  //return the InverseDFT of the product as a vector of integers
  iDft = new int[2*K+1];
  for (Ai = 0; Ai < 2*K +1; Ai++)
  {
     iDft[Ai] = 0;
  }

  DFT(prodDftADftB,2*K,temp,1,0,iDft); // inverseDFT with no zero flag

  cout << "*** Done with inverse DFT of product" << endl;   
  convertIntVecToIntRep(iDft, L, 2*K, r);
  cout << "*** Done with inverse DFT of product" << endl;

  // free the memory
  delete [] A;
  delete [] dftA;

  if (!xysame)
  {
    delete [] B;
    delete [] dftB;
  }

  delete [] temp;

  delete [] prodDftADftB;
  delete [] iDft;
}


IntRep* multiply(const IntRep* x, const IntRep* y, IntRep* r)
{
  //FS have some lower bound: if it's bigger than that, use the DFT
  static const int CORE_MULT_THRESHOLD=100; //FS
  //FS 
  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 (rl > CORE_MULT_THRESHOLD) //FS 
  //    return (multiplyFFT(x, y, r)); //FS 
  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 = static_cast<unsigned short>(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) my_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 = static_cast<unsigned short>((I_RADIX / (1 + y->s[yl - 1])));
    if (prescale != 1 || y == q)
    {
      yy = multiply(y, static_cast<long>((prescale & I_MAXNUM)), yy);
      r = multiply(x, static_cast<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 && !STATIC_IntRep(yy)) delete yy;
    if (!STATIC_IntRep(r)) 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) my_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 = static_cast<unsigned short>((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, static_cast<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 (!STATIC_IntRep(r)) 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) my_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 = static_cast<unsigned short>((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, static_cast<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);
    if (!STATIC_IntRep(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)
    my_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 = static_cast<unsigned short>((I_RADIX / (1 + y->s[yl - 1])));
    if (prescale != 1 || y == q || y == r)
    {
      yy = multiply(y, static_cast<long>((prescale & I_MAXNUM)), yy);
      r = multiply(x, static_cast<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 && !STATIC_IntRep(yy)) 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) my_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 = static_cast<unsigned short>((I_RADIX / (1 + y->s[yl - 1])));
    if (prescale != 1 || y == r)
    {
      yy = multiply(y, static_cast<long>((prescale & I_MAXNUM)), yy);
      r = multiply(x, static_cast<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 && !STATIC_IntRep(yy)) 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) my_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 = static_cast<unsigned short>((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, static_cast<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 = static_cast<int>(ay / I_SHIFT);
  int sw = static_cast<int>(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))
  if ( 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 = static_cast<int>((unsigned long)b / I_SHIFT);
    int sw = static_cast<int>((unsigned long)b % I_SHIFT);
    int xl = x.rep ? x.rep->len : 0;
    if (xl <= bw)
      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)
    {
      if (x.rep == 0)
	x.rep = &_ZeroRep;
      else
	{
	  int bw = static_cast<int>((unsigned long)b / I_SHIFT);
	  int sw = static_cast<int>((unsigned long)b % I_SHIFT);
	  if (x.rep->len > bw)
	    x.rep->s[bw] &= ~(1 << sw);
	}
    Icheck(x.rep);
  }
}

int testbit(const Integer& x, long b)
{
  if (x.rep != 0 && b >= 0)
  {
    int bw = static_cast<int>((unsigned long)b / I_SHIFT);
    int sw = static_cast<int>((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 (int 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);
    }
  }
  if (!STATIC_IntRep(t)) delete t;
  if (!STATIC_IntRep(v)) delete v;
  if (k != 0) u = lshift(u, k, u);
  return u;
}



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

  int l = (xl - 1) * I_SHIFT - 1; // this number can be negative
  unsigned short a = x->s[xl-1];

  while (a != 0)
  {
    a = a >> 1;
    ++l;
  }
#ifdef DEBUG
  assert(l >= 0);
#endif
  return l;
}

unsigned long ceilLog(const IntRep* x) {
  // chen li 07/99
  nonnil(x);
  unsigned short xl = x->len;
  if (xl == 0) 
    return 0;

  int low_bits = (xl - 1) * I_SHIFT - 1; // can be negative.
  unsigned long ceil_lg = 0;
  unsigned short a = x->s[xl-1];

  while (a != 0)
  {
    a = a >> 1;
    ++ceil_lg;
  }

  if ((ceil_lg > 0) && 
      ((1 << (ceil_lg - 1)) != a)) {
    ceil_lg++;
  } else {
    for (unsigned int k = 0; k < xl-1; k++)
      if (x->s[k] != 0) {
	ceil_lg++;
	break;
      }
  }
#ifdef DEBUG
  assert((low_bits >= 0) || (- low_bits <= ceil_lg)); // check
#endif
  return (ceil_lg + low_bits); // shoud alwasy be non-negative!
}
  
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);
    }
    if (!STATIC_IntRep(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(_G_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, static_cast<int>(sl * (log(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;
}


// Chen Li: AllocRing is not supported by latest g++ anymore.
// extern AllocRing _libgxx_fmtq;

char* Itoa(const IntRep* x, int base, int width)
{
  int fmtlen = static_cast<int>((x->len + 1) * I_SHIFT / log(base) + 4 + width);
  //  char* fmtbase = (char *) _libgxx_fmtq.alloc(fmtlen);
  char *fmtbase = new char[fmtlen];
  char* f = cvtItoa(x, fmtbase, fmtlen, base, 0, width, 0, ' ', 'X', 0);
  return f;
}

ostream& operator << (ostream& s, const Integer& y)
{
#ifdef _OLD_STREAMS
  return s << Itoa(y.rep);
#else
  // Chen Li: modified it to fit in new Sun's CC compiler
  // change:
  //  if (s.opfx())
  // to:
  if (s.good())
    {
      int base = (s.flags() & ios::oct) ? 8 : (s.flags() & ios::hex) ? 16 : 10;
      int width = s.width();
      y.printon(s, base, width);
    }
  return s;
#endif
}

void Integer::printon(ostream& s, int base /* =10 */, int width /* =0 */) const
{
  int align_right = !(s.flags() & ios::left);
  int showpos = s.flags() & ios::showpos;
  int showbase = s.flags() & ios::showbase;
  char fillchar = s.fill();
  char Xcase = (s.flags() & ios::uppercase)? 'X' : 'x';
  const IntRep* x = rep;
  int fmtlen = static_cast<int>((x->len + 1) * I_SHIFT / log(base) + 4 + width);
  char* fmtbase = new char[fmtlen];
  char* f = cvtItoa(x, fmtbase, fmtlen, base, showbase, width, align_right,
		    fillchar, Xcase, showpos);
  s.write(f, fmtlen);
  delete [] fmtbase;
}

char*  cvtItoa(const IntRep* x, char* fmt, int& fmtlen, int base, int showbase,
              int width, int align_right, char fillchar, char Xcase, 
              int showpos)
{
  char* e = fmt + fmtlen - 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;
        }
	if (!STATIC_IntRep(z)) 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 (base == 8 && showbase) 
    *--s = '0';
  else if (base == 16 && showbase) 
  { 
    *--s = Xcase; 
    *--s = '0'; 
  }
  if (x->sgn == I_NEGATIVE) *--s = '-';
  else if (showpos) *--s = '+';
  int w = e - s - 1;
  if (!align_right || w >= width)
  {
    while (w++ < width)
      *--s = fillchar;
    fmtlen = e - s - 1;
    return s;
  }
  else
  {
    char* p = fmt;
    for (char* t = s; *t != 0; ++t, ++p) *p = *t;
    while (w++ < width) *p++ = fillchar;
    *p = 0;
    fmtlen = p - fmt;
    return fmt;
  }
}

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& stream, Integer& val)
{
  // Chen Li: modified to fit in new Sun's CC compiler. change:
  //  if (!stream.ipfx(0))
  // to:
  if (!stream.good())
    return stream;
  int sign = ' ';
  register streambuf* sb = stream.rdbuf();
  int base = 10;
  int ndigits = 0;
  register int ch = sb->sbumpc();
  while (ch != EOF && isspace(ch))
    ch = sb->sbumpc();
  if (ch == '+' || ch == '-')
    {
      sign = ch;
      ch = sb->sbumpc();
      while (ch != EOF && isspace(ch))
	ch = sb->sbumpc();
    }
  if (ch == EOF) goto eof_fail;
  if (!(stream.flags() & ios::basefield))
    {
      if (ch == '0')
	{
	  ch = sb->sbumpc();
	  if (ch == EOF) { }
	  else if (ch == 'x' || ch == 'X')
	    {
	      base = 16;
	      ch = sb->sbumpc();
	      if (ch == EOF) goto eof_fail;
	    }
	  else
	    {
	      sb->sputbackc(ch);
	      base = 8;
	      ch = '0';
	    }
	}
    }
  else if ((stream.flags() & ios::basefield) == ios::hex)
    base = 16;
  else if ((stream.flags() & ios::basefield) == ios::oct)
    base = 8;

  val.rep = Icopy_zero(val.rep);

  for (;;)
    {
      if (ch == EOF)
	break;
      int digit;
      if (ch >= '0' && ch <= '9')
	digit = ch - '0';
      else if (ch >= 'A' && ch <= 'F')
	digit = ch - 'A' + 10;
      else if (ch >= 'a' && ch <= 'f')
	digit = ch - 'a' + 10;
      else
	digit = 999;
      if (digit >= base)
	{
	  sb->sputbackc(ch);
	  if (ndigits == 0)
	    goto fail;
	  else
	    goto done;
	}
      ndigits++;
      switch (base)
	{
	case 8:
	  val <<= 3;
	  break;
	case 16:
	  val <<= 4;
	  break;
	default:
	  val *= base;
	  break;
	}
      val += digit;
      ch = sb->sbumpc();
    }
 fail:
  // Chen Li: modify to fonform with Standard C++
  //  stream.set(ios::failbit);
  stream.clear(ios::failbit);
 done:
  if (sign == '-')
    val.negate();
  return stream;
 eof_fail:
  // Chen Li: modify to conform with standard C++
  //  stream.set(ios::failbit|ios::eofbit);
  stream.clear(ios::failbit|ios::eofbit);
  return stream;
} // operator >>


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

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

long lg(const Integer& x)
{
  return lg(x.rep);
}

unsigned long ceilLog(const Integer& x)
{
  return ceilLog(x.rep);
}

double ratio(const Integer& x, const Integer& y) {
  return (x.as_double() / y.as_double());
}