core/vnl/vnl_bignum.cxx

Go to the documentation of this file.

// This is core/vnl/vnl_bignum.cxx
#include "vnl_bignum.h"
//:
// \file

#include <vcl_cctype.h>    // Include character macros
#include <vcl_cstdlib.h>   // for vcl_atol
#include <vcl_cstring.h>   // for vcl_strlen
#include <vcl_cmath.h>     // for vcl_fmod
#include <vcl_algorithm.h> // for vcl_copy
#include <vcl_vector.h>
#include <vcl_cassert.h>
#include <vcl_iostream.h>
#include <vcl_limits.h>
#include <vnl/vnl_math.h> // for vnl_math_isfinite(double)

typedef unsigned short Counter;
typedef unsigned short Data;

//: Creates a zero vnl_bignum.

vnl_bignum::vnl_bignum ()
: count(0), sign(1), data(0)
{
}

//: Creates a vnl_bignum from a long integer.

vnl_bignum::vnl_bignum (long l)
: count(0), sign(1), data(0)
{
  if (l < 0) {                  // Get correct sign
    l = -l;                     // Get absolute value of l
    this->sign = -1;
  }
  Data buf[sizeof(l)];          // Temp buffer to store l in
  Counter i = 0;                // buffer index
  while (l) {                   // While more bits in l
    assert(i < sizeof(l));      // no more buffer space
    buf[i] = Data(l);           // Peel off lower order bits
    l >>= 16;   // Shift next bits into place
    ++i;
  }
  if (i > 0)
    this->data = new Data[this->count=i]; // Allocate permanent data

  while (i--)     // Save buffer into perm. data
    this->data[i] = buf[i];
}

//: Creates a vnl_bignum from an integer.

vnl_bignum::vnl_bignum (int l)
: count(0), sign(1), data(0)
{
  if (l < 0) {                 // Get correct sign
    l = -l;                     // Get absolute value of l
    this->sign = -1;
  }
  Data buf[sizeof(l)];          // Temp buffer to store l in
  Counter i = 0;                // buffer index
  while (l) {                   // While more bits in l
    assert(i < sizeof(l));      // no more buffer space
    buf[i] = Data(l);           // Peel off lower order bits
    l >>= 16;   // Shift next bits into place
    i++;
  }
  if (i > 0)
    this->data = new Data[this->count=i]; // Allocate permanent data

  while (i--)     // Save buffer into perm. data
    this->data[i] = buf[i];
}

//: Creates a vnl_bignum from an unsigned long integer.

vnl_bignum::vnl_bignum (unsigned long l)
: count(0), sign(1), data(0)
{
  Data buf[sizeof(l)];          // Temp buffer to store l in
  Counter i = 0;                // buffer index
  while (l) {                   // While more bits in l
    assert(i < sizeof(l));      // no more buffer space
    buf[i] = Data(l);           // Peel off lower order bits
    l >>= 16;   // Shift next bits into place
    i++;
  }
  if (i > 0)
    this->data = new Data[this->count=i]; // Allocate permanent data

  while (i--)     // Save buffer into perm. data
    this->data[i] = buf[i];
}

//: Creates a vnl_bignum from an unsigned integer.

vnl_bignum::vnl_bignum (unsigned int l)
: count(0), sign(1), data(0)
{
  Data buf[sizeof(l)];          // Temp buffer to store l in
  Counter i = 0;                // buffer index
  while (l) {                   // While more bits in l
    assert(i < sizeof(l));      // no more buffer space
    buf[i] = Data(l);           // Peel off lower order bits
    l >>= 16;   // Shift next bits into place
    i++;
  }
  if (i > 0)
    this->data = new Data[this->count=i]; // Allocate permanent data

  while (i--)     // Save buffer into perm. data
    this->data[i] = buf[i];
}

//: Creates a vnl_bignum from a single-precision floating point number.

vnl_bignum::vnl_bignum (float f)
: count(0), sign(1), data(0)
{
  double d = f;
  if (d < 0.0) {                // Get sign of d
    d = -d;                     // Get absolute value of d
    this->sign = -1;
  }

  if (!vnl_math_isfinite(d)) {
    // Infinity is represented as: count=1, data[0]=0.
    // This is an otherwise unused representation, since 0 is represented as count=0.
    this->count = 1;
    this->data = new Data[1];
    this->data[0] = 0;
  } else if (d >= 1.0) {
    // Note: 0x10000L == 1 >> 16: the (assumed) size of unsigned short is 16 bits.
    vcl_vector<Data> buf;
    while (d >= 1.0) {
      buf.push_back( Data(vcl_fmod(d,0x10000L)) );  // Get next data "digit" from d
      d /= 0x10000L;                                // Shift d right 1 data "digit"
    }
    // Allocate and copy into permanent buffer
    this->data = buf.size()>0 ? new Data[buf.size()] : 0;
    this->count = buf.size();
    vcl_copy( buf.begin(), buf.end(), data );
  }
}

