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_cstdlib.h>   // for vcl_atol
00007 #include <vcl_cstring.h>   // for vcl_strlen
00008 #include <vcl_cmath.h>     // for vcl_fmod
00009 #include <vcl_algorithm.h> // for vcl_copy
00010 #include <vcl_vector.h>
00011 #include <vcl_cassert.h>
00012 #include <vcl_iostream.h>
00013 #include <vcl_limits.h>
00014 #include <vnl/vnl_math.h> // for vnl_math_isfinite(double)
00015 
00016 typedef unsigned short Counter;
00017 typedef unsigned short Data;
00018 
00019 //: Creates a zero vnl_bignum.
00020 
00021 vnl_bignum::vnl_bignum()
00022 : count(0), sign(1), data(0)
00023 {
00024 }
00025 
00026 //: Creates a vnl_bignum from a long integer.
00027 
00028 vnl_bignum::vnl_bignum(long l)
00029 : count(0), sign(1), data(0)
00030 {
00031   if (l < 0) {                  // Get correct sign
00032     l = -l;                     // Get absolute value of l
00033     this->sign = -1;
00034   }
00035   Data buf[sizeof(l)];          // Temp buffer to store l in
00036   Counter i = 0;                // buffer index
00037   while (l) {                   // While more bits in l
00038     assert(i < sizeof(l));      // no more buffer space
00039     buf[i] = Data(l);           // Peel off lower order bits
00040     l >>= 16;   // Shift next bits into place
00041     ++i;
00042   }
00043   if (i > 0)
00044     this->data = new Data[this->count=i]; // Allocate permanent data
00045 
00046   while (i--)     // Save buffer into perm. data
00047     this->data[i] = buf[i];
00048 }
00049 
00050 //: Creates a vnl_bignum from an integer.
00051 
00052 vnl_bignum::vnl_bignum(int l)
00053 : count(0), sign(1), data(0)
00054 {
00055   if (l < 0) {                  // Get correct sign
00056     l = -l;                     // Get absolute value of l
00057     this->sign = -1;
00058   }
00059   Data buf[sizeof(l)];          // Temp buffer to store l in
00060   Counter i = 0;                // buffer index
00061   while (l) {                   // While more bits in l
00062     assert(i < sizeof(l));      // no more buffer space
00063     buf[i] = Data(l);           // Peel off lower order bits
00064     l >>= 16;                   // Shift next bits into place
00065     i++;
00066   }
00067   if (i > 0)
00068     this->data = new Data[this->count=i]; // Allocate permanent data
00069 
00070   while (i--)     // Save buffer into perm. data
00071     this->data[i] = buf[i];
00072 }
00073 
00074 //: Creates a vnl_bignum from an unsigned long integer.
00075 
00076 vnl_bignum::vnl_bignum(unsigned long l)
00077 : count(0), sign(1), data(0)
00078 {
00079   Data buf[sizeof(l)];          // Temp buffer to store l in
00080   Counter i = 0;                // buffer index
00081   while (l) {                   // While more bits in l
00082     assert(i < sizeof(l));      // no more buffer space
00083     buf[i] = Data(l);           // Peel off lower order bits
00084     l >>= 16;   // Shift next bits into place
00085     i++;
00086   }
00087   if (i > 0)
00088     this->data = new Data[this->count=i]; // Allocate permanent data
00089 
00090   while (i--)     // Save buffer into perm. data
00091     this->data[i] = buf[i];
00092 }
00093 
00094 //: Creates a vnl_bignum from an unsigned integer.
00095 
00096 vnl_bignum::vnl_bignum(unsigned int l)
00097 : count(0), sign(1), data(0)
00098 {
00099   Data buf[sizeof(l)];          // Temp buffer to store l in
00100   Counter i = 0;                // buffer index
00101   while (l) {                   // While more bits in l
00102     assert(i < sizeof(l));      // no more buffer space
00103     buf[i] = Data(l);           // Peel off lower order bits
00104     l >>= 16;   // Shift next bits into place
00105     i++;
00106   }
00107   if (i > 0)
00108     this->data = new Data[this->count=i]; // Allocate permanent data
00109 
00110   while (i--)     // Save buffer into perm. data
00111     this->data[i] = buf[i];
00112 }
00113 
00114 //: Creates a vnl_bignum from a single-precision floating point number.
00115 
00116 vnl_bignum::vnl_bignum(float f)
00117 : count(0), sign(1), data(0)
00118 {
00119   double d = f;
00120   if (d < 0.0) {                // Get sign of d
00121     d = -d;                     // Get absolute value of d
00122     this->sign = -1;
00123   }
00124 
00125   if (!vnl_math_isfinite(d)) {
00126     // Infinity is represented as: count=1, data[0]=0.
00127     // This is an otherwise unused representation, since 0 is represented as count=0.
00128     this->count = 1;
00129     this->data = new Data[1];
00130     this->data[0] = 0;
00131   }
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 = (unsigned short)(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   }
00163   else if (d >= 1.0) {
00164     // Note: 0x10000L == 1 >> 16: the (assumed) size of unsigned short is 16 bits.
00165     vcl_vector<Data> buf;
00166     while (d >= 1.0) {
00167       buf.push_back( Data(vcl_fmod(d,0x10000L)) );  // Get next data "digit" from d
00168       d /= 0x10000L;                                // Shift d right 1 data "digit"
00169     }
00170     // Allocate and copy into permanent buffer
00171     this->data = buf.size()>0 ? new Data[buf.size()] : 0;
00172     this->count = (unsigned short)(buf.size());
00173     vcl_copy( buf.begin(), buf.end(), data );
00174   }
00175 }
00176 
00177 //: Creates a vnl_bignum from a "long double" floating point number.
00178 
00179 vnl_bignum::vnl_bignum(long double d)
00180 : count(0), sign(1), data(0)
00181 {
00182   if (d < 0.0) {                // Get sign of d
00183     d = -d;                     // Get absolute value of d
00184     this->sign = -1;
00185   }
00186 
00187   if (!vnl_math_isfinite(d)) {
00188     // Infinity is represented as: count=1, data[0]=0.
00189     // This is an otherwise unused representation, since 0 is represented as count=0.
00190     this->count = 1;
00191     this->data = new Data[1];
00192     this->data[0] = 0;
00193   }
00194   else if (d >= 1.0) {
00195     // Note: 0x10000L == 1 >> 16: the (assumed) size of unsigned short is 16 bits.
00196     vcl_vector<Data> buf;
00197     while (d >= 1.0) {
00198       buf.push_back( Data(vcl_fmod(d,0x10000L)) );  // Get next data "digit" from d
00199       d /= 0x10000L;                                // Shift d right 1 data "digit"
00200     }
00201     // Allocate and copy into permanent buffer
00202     this->data = (buf.size()>0 ? new Data[buf.size()] : 0);
00203     this->count = (unsigned short)(buf.size());
00204     vcl_copy( buf.begin(), buf.end(), data );
00205   }
00206 }
00207 
00208 #if 0 // Old, original Texas Instruments implementation - PVr
00209       // Left here for ducumentation purposes only!
00210 
00211 static bool is_decimal(const char *s)
00212 {
00213   if (*s == '+' || *s == '-') ++s;
00214   if (*s < '1' || *s > '9') return false;
00215   while (*s >= '0' && *s <= '9') ++s;
00216   if (*s == 'l' || *s == 'L') ++s;
00217   return *s == '\0';
00218 }
00219 
00220 static bool is_exponential(const char *s)
00221 {
00222   if (*s == '+' || *s == '-') ++s;
00223   if (*s < '1' || *s > '9') return false;
00224   while (*s >= '0' && *s <= '9') ++s;
00225   if (*s != 'e' && *s != 'E') return false;
00226   ++s;
00227   if (*s < '1' || *s > '9') return false;
00228   while (*s >= '0' && *s <= '9') ++s;
00229   return *s == '\0';
00230 }
00231 
00232 static bool is_hexadecimal(const char *s)
00233 {
00234   if (*s == '+' || *s == '-') ++s;
00235   if (*s != '0') return false;
00236   ++s;
00237   if (*s != 'x' && *s != 'X') return false;
00238   ++s;
00239   if ((*s < '0' || *s > '9') &&
00240       (*s < 'a' || *s > 'f') &&
00241       (*s < 'A' || *s > 'F')) return false;
00242   while ((*s >= '0' && *s <= '9') ||
00243          (*s >= 'a' && *s <= 'f') ||
00244          (*s >= 'A' && *s <= 'F')) ++s;
00245   if (*s == 'l' || *s == 'L') ++s;
00246   return *s == '\0';
00247 }
00248 
00249 static bool is_octal(const char *s)
00250 {
00251   if (*s == '+' || *s == '-') ++s;
00252   if (*s != '0') return false;
00253   while (*s >= '0' && *s <= '7') ++s;
00254   if (*s == 'l' || *s == 'L') ++s;
00255   return *s == '\0';
00256 }
00257 
00258 #else // new implementation, also to be used for operator>> - PVr
00259 
00260 static char rt[4096];
00261 static int rt_pos = 0;
00262 
00263 static char next(const char*& s, vcl_istream** is)
00264 {
00265   if (!is || *s) { char c = *s; if (c) ++rt_pos, ++s; return c; }
00266   if (rt_pos == 4096) return '\0';
00267   (*is)->get(rt[rt_pos]); // read a single byte from istream
00268   if (*s) ++s; // in case s == rt+rt_pos
00269   rt[++rt_pos] = '\0'; return rt[rt_pos-1];
00270 }
00271 
00272 static bool is_decimal(const char* s, vcl_istream** is = 0)
00273 {
00274   rt_pos = 0;
00275   char c = next(s,is);
00276   while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
00277   if (c == '+' || c == '-') c = next(s,is);
00278   if (c < '1' || c > '9') return false;
00279   while (c >= '0' && c <= '9') c = next(s,is);
00280   if (c == 'l' || c == 'L') c = next(s,is);
00281   if (rt_pos > 0) rt[++rt_pos] = '\0';
00282   return is ? true : c == '\0';
00283 }
00284 
00285 static bool is_exponential(const char* s, vcl_istream** is = 0)
00286 {
00287   rt_pos = 0;
00288   char c = next(s,is);
00289   while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
00290   if (c == '+' || c == '-') c = next(s,is);
00291   if (c < '1' || c > '9') return false;
00292   while (c >= '0' && c <= '9') c = next(s,is);
00293   if (c != 'e' && c != 'E') return false;
00294   c = next(s,is);
00295   if (c == '+') c = next(s,is); // no negative exponent!
00296   if (c < '0' || c > '9') return false;
00297   while (c >= '0' && c <= '9') c = next(s,is);
00298   if (rt_pos > 0) rt[++rt_pos] = '\0';
00299   return is ? true : c == '\0';
00300 }
00301 
00302 static bool is_hexadecimal(const char* s, vcl_istream** is = 0)
00303 {
00304   rt_pos = 0;
00305   char c = next(s,is);
00306   while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
00307   if (c == '+' || c == '-') c = next(s,is);
00308   if (c != '0') return false;
00309   c = next(s,is);
00310   if (c != 'x' && c != 'X') return false;
00311   c = next(s,is);
00312   if ((c < '0' || c > '9') &&
00313       (c < 'a' || c > 'f') &&
00314       (c < 'A' || c > 'F')) return false;
00315   while ((c >= '0' && c <= '9') ||
00316          (c >= 'a' && c <= 'f') ||
00317          (c >= 'A' && c <= 'F')) c = next(s,is);
00318   if (c == 'l' || c == 'L') c = next(s,is);
00319   if (rt_pos > 0) rt[++rt_pos] = '\0';
00320   return is ? true : c == '\0';
00321 }
00322 
00323 static bool is_octal(const char* s, vcl_istream** is = 0)
00324 {
00325   rt_pos = 0;
00326   char c = next(s,is);
00327   while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
00328   if (c == '+' || c == '-') c = next(s,is);
00329   if (c != '0') return false;
00330   while (c >= '0' && c <= '7') c = next(s,is);
00331   if (c == 'l' || c == 'L') c = next(s,is);
00332   if (rt_pos > 0) rt[++rt_pos] = '\0';
00333   return is ? true : c == '\0';
00334 }
00335 
00336 static bool is_plus_inf(const char* s, vcl_istream** is = 0)
00337 {
00338   rt_pos = 0;
00339   char c = next(s,is);
00340   while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
00341   if (c == '+') c = next(s,is);
00342   if (c != 'I') return false; c = next(s,is);
00343   if (c != 'n') return false; c = next(s,is);
00344   if (c != 'f') return false; c = next(s,is);
00345   if (c == 'i') c = next(s,is);
00346   if (c == 'n') c = next(s,is);
00347   if (c == 'i') c = next(s,is);
00348   if (c == 't') c = next(s,is);
00349   if (c == 'y') c = next(s,is);
00350   if (rt_pos > 0) rt[++rt_pos] = '\0';
00351   return is ? true : c == '\0';
00352 }
00353 
00354 static bool is_minus_inf(const char* s, vcl_istream** is = 0)
00355 {
00356   rt_pos = 0;
00357   char c = next(s,is);
00358   while (c == ' ' || c == '\t' || c == '\n' || c == '\r') c = next(s,is);
00359   if (c != '-') return false; c = next(s,is);
00360   if (c != 'I') return false; c = next(s,is);
00361   if (c != 'n') return false; c = next(s,is);
00362   if (c != 'f') return false; c = next(s,is);
00363   if (c == 'i') c = next(s,is);
00364   if (c == 'n') c = next(s,is);
00365   if (c == 'i') c = next(s,is);
00366   if (c == 't') c = next(s,is);
00367   if (c == 'y') c = next(s,is);
00368   if (rt_pos > 0) rt[++rt_pos] = '\0';
00369   return is ? true : c == '\0';
00370 }
00371 
00372 #endif // new implementation - PVr
00373 
00374 //: Creates a vnl_bignum from the character string representation.
00375 
00376 vnl_bignum::vnl_bignum(const char *s)
00377 : count(0), sign(1), data(0)
00378 {
00379   // decimal:     "^ *[-+]?[1-9][0-9]*$"
00380   // exponential: "^ *[-+]?[1-9][0-9]*[eE][+]?[0-9]+$"
00381   // hexadecimal: "^ *[-+]?0[xX][0-9a-fA-F]+$"
00382   // octal:       "^ *[-+]?0[0-7]*$"
00383   // infinity:    "^ *[-+]?Inf(inity)?$"
00384 
00385   if (is_plus_inf(s))
00386     count=1,data=new Data[1],data[0]=0;
00387   else if (is_minus_inf(s))
00388     sign=-1,count=1,data=new Data[1],data[0]=0;
00389   else if (is_decimal(s))               // If string is decimal
00390     this->dtoBigNum(s);                 // convert decimal to vnl_bignum
00391   else if (is_exponential(s))           // If string is exponential
00392     this->exptoBigNum(s);               // convert exp. to vnl_bignum
00393   else if (is_hexadecimal(s))           // If string is hex,
00394     this->xtoBigNum(s);                 // convert hex to vnl_bignum
00395   else if (is_octal(s))                 // If string is octal
00396     this->otoBigNum(s);                 // convert octal to vnl_bignum
00397   else {                                // Otherwise
00398     vcl_cerr << "Cannot convert string " << s << " to vnl_bignum\n";
00399   }
00400 }
00401 
00402 //: Reads a vnl_bignum from a stream
00403 
00404 vcl_istream& operator>>(vcl_istream& is, vnl_bignum& x)
00405 {
00406   // decimal:     "^ *[-+]?[1-9][0-9]*$"
00407   // exponential: "^ *[-+]?[1-9][0-9]*[eE][+]?[0-9]+$"
00408   // hexadecimal: "^ *[-+]?0[xX][0-9a-fA-F]+$"
00409   // octal:       "^ *[-+]?0[0-7]*$"
00410   vcl_istream* isp = &is;
00411   rt[0] = '\0';
00412 
00413   x = 0L;
00414   if (is_plus_inf(rt,&isp))
00415     x.sign=1,x.count=1,x.data=new Data[1],x.data[0]=0;
00416   else if (is_minus_inf(rt,&isp))
00417     x.sign=-1,x.count=1,x.data=new Data[1],x.data[0]=0;
00418   else if (is_exponential(rt,&isp))     // If input stream string is exponential
00419     x.exptoBigNum(rt);                  // convert exp. to vnl_bignum
00420   else if (is_decimal(rt,&isp))         // If string is decimal
00421     x.dtoBigNum(rt);                    // convert decimal to vnl_bignum
00422   else if (is_hexadecimal(rt,&isp))     // If string is hex,
00423     x.xtoBigNum(rt);                    // convert hex to vnl_bignum
00424   else if (is_octal(rt,&isp))           // If string is octal
00425     x.otoBigNum(rt);                    // convert octal to vnl_bignum
00426   else {                                // Otherwise
00427     vcl_cerr << "Cannot convert string " << rt << " to vnl_bignum\n";
00428   }
00429   return is; // FIXME - should probably push back read characters to istream
00430 }
00431 
00432 //: Copies the contents of vnl_bignum b.
00433 
00434 vnl_bignum::vnl_bignum(const vnl_bignum& b)
00435 : count(b.count), sign(b.sign)
00436 {
00437   if (b.data) {
00438     this->data = new Data[b.count];            // Allocate data
00439     for (Counter i = 0; i < this->count; ++i)  // Copy b data
00440       this->data[i] = b.data[i];
00441   }
00442   else {
00443     this->data = 0;
00444   }
00445 }
00446 
00447 //: Frees space for vnl_bignum.
00448 
00449 vnl_bignum::~vnl_bignum()
00450 {
00451   delete [] this->data; this->count = 0;        // Delete any allocated data
00452 }
00453 
00454 //: Copies rhs vnl_bignum to lhs vnl_bignum.
00455 
00456 vnl_bignum& vnl_bignum::operator=(const vnl_bignum& rhs)
00457 {
00458   if (this != &rhs) {                           // Avoid self-assignment
00459     delete[] this->data;                        // Delete existing data
00460     this->count = rhs.count;                    // Copy rhs's count
00461     this->data = rhs.data ? new Data[rhs.count] : 0; // Allocate data if necessary
00462     for (Counter i = 0; i < rhs.count; ++i)     // Copy rhs's data
00463       this->data[i] = rhs.data[i];
00464     this->sign = rhs.sign;                      // Copy rhs's sign
00465   }
00466   return *this;                                 // Return reference
00467 }
00468 
00469 //: Returns the negation of a vnl_bignum.
00470 
00471 vnl_bignum vnl_bignum::operator-() const
00472 {
00473   vnl_bignum neg(*this);
00474   if (neg.count)                // So long as this is non-zero
00475     neg.sign = -neg.sign;       // Flip its sign
00476   return neg;
00477 }
00478 
00479 //: Prefix increment. Increments a vnl_bignum by 1, and returns it.
00480 
00481 vnl_bignum& vnl_bignum::operator++()
00482 {
00483   if (this->is_infinity()) return *this;
00484   if (this->count==0)
00485   {
00486     this->resize(1);
00487     this->data[0] = 1;
00488     this->sign = +1;
00489     return *this;
00490   }
00491 
00492   if (this->sign > 0) increment(*this);
00493   else decrement(*this);
00494 
00495   return *this;
00496 }
00497 
00498 //: Prefix decrement. Decrements a vnl_bignum by 1, and returns it.
00499 
00500 vnl_bignum& vnl_bignum::operator--()
00501 {
00502   if (this->is_infinity()) return *this;
00503   if (this->count==0)
00504   {
00505     this->resize(1);
00506     this->data[0] = 1;
00507     this->sign = -1;
00508     return *this;
00509   }
00510 
00511   if (this->sign < 0) increment(*this);
00512   else decrement(*this);
00513 
00514   return *this;
00515 }
00516 
00517 //: Adds two vnl_bignums, and returns new sum.
00518 
00519 vnl_bignum vnl_bignum::operator+(const vnl_bignum& b) const
00520 {
00521   // Infinity arithmetic:
00522   assert (! b.is_minus_infinity() || ! this->is_plus_infinity() ); // +Inf-Inf
00523   assert (! b.is_plus_infinity() || ! this->is_minus_infinity() ); // -Inf+Inf
00524   if (b.is_infinity()) { return b; }
00525   if (this->is_infinity()) { return *this; }
00526 
00527   vnl_bignum sum;                       // Init sum to zero
00528   if (this->sign == b.sign) {           // If both have same sign
00529     add(*this,b,sum);                   //   Do simple addition
00530     sum.sign = this->sign;              // Attach proper sign
00531   }
00532   else {                                // Else different signs
00533     int mag = magnitude_cmp(*this,b);   // Determine relative sizes
00534     if (mag > 0) {                      // If abs(*this) > abs(b)
00535       subtract(*this,b,sum);            //   sum = *this - b
00536       sum.sign = this->sign;            // Sign of sum follows *this
00537     }
00538     else if (mag < 0) {                 // Else if abs(*this) < abs(b)
00539       subtract(b,*this,sum);            //   sum = b - *this
00540       sum.sign = b.sign;                // Sign of sum follows b
00541     }                                   // (Else abs(*this) == abs(b)
00542   }                                     //   so sum must be zero)
00543   return sum;                           // shallow swap on return
00544 }
00545 
00546 //: Multiplies this with a vnl_bignum
00547 
00548 vnl_bignum& vnl_bignum::operator*=(const vnl_bignum& b)
00549 {
00550   // Infinity arithmetic:
00551   assert (! b.is_infinity() || this->count != 0 ); // multiplication 0*Inf
00552   assert (! this->is_infinity() || b.count != 0 ); // multiplication Inf*0
00553   if (b.is_infinity()) return (*this) = (this->sign<0 ? -b : b);
00554   if (this->is_infinity()) return (*this) = (b.sign<0 ? -(*this) : *this);
00555 
00556   if (b.count == 0 || this->count == 0)
00557     return (*this)=0L;
00558   vnl_bignum prod;
00559   prod.resize(this->count + b.count);           //   allocate data for product
00560   for (Counter i = 0; i < b.count; i++)         //   multiply each b "digit"
00561     multiply_aux(*this, b.data[i], prod, i);    //   times b1 and add to total
00562   prod.sign = this->sign * b.sign;              //   determine correct sign
00563   prod.trim();                                  //   trim excess data and ret.
00564   return (*this)=prod;
00565 }
00566 
00567 //: Divides this by a vnl_bignum
00568 
00569 vnl_bignum& vnl_bignum::operator/=(const vnl_bignum& b)
00570 {
00571   // Infinity arithmetic:
00572   assert (! b.is_infinity() || ! this->is_infinity() ); // division Inf/Inf
00573   if (b.is_infinity()) return (*this)=0L;
00574   if (this->is_infinity()) return (*this) = (b.sign<0 ? -(*this) : *this);
00575   assert (b.count!=0 || this->count != 0); // division 0/0
00576   if (b.count == 0)
00577     return (*this) = (this->sign < 0 ? vnl_bignum("-Inf") : vnl_bignum("+Inf"));
00578 
00579   vnl_bignum quot, r;          // Quotient and remainder
00580   divide(*this,b,quot,r);      // Call divide fn
00581   return (*this) = quot;
00582 }
00583 
00584 //: Divides this by a vnl_bignum and replaces this by remainder.
00585 
00586 vnl_bignum& vnl_bignum::operator%=(const vnl_bignum& b)
00587 {
00588   // Infinity arithmetic:
00589   assert (! b.is_infinity() || ! this->is_infinity() ); // division Inf/Inf
00590   if (b.is_infinity()) return *this;                    // remainder of x/Inf is x.
00591   if (this->is_infinity()) return (*this) = 0L;         // convention: remainder is 0
00592   assert (b.count!=0 || this->count != 0);              // division 0/0
00593   if (b.count == 0) return (*this) = 0L;                // convention: remainder is 0
00594 
00595   vnl_bignum remain, q;        // Quotient and remainder
00596   divide(*this,b,q,remain);    // divide by b and save remainder
00597   return (*this) = remain;     // shallow swap on return
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 //: Shifts bignum to the right l digits.
00616 
00617 vnl_bignum vnl_bignum::operator>>(int l) const
00618 {
00619   // Infinity arithmetic:
00620   if (this->is_infinity()) return *this;
00621 
00622   if (l == 0 || *this == 0L)            // if either arg is zero
00623     return *this;
00624   if (l < 0)                            // if shift amt is negative
00625     return left_shift(*this,-l);        //   do an actual left shift
00626   else                                  // else
00627     return right_shift(*this,l);        //   do a right shift
00628 }
00629 
00630 //: Two vnl_bignums are equal if and only if they have the same integer representation.
00631 
00632 bool vnl_bignum::operator==(const vnl_bignum& rhs) const
00633 {
00634   if (this != &rhs) {                           // Check address
00635     if (this->sign != rhs.sign) return false;   // Different sign implies !=
00636     if (this->count != rhs.count) return false; // Different size implies !=
00637     for (Counter i = 0; i < this->count; i++)   // Each data element the same?
00638       if (this->data[i] != rhs.data[i]) return false; // No. Return !=
00639   }
00640   return true;                                    // Yes. Return ==
00641 }
00642 
00643 //: Compares two vnl_bignums.
00644 
00645 bool vnl_bignum::operator<(const vnl_bignum& rhs) const
00646 {
00647   if (this->sign < rhs.sign) return true;       // Different signs?
00648   if (this->sign > rhs.sign) return false;
00649   if (this->sign == 1)                          // Both signs == 1
00650     return magnitude_cmp(*this,rhs) < 0;        // this must be smaller
00651   else                                          // Both signs == -1
00652     return magnitude_cmp(*this,rhs) > 0;        // this must be larger
00653 }
00654 
00655 //: Formatted output for bignum.
00656 
00657 vcl_ostream& operator<<(vcl_ostream& os, const vnl_bignum& b)
00658 {
00659   vnl_bignum d = b;                     // Copy the input vnl_bignum
00660   if (d.sign == -1) {                   // If it's negative
00661     os << '-';                          //   Output leading minus sign
00662     d.sign = 1;                         //   Make d positive for divide
00663   }
00664   if (d.is_infinity()) return os<<"Inf";
00665   vnl_bignum q,r;                       // Temp quotient and remainder
00666   char *cbuf = new char[5 * (b.count+1)];   // Temp character buffer
00667   Counter i = 0;
00668   do {                                  // repeat:
00669     divide(d,10L,q,r);                  //   Divide vnl_bignum by ten
00670     cbuf[i++] = char(long(r) + '0');    //   Get one's digit
00671     d = q;                              //   Then discard one's digit
00672     q = r = 0L;                         //   Prep for next divide
00673   } while (d != 0L);                    // until no more one's digits
00674   do {                                  // repeat;
00675     os << cbuf[--i];                    //   output char buf in reverse
00676   } while (i);                          // until no more chars
00677   delete [] cbuf;                       // delete temp char buf
00678   return os;                            // return output stream
00679 }
00680 
00681 //: Convert the number to a decimal representation in a string.
00682 
00683 vcl_string& vnl_bignum_to_string(vcl_string& s, const vnl_bignum& b)
00684 {
00685   s.erase();
00686   vcl_string::size_type insert_point = 0; // keep record of location of first number.
00687 
00688   vnl_bignum d = b;                     // Copy the input vnl_bignum
00689   if (d.sign == -1) {                   // If it's negative
00690     s.insert(insert_point,"-");         //   Output leading minus sign
00691     d.sign = 1;                         //   Make d positive for divide
00692     ++insert_point;                     // keep record of location of first number.
00693   }
00694   if (d.is_infinity()) return s+="Inf";
00695   vnl_bignum q,r;                       // Temp quotient and remainder
00696   do {                                  // repeat:
00697     divide(d,10L,q,r);                  //   Divide vnl_bignum by ten
00698     s.insert(insert_point, 1, char('0'+long(r))); //   Get one's digit, and insert it at head.
00699     d = q;                              //   Then discard one's digit
00700     q = r = 0L;                         //   Prep for next divide
00701   } while (d != 0L);                    // until no more one's digits
00702   return s;
00703 }
00704 
00705 //: Convert the number from a decimal representation in a string.
00706 
00707 vnl_bignum& vnl_bignum_from_string(vnl_bignum& b, const vcl_string& s)
00708 {
00709   // decimal:     "^ *[-+]?[1-9][0-9]*$"
00710   // Infinity:    "^ *[-+]?Inf(inity)?$"
00711 
00712   if (is_plus_inf(s.c_str()))
00713     b=vnl_bignum("+Inf");
00714   else if (is_minus_inf(s.c_str()))
00715     b=vnl_bignum("-Inf");
00716   else
00717     b.dtoBigNum(s.c_str());             // convert decimal to vnl_bignum
00718   return b;
00719 }
00720 
00721 //: Implicit conversion from a vnl_bignum to a short.
00722 
00723 vnl_bignum::operator short() const
00724 {
00725   int j = this->operator int();
00726   return (short)j;
00727 }
00728 
00729 //: Implicit conversion from a vnl_bignum to an int.
00730 
00731 vnl_bignum::operator int() const
00732 {
00733   int j = 0;
00734   for (Counter i = this->count; i > 0; )
00735     j = int(j*0x10000 + this->data[--i]);
00736   return (this->sign < 0) ? -j : j;
00737 }
00738 
00739 //: Implicit conversion from a vnl_bignum to a long.
00740 
00741 vnl_bignum::operator long() const
00742 {
00743   long l = 0;
00744   for (Counter i = this->count; i > 0; )
00745     l = l*0x10000L + this->data[--i];
00746   return (this->sign < 0) ? -l : l;
00747 }
00748 
00749 //: Implicit conversion from a vnl_bignum to a float.
00750 
00751 vnl_bignum::operator float() const
00752 {
00753   float f = 0.0f;
00754   for (Counter i = this->count; i > 0; )
00755     f = f*0x10000 + this->data[--i];
00756   if (this->is_infinity()) f = vcl_numeric_limits<float>::infinity();
00757   return (this->sign < 0) ? -f : f;
00758 }
00759 
00760 //: Implicit conversion from a vnl_bignum to a double.
00761 
00762 vnl_bignum::operator double() const
00763 {
00764   double d = 0.0;
00765   for (Counter i = this->count; i > 0; )
00766     d = d*0x10000 + this->data[--i];
00767   if (this->is_infinity()) d = vcl_numeric_limits<double>::infinity();
00768   return (this->sign < 0) ? -d : d;
00769 }
00770 
00771 //: Implicit conversion from a vnl_bignum to a long double.
00772 
00773 vnl_bignum::operator long double() const
00774 {
00775   long double d = 0.0;
00776   for (Counter i = this->count; i > 0; )
00777     d = d*0x10000 + this->data[--i];
00778   if (this->is_infinity()) d = vcl_numeric_limits<long double>::infinity();
00779   return (this->sign < 0) ? -d : d;
00780 }
00781 
00782 //: dump the contents of a vnl_bignum to a stream, default cout.
00783 
00784 void vnl_bignum::dump(vcl_ostream& os) const
00785 {
00786   os << "{count=" << this->count       // output count field
00787      << ", sign=" << this->sign        // output sign field
00788      << ", data=" << this->data        // output data pointer
00789      << ", value=" << *this
00790      << ", {";
00791   // format string == "%04X%s" or "%02X%s", etc.
00792   //  static char format_str[10] =
00793   //    {'%','0',char(2*2 + '0'),'X','%','s'};
00794   //  format_str[2] = char(2*2 + '0');
00795   if (this->count > 0) { // output data array
00796     os << vcl_hex << (this->data[this->count-1]);
00797     for (Counter i = this->count - 1; i > 0; --i) {
00798       os  << ',';
00799       if (this->data[i-1] < 0x10) os << '0';
00800       if (this->data[i-1] < 0x100) os << '0';
00801       if (this->data[i-1] < 0x1000) os << '0';
00802       os  << this->data[i-1];
00803     }
00804     os << vcl_dec;
00805   }
00806   os << "}}\n";                         // close brackets
00807 }
00808 
00809 //: Converts decimal string to a vnl_bignum.
00810 
00811 int vnl_bignum::dtoBigNum(const char *s)
00812 {
00813   this->resize(0); this->sign = 1;      // Reset number to 0.
00814   Counter len = 0;                      // No chars converted yet
00815   vnl_bignum sum;
00816   while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s; // skip whitespace
00817   if (*s == '-' || *s == '+') ++len;    // Skip over leading +,-
00818   while (s[len]>='0' && s[len]<='9') {  // While current char is digit
00819     *this *= vnl_bignum(10L),           // Shift vnl_bignum left a decimal
00820     add(*this,vnl_bignum(long(s[len++]-'0')),sum), // digit and add new digit
00821     *this = sum;
00822   }
00823   if (*s == '-') this->sign = -1;       // If s had leading -, note it
00824   return len;                           // Return # of chars processed
00825 }
00826 
00827 //: convert exponential string to a vnl_bignum
00828 
00829 void vnl_bignum::exptoBigNum(const char *s)
00830 {
00831   while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s; // skip whitespace
00832   Counter pos = Counter(this->dtoBigNum(s) + 1); // Convert the base, skip [eE]
00833   long pow = vcl_atol(s + pos);         // Convert the exponent to long
00834   while (pow-- > 0)                     // Raise vnl_bignum to the given
00835     *this = (*this) * 10L;              // power
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 = Counter(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 //: convert octal string to vnl_bignum
00866 
00867 void vnl_bignum::otoBigNum(const char *s)
00868 {
00869   this->resize(0); sign = 1;           // Reset number to 0.
00870   while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s; // skip whitespace
00871   Counter size = Counter(vcl_strlen(s));
00872   Counter len = 0;                      // No chars converted yet
00873   while (len < size) {                  // While there are more chars
00874     (*this) = ((*this) * 8L) +          // Shift vnl_bignum left 1 oct dig.
00875       vnl_bignum(long(s[len++] - '0')); // Add next character value
00876   }
00877 }
00878 
00879 //: change the data allotment for a vnl_bignum
00880 
00881 void vnl_bignum::resize(short new_count)
00882 {
00883   assert(new_count >= 0);
00884   if (new_count == this->count) return;
00885   Data *new_data = (new_count > 0 ? new Data[new_count] : 0); // Allocate data if necessary
00886 
00887   if (this->count <= new_count) {       // Copy old data into new
00888     short i = 0;
00889     for (; i < this->count; i++)
00890       new_data[i] = this->data[i];
00891     for (; i < new_count; i++)
00892       new_data[i] = 0;
00893   }
00894   else {
00895     for (short i = 0; i < new_count; i++)
00896       new_data[i] = this->data[i];
00897   }
00898 
00899   delete [] this->data;                 // Get rid of old data
00900   this->data = new_data;                // Point to new data
00901   this->count = new_count;              // Save new count
00902 }
00903 
00904 //: trim non-infinite vnl_bignum of excess data allotment
00905 
00906 vnl_bignum& vnl_bignum::trim()
00907 {
00908   Counter i = this->count;
00909   for (; i > 0; i--)                    // Skip over high-order words
00910     if (this->data[i - 1] != 0) break;  //   that are zero
00911   if (i < this->count) {                // If there are some such words
00912     this->count = i;                    // Update the count
00913     Data *new_data = (i > 0 ? new Data[i] : 0); // Allocate data if necessary
00914     for (; i > 0; i--)                  // Copy old data into new
00915       new_data[i - 1] = this->data[i - 1];
00916     delete [] this->data;               // Delete old data
00917     this->data = new_data;              // Point to new data
00918   }
00919   return *this;                         // return reference to vnl_bignum
00920 }
00921 
00922 //: add two non-infinite vnl_bignum values and save their sum
00923 
00924 void add(const vnl_bignum& b1, const vnl_bignum& b2, vnl_bignum& sum)
00925 {
00926   const vnl_bignum *bmax, *bmin;        // Determine which of the two
00927   if (b1.count >= b2.count) {           //   addends has the most
00928     bmax = &b1;                         //   data.
00929     bmin = &b2;
00930   }
00931   else {
00932     bmax = &b2;
00933     bmin = &b1;
00934   }
00935   sum.resize(bmax->count);              // Allocate data for their sum
00936   unsigned long temp, carry = 0;
00937   Counter i = 0;
00938   while (i < bmin->count) {             // Add, element by element.
00939     // Add both elements and carry
00940     temp = (unsigned long)b1.data[i] + (unsigned long)b2.data[i] + carry;
00941     carry = temp/0x10000L;              // keep track of the carry
00942     sum.data[i] = Data(temp);           // store sum
00943     i++;                                // go to next element
00944   }
00945   while (i < bmax->count) {             // bmin has no more elements
00946     temp = bmax->data[i] + carry;       // propagate the carry through
00947     carry = temp/0x10000L;              // the rest of bmax's elements
00948     sum.data[i] = Data(temp);           // store sum
00949     i++;
00950   }
00951   if (carry) {                          // if carry left over
00952     sum.resize(bmax->count + 1);        //   allocate another word
00953     sum.data[bmax->count] = 1;          //   save the carry in it
00954   }
00955 }
00956 
00957 //: Add 1 to bnum (unsigned, non-infinite)
00958 
00959 void increment(vnl_bignum& bnum)
00960 {
00961   Counter i = 0;
00962   unsigned long carry = 1;
00963   while (i < bnum.count && carry) {             // increment, element by element.
00964     unsigned long temp = (unsigned long)bnum.data[i] + carry;
00965     carry = temp/0x10000L;
00966     bnum.data[i] = (Data)temp;
00967     ++i;
00968   }
00969   if (carry)
00970   {
00971     bnum.resize(bnum.count+1);
00972     bnum.data[bnum.count-1] = 1;
00973   }
00974 }
00975 
00976 //: subtract bmin from bmax (unsigned, non-infinite), result in diff
00977 
00978 void subtract(const vnl_bignum& bmax, const vnl_bignum& bmin, vnl_bignum& diff)
00979 {
00980   diff.resize(bmax.count);                      // Allocate data for difference
00981   unsigned long temp;
00982   int borrow = 0;
00983   Counter i = 0;
00984   for (; i < bmin.count; i++) {                 // Subtract word by word.
00985     temp = (unsigned long)bmax.data[i] + 0x10000L - borrow; // Add radix to bmax's data
00986     temp -= (unsigned long)bmin.data[i];        // Subtract off bmin's data
00987     borrow = (temp/0x10000L == 0);              // Did we have to borrow?
00988     diff.data[i] = (Data) temp;                 // Reduce modulo radix and save
00989   }
00990   for (; i < bmax.count; i++) {                 // No more data for bmin
00991     temp = (unsigned long)bmax.data[i] + 0x10000L - borrow; // Propagate the borrow through
00992     borrow = (temp/0x10000L == 0);              //   rest of bmax's data
00993     diff.data[i] = (Data) temp;
00994   }
00995   diff.trim();                                  // Done. Now trim excess data
00996 }
00997 
00998 //: Subtract 1 from bnum (unsigned, non-infinite, non-zero)
00999 
01000 void decrement(vnl_bignum& bnum)
01001 {
01002   Counter i = 0;
01003   unsigned long borrow = 1;
01004   while (i < bnum.count && borrow) {            // decrement, element by element.
01005     unsigned long temp = (unsigned long)bnum.data[i] + 0x10000L - borrow;
01006     borrow = (temp/0x10000L == 0);              // Did we have to borrow?
01007     bnum.data[i] = (Data)temp;                  // Reduce modulo radix and save
01008     ++i;
01009   }
01010   bnum.trim();                                  // Done. Now trim excess data
01011   if (bnum.count==0) bnum.sign=+1;              // Re-establish sign invariant
01012 }
01013 
01014 //: compare absolute values of two vnl_bignums
01015 // Outputs:  result of comparison:  -1 if abs(b1) < abs(b2)
01016 //                                   0 if abs(b1) == abs(b2)
01017 //                                  +1 if abs(b1) > abs(b2)
01018 
01019 int magnitude_cmp(const vnl_bignum& b1, const vnl_bignum& b2)
01020 {
01021   if (b1.is_infinity()) return b2.is_infinity() ? 0 : 1;
01022   if (b2.is_infinity()) return -1;
01023   if (b1.count > b2.count) return 1;    // If one has more data than
01024   if (b2.count > b1.count) return -1;   //   the other, it wins
01025   Counter i = b1.count;                 // Else same number of elements
01026   while (i > 0) {                       // Do lexicographic comparison
01027     if (b1.data[i - 1] > b2.data[i - 1])
01028       return 1;
01029     else if (b1.data[i - 1] < b2.data[i - 1])
01030       return -1;
01031     i--;
01032   }                                     // No data, or all elements same
01033   return 0;                             //  so must be equal
01034 }
01035 
01036 //: multiply a non-infinite vnl_bignum by a "single digit"
01037 // - Inputs:  vnl_bignum reference, single word multiplier, reference to the product,
01038 //            and index of starting storage location to use in product
01039 
01040 void multiply_aux(const vnl_bignum& b, Data d, vnl_bignum& prod, Counter i)
01041 {
01042   // this function works just like normal multiplication by hand, in that the
01043   // top number is multiplied by the first digit of the bottom number, then the
01044   // second digit, and so on.  The only difference is that instead of doing all
01045   // of the multiplication before adding the rows, addition is done
01046   // concurrently.
01047   if (i == 0) {                         // if index is zero
01048     Counter j = 0;                      //   then zero out all of
01049     while (j < prod.count)              //   prod's data elements
01050       prod.data[j++] = 0;
01051   }
01052   if (d != 0) {                         // if d == 0, nothing to do
01053     unsigned long temp;
01054     Data carry = 0;
01055 
01056     Counter j = 0;
01057     for (; j < b.count; j++) {
01058       // for each of b's data elements, multiply times d and add running product
01059       temp = (unsigned long)b.data[j] * (unsigned long)d
01060            + (unsigned long)prod.data[i + j] + carry;
01061       prod.data[i + j] = Data(temp % 0x10000L); //   store result in product
01062       carry = Data(temp/0x10000L);              //   keep track of carry
01063     }
01064     if (i+j < prod.count)
01065       prod.data[i + j] = carry;                 // Done.  Store the final carry
01066   }
01067 }
01068 
01069 //: normalize two vnl_bignums
01070 // (Refer to Knuth, V.2, Section 4.3.1, Algorithm D for details.
01071 //  A digit here is one data element in the radix 2**2.)
01072 // - Inputs:  references to two vnl_bignums b1, b2
01073 // - Outputs: their normalized counterparts u and v,
01074 //            and the integral normalization factor used
01075 
01076 Data normalize(const vnl_bignum& b1, const vnl_bignum& b2, vnl_bignum& u, vnl_bignum& v)
01077 {
01078   Data d = Data(0x10000L/((unsigned long)(b2.data[b2.count - 1]) + 1L)); // Calculate normalization factor.
01079   u.resize(b1.count + 1);                       // Get data for u (plus extra)
01080   v.resize(b2.count);                           // Get data for v
01081   u.data[b1.count] = 0;                         // Set u's leading digit to 0
01082   multiply_aux(b1,d,u,0);                       // u = b1 * d
01083   multiply_aux(b2,d,v,0);                       // v = b2 * d
01084   return d;                                     // return normalization factor
01085 }
01086 
01087 //: divide a vnl_bignum by a "single digit"
01088 // (Refer to Knuth, V.2, Section 4.3.2, exercise 16 for details.
01089 //  A digit here is one data element in the radix 2**2.)
01090 // - Inputs:  reference to vnl_bignum dividend, single digit divisor d, vnl_bignum
01091 //            quotient, and single digit remainder r
01092 
01093 void divide_aux(const vnl_bignum& b1, Data d, vnl_bignum& q, Data& r)
01094 {
01095   r = 0;                                        // init remainder to zero
01096   unsigned long temp;
01097   for (Counter j = b1.count; j > 0; j--) {
01098     temp = (unsigned long)r*0x10000L + (unsigned long)b1.data[j - 1]; // get remainder, append next
01099     if (j < 1 + q.count)
01100       q.data[j - 1] = Data(temp/d);             //   digit, then divide
01101     r = Data(temp % d);                         // calculate new remainder
01102   }
01103 }
01104 
01105 //: estimate next dividend
01106 // (Refer to Knuth, V.2, Section 4.3.1, Algorithm D for details.
01107 //  This function estimates how many times v goes into u, starting at u's
01108 //  jth digit.  A digit here is actually a data element, thought of as
01109 //  being in the radix 2**2.)
01110 // - Inputs:  reference to vnl_bignum dividend and divisor and starting digit j
01111 // - Outputs: estimated number of times v goes into u
01112 
01113 Data estimate_q_hat(const vnl_bignum& u, const vnl_bignum& v, Counter j)
01114 {
01115   Data q_hat,
01116        v1 = v.data[v.count - 1],                // localize frequent data
01117        v2 = v.data[v.count - 2],
01118        u0 = u.data[u.count - 1 - j],
01119        u1 = u.data[u.count - 2 - j],
01120        u2 = u.data[u.count - 3 - j];
01121 
01122   // Initial Knuth estimate, usually correct
01123   q_hat = (u0 == v1 ? Data(0xffff) : Data(((unsigned long)u0 * 0x10000L + (unsigned long)u1) / (unsigned long)v1));
01124 
01125   // high speed test to determine most of the cases in which initial
01126   // estimate is too large.  Eliminates most cases in which q_hat is one too
01127   // large.  Eliminates all cases in which q_hat is two too large.  The test
01128   // looks hairy because we have to watch out for overflow.  In the book, this
01129   // test is the simple inequality:
01130   //     v2*q_hat > (u0*0x10000L + u1 - q_hat*v1)*0x10000L + u2.
01131   // If the inequality is true, decrease q_hat by 1.  If inequality is still
01132   // true, decrease q_hat again.
01133   unsigned long lhs, rhs;               // lhs, rhs of Knuth inequality
01134   for (Counter i = 0; i < 2; i++) {     // loop at most twice
01135     lhs = (unsigned long)v2 * (unsigned long)q_hat;       // Calculate left-hand side of ineq.
01136     rhs = (unsigned long)u0 * 0x10000L +(unsigned long)u1;// Calculate part of right-hand side
01137     rhs -= ((unsigned long)q_hat * (unsigned long)v1);    // Now subtract off part
01138 
01139     if ( rhs >= 0x10000L )              // if multiplication with 0x10000L causes overflow
01140       break;                            //   then rhs > lhs, so test fails
01141     rhs *= 0x10000L;                    // No overflow:  ok to multiply
01142 
01143     if (rhs > rhs + (unsigned long)u2)  // if addition yields overflow
01144       break;                            //   then rhs > lhs, so test fails
01145     rhs += u2;                          // No overflow: ok to add.
01146     if (lhs <= rhs)                     // if lhs <= rhs
01147       break;                            //   test fails
01148     q_hat--;                            // Test passes:  decrement q_hat
01149   }                                     // Loop again
01150   return q_hat;                         // Return estimate
01151 }
01152 
01153 //: calculate u - v*q_hat
01154 // (Refer to Knuth, V. 2, Section 4.3.1, Algorithm D for details.
01155 //  A digit here is a data element, thought of as being in the radix 2**2.)
01156 // - Inputs:  reference to vnl_bignum dividend, divisor, estimated result, and index
01157 //            into jth digit of dividend
01158 // - Outputs: true number of times v goes into u
01159 
01160 Data multiply_subtract(vnl_bignum& u, const vnl_bignum& v, Data q_hat, Counter j)
01161 {
01162   // At this point it has been estimated that v goes into the jth and higher
01163   // digits of u about q_hat times, and in fact that q_hat is either the
01164   // correct number of times or one too large.
01165 
01166   if (q_hat == 0) return q_hat;         // if q_hat 0, nothing to do
01167   vnl_bignum rslt;                      // create a temporary vnl_bignum
01168   Counter tmpcnt;
01169   rslt.resize(v.count + 1u);            // allocate data for it
01170 
01171   // simultaneous computation of u - v*q_hat
01172   unsigned long prod, diff;
01173   Data carry = 0, borrow = 0;
01174   Counter i = 0;
01175   for (; i < v.count; ++i) {
01176     // for each digit of v, multiply it by q_hat and subtract the result
01177     prod = (unsigned long)v.data[i] * (unsigned long)q_hat + carry;
01178     diff = (unsigned long)u.data[u.count - v.count - 1 - j + i] + (0x10000L - (unsigned long)borrow);
01179     diff -= (unsigned long)Data(prod);  //   form proper digit of u
01180     rslt.data[i] = Data(diff);          //   save the result
01181     borrow = (diff/0x10000L == 0) ? 1 : 0; //   keep track of any borrows
01182     carry = Data(prod/0x10000L);        //   keep track of carries
01183   }
01184   tmpcnt = Counter(u.count - v.count + i - j - 1);
01185   diff = (unsigned long)u.data[tmpcnt]  //  special case for the last digit
01186          + (0x10000L - (unsigned long)borrow);
01187   diff -= (unsigned long)carry;
01188   rslt.data[i] = Data(diff);
01189   borrow = (diff/0x10000L == 0) ? 1 : 0;
01190 
01191   // A leftover borrow indicates that u - v*q_hat is negative, i.e., that
01192   // q_hat was one too large.  So to get correct result, decrement q_hat and
01193   // add back one multiple of v
01194   if (borrow) {
01195     q_hat--;
01196     carry = 0;
01197     unsigned long sum;
01198     for (i = 0; i < v.count; ++i) {
01199       sum = (unsigned long)rslt.data[i] + (unsigned long)v.data[i] + carry;
01200       carry = Data(sum/0x10000L);
01201       u.data[u.count - v.count + i - 1 - j] = Data(sum);
01202     }
01203     u.data[u.count - v.count + i - 1 - j] = rslt.data[i] + carry;
01204   }
01205   else {                                // otherwise, result is ok
01206     for (i = 0; i < rslt.count; ++i)    // store result back into u
01207       u.data[u.count - v.count + i - 1 - j] = rslt.data[i];
01208   }
01209   return q_hat;                         // return corrected q_hat
01210 }
01211 
01212 //: divide b2 into b1, getting quotient q and remainder r.
01213 // (Refer to Knuth, V.2, Section 4.3.1, Algorithm D for details.
01214 //  This function implements Algorithm D.)
01215 // - Inputs: references to a vnl_bignum dividend b1, divisor b2, quotient q, and
01216 //           remainder r.
01217 
01218 void divide(const vnl_bignum& b1, const vnl_bignum& b2, vnl_bignum& q, vnl_bignum& r)
01219 {
01220   // Note that q or r may *not* be identical to either b1 or b2 !
01221   assert(&b1 != &q && &b2 != &q && &b1 != &r && &b2 != &r);
01222   q = r = 0L;
01223   if (b1 == 0L)                      // If divisor is zero
01224     return;                          //   return zero quotient and remainder
01225   int mag = magnitude_cmp(b1,b2);    // Compare magnitudes
01226   if (mag < 0)                       // if abs(b1) < abs(b2)
01227     r = b1;                          //   return zero quotient, b1 remainder
01228   else if (mag == 0)                 // if abs(b1) == abs(b2)
01229     q = 1L;                          //   quotient is 1, remainder is 0
01230   else {                             // otherwise abs(b1) > abs(b2), so divide
01231     q.resize(b1.count + 1u - b2.count); // Allocate quotient
01232     r.resize(b2.count);              // Allocate remainder
01233     if (b2.count == 1) {                        // Single digit divisor?
01234       divide_aux(b1,b2.data[0],q,r.data[0]);    // Do single digit divide
01235     }
01236     else {                                      // Else full-blown divide
01237       vnl_bignum u,v;
01238 #ifdef DEBUG
01239       vcl_cerr << "\nvnl_bignum::divide: b1 ="; b1.dump(vcl_cerr);
01240       vcl_cerr << "vnl_bignum::divide: b2 ="; b2.dump(vcl_cerr);
01241 #endif
01242       Data d = normalize(b1,b2,u,v);            // Set u = b1*d, v = b2*d
01243 #ifdef DEBUG
01244       vcl_cerr << "vnl_bignum::divide: d = " << d << vcl_hex << " (0x" << d << ")\n" << vcl_dec;
01245       vcl_cerr << "vnl_bignum::divide: u = "; u.dump(vcl_cerr);
01246       vcl_cerr << "vnl_bignum::divide: v = "; v.dump(vcl_cerr);
01247 #endif
01248       Counter j = 0;
01249       while (j <= b1.count - b2.count) {        // Main division loop
01250         Data q_hat = estimate_q_hat(u,v,j);     // Estimate # times v divides u
01251         q.data[q.count - 1 - j] =               // Do division, get true answ.
01252           multiply_subtract(u,v,q_hat,j);
01253         j++;
01254 #ifdef DEBUG
01255         vcl_cerr << "vnl_bignum::divide: q_hat = " << q_hat << vcl_hex << " (0x" << q_hat << ")\n" << vcl_dec;
01256         vcl_cerr << "vnl_bignum::divide: u = "; u.dump(vcl_cerr);
01257 #endif
01258       }
01259       static Data dufus;                // dummy variable
01260       divide_aux(u,d,r,dufus);          // Unnormalize u for remainder
01261 
01262 #ifdef DEBUG
01263       vcl_cerr << "vnl_bignum::divide: q = "; q.dump(vcl_cerr);
01264       vcl_cerr << "vnl_bignum::divide: r = "; r.dump(vcl_cerr);
01265 #endif
01266     }
01267     q.trim();                           // Trim leading zeros of quot.
01268     r.trim();                           // Trim leading zeros of rem.
01269   }
01270   q.sign = r.sign = b1.sign * b2.sign;  // Calculate signs
01271 }
01272 
01273 //: left shift (arithmetic) non-infinite vnl_bignum by positive number.
01274 // - Inputs:  reference to vnl_bignum, positive shift value
01275 // - Outputs: vnl_bignum, multiplied by the corresponding power of two
01276 
01277 vnl_bignum left_shift(const vnl_bignum& b1, int l)
01278 {
01279   // to carry out this arithmetic left shift, we cheat.  Instead of physically
01280   // shifting the data array l bits to the left, we shift just enough to get
01281   // the correct word alignment, and then pad the array on the right with as
01282   // many zeros as we need.
01283   vnl_bignum rslt;                      // result of shift
01284   rslt.sign = b1.sign;                  // result follows sign of input
01285   Counter growth = Counter(l / 16);     // # of words rslt will grow by
01286   Data shift = Data(l % 16);            // amount to actually shift
01287   Data rshift = Data(16 - shift);       // amount to shift next word by
01288   Data carry = Data(                    // value that will be shifted
01289     b1.data[b1.count - 1] >> (16 - shift));// out end of current array
01290   rslt.resize(b1.count + growth + (carry ? 1u : 0u)); // allocate new data array
01291   Counter i = 0;
01292   while (i < growth)                            // zero out padded elements
01293     rslt.data[i++] = 0;
01294   rslt.data[i++] = b1.data[0] << shift;         // shift first non-zero element
01295   while (i < rslt.count - 1) {                  // for remaining data words
01296     rslt.data[i] = (b1.data[i - growth] << shift) + // shift current data word
01297       (b1.data[i - 1 - growth] >> rshift);          // propagate adjacent
01298     i++;                                        // carry into current word
01299   }
01300   if (i < rslt.count) {
01301     if (carry)                                  // if last word had overflow
01302       rslt.data[i] = carry;                     //   store it new data
01303     else                                        // otherwise,
01304       rslt.data[i] = (b1.data[i - growth] << shift) // do like the rest
01305                    + (b1.data[i - 1 - growth] >> rshift);
01306   }
01307   vnl_bignum& result = *((vnl_bignum*) &rslt);// same physical object
01308   return result;                              // shallow swap on return
01309 }
01310 
01311 //: right shift (arithmetic) non-infinite vnl_bignum by positive number.
01312 // - Inputs:  reference to vnl_bignum, positive shift value
01313 // - Outputs: vnl_bignum, divided by the corresponding power of two
01314 
01315 vnl_bignum right_shift(const vnl_bignum& b1, int l)
01316 {
01317   vnl_bignum rslt;                              // result of shift
01318   Counter shrinkage = Counter(l / 16);          // # of words rslt will shrink
01319   Data shift = Data(l % 16);                    // amount to actually shift
01320   Data dregs = Data(b1.data[b1.count-1] >> shift);// high end data to save
01321   if (shrinkage + (dregs == 0) < b1.count) {    // if not all data shifted out
01322     rslt.sign = b1.sign;                        // rslt follows sign of input
01323     rslt.resize(b1.count - shrinkage - (dregs == 0 ? 1 : 0)); // allocate new data
01324     Data lshift = Data(16 - shift);             // amount to shift high word
01325     Counter i = 0;
01326     while (i < rslt.count - 1) {                // shift current word
01327       rslt.data[i] = (b1.data[i + shrinkage] >> shift) + // propagate adjacent
01328         (b1.data[i + shrinkage + 1u] << lshift); // word into current word
01329       i++;
01330     }
01331     if (dregs)                                  // don't lose dregs
01332       rslt.data[i] = dregs;
01333     else {
01334       rslt.data[i] = (b1.data[i + shrinkage] >> shift) +
01335         (b1.data[i + shrinkage + 1u] << lshift);
01336     }
01337   }
01338   vnl_bignum& result = *((vnl_bignum*) &rslt);  // same physical object
01339   return result;                                // shallow swap on return
01340 }