integer.cc

早期freebsd实现

    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);
    }
  }
  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);
  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);
    }
    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(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 fmtlen = (x->len + 1) * I_SHIFT / lg(base) + 4 + width;
  char* fmtbase = (char *) _libgxx_fmtq.alloc(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
  if (s.opfx())
    {
      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 = (x->len + 1) * I_SHIFT / lg(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;
    int gap = s - p;
    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& s, Integer& y)
{
#ifdef _OLD_STREAMS
  if (!s.good())
    return s;
#else
  if (!s.ipfx(0))
  {
    s.set(ios::failbit);
    return s;
  }
#endif
  s >> ws;
  if (!s.good())
  {
    s.set(_fail);
    return s;
  }
  
#ifdef _OLD_STREAMS
  int know_base = 0;
  int base = 10;
#else
  int know_base = s.flags() & (ios::oct | ios::hex | ios::dec);
  int base = (s.flags() & ios::oct) ? 8 : (s.flags() & ios::hex) ? 16 : 10;
#endif

  int got_one = 0;
  char sgn = 0;
  char ch;
  y.rep = Icopy_zero(y.rep);

  while (s.get(ch))
  {
    if (ch == '-')
    {
      if (sgn == 0)
        sgn = '-';
      else
        break;
    }
    else if (!know_base & !got_one && ch == '0')
      base = 8, got_one = 1;
    else if (!know_base & !got_one && base == 8 && (ch == 'X' || ch == 'x'))
      base = 16;
    else if (base == 8)
    {
      if (ch >= '0' && ch <= '7')
      {
        long digit = ch - '0';
        y <<= 3;
        y += digit;
        got_one = 1;
      }
      else
        break;
    }
    else if (base == 16)
    {
      long 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 = base;
      if (digit < base)
      {
        y <<= 4;
        y += digit;
        got_one = 1;
      }
      else
        break;
    }
    else if (base == 10)
    {
      if (ch >= '0' && ch <= '9')
      {
        long digit = ch - '0';
        y *= 10;
        y += digit;
        got_one = 1;
      }
      else
        break;
    }
    else
      abort(); // can't happen for now
  }
  if (s.good())
    s.putback(ch);
  if (!got_one)
    s.set(_fail);

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

  return s;
}

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
{
  (*lib_error_handler)("Integer", msg);
}