//: Creates a vnl_bignum from a double floating point number.

vnl_bignum::vnl_bignum (double d)
: count(0), sign(1), data(0)
{
  if (d < 0.0) {                // Get sign of d
    d = -d;                     // Get absolute value of d
    this->sign = -1;
  }

  if (!vnl_math_isfinite(d)) {
    // Infinity is represented as: count=1, data[0]=0.
    // This is an otherwise unused representation, since 0 is represented as count=0.
    this->count = 1;
    this->data = new Data[1];
    this->data[0] = 0;
  } else if (d >= 1.0) {
    // Note: 0x10000L == 1 >> 16: the (assumed) size of unsigned short is 16 bits.
    vcl_vector<Data> buf;
    while (d >= 1.0) {
      buf.push_back( Data(vcl_fmod(d,0x10000L)) );  // Get next data "digit" from d
      d /= 0x10000L;                                // Shift d right 1 data "digit"
    }
    // Allocate and copy into permanent buffer
    this->data = buf.size()>0 ? new Data[buf.size()] : 0;
    this->count = buf.size();
    vcl_copy( buf.begin(), buf.end(), data );
  }
}

//: Creates a vnl_bignum from a "long double" floating point number.

vnl_bignum::vnl_bignum (long double d)
: count(0), sign(1), data(0)
{
  if (d < 0.0) {                // Get sign of d
    d = -d;                     // Get absolute value of d
    this->sign = -1;
  }

  if (!vnl_math_isfinite(d)) {
    // Infinity is represented as: count=1, data[0]=0.
    // This is an otherwise unused representation, since 0 is represented as count=0.
    this->count = 1;
    this->data = new Data[1];
    this->data[0] = 0;
  } else if (d >= 1.0) {
    // Note: 0x10000L == 1 >> 16: the (assumed) size of unsigned short is 16 bits.
    vcl_vector<Data> buf;
    while (d >= 1.0) {
      buf.push_back( Data(vcl_fmod(d,0x10000L)) );  // Get next data "digit" from d
      d /= 0x10000L;                                // Shift d right 1 data "digit"
    }
    // Allocate and copy into permanent buffer
    this->data = (buf.size()>0 ? new Data[buf.size()] : 0);
    this->count = buf.size();
    vcl_copy( buf.begin(), buf.end(), data );
  }
}


#if 0 // old, original Texas Instruments implementation - PVr

static bool is_decimal(const char *s)
{
  if (*s == '+' || *s == '-') ++s;
  if (*s < '1' || *s > '9') return false;
  while (*s >= '0' && *s <= '9') ++s;
  if (*s == 'l' || *s == 'L') ++s;
  return *s == '\0';
}

static bool is_exponential(const char *s)
{
  if (*s == '+' || *s == '-') ++s;
  if (*s < '1' || *s > '9') return false;
  while (*s >= '0' && *s <= '9') ++s;
  if (*s != 'e' && *s != 'E') return false;
  ++s;
  if (*s < '1' || *s > '9') return false;
  while (*s >= '0' && *s <= '9') ++s;
  return *s == '\0';
}

static bool is_hexadecimal(const char *s)
{
  if (*s == '+' || *s == '-') ++s;
  if (*s != '0') return false;
  ++s;
  if (*s != 'x' && *s != 'X') return false;
  ++s;
  if ((*s < '0' || *s > '9') &&
      (*s < 'a' || *s > 'f') &&
      (*s < 'A' || *s > 'F')) return false;
  while ((*s >= '0' && *s <= '9') ||
         (*s >= 'a' && *s <= 'f') ||
         (*s >= 'A' && *s <= 'F')) ++s;
  if (*s == 'l' || *s == 'L') ++s;
  return *s == '\0';
}

static bool is_octal(const char *s)
{
  if (*s == '+' || *s == '-') ++s;
  if (*s != '0') return false;
  while (*s >= '0' && *s <= '7') ++s;
  if (*s == 'l' || *s == 'L') ++s;
  return *s == '\0';
}

