core/vnl/vnl_bignum.cxx

Go to the documentation of this file.
00001 // This is core/vnl/vnl_bignum.cxx
00002 #include "vnl_bignum.h"
00003 //:
00004 // \file
00005 
00006 #include <vcl_cctype.h>    // Include character macros
00007 #include <vcl_cstdlib.h>   // for vcl_atol
00008 #include <vcl_cstring.h>   // for vcl_strlen
00009 #include <vcl_cmath.h>     // for vcl_fmod
00010 #include <vcl_algorithm.h> // for vcl_copy
00011 #include <vcl_vector.h>
00012 #include <vcl_cassert.h>
00013 #include <vcl_iostream.h>
00014 #include <vcl_limits.h>
00015 #include <vnl/vnl_math.h> // for vnl_math_isfinite(double)
00016 
00017 typedef unsigned short Counter;
00018 typedef unsigned short Data;
00019 
00020 //: Creates a zero vnl_bignum.
00021 
00022 vnl_bignum::vnl_bignum ()
00023 : count(0), sign(1), data(0)
00024 {
00025 }
00026 
00027 //: Creates a vnl_bignum from a long integer.
00028 
00029 vnl_bignum::vnl_bignum (long l)
00030 : count(0), sign(1), data(0)
00031 {
00032   if (l < 0) {                  // Get correct sign
00033     l = -l;                     // Get absolute value of l
00034     this->sign = -1;
00035   }
00036   Data buf[sizeof(l)];          // Temp buffer to store l in
00037   Counter i = 0;                // buffer index
00038   while (l) {                   // While more bits in l
00039     assert(i < sizeof(l));      // no more buffer space
00040     buf[i] = Data(l);           // Peel off lower order bits
00041     l >>= 16;   // Shift next bits into place
00042     ++i;
00043   }
00044   if (i > 0)
00045     this->data = new Data[this->count=i]; // Allocate permanent data
00046 
00047   while (i--)     // Save buffer into perm. data
00048     this->data[i] = buf[i];
00049 }
00050 
00051 //: Creates a vnl_bignum from an integer.
00052 
00053 vnl_bignum::vnl_bignum (int l)
00054 : count(0), sign(1), data(0)
00055 {
00056   if (l < 0) {                 // Get correct sign
00057     l = -l;                     // Get absolute value of l
00058     this->sign = -1;
00059   }
00060   Data buf[sizeof(l)];          // Temp buffer to store l in
00061   Counter i = 0;                // buffer index
00062   while (l) {                   // While more bits in l
00063     assert(i < sizeof(l));      // no more buffer space
00064     buf[i] = Data(l);           // Peel off lower order bits
00065     l >>= 16;   // Shift next bits into place
00066     i++;
00067   }
00068   if (i > 0)
00069     this->data = new Data[this->count=i]; // Allocate permanent data
00070 
00071   while (i--)     // Save buffer into perm. data
00072     this->data[i] = buf[i];
00073 }
00074 
00075 //: Creates a vnl_bignum from an unsigned long integer.
00076 
00077 vnl_bignum::vnl_bignum (unsigned long l)
00078 : count(0), sign(1), data(0)
00079 {
00080   Data buf[sizeof(l)];          // Temp buffer to store l in
00081   Counter i = 0;                // buffer index
00082   while (l) {                   // While more bits in l
00083     assert(i < sizeof(l));      // no more buffer space
00084     buf[i] = Data(l);           // Peel off lower order bits
00085     l >>= 16;   // Shift next bits into place
00086     i++;
00087   }
00088   if (i > 0)
00089     this->data = new Data[this->count=i]; // Allocate permanent data
00090 
00091   while (i--)     // Save buffer into perm. data
00092     this->data[i] = buf[i];
00093 }
00094 
00095 //: Creates a vnl_bignum from an unsigned integer.
00096 
00097 vnl_bignum::vnl_bignum (unsigned int l)
00098 : count(0), sign(1), data(0)
00099 {
00100   Data buf[sizeof(l)];          // Temp buffer to store l in
00101   Counter i = 0;                // buffer index
00102   while (l) {                   // While more bits in l
00103     assert(i < sizeof(l));      // no more buffer space
00104     buf[i] = Data(l);           // Peel off lower order bits
00105     l >>= 16;   // Shift next bits into place
00106     i++;
00107   }
00108   if (i > 0)
00109     this->data = new Data[this->count=i]; // Allocate permanent data
00110 
00111   while (i--)     // Save buffer into perm. data
00112     this->data[i] = buf[i];
00113 }
00114 
00115 //: Creates a vnl_bignum from a single-precision floating point number.
00116 
00117 vnl_bignum::vnl_bignum (float f)
00118 : count(0), sign(1), data(0)
00119 {
00120   double d = f;
00121   if (d < 0.0) {                // Get sign of d
00122     d = -d;                     // Get absolute value of d
00123     this->sign = -1;
00124   }
00125 
00126   if (!vnl_math_isfinite(d)) {
00127     // Infinity is represented as: count=1, data[0]=0.
00128     // This is an otherwise unused representation, since 0 is represented as count=0.
00129     this->count = 1;
00130     this->data = new Data[1];
00131     this->data[0] = 0;
00132   } else if (d >= 1.0) {
00133     // Note: 0x10000L == 1 >> 16: the (assumed) size of unsigned short is 16 bits.
00134     vcl_vector<Data> buf;
00135     while (d >= 1.0) {
00136       buf.push_back( Data(vcl_fmod(d,0x10000L)) );  // Get next data "digit" from d
00137       d /= 0x10000L;                                // Shift d right 1 data "digit"
00138     }
00139     // Allocate and copy into permanent buffer
00140     this->data = buf.size()>0 ? new Data[buf.size()] : 0;
00141     this->count = buf.size();
00142     vcl_copy( buf.begin(), buf.end(), data );
00143   }
00144 }
00145 
00146 //: Creates a vnl_bignum from a double floating point number.
00147 
00148 vnl_bignum::vnl_bignum (double d)
00149 : count(0), sign(1), data(0)
00150 {
00151   if (d < 0.0) {                // Get sign of d
00152     d = -d;                     // Get absolute value of d
00153     this->sign = -1;
00154   }
00155 
00156   if (!vnl_math_isfinite(d)) {
00157     // Infinity is represented as: count=1, data[0]=0.
00158     // This is an otherwise unused representation, since 0 is represented as count=0.
00159     this->count = 1;
00160     this->data = new Data[1];
00161     this->data[0] = 0;
00162   } else if (d >= 1.0) {
00163     // Note: 0x10000L == 1 >> 16: the (assumed) size of unsigned short is 16 bits.
00164     vcl_vector<Data> buf;
00165     while (d >= 1.0) {
00166       buf.push_back( Data(vcl_fmod(d,0x10000L)) );  // Get next data "digit" from d
00167       d /= 0x10000L;                                // Shift d right 1 data "digit"
00168     }
00169     // Allocate and copy into permanent buffer
00170     this->data = buf.size()>0 ? new Data[buf.size()] : 0;
00171     this->count = buf.size();
00172     vcl_copy( buf.begin(), buf.end(), data );
00173   }
00174 }
00175 
00176 //: Creates a vnl_bignum from a "long double" floating point number.
00177 
00178 vnl_bignum::vnl_bignum (long double d)
00179 : count(0), sign(1), data(0)
00180 {
00181   if (d < 0.0) {                // Get sign of d
00182     d = -d;                     // Get absolute value of d
00183     this->sign = -1;
00184   }
00185 
00186   if (!vnl_math_isfinite(d)) {
00187     // Infinity is represented as: count=1, data[0]=0.
00188     // This is an otherwise unused representation, since 0 is represented as count=0.
00189     this->count = 1;
00190     this->data = new Data[1];
00191     this->data[0] = 0;
00192   } else if (d >= 1.0) {
00193     // Note: 0x10000L == 1 >> 16: the (assumed) size of unsigned short is 16 bits.
00194     vcl_vector<Data> buf;
00195     while (d >= 1.0) {
00196       buf.push_back( Data(vcl_fmod(d,0x10000L)) );  // Get next data "digit" from d
00197       d /= 0x10000L;                                // Shift d right 1 data "digit"
00198     }
00199     // Allocate and copy into permanent buffer
00200     this->data = (buf.size()>0 ? new Data[buf.size()] : 0);
00201     this->count = buf.size();
00202     vcl_copy( buf.begin(), buf.end(), data );
00203   }
00204 }
00205 
00206 
00207 #if 0 // old, original Texas Instruments implementation - PVr
00208 
00209 static bool is_decimal(const char *s)
00210 {
00211   if (*s == '+' || *s == '-') ++s;
00212   if (*s < '1' || *s > '9') return false;
00213   while (*s >= '0' && *s <= '9') ++s;
00214   if (*s == 'l' || *s == 'L') ++s;
00215   return *s == '\0';
00216 }
00217 
00218 static bool is_exponential(const char *s)
00219 {
00220   if (*s == '+' || *s == '-') ++s;
00221   if (*s < '1' || *s > '9') return false;
00222   while (*s >= '0' && *s <= '9') ++s;
00223   if (*s != 'e' && *s != 'E') return false;
00224   ++s;
00225   if (*s < '1' || *s > '9') return false;
00226   while (*s >= '0' && *s <= '9') ++s;
00227   return *s == '\0';
00228 }
00229 
00230 static bool is_hexadecimal(const char *s)
00231 {
00232   if (*s == '+' || *s == '-') ++s;
00233   if (*s != '0') return false;
00234   ++s;
00235   if (*s != 'x' && *s != 'X') return false;
00236   ++s;
00237   if ((*s < '0' || *s > '9') &&
00238       (*s < 'a' || *s > 'f') &&
00239       (*s < 'A' || *s > 'F')) return false;
00240   while ((*s >= '0' && *s <= '9') ||
00241          (*s >= 'a' && *s <= 'f') ||
00242          (*s >= 'A' && *s <= 'F')) ++s;
00243   if (*s == 'l' || *s == 'L') ++s;
00244   return *s == '\0';
00245 }
00246 
00247 static bool is_octal(const char *s)
00248 {
00249   if (*s == '+' || *s == '-') ++s;
00250   if (*s != '0') return false;
00251   while (*s >= '0' && *s <= '7') ++s;
00252   if (*s == 'l' || *s == 'L') ++s;
00253   return *s == '\0';
00254 }
00255 
00256 #else // new implementation, also to be used for operator>> - PVr
00257 
00258 static char rt[4096];
00259 static int rt_pos = 0;
00260 
00261 static char next(const char*& s, vcl_istream** is)
00262 {
00263   if (!is || *s) { char c = *s; if (c) ++rt_pos, ++s; return c; }
00264   if (rt_pos == 4096) return '\0';
00265   (*is)->get(rt[rt_pos]); // read a single byte from istream
00266   if (*s) ++s; // in case s == rt+rt_pos
00267   rt[++rt_pos] = '\0'; return rt[rt_pos-1];
00268 }
00269 
00270 static bool is_decimal(const char* s, vcl_istream** is = 0)
00271 {
00272   rt_pos = 0;
00273   char c = next(s,is);
00274   while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
00275   if (c == '+' || c == '-') c = next(s,is);
00276   if (c < '1' || c > '9') return false;
00277   while (c >= '0' && c <= '9') c = next(s,is);
00278   if (c == 'l' || c == 'L') c = next(s,is);
00279   if (rt_pos > 0) rt[++rt_pos] = '\0';
00280   return is ? true : c == '\0';
00281 }
00282 
00283 static bool is_exponential(const char* s, vcl_istream** is = 0)
00284 {
00285   rt_pos = 0;
00286   char c = next(s,is);
00287   while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
00288   if (c == '+' || c == '-') c = next(s,is);
00289   if (c < '1' || c > '9') return false;
00290   while (c >= '0' && c <= '9') c = next(s,is);
00291   if (c != 'e' && c != 'E') return false;
00292   c = next(s,is);
00293   if (c == '+') c = next(s,is); // no negative exponent!
00294   if (c < '0' || c > '9') return false;
00295   while (c >= '0' && c <= '9') c = next(s,is);
00296   if (rt_pos > 0) rt[++rt_pos] = '\0';
00297   return is ? true : c == '\0';
00298 }
00299 
00300 static bool is_hexadecimal(const char* s, vcl_istream** is = 0)
00301 {
00302   rt_pos = 0;
00303   char c = next(s,is);
00304   while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
00305   if (c == '+' || c == '-') c = next(s,is);
00306   if (c != '0') return false;
00307   c = next(s,is);
00308   if (c != 'x' && c != 'X') return false;
00309   c = next(s,is);
00310   if ((c < '0' || c > '9') &&
00311       (c < 'a' || c > 'f') &&
00312       (c < 'A' || c > 'F')) return false;
00313   while ((c >= '0' && c <= '9') ||
00314          (c >= 'a' && c <= 'f') ||
00315          (c >= 'A' && c <= 'F')) c = next(s,is);
00316   if (c == 'l' || c == 'L') c = next(s,is);
00317   if (rt_pos > 0) rt[++rt_pos] = '\0';
00318   return is ? true : c == '\0';
00319 }
00320 
00321 static bool is_octal(const char* s, vcl_istream** is = 0)
00322 {
00323   rt_pos = 0;
00324   char c = next(s,is);
00325   while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
00326   if (c == '+' || c == '-') c = next(s,is);
00327   if (c != '0') return false;
00328   while (c >= '0' && c <= '7') c = next(s,is);
00329   if (c == 'l' || c == 'L') c = next(s,is);
00330   if (rt_pos > 0) rt[++rt_pos] = '\0';
00331   return is ? true : c == '\0';
00332 }
00333 
00334 static bool is_plus_inf(const char* s, vcl_istream** is = 0)
00335 {
00336   rt_pos = 0;
00337   char c = next(s,is);
00338   while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
00339   if (c == '+') c = next(s,is);
00340   if (c != 'I') return false; c = next(s,is);
00341   if (c != 'n') return false; c = next(s,is);
00342   if (c != 'f') return false; c = next(s,is);
00343   if (c == 'i') c = next(s,is);
00344   if (c == 'n') c = next(s,is);
00345   if (c == 'i') c = next(s,is);
00346   if (c == 't') c = next(s,is);
00347   if (c == 'y') c = next(s,is);
00348   if (rt_pos > 0) rt[++rt_pos] = '\0';
00349   return is ? true : c == '\0';
00350 }
00351 
00352 static bool is_minus_inf(const char* s, vcl_istream** is = 0)
00353 {
00354   rt_pos = 0;
00355   char c = next(s,is);
00356   while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
00357   if (c != '-') return false; c = next(s,is);
00358   if (c != 'I') return false; c = next(s,is);
00359   if (c != 'n') return false; c = next(s,is);
00360   if (c != 'f') return false; c = next(s,is);
00361   if (c == 'i') c = next(s,is);
00362   if (c == 'n') c = next(s,is);
00363   if (c == 'i') c = next(s,is);
00364   if (c == 't') c = next(s,is);
00365   if (c == 'y') c = next(s,is);
00366   if (rt_pos > 0) rt[++rt_pos] = '\0';
00367   return is ? true : c == '\0';
00368 }
00369 
00370 #endif // new implementation - PVr
00371 
00372 //: Creates a vnl_bignum from the character string representation.
00373 
00374 vnl_bignum::vnl_bignum (const char *s)
00375 : count(0), sign(1), data(0)
00376 {
00377   // decimal:     "^ *[-+]?[1-9][0-9]*$"
00378   // exponential: "^ *[-+]?[1-9][0-9]*[eE][+]?[0-9]+$"
00379   // hexadecimal: "^ *[-+]?0[xX][0-9a-fA-F]+$"
00380   // octal:       "^ *[-+]?0[0-7]*$"
00381   // infinity:    "^ *[-+]?Inf(inity)?$"
00382 
00383   if (is_plus_inf(s))
00384     sign=1,count=1,data=new Data[1],data[0]=0;
00385   else if (is_minus_inf(s))
00386     sign=-1,count=1,data=new Data[1],data[0]=0;
00387   else if (is_decimal(s))               // If string is decimal
00388     this->dtoBigNum(s);                 // convert decimal to vnl_bignum
00389   else if (is_exponential(s))           // If string is exponential
00390     this->exptoBigNum(s);               // convert exp. to vnl_bignum
00391   else if (is_hexadecimal(s))           // If string is hex,
00392     this->xtoBigNum(s);                 // convert hex to vnl_bignum
00393   else if (is_octal(s))                 // If string is octal
00394     this->otoBigNum(s);                 // convert octal to vnl_bignum
00395   else {                                // Otherwise
00396     vcl_cerr << "Cannot convert string " << s << " to vnl_bignum\n";
00397   }
00398 }
00399 
00400 //: Reads a vnl_bignum from a stream
00401 
00402 vcl_istream& operator>> (vcl_istream& is, vnl_bignum& x)
00403 {
00404   // decimal:     "^ *[-+]?[1-9][0-9]*$"
00405   // exponential: "^ *[-+]?[1-9][0-9]*[eE][+]?[0-9]+$"
00406   // hexadecimal: "^ *[-+]?0[xX][0-9a-fA-F]+$"
00407   // octal:       "^ *[-+]?0[0-7]*$"
00408   vcl_istream* isp = &is;
00409   rt[0] = '\0';
00410 
00411 //delete[] x.data;                      // Delete existing data
00412   if (is_plus_inf(rt,&isp))
00413     x.sign=1,x.count=1,x.data=new Data[1],x.data[0]=0;
00414   else if (is_minus_inf(rt,&isp))
00415     x.sign=-1,x.count=1,x.data=new Data[1],x.data[0]=0;
00416   if (is_exponential(rt,&isp))          // If input stream string is exponential
00417     x.exptoBigNum(rt);                  // convert exp. to vnl_bignum
00418   else if (is_decimal(rt,&isp))         // If string is decimal
00419     x.dtoBigNum(rt);                    // convert decimal to vnl_bignum
00420   else if (is_hexadecimal(rt,&isp))     // If string is hex,
00421     x.xtoBigNum(rt);                    // convert hex to vnl_bignum
00422   else if (is_octal(rt,&isp))           // If string is octal
00423     x.otoBigNum(rt);                    // convert octal to vnl_bignum
00424   else {                                // Otherwise
00425     vcl_cerr << "Cannot convert string " << rt << " to vnl_bignum\n";
00426     x = 0L;
00427   }
00428   return is; // FIXME - should probably push back read characters to istream
00429 }
00430 
00431 //: Copies the contents of vnl_bignum b.
00432 
00433 vnl_bignum::vnl_bignum (const vnl_bignum& b)
00434 : count(b.count), sign(b.sign)
00435 {
00436   this->data = b.data ? new Data[b.count] : 0;  // Allocate data if necessary
00437   for (Counter i = 0; i < this->count; ++i)     // Copy b data
00438     this->data[i] = b.data[i];
00439 }
00440 
00441 
00442 //: Frees space for vnl_bignum.
00443 
00444 vnl_bignum::~vnl_bignum ()
00445 {
00446   delete [] this->data; this->count = 0;        // Delete any allocated data
00447 }
00448 
00449 //: Copies rhs vnl_bignum to lhs vnl_bignum.
00450 
00451 vnl_bignum& vnl_bignum::operator= (const vnl_bignum& rhs)
00452 {
00453   if (this != &rhs) {                           // Avoid self-assignment
00454     delete[] this->data;                        // Delete existing data
00455     this->count = rhs.count;                    // Copy rhs's count
00456     this->data = rhs.data ? new Data[rhs.count] : 0; // Allocate data if necessary
00457     for (Counter i = 0; i < rhs.count; ++i)     // Copy rhs's data
00458       this->data[i] = rhs.data[i];
00459     this->sign = rhs.sign;                      // Copy rhs's sign
00460   }
00461   return *this;                                 // Return reference
00462 }
00463 
00464 //: Returns the negation of a vnl_bignum.
00465 
00466 vnl_bignum vnl_bignum::operator- () const
00467 {
00468   vnl_bignum neg(*this);
00469   if (neg.count)                // So long as this is non-zero
00470     neg.sign = -neg.sign;       // Flip its sign
00471   return neg;
00472 }
00473 
00474 
00475 //: Prefix increment. Increments a vnl_bignum by 1, and returns it.
00476 
00477 vnl_bignum& vnl_bignum::operator++ ()
00478 {
00479   if (this->is_infinity()) return *this;
00480   if (this->count==0)
00481   {
00482     this->resize(1);
00483     this->data[0] = 1;
00484     this->sign = +1;
00485     return *this;
00486   }
00487 
00488   if (this->sign > 0) increment(*this);
00489   else decrement(*this);
00490 
00491   return *this;
00492 }
00493 
00494 
00495 //: Prefix decrement. Decrements a vnl_bignum by 1, and returns it.
00496 
00497 vnl_bignum& vnl_bignum::operator-- ()
00498 {
00499   if (this->is_infinity()) return *this;
00500   if (this->count==0)
00501   {
00502     this->resize(1);
00503     this->data[0] = 1;
00504     this->sign = -1;
00505     return *this;
00506   }
00507 
00508   if (this->sign < 0) increment(*this);
00509   else decrement(*this);
00510 
00511   return *this;
00512 }
00513 
00514 //: Adds two vnl_bignums, and returns new sum.
00515 
00516 vnl_bignum vnl_bignum::operator+(const vnl_bignum& b) const
00517 {
00518   // Infinity arithmetic:
00519   assert (! b.is_minus_infinity() || ! this->is_plus_infinity() ); // +Inf-Inf
00520   assert (! b.is_plus_infinity() || ! this->is_minus_infinity() ); // -Inf+Inf
00521   if (b.is_infinity()) { return b; }
00522   if (this->is_infinity()) { return *this; }
00523 
00524   vnl_bignum sum;                       // Init sum to zero
00525   if (this->sign == b.sign) {           // If both have same sign
00526     add(*this,b,sum);                   //   Do simple addition
00527     sum.sign = this->sign;              // Attach proper sign
00528   }
00529   else {                                // Else different signs
00530     int mag = magnitude_cmp(*this,b);   // Determine relative sizes
00531     if (mag > 0) {                      // If abs(*this) > abs(b)
00532       subtract(*this,b,sum);            //   sum = *this - b
00533       sum.sign = this->sign;            // Sign of sum follows *this
00534     }
00535     else if (mag < 0) {                 // Else if abs(*this) < abs(b)
00536       subtract(b,*this,sum);            //   sum = b - *this
00537       sum.sign = b.sign;                // Sign of sum follows b
00538     }                                   // (Else abs(*this) == abs(b)
00539   }                                     //   so sum must be zero)
00540   return sum;                           // shallow swap on return
00541 }
00542 
00543 
00544 //: Multiplies this with a vnl_bignum
00545 
00546 vnl_bignum& vnl_bignum::operator*= (const vnl_bignum& b)
00547 {
00548   // Infinity arithmetic:
00549   assert (! b.is_infinity() || this->count != 0 ); // multiplication 0*Inf
00550   assert (! this->is_infinity() || b.count != 0 ); // multiplication Inf*0
00551   if (b.is_infinity()) return (*this) = (this->sign<0 ? -b : b);
00552   if (this->is_infinity()) return (*this) = (b.sign<0 ? -(*this) : *this);
00553 
00554   if (b.count == 0 || this->count == 0)
00555     return (*this)=0L;
00556   vnl_bignum prod;
00557   prod.data = new Data[prod.count = this->count + b.count]; // allocate data for product
00558   for (Counter i = 0; i < b.count; i++)         //   multiply each b "digit"
00559     multiply_aux(*this, b.data[i], prod, i);    //   times b1 and add to total
00560   prod.sign = this->sign * b.sign;              //   determine correct sign
00561   prod.trim();                                  //   trim excess data and ret.
00562   return (*this)=prod;
00563 }
00564 
00565 
00566 //: Divides this by a vnl_bignum
00567 
00568 vnl_bignum& vnl_bignum::operator/= (const vnl_bignum& b)
00569 {
00570   // Infinity arithmetic:
00571   assert (! b.is_infinity() || ! this->is_infinity() ); // division Inf/Inf
00572   if (b.is_infinity()) return (*this)=0L;
00573   if (this->is_infinity()) return (*this) = (b.sign<0 ? -(*this) : *this);
00574   assert (b.count!=0 || this->count != 0); // division 0/0
00575   if (b.count == 0)
00576     return (*this) = (this->sign < 0 ? vnl_bignum("-Inf") : vnl_bignum("+Inf"));
00577 
00578   vnl_bignum quot, r;          // Quotient and remainder
00579   divide(*this,b,quot,r);      // Call divide fn
00580   return (*this) = quot;
00581 }
00582 
00583 //: Divides this by a vnl_bignum and replaces this by remainder.
00584 
00585 vnl_bignum& vnl_bignum::operator%= (const vnl_bignum& b)
00586 {
00587   // Infinity arithmetic:
00588   assert (! b.is_infinity() || ! this->is_infinity() ); // division Inf/Inf
00589   if (b.is_infinity()) return *this;                    // remainder of x/Inf is x.
00590   if (this->is_infinity()) return (*this) = 0L;         // convention: remainder is 0
00591   assert (b.count!=0 || this->count != 0);              // division 0/0
00592   if (b.count == 0) return (*this) = 0L;                // convention: remainder is 0
00593 
00594   vnl_bignum remain, q;        // Quotient and remainder
00595   divide(*this,b,q,remain);    // divide by b and save remainder
00596   return (*this) = remain;     // shallow swap on return
00597 }
00598 
00599 
00600 //: Shifts bignum to the left l digits.
00601 
00602 vnl_bignum vnl_bignum::operator<< (int l) const
00603 {
00604   // Infinity arithmetic:
00605   if (this->is_infinity()) return *this;
00606 
00607   if (l == 0 || *this == 0L)            // if either arg is zero
00608     return *this;
00609   if (l < 0)                            // if shift amt is negative
00610     return right_shift(*this,-l);       //   do an actual right shift
00611   else                                  // otherwise
00612     return left_shift(*this,l);         //   do a left shift
00613 }
00614 
00615 
00616 //: Shifts bignum to the right l digits.
00617 
00618 vnl_bignum vnl_bignum::operator>> (int l) const
00619 {
00620   // Infinity arithmetic:
00621   if (this->is_infinity()) return *this;
00622 
00623   if (l == 0 || *this == 0L)            // if either arg is zero
00624     return *this;
00625   if (l < 0)                            // if shift amt is negative
00626     return left_shift(*this,-l);        //   do an actual left shift
00627   else                                  // else
00628     return right_shift(*this,l);        //   do a right shift
00629 }
00630 
00631 
00632 //: Two vnl_bignums are equal if and only if they have the same integer representation.
00633 
00634 bool vnl_bignum::operator== (const vnl_bignum& rhs) const
00635 {
00636   if (this != &rhs) {                           // Check address
00637     if (this->sign != rhs.sign) return false;   // Different sign implies !=
00638     if (this->count != rhs.count) return false; // Different size implies !=
00639     for (Counter i = 0; i < this->count; i++)   // Each data element the same?
00640       if (this->data[i] != rhs.data[i]) return false; // No. Return !=
00641   }
00642   return true;                                    // Yes. Return ==
00643 }
00644 
00645 
00646 //: Compares two vnl_bignums.
00647 
00648 bool vnl_bignum::operator< (const vnl_bignum& rhs) const
00649 {
00650   if (this->sign < rhs.sign) return true;       // Different signs?
00651   if (this->sign > rhs.sign) return false;
00652   if (this->sign == 1)                          // Both signs == 1
00653     return magnitude_cmp(*this,rhs) < 0;        // this must be smaller
00654   else                                          // Both signs == -1
00655     return magnitude_cmp(*this,rhs) > 0;        // this must be larger
00656 }
00657 
00658 
00659 //: Formatted output for bignum.
00660 
00661 vcl_ostream& operator<< (vcl_ostream& os, const vnl_bignum& b)
00662 {
00663   vnl_bignum d = b;                     // Copy the input vnl_bignum
00664   if (d.sign == -1) {                   // If it's negative
00665     os << '-';                          //   Output leading minus sign
00666     d.sign = 1;                         //   Make d positive for divide
00667   }
00668   if (d.is_infinity()) return os<<"Inf";
00669   vnl_bignum q,r;                       // Temp quotient and remainder
00670   char *cbuf = new char[5 * (b.count+1)];   // Temp character buffer
00671   Counter i = 0;
00672   do {                                  // repeat:
00673     divide(d,10L,q,r);                  //   Divide vnl_bignum by ten
00674     cbuf[i++] = char(long(r) + '0');    //   Get one's digit
00675     d = q;                              //   Then discard one's digit
00676     q = r = 0L;                         //   Prep for next divide
00677   } while (d != 0L);                    // until no more one's digits
00678   do {                                  // repeat;
00679     os << cbuf[--i];                    //   output char buf in reverse
00680   } while (i);                          // until no more chars
00681   delete [] cbuf;                       // delete temp char buf
00682   return os;                            // return output stream
00683 }
00684 
00685 //: Convert the number to a decimal representation in a string.
00686 vcl_string& vnl_bignum_to_string (vcl_string& s, const vnl_bignum& b)
00687 {
00688   s.erase();
00689   vcl_string::size_type insert_point = 0; // keep record of location of first number.
00690 
00691   vnl_bignum d = b;                     // Copy the input vnl_bignum
00692   if (d.sign == -1) {                   // If it's negative
00693     s.insert(insert_point,"-");         //   Output leading minus sign
00694     d.sign = 1;                         //   Make d positive for divide
00695     ++insert_point;                     // keep record of location of first number.
00696   }
00697   if (d.is_infinity()) return s+="Inf";
00698   vnl_bignum q,r;                       // Temp quotient and remainder
00699   do {                                  // repeat:
00700     divide(d,10L,q,r);                  //   Divide vnl_bignum by ten
00701     s.insert(insert_point, 1, char('0'+long(r))); //   Get one's digit, and insert it at head.
00702     d = q;                              //   Then discard one's digit
00703     q = r = 0L;                         //   Prep for next divide
00704   } while (d != 0L);                    // until no more one's digits
00705   return s;
00706 }
00707 
00708 //: Convert the number from a decimal representation in a string.
00709 vnl_bignum& vnl_bignum_from_string (vnl_bignum& b, const vcl_string& s)
00710 {
00711   // decimal:     "^ *[-+]?[1-9][0-9]*$"
00712   // Infinity:    "^ *[-+]?Inf(inity)?$"
00713 
00714   if (is_plus_inf(s.c_str()))
00715     b=vnl_bignum("+Inf");
00716   else if (is_minus_inf(s.c_str()))
00717     b=vnl_bignum("-Inf");
00718   else
00719     b.dtoBigNum(s.c_str());             // convert decimal to vnl_bignum
00720   return b;
00721 }
00722 
00723 
00724 //: Implicit conversion from a vnl_bignum to a short.
00725 vnl_bignum::operator short () const
00726 {
00727   short s = 0;
00728   for (Counter i = this->count; i > 0; )
00729     s = short(s*0x10000 + this->data[--i]);
00730   return this->sign*s;
00731 }
00732 
00733 
00734 //: Implicit conversion from a vnl_bignum to an int.
00735 vnl_bignum::operator int () const
00736 {
00737   int j = 0;
00738   for (Counter i = this->count; i > 0; )
00739     j = int(j*0x10000 + this->data[--i]);
00740   return this->sign*j;
00741 }
00742 
00743 
00744 //: Implicit conversion from a vnl_bignum to a long.
00745 vnl_bignum::operator long () const
00746 {
00747   long l = 0;
00748   for (Counter i = this->count; i > 0; )
00749     l = l*0x10000L + this->data[--i];
00750   return this->sign*l;
00751 }
00752 
00753 
00754 //: Implicit conversion from a vnl_bignum to a float.
00755 vnl_bignum::operator float () const
00756 {
00757   float f = 0.0f;
00758   for (Counter i = this->count; i > 0; )
00759     f = f*0x10000 + this->data[--i];
00760   if (this->is_infinity()) f = vcl_numeric_limits<float>::infinity();
00761   return this->sign*f;
00762 }
00763 
00764 
00765 //: Implicit conversion from a vnl_bignum to a double.
00766 
00767 vnl_bignum::operator double () const
00768 {
00769   double d = 0.0;
00770   for (Counter i = this->count; i > 0; )
00771     d = d*0x10000 + this->data[--i];
00772   if (this->is_infinity()) d = vcl_numeric_limits<double>::infinity();
00773   return this->sign*d;
00774 }
00775 
00776 //: Implicit conversion from a vnl_bignum to a long double.
00777 
00778 vnl_bignum::operator long double () const
00779 {
00780   long double d = 0.0;
00781   for (Counter i = this->count; i > 0; )
00782     d = d*0x10000 + this->data[--i];
00783   if (this->is_infinity()) d = vcl_numeric_limits<long double>::infinity();
00784   return this->sign*d;
00785 }
00786 
00787 //: dump the contents of a vnl_bignum to a stream, default cout.
00788 
00789 void vnl_bignum::dump (vcl_ostream& os) const
00790 {
00791   os << "{count=" << this->count       // output count field
00792      << ", sign=" << this->sign        // output sign field
00793      << ", data=" << this->data        // output data pointer
00794      << ", value=" << *this
00795      << ", {";
00796   // format string == "%04X%s" or "%02X%s", etc.
00797   //  static char format_str[10] =
00798   //    {'%','0',char(2*2 + '0'),'X','%','s'};
00799   //  format_str[2] = char(2*2 + '0');
00800   if (this->count > 0) { // output data array
00801     os << vcl_hex;
00802     for (Counter i = this->count; i > 1; --i)
00803       os << (this->data[i - 1]) << ',';
00804     os << (this->data[0]) << vcl_dec;
00805   }
00806   os << "}}\n";                         // close brackets
00807 }
00808 
00809 
00810 //: Converts decimal string to a vnl_bignum.
00811 
00812 int vnl_bignum::dtoBigNum (const char *s)
00813 {
00814   this->resize(0); sign = 1;            // Reset number to 0.
00815   Counter len = 0;                      // No chars converted yet
00816   while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s; // skip whitespace
00817   if (s[0] == '-' || s[0] == '+') len++;// Skip over leading +,-
00818   while (vcl_isdigit(s[len])) {         // If current char is digit
00819     (*this) = ((*this) * 10L) +         // Shift vnl_bignum left a decimal
00820       vnl_bignum(long(s[len++] - '0')); // digit and add new digit
00821   }
00822   if (s[0] == '-') this->sign = -1;     // If s had leading -, note it
00823   return len;                           // Return # of chars processed
00824 }
00825 
00826 //: convert exponential string to a vnl_bignum
00827 
00828 void vnl_bignum::exptoBigNum (const char *s)
00829 {
00830   while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s; // skip whitespace
00831   Counter pos = this->dtoBigNum(s) + 1; // Convert the base, skip [eE]
00832   long pow = vcl_atol(s + pos);         // Convert the exponent to long
00833   while (pow-- > 0)                     // Raise vnl_bignum to the given
00834     *this = (*this) * 10L;              // power
00835 }
00836 
00837 
00838 //: convert hex character to integer hex value (ASCII or EBCDIC)
00839 // - Inputs:  character representation of a hex number
00840 // - Outputs: integer value of the hex number
00841 
00842 unsigned int ctox (int c)
00843 {
00844   if ('0' <= c && c <= '9')
00845     return c - '0';
00846   if ('a' <= c && c <= 'f')
00847     return c - 'a' + 10;
00848   return c - 'A' + 10;
00849 }
00850 
00851 //: convert hex string to vnl_bignum
00852 
00853 void vnl_bignum::xtoBigNum (const char *s)
00854 {
00855   this->resize(0); sign = 1;            // Reset number to 0.
00856   while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s; // skip whitespace
00857   Counter size = vcl_strlen(s);
00858   Counter len = 2;                      // skip leading "0x"
00859   while (len < size) {                  // While there are more chars
00860     (*this) = ((*this) * 16L) +         // Shift vnl_bignum left one hex
00861       vnl_bignum(long(ctox(s[len++]))); //   digit and add next digit
00862   }
00863 }
00864 
00865 
00866 //: convert octal string to vnl_bignum
00867 
00868 void vnl_bignum::otoBigNum (const char *s)
00869 {
00870   this->resize(0); sign = 1;           // Reset number to 0.
00871   while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s; // skip whitespace
00872   Counter size = vcl_strlen(s);
00873   Counter len = 0;                      // No chars converted yet
00874   while (len < size) {                  // While there are more chars
00875     (*this) = ((*this) * 8L) +          // Shift vnl_bignum left 1 oct dig.
00876       vnl_bignum(long(s[len++] - '0')); // Add next character value
00877   }
00878 }
00879 
00880 //: change the data allotment for a vnl_bignum
00881 
00882 void vnl_bignum::resize (short new_count)
00883 {
00884   assert(new_count >= 0);
00885   if (new_count == this->count) return;
00886   Data *new_data = (new_count > 0 ? new Data[new_count] : 0); // Allocate data if necessary
00887 
00888   if (this->count <= new_count) {       // Copy old data into new
00889     short i = 0;
00890     for (; i < this->count; i++)
00891       new_data[i] = this->data[i];
00892     for (; i < new_count; i++)
00893       new_data[i] = 0;
00894   }
00895   else {
00896     for (short i = 0; i < new_count; i++)
00897       new_data[i] = this->data[i];
00898   }
00899 
00900   delete [] this->data;                 // Get rid of old data
00901   this->data = new_data;                // Point to new data
00902   this->count = new_count;              // Save new count
00903 }
00904 
00905 
00906 //: trim non-infinite vnl_bignum of excess data allotment
00907 
00908 vnl_bignum& vnl_bignum::trim ()
00909 {
00910   Counter i = this->count;
00911   for (; i > 0; i--)                    // Skip over high-order words
00912     if (this->data[i - 1] != 0) break;  //   that are zero
00913   if (i < this->count) {                // If there are some such words
00914     this->count = i;                    // Update the count
00915     Data *new_data = (i > 0 ? new Data[i] : 0); // Allocate data if necessary
00916     for (; i > 0; i--)                  // Copy old data into new
00917       new_data[i - 1] = this->data[i - 1];
00918     delete [] this->data;               // Delete old data
00919     this->data = new_data;              // Point to new data
00920   }
00921   return *this;                         // return reference to vnl_bignum
00922 }
00923 
00924 
00925 //: add two non-infinite vnl_bignum values and save their sum
00926 
00927 void add (const vnl_bignum& b1, const vnl_bignum& b2, vnl_bignum& sum)
00928 {
00929   const vnl_bignum *bmax, *bmin;        // Determine which of the two
00930   if (b1.count >= b2.count) {           //   addends has the most
00931     bmax = &b1;                         //   data.
00932     bmin = &b2;
00933   }
00934   else {
00935     bmax = &b2;
00936     bmin = &b1;
00937   }
00938 //delete[] sum.data;                    // Delete existing data
00939   sum.data = (sum.count = bmax->count) > 0 ?   // Allocate data for their sum
00940              new Data[sum.count] : 0;
00941   unsigned long temp, carry = 0;
00942   Counter i = 0;
00943   while (i < bmin->count) {             // Add, element by element.
00944     // Add both elements and carry
00945     temp = (unsigned long)b1.data[i] + (unsigned long)b2.data[i] + carry;
00946     carry = temp/0x10000L;              // keep track of the carry
00947     sum.data[i] = Data(temp);           // store sum
00948     i++;                                // go to next element
00949   }
00950   while (i < bmax->count) {             // bmin has no more elements
00951     temp = bmax->data[i] + carry;       // propagate the carry through
00952     carry = temp/0x10000L;              // the rest of bmax's elements
00953     sum.data[i] = Data(temp);           // store sum
00954     i++;
00955   }
00956   if (carry) {                          // if carry left over
00957     sum.resize(bmax->count + 1);        //   allocate another word
00958     sum.data[bmax->count] = 1;          //   save the carry in it
00959   }
00960 }
00961 
00962 //: Add 1 to bnum (unsigned, non-infinite)
00963 void increment (vnl_bignum& bnum)
00964 {
00965   Counter i = 0;
00966   unsigned long carry = 1;
00967   while (i < bnum.count && carry) {             // increment, element by element.
00968     unsigned long temp = (unsigned long)bnum.data[i] + carry;
00969     carry = temp/0x10000L;
00970     bnum.data[i] = (Data)temp;
00971     ++i;
00972   }
00973   if (carry)
00974   {
00975     bnum.resize(bnum.count+1);
00976     bnum.data[bnum.count-1] = 1;
00977   }
00978 }
00979 
00980 
00981 //: subtract bmin from bmax (unsigned, non-infinite), result in diff
00982 
00983 void subtract (const vnl_bignum& bmax, <