integer.cc

Gambit 是一个游戏库理论软件

//
// $Source: /home/gambit/CVS/gambit/sources/math/integer.cc,v $
// $Date: 2002/08/26 05:50:04 $
// $Revision: 1.7 $
//
// DESCRIPTION;
// Implementation of an arbitrary-length integer class
//
// This file is part of Gambit
// Modifications copyright (c) 2002, The Gambit Project
//
// The original copyright and license are included below.

/* 
Copyright (C) 1988 Free Software Foundation
    written by Doug Lea ([e-mail])

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, 675 Mass Ave, Cambridge, MA 02139, 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
#include "base/gstream.h"
#include "math/integer.h"
#include "gnulib.h"
#include <ctype.h>
#include <float.h>
#include <limits.h>
#include <math.h>
//#include <new.h>
#include <assert.h>

long lg(unsigned long x)
{
  long l = 0;
  while (x > 1)
  {
    x = x >> 1;
    ++l;
  }
  return l;
}


#ifndef HUGE_VAL
#ifdef HUGE
#define HUGE_VAL HUGE
#else
#define HUGE_VAL DBL_MAX
#endif
#endif

/*
 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         (sizeof(short) * CHAR_BIT)
#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  ((unsigned)(((sizeof(long) + sizeof(short) - 1) / sizeof(short))))
#define CHAR_PER_LONG   ((unsigned)sizeof(long))

/*
  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

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


// utilities to extract and transfer bits

// get low bits

inline static unsigned short extract(unsigned long x)
{
  return (unsigned short) (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

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

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

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

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)
{
  assert(rep != 0);
}

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

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;
  assert((unsigned long) allocsiz < MAXIntRep_SIZE * sizeof(short));
    
  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 && !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];
#ifndef __BCC55__
    // This condition is always false under BCC55
    if (SHORT_PER_LONG > 2) // normally optimized out
    {
		for (int i = rep->len - 2; i >= 0; --i)
		  a = up(a) | rep->s[i];
    }
#endif  // __BCC55__
    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) ? LONG_MAX : LONG_MIN;
    else
    {
      a = up(a) | rep->s[SHORT_PER_LONG - 2];
#ifndef __BCC55__
      // This condition is always false under BCC55
      if (SHORT_PER_LONG > 2)
      {
		  for (int i = SHORT_PER_LONG - 3; i >= 0; --i)
			 a = up(a) | rep->s[i];
      }
#endif  // __BCC55__
      return (rep->sgn == I_POSITIVE)? a : -((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;
}

// convert to a double 

double Itodouble(const IntRep* rep)
{ 
  double d = 0.0;
  double bound = DBL_MAX / 2.0;
  for (int i = rep->len - 1; i >= 0; --i)
  {
	 unsigned short a = (unsigned short) (I_RADIX >> 1);
	 while (a != 0)
    {
      if (d >= bound)
        return (rep->sgn == I_NEGATIVE) ? -HUGE_VAL : HUGE_VAL;
      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 = DBL_MAX / 2.0;
  for (int i = rep->len - 1; i >= 0; --i)
  {
	 unsigned short a = (unsigned short) (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 gInteger& num, const gInteger& den)
{
  gInteger q, r;
  divide(num, den, q, r);
  double d1 = q.as_double();
 
  if (d1 >= DBL_MAX || d1 <= -DBL_MAX || sign(r) == 0)
    return d1;
  else      // use as much precision as available for fractional part
  {
    double  d2 = 0.0;
    double  d3 = 0.0; 
    int cont = 1;
    for (int i = den.rep->len - 1; i >= 0 && cont; --i)
    {
		unsigned short a = (unsigned short) (I_RADIX >> 1);
      while (a != 0)
      {
        if (d2 + 1.0 == d2) // out of precision when we get here
        {
          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);