#else // new implementation, also to be used for operator>> - PVr

static char rt[4096];
static int rt_pos = 0;

static char next(const char*& s, vcl_istream** is)
{
  if (!is || *s) { char c = *s; if (c) ++rt_pos, ++s; return c; }
  if (rt_pos == 4096) return '\0';
  (*is)->get(rt[rt_pos]); // read a single byte from istream
  if (*s) ++s; // in case s == rt+rt_pos
  rt[++rt_pos] = '\0'; return rt[rt_pos-1];
}

static bool is_decimal(const char* s, vcl_istream** is = 0)
{
  rt_pos = 0;
  char c = next(s,is);
  while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
  if (c == '+' || c == '-') c = next(s,is);
  if (c < '1' || c > '9') return false;
  while (c >= '0' && c <= '9') c = next(s,is);
  if (c == 'l' || c == 'L') c = next(s,is);
  if (rt_pos > 0) rt[++rt_pos] = '\0';
  return is ? true : c == '\0';
}

static bool is_exponential(const char* s, vcl_istream** is = 0)
{
  rt_pos = 0;
  char c = next(s,is);
  while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
  if (c == '+' || c == '-') c = next(s,is);
  if (c < '1' || c > '9') return false;
  while (c >= '0' && c <= '9') c = next(s,is);
  if (c != 'e' && c != 'E') return false;
  c = next(s,is);
  if (c == '+') c = next(s,is); // no negative exponent!
  if (c < '0' || c > '9') return false;
  while (c >= '0' && c <= '9') c = next(s,is);
  if (rt_pos > 0) rt[++rt_pos] = '\0';
  return is ? true : c == '\0';
}

static bool is_hexadecimal(const char* s, vcl_istream** is = 0)
{
  rt_pos = 0;
  char c = next(s,is);
  while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
  if (c == '+' || c == '-') c = next(s,is);
  if (c != '0') return false;
  c = next(s,is);
  if (c != 'x' && c != 'X') return false;
  c = next(s,is);
  if ((c < '0' || c > '9') &&
      (c < 'a' || c > 'f') &&
      (c < 'A' || c > 'F')) return false;
  while ((c >= '0' && c <= '9') ||
         (c >= 'a' && c <= 'f') ||
         (c >= 'A' && c <= 'F')) c = next(s,is);
  if (c == 'l' || c == 'L') c = next(s,is);
  if (rt_pos > 0) rt[++rt_pos] = '\0';
  return is ? true : c == '\0';
}

static bool is_octal(const char* s, vcl_istream** is = 0)
{
  rt_pos = 0;
  char c = next(s,is);
  while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
  if (c == '+' || c == '-') c = next(s,is);
  if (c != '0') return false;
  while (c >= '0' && c <= '7') c = next(s,is);
  if (c == 'l' || c == 'L') c = next(s,is);
  if (rt_pos > 0) rt[++rt_pos] = '\0';
  return is ? true : c == '\0';
}

static bool is_plus_inf(const char* s, vcl_istream** is = 0)
{
  rt_pos = 0;
  char c = next(s,is);
  while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
  if (c == '+') c = next(s,is);
  if (c != 'I') return false; c = next(s,is);
  if (c != 'n') return false; c = next(s,is);
  if (c != 'f') return false; c = next(s,is);
  if (c == 'i') c = next(s,is);
  if (c == 'n') c = next(s,is);
  if (c == 'i') c = next(s,is);
  if (c == 't') c = next(s,is);
  if (c == 'y') c = next(s,is);
  if (rt_pos > 0) rt[++rt_pos] = '\0';
  return is ? true : c == '\0';
}

static bool is_minus_inf(const char* s, vcl_istream** is = 0)
{
  rt_pos = 0;
  char c = next(s,is);
  while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
  if (c != '-') return false; c = next(s,is);
  if (c != 'I') return false; c = next(s,is);
  if (c != 'n') return false; c = next(s,is);
  if (c != 'f') return false; c = next(s,is);
  if (c == 'i') c = next(s,is);
  if (c == 'n') c = next(s,is);
  if (c == 'i') c = next(s,is);
  if (c == 't') c = next(s,is);
  if (c == 'y') c = next(s,is);
  if (rt_pos > 0) rt[++rt_pos] = '\0';
  return is ? true : c == '\0';
}

#endif // new implementation - PVr

//: Creates a vnl_bignum from the character string representation.

vnl_bignum::vnl_bignum (const char *s)
: count(0), sign(1), data(0)
{
  // decimal:     "^ *[-+]?[1-9][0-9]*$"
  // exponential: "^ *[-+]?[1-9][0-9]*[eE][+]?[0-9]+$"
  // hexadecimal: "^ *[-+]?0[xX][0-9a-fA-F]+$"
  // octal:       "^ *[-+]?0[0-7]*$"
  // infinity:    "^ *[-+]?Inf(inity)?$"

  if (is_plus_inf(s))
    sign=1,count=1,data=new Data[1],data[0]=0;
  else if (is_minus_inf(s))
    sign=-1,count=1,data=new Data[1],data[0]=0;
  else if (is_decimal(s))               // If string is decimal
    this->dtoBigNum(s);                 // convert decimal to vnl_bignum
  else if (is_exponential(s))           // If string is exponential
    this->exptoBigNum(s);               // convert exp. to vnl_bignum
  else if (is_hexadecimal(s))           // If string is hex,
    this->xtoBigNum(s);                 // convert hex to vnl_bignum
  else if (is_octal(s))                 // If string is octal
    this->otoBigNum(s);                 // convert octal to vnl_bignum
  else {                                // Otherwise
    vcl_cerr << "Cannot convert string " << s << " to vnl_bignum\n";
  }
}

//: Reads a vnl_bignum from a stream

vcl_istream& operator>> (vcl_istream& is, vnl_bignum& x)
{
  // decimal:     "^ *[-+]?[1-9][0-9]*$"
  // exponential: "^ *[-+]?[1-9][0-9]*[eE][+]?[0-9]+$"
  // hexadecimal: "^ *[-+]?0[xX][0-9a-fA-F]+$"
  // octal:       "^ *[-+]?0[0-7]*$"
  vcl_istream* isp = &is;
  rt[0] = '\0';

//delete[] x.data;                      // Delete existing data
  if (is_plus_inf(rt,&isp))
    x.sign=1,x.count=1,x.data=new Data[1],x.data[0]=0;
  else if (is_minus_inf(rt,&isp))
    x.sign=-1,x.count=1,x.data=new Data[1],x.data[0]=0;
  if (is_exponential(rt,&isp))          // If input stream string is exponential
    x.exptoBigNum(rt);                  // convert exp. to vnl_bignum
  else if (is_decimal(rt,&isp))         // If string is decimal
    x.dtoBigNum(rt);                    // convert decimal to vnl_bignum
  else if (is_hexadecimal(rt,&isp))     // If string is hex,
    x.xtoBigNum(rt);                    // convert hex to vnl_bignum
  else if (is_octal(rt,&isp))           // If string is octal
    x.otoBigNum(rt);                    // convert octal to vnl_bignum
  else {                                // Otherwise
    vcl_cerr << "Cannot convert string " << rt << " to vnl_bignum\n";
    x = 0L;
  }
  return is; // FIXME - should probably push back read characters to istream
}

//: Copies the contents of vnl_bignum b.

vnl_bignum::vnl_bignum (const vnl_bignum& b)
: count(b.count), sign(b.sign)
{
  this->data = b.data ? new Data[b.count] : 0;  // Allocate data if necessary
  for (Counter i = 0; i < this->count; ++i)     // Copy b data
    this->data[i] = b.data[i];
}


//: Frees space for vnl_bignum.

vnl_bignum::~vnl_bignum ()
{
  delete [] this->data; this->count = 0;        // Delete any allocated data
}

//: Copies rhs vnl_bignum to lhs vnl_bignum.

vnl_bignum& vnl_bignum::operator= (const vnl_bignum& rhs)
{
  if (this != &rhs) {                           // Avoid self-assignment
    delete[] this->data;                        // Delete existing data
    this->count = rhs.count;                    // Copy rhs's count
    this->data = rhs.data ? new Data[rhs.count] : 0; // Allocate data if necessary
    for (Counter i = 0; i < rhs.count; ++i)     // Copy rhs's data
      this->data[i] = rhs.data[i];
    this->sign = rhs.sign;                      // Copy rhs's sign
  }
  return *this;                                 // Return reference
}

//: Returns the negation of a vnl_bignum.

vnl_bignum vnl_bignum::operator- () const
{
  vnl_bignum neg(*this);
  if (neg.count)                // So long as this is non-zero
    neg.sign = -neg.sign;       // Flip its sign
  return neg;
}


//: Prefix increment. Increments a vnl_bignum by 1, and returns it.

vnl_bignum& vnl_bignum::operator++ ()
{
  if (this->is_infinity()) return *this;
  if (this->count==0)
  {
    this->resize(1);
    this->data[0] = 1;
    this->sign = +1;
    return *this;
  }

  if (this->sign > 0) increment(*this);
  else decrement(*this);

  return *this;
}


//: Prefix decrement. Decrements a vnl_bignum by 1, and returns it.

vnl_bignum& vnl_bignum::operator-- ()
{
  if (this->is_infinity()) return *this;
  if (this->count==0)
  {
    this->resize(1);
    this->data[0] = 1;
    this->sign = -1;
    return *this;
  }

  if (this->sign < 0) increment(*this);
  else decrement(*this);

  return *this;
}

//: Adds two vnl_bignums, and returns new sum.

vnl_bignum vnl_bignum::operator+(const vnl_bignum& b) const
{
  // Infinity arithmetic:
  assert (! b.is_minus_infinity() || ! this->is_plus_infinity() ); // +Inf-Inf
  assert (! b.is_plus_infinity() || ! this->is_minus_infinity() ); // -Inf+Inf
  if (b.is_infinity()) { return b; }
  if (this->is_infinity()) { return *this; }

  vnl_bignum sum;                       // Init sum to zero
  if (this->sign == b.sign) {           // If both have same sign
    add(*this,b,sum);                   //   Do simple addition
    sum.sign = this->sign;              // Attach proper sign
  }
  else {                                // Else different signs
    int mag = magnitude_cmp(*this,b);   // Determine relative sizes
    if (mag > 0) {                      // If abs(*this) > abs(b)
      subtract(*this,b,sum);            //   sum = *this - b
      sum.sign = this->sign;            // Sign of sum follows *this
    }
    else if (mag < 0) {                 // Else if abs(*this) < abs(b)
      subtract(b,*this,sum);            //   sum = b - *this
      sum.sign = b.sign;                // Sign of sum follows b
    }                                   // (Else abs(*this) == abs(b)
  }                                     //   so sum must be zero)
  return sum;                           // shallow swap on return
}


//: Multiplies this with a vnl_bignum

vnl_bignum& vnl_bignum::operator*= (const vnl_bignum& b)
{
  // Infinity arithmetic:
  assert (! b.is_infinity() || this->count != 0 ); // multiplication 0*Inf
  assert (! this->is_infinity() || b.count != 0 ); // multiplication Inf*0
  if (b.is_infinity()) return (*this) = (this->sign<0 ? -b : b);
  if (this->is_infinity()) return (*this) = (b.sign<0 ? -(*this) : *this);

  if (b.count == 0 || this->count == 0)
    return (*this)=0L;
  vnl_bignum prod;
  prod.data = new Data[prod.count = this->count + b.count]; // allocate data for product
  for (Counter i = 0; i < b.count; i++)         //   multiply each b "digit"
    multiply_aux(*this, b.data[i], prod, i);    //   times b1 and add to total
  prod.sign = this->sign * b.sign;              //   determine correct sign
  prod.trim();                                  //   trim excess data and ret.
  return (*this)=prod;
}


//: Divides this by a vnl_bignum

vnl_bignum& vnl_bignum::operator/= (const vnl_bignum& b)
{
  // Infinity arithmetic:
  assert (! b.is_infinity() || ! this->is_infinity() ); // division Inf/Inf
  if (b.is_infinity()) return (*this)=0L;
  if (this->is_infinity()) return (*this) = (b.sign<0 ? -(*this) : *this);
  assert (b.count!=0 || this->count != 0); // division 0/0
  if (b.count == 0)
    return (*this) = (this->sign < 0 ? vnl_bignum("-Inf") : vnl_bignum("+Inf"));

  vnl_bignum quot, r;          // Quotient and remainder
  divide(*this,b,quot,r);      // Call divide fn
  return (*this) = quot;
}

//: Divides this by a vnl_bignum and replaces this by remainder.

vnl_bignum& vnl_bignum::operator%= (const vnl_bignum& b)
{
  // Infinity arithmetic:
  assert (! b.is_infinity() || ! this->is_infinity() ); // division Inf/Inf
  if (b.is_infinity()) return *this;                    // remainder of x/Inf is x.
  if (this->is_infinity()) return (*this) = 0L;         // convention: remainder is 0
  assert (b.count!=0 || this->count != 0);              // division 0/0
  if (b.count == 0) return (*this) = 0L;                // convention: remainder is 0

  vnl_bignum remain, q;        // Quotient and remainder
  divide(*this,b,q,remain);    // divide by b and save remainder
  return (*this) = remain;     // shallow swap on return
}


//: Shifts bignum to the left l digits.

vnl_bignum vnl_bignum::operator<< (int l) const
{
  // Infinity arithmetic:
  if (this->is_infinity()) return *this;

  if (l == 0 || *this == 0L)            // if either arg is zero
    return *this;
  if (l < 0)                            // if shift amt is negative
    return right_shift(*this,-l);       //   do an actual right shift
  else                                  // otherwise
    return left_shift(*this,l);         //   do a left shift
}


//: Shifts bignum to the right l digits.

vnl_bignum vnl_bignum::operator>> (int l) const
{
  // Infinity arithmetic:
  if (this->is_infinity()) return *this;

  if (l == 0 || *this == 0L)            // if either arg is zero
    return *this;
  if (l < 0)                            // if shift amt is negative
    return left_shift(*this,-l);        //   do an actual left shift
  else                                  // else
    return right_shift(*this,l);        //   do a right shift
}


//: Two vnl_bignums are equal if and only if they have the same integer representation.

bool vnl_bignum::operator== (const vnl_bignum& rhs) const
{
  if (this != &rhs) {                           // Check address
    if (this->sign != rhs.sign) return false;   // Different sign implies !=
    if (this->count != rhs.count) return false; // Different size implies !=
    for (Counter i = 0; i < this->count; i++)   // Each data element the same?
      if (this->data[i] != rhs.data[i]) return false; // No. Return !=
  }
  return true;                                    // Yes. Return ==
}


//: Compares two vnl_bignums.

bool vnl_bignum::operator< (const vnl_bignum& rhs) const
{
  if (this->sign < rhs.sign) return true;       // Different signs?
  if (this->sign > rhs.sign) return false;
  if (this->sign == 1)                          // Both signs == 1
    return magnitude_cmp(*this,rhs) < 0;        // this must be smaller
  else                                          // Both signs == -1
    return magnitude_cmp(*this,rhs) > 0;        // this must be larger
}


//: Formatted output for bignum.

vcl_ostream& operator<< (vcl_ostream& os, const vnl_bignum& b)
{
  vnl_bignum d = b;                     // Copy the input vnl_bignum
  if (d.sign == -1) {                   // If it's negative
    os << '-';                          //   Output leading minus sign
    d.sign = 1;                         //   Make d positive for divide
  }
  if (d.is_infinity()) return os<<"Inf";
  vnl_bignum q,r;                       // Temp quotient and remainder
  char *cbuf = new char[5 * (b.count+1)];   // Temp character buffer
  Counter i = 0;
  do {                                  // repeat:
    divide(d,10L,q,r);                  //   Divide vnl_bignum by ten
    cbuf[i++] = char(long(r) + '0');    //   Get one's digit
    d = q;                              //   Then discard one's digit
    q = r = 0L;                         //   Prep for next divide
  } while (d != 0L);                    // until no more one's digits
  do {                                  // repeat;
    os << cbuf[--i];                    //   output char buf in reverse
  } while (i);                          // until no more chars
  delete [] cbuf;                       // delete temp char buf
  return os;                            // return output stream
}

//: Convert the number to a decimal representation in a string.
vcl_string& vnl_bignum_to_string (vcl_string& s, const vnl_bignum& b)
{
  s.erase();
  vcl_string::size_type insert_point = 0; // keep record of location of first number.

  vnl_bignum d = b;                     // Copy the input vnl_bignum
  if (d.sign == -1) {                   // If it's negative
    s.insert(insert_point,"-");         //   Output leading minus sign
    d.sign = 1;                         //   Make d positive for divide
    ++insert_point;                     // keep record of location of first number.
  }
  if (d.is_infinity()) return s+="Inf";
  vnl_bignum q,r;                       // Temp quotient and remainder
  do {                                  // repeat:
    divide(d,10L,q,r);                  //   Divide vnl_bignum by ten
    s.insert(insert_point, 1, char('0'+long(r))); //   Get one's digit, and insert it at head.
    d = q;                              //   Then discard one's digit
    q = r = 0L;                         //   Prep for next divide
  } while (d != 0L);                    // until no more one's digits
  return s;
}

//: Convert the number from a decimal representation in a string.
vnl_bignum& vnl_bignum_from_string (vnl_bignum& b, const vcl_string& s)
{
  // decimal:     "^ *[-+]?[1-9][0-9]*$"
  // Infinity:    "^ *[-+]?Inf(inity)?$"

  if (is_plus_inf(s.c_str()))
    b=vnl_bignum("+Inf");
  else if (is_minus_inf(s.c_str()))
    b=vnl_bignum("-Inf");
  else
    b.dtoBigNum(s.c_str());             // convert decimal to vnl_bignum
  return b;
}


//: Implicit conversion from a vnl_bignum to a short.
vnl_bignum::operator short () const
{
  short s = 0;
  for (Counter i = this->count; i > 0; )
    s = short(s*0x10000 + this->data[--i]);
  return this->sign*s;
}


//: Implicit conversion from a vnl_bignum to an int.
vnl_bignum::operator int () const
{
  int j = 0;
  for (Counter i = this->count; i > 0; )
    j = int(j*0x10000 + this->data[--i]);
  return this->sign*j;
}


//: Implicit conversion from a vnl_bignum to a long.
vnl_bignum::operator long () const
{
  long l = 0;
  for (Counter i = this->count; i > 0; )
    l = l*0x10000L + this->data[--i];
  return this->sign*l;
}


//: Implicit conversion from a vnl_bignum to a float.
vnl_bignum::operator float () const
{
  float f = 0.0f;
  for (Counter i = this->count; i > 0; )
    f = f*0x10000 + this->data[--i];
  if (this->is_infinity()) f = vcl_numeric_limits<float>::infinity();
  return this->sign*f;
}


//: Implicit conversion from a vnl_bignum to a double.

vnl_bignum::operator double () const
{
  double d = 0.0;
  for (Counter i = this->count; i > 0; )
    d = d*0x10000 + this->data[--i];
  if (this->is_infinity()) d = vcl_numeric_limits<double>::infinity();
  return this->sign*d;
}

//: Implicit conversion from a vnl_bignum to a long double.

vnl_bignum::operator long double () const
{
  long double d = 0.0;
  for (Counter i = this->count; i > 0; )
    d = d*0x10000 + this->data[--i];
  if (this->is_infinity()) d = vcl_numeric_limits<long double>::infinity();
  return this->sign*d;
}

//: dump the contents of a vnl_bignum to a stream, default cout.

void vnl_bignum::dump (vcl_ostream& os) const
{
  os << "{count=" << this->count       // output count field
     << ", sign=" << this->sign        // output sign field
     << ", data=" << this->data        // output data pointer
     << ", value=" << *this
     << ", {";
  // format string == "%04X%s" or "%02X%s", etc.
  //  static char format_str[10] =
  //    {'%','0',char(2*2 + '0'),'X','%','s'};
  //  format_str[2] = char(2*2 + '0');
  if (this->count > 0) { // output data array
    os << vcl_hex;
    for (Counter i = this->count; i > 1; --i)
      os << (this->data[i - 1]) << ',';
    os << (this->data[0]) << vcl_dec;
  }
  os << "}}\n";                         // close brackets
}


//: Converts decimal string to a vnl_bignum.

int vnl_bignum::dtoBigNum (const char *s)
{
  this->resize(0); sign = 1;            // Reset number to 0.
  Counter len = 0;                      // No chars converted yet
  while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s; // skip whitespace
  if (s[0] == '-' || s[0] == '+') len++;// Skip over leading +,-
  while (vcl_isdigit(s[len])) {         // If current char is digit
    (*this) = ((*this) * 10L) +         // Shift vnl_bignum left a decimal
      vnl_bignum(long(s[len++] - '0')); // digit and add new digit
  }
  if (s[0] == '-') this->sign = -1;     // If s had leading -, note it
  return len;                           // Return # of chars processed
}

//: convert exponential string to a vnl_bignum

void vnl_bignum::exptoBigNum (const char *s)
{
  while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s; // skip whitespace
  Counter pos = this->dtoBigNum(s) + 1; // Convert the base, skip [eE]
  long pow = vcl_atol(s + pos);         // Convert the exponent to long
  while (pow-- > 0)                     // Raise vnl_bignum to the given
    *this = (*this) * 10L;              // power
}


//: convert hex character to integer hex value (ASCII or EBCDIC)
// - Inputs:  character representation of a hex number
// - Outputs: integer value of the hex number

unsigned int ctox (int c)
{
  if ('0' <= c && c <= '9')
    return c - '0';
  if ('a' <= c && c <= 'f')
    return c - 'a' + 10;
  return c - 'A' + 10;
}

//: convert hex string to vnl_bignum

void vnl_bignum::xtoBigNum (const char *s)
{
  this->resize(0); sign = 1;            // Reset number to 0.
  while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s; // skip whitespace
  Counter size = vcl_strlen(s);
  Counter len = 2;                      // skip leading "0x"
  while (len < size) {                  // While there are more chars
    (*this) = ((*this) * 16L) +         // Shift vnl_bignum left one hex
      vnl_bignum(long(ctox(s[len++]))); //   digit and add next digit
  }
}


//: convert octal string to vnl_bignum

void vnl_bignum::otoBigNum (const char *s)
{
  this->resize(0); sign = 1;           // Reset number to 0.
  while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s; // skip whitespace
  Counter size = vcl_strlen(s);
  Counter len = 0;                      // No chars converted yet
  while (len < size) {                  // While there are more chars
    (*this) = ((*this) * 8L) +          // Shift vnl_bignum left 1 oct dig.
      vnl_bignum(long(s[len++] - '0')); // Add next character value
  }
}

//: change the data allotment for a vnl_bignum

void vnl_bignum::resize (short new_count)
{
  assert(new_count >= 0);
  if (new_count == this->count) return;
  Data *new_data = (new_count > 0 ? new Data[new_count] : 0); // Allocate data if necessary

  if (this->count <= new_count) {       // Copy old data into new
    short i = 0;
    for (; i < this->count; i++)
      new_data[i] = this->data[i];
    for (; i < new_count; i++)
      new_data[i] = 0;
  }
  else {
    for (short i = 0; i < new_count; i++)
      new_data[i] = this->data[i];
  }

  delete [] this->data;                 // Get rid of old data
  this->data = new_data;                // Point to new data
  this->count = new_count;              // Save new count
}


//: trim non-infinite vnl_bignum of excess data allotment

vnl_bignum& vnl_bignum::trim ()
{
  Counter i = this->count;
  for (; i > 0; i--)                    // Skip over high-order words
    if (this->data[i - 1] != 0) break;  //   that are zero
  if (i < this->count) {                // If there are some such words
    this->count = i;                    // Update the count
    Data *new_data = (i > 0 ? new Data[i] : 0); // Allocate data if necessary
    for (; i > 0; i--)                  // Copy old data into new
      new_data[i - 1] = this->data[i - 1];
    delete [] this->data;               // Delete old data
    this->data = new_data;              // Point to new data
  }
  return *this;                         // return reference to vnl_bignum
}


//: add two non-infinite vnl_bignum values and save their sum

void add (const vnl_bignum& b1, const vnl_bignum& b2, vnl_bignum& sum)
{
  const vnl_bignum *bmax, *bmin;        // Determine which of the two
  if (b1.count >= b2.count) {           //   addends has the most
    bmax = &b1;                         //   data.
    bmin = &b2;
  }
  else {
    bmax = &b2;
    bmin = &b1;
  }
//delete[] sum.data;                    // Delete existing data
  sum.data = (sum.count = bmax->count) > 0 ?   // Allocate data for their sum
             new Data[sum.count] : 0;
  unsigned long temp, carry = 0;
  Counter i = 0;
  while (i < bmin->count) {             // Add, element by element.
    // Add both elements and carry
    temp = (unsigned long)b1.data[i] + (unsigned long)b2.data[i] + carry;
    carry = temp/0x10000L;              // keep track of the carry
    sum.data[i] = Data(temp);           // store sum
    i++;                                // go to next element
  }
  while (i < bmax->count) {             // bmin has no more elements
    temp = bmax->data[i] + carry;       // propagate the carry through
    carry = temp/0x10000L;              // the rest of bmax's elements
    sum.data[i] = Data(temp);           // store sum
    i++;
  }
  if (carry) {                          // if carry left over
    sum.resize(bmax->count + 1);        //   allocate another word
    sum.data[bmax->count] = 1;          //   save the carry in it
  }
}

//: Add 1 to bnum (unsigned, non-infinite)
void increment (vnl_bignum& bnum)
{
  Counter i = 0;
  unsigned long carry = 1;
  while (i < bnum.count && carry) {             // increment, element by element.
    unsigned long temp = (unsigned long)bnum.data[i] + carry;
    carry = temp/0x10000L;
    bnum.data[i] = (Data)temp;
    ++i;
  }
  if (carry)
  {
    bnum.resize(bnum.count+1);
    bnum.data[bnum.count-1] = 1;
  }
}


//: subtract bmin from bmax (unsigned, non-infinite), result in diff

void subtract (const vnl_bignum& bmax, <