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 
00209 #if 0 // old, original Texas Instruments implementation - PVr
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   this->data = b.data ? new Data[b.count] : 0;  // Allocate data if necessary
00438   for (Counter i = 0; i < this->count; ++i)     // Copy b data
00439     this->data[i] = b.data[i];
00440 }
00441 
00442 
00443 //: Frees space for vnl_bignum.
00444 
00445 vnl_bignum::~vnl_bignum()
00446 {
00447   delete [] this->data; this->count = 0;        // Delete any allocated data
00448 }
00449 
00450 //: Copies rhs vnl_bignum to lhs vnl_bignum.
00451 
00452 vnl_bignum& vnl_bignum::operator=(const vnl_bignum& rhs)
00453 {
00454   if (this != &rhs) {                           // Avoid self-assignment
00455     delete[] this->data;                        // Delete existing data
00456     this->count = rhs.count;                    // Copy rhs's count
00457     this->data = rhs.data ? new Data[rhs.count] : 0; // Allocate data if necessary
00458     for (Counter i = 0; i < rhs.count; ++i)     // Copy rhs's data
00459       this->data[i] = rhs.data[i];
00460     this->sign = rhs.sign;                      // Copy rhs's sign
00461   }
00462   return *this;                                 // Return reference
00463 }
00464 
00465 //: Returns the negation of a vnl_bignum.
00466 
00467 vnl_bignum vnl_bignum::operator-() const
00468 {
00469   vnl_bignum neg(*this);
00470   if (neg.count)                // So long as this is non-zero
00471     neg.sign = -neg.sign;       // Flip its sign
00472   return neg;
00473 }
00474 
00475 
00476 //: Prefix increment. Increments a vnl_bignum by 1, and returns it.
00477 
00478 vnl_bignum& vnl_bignum::operator++()
00479 {
00480   if (this->is_infinity()) return *this;
00481   if (this->count==0)
00482   {
00483     this->resize(1);
00484     this->data[0] = 1;
00485     this->sign = +1;
00486     return *this;
00487   }
00488 
00489   if (this->sign > 0) increment(*this);
00490   else decrement(*this);
00491 
00492   return *this;
00493 }
00494 
00495 
00496 //: Prefix decrement. Decrements a vnl_bignum by 1, and returns it.
00497 
00498 vnl_bignum& vnl_bignum::operator--()
00499 {
00500   if (this->is_infinity()) return *this;
00501   if (this->count==0)
00502   {
00503     this->resize(1);
00504     this->data[0] = 1;
00505     this->sign = -1;
00506     return *this;
00507   }
00508 
00509   if (this->sign < 0) increment(*this);
00510   else decrement(*this);
00511 
00512   return *this;
00513 }
00514 
00515 //: Adds two vnl_bignums, and returns new sum.
00516 
00517 vnl_bignum vnl_bignum::operator+(const vnl_bignum& b) const
00518 {
00519   // Infinity arithmetic:
00520   assert (! b.is_minus_infinity() || ! this->is_plus_infinity() ); // +Inf-Inf
00521   assert (! b.is_plus_infinity() || ! this->is_minus_infinity() ); // -Inf+Inf
00522   if (b.is_infinity()) { return b; }
00523   if (this->is_infinity()) { return *this; }
00524 
00525   vnl_bignum sum;                       // Init sum to zero
00526   if (this->sign == b.sign) {           // If both have same sign
00527     add(*this,b,sum);                   //   Do simple addition
00528     sum.sign = this->sign;              // Attach proper sign
00529   }
00530   else {                                // Else different signs
00531     int mag = magnitude_cmp(*this,b);   // Determine relative sizes
00532     if (mag > 0) {                      // If abs(*this) > abs(b)
00533       subtract(*this,b,sum);            //   sum = *this - b
00534       sum.sign = this->sign;            // Sign of sum follows *this
00535     }
00536     else if (mag < 0) {                 // Else if abs(*this) < abs(b)
00537       subtract(b,*this,sum);            //   sum = b - *this
00538       sum.sign = b.sign;                // Sign of sum follows b
00539     }                                   // (Else abs(*this) == abs(b)
00540   }                                     //   so sum must be zero)
00541   return sum;                           // shallow swap on return
00542 }
00543 
00544 
00545 //: Multiplies this with a vnl_bignum
00546 
00547 vnl_bignum& vnl_bignum::operator*=(const vnl_bignum& b)
00548 {
00549   // Infinity arithmetic:
00550   assert (! b.is_infinity() || this->count != 0 ); // multiplication 0*Inf
00551   assert (! this->is_infinity() || b.count != 0 ); // multiplication Inf*0
00552   if (b.is_infinity()) return (*this) = (this->sign<0 ? -b : b);
00553   if (this->is_infinity()) return (*this) = (b.sign<0 ? -(*this) : *this);
00554 
00555   if (b.count == 0 || this->count == 0)
00556     return (*this)=0L;
00557   vnl_bignum prod;
00558   prod.resize(this->count + b.count);           //   allocate data for product
00559   for (Counter i = 0; i < b.count; i++)         //   multiply each b "digit"
00560     multiply_aux(*this, b.data[i], prod, i);    //   times b1 and add to total
00561   prod.sign = this->sign * b.sign;              //   determine correct sign
00562   prod.trim();                                  //   trim excess data and ret.
00563   return (*this)=prod;
00564 }
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 
00601 //: Shifts bignum to the left l digits.
00602 
00603 vnl_bignum vnl_bignum::operator<<(int l) const
00604 {
00605   // Infinity arithmetic:
00606   if (this->is_infinity()) return *this;
00607 
00608   if (l == 0 || *this == 0L)            // if either arg is zero
00609     return *this;
00610   if (l < 0)                            // if shift amt is negative
00611     return right_shift(*this,-l);       //   do an actual right shift
00612   else                                  // otherwise
00613     return left_shift(*this,l);         //   do a left shift
00614 }
00615 
00616 
00617 //: Shifts bignum to the right l digits.
00618 
00619 vnl_bignum vnl_bignum::operator>>(int l) const
00620 {
00621   // Infinity arithmetic:
00622   if (this->is_infinity()) return *this;
00623 
00624   if (l == 0 || *this == 0L)            // if either arg is zero
00625     return *this;
00626   if (l < 0)                            // if shift amt is negative
00627     return left_shift(*this,-l);        //   do an actual left shift
00628   else                                  // else
00629     return right_shift(*this,l);        //   do a right shift
00630 }
00631 
00632 
00633 //: Two vnl_bignums are equal if and only if they have the same integer representation.
00634 
00635 bool vnl_bignum::operator==(const vnl_bignum& rhs) const
00636 {
00637   if (this != &rhs) {                           // Check address
00638     if (this->sign != rhs.sign) return false;   // Different sign implies !=
00639     if (this->count != rhs.count) return false; // Different size implies !=
00640     for (Counter i = 0; i < this->count; i++)   // Each data element the same?
00641       if (this->data[i] != rhs.data[i]) return false; // No. Return !=
00642   }
00643   return true;                                    // Yes. Return ==
00644 }
00645 
00646 
00647 //: Compares two vnl_bignums.
00648 
00649 bool vnl_bignum::operator<(const vnl_bignum& rhs) const
00650 {
00651   if (this->sign < rhs.sign) return true;       // Different signs?
00652   if (this->sign > rhs.sign) return false;
00653   if (this->sign == 1)                          // Both signs == 1
00654     return magnitude_cmp(*this,rhs) < 0;        // this must be smaller
00655   else                                          // Both signs == -1
00656     return magnitude_cmp(*this,rhs) > 0;        // this must be larger
00657 }
00658 
00659 
00660 //: Formatted output for bignum.
00661 
00662 vcl_ostream& operator<<(vcl_ostream& os, const vnl_bignum& b)
00663 {
00664   vnl_bignum d = b;                     // Copy the input vnl_bignum
00665   if (d.sign == -1) {                   // If it's negative
00666     os << '-';                          //   Output leading minus sign
00667     d.sign = 1;                         //   Make d positive for divide
00668   }
00669   if (d.is_infinity()) return os<<"Inf";
00670   vnl_bignum q,r;                       // Temp quotient and remainder
00671   char *cbuf = new char[5 * (b.count+1)];   // Temp character buffer
00672   Counter i = 0;
00673   do {                                  // repeat:
00674     divide(d,10L,q,r);                  //   Divide vnl_bignum by ten
00675     cbuf[i++] = char(long(r) + '0');    //   Get one's digit
00676     d = q;                              //   Then discard one's digit
00677     q = r = 0L;                         //   Prep for next divide
00678   } while (d != 0L);                    // until no more one's digits
00679   do {                                  // repeat;
00680     os << cbuf[--i];                    //   output char buf in reverse
00681   } while (i);                          // until no more chars
00682   delete [] cbuf;                       // delete temp char buf
00683   return os;                            // return output stream
00684 }
00685 
00686 //: Convert the number to a decimal representation in a string.
00687 vcl_string& vnl_bignum_to_string(vcl_string& s, const vnl_bignum& b)
00688 {
00689   s.erase();
00690   vcl_string::size_type insert_point = 0; // keep record of location of first number.
00691 
00692   vnl_bignum d = b;                     // Copy the input vnl_bignum
00693   if (d.sign == -1) {                   // If it's negative
00694     s.insert(insert_point,"-");         //   Output leading minus sign
00695     d.sign = 1;                         //   Make d positive for divide
00696     ++insert_point;                     // keep record of location of first number.
00697   }
00698   if (d.is_infinity()) return s+="Inf";
00699   vnl_bignum q,r;                       // Temp quotient and remainder
00700   do {                                  // repeat:
00701     divide(d,10L,q,r);                  //   Divide vnl_bignum by ten
00702     s.insert(insert_point, 1, char('0'+long(r))); //   Get one's digit, and insert it at head.
00703     d = q;                              //   Then discard one's digit
00704     q = r = 0L;                         //   Prep for next divide
00705   } while (d != 0L);                    // until no more one's digits
00706   return s;
00707 }
00708 
00709 //: Convert the number from a decimal representation in a string.
00710 vnl_bignum& vnl_bignum_from_string(vnl_bignum& b, const vcl_string& s)
00711 {
00712   // decimal:     "^ *[-+]?[1-9][0-9]*$"
00713   // Infinity:    "^ *[-+]?Inf(inity)?$"
00714 
00715   if (is_plus_inf(s.c_str()))
00716     b=vnl_bignum("+Inf");
00717   else if (is_minus_inf(s.c_str()))
00718     b=vnl_bignum("-Inf");
00719   else
00720     b.dtoBigNum(s.c_str());             // convert decimal to vnl_bignum
00721   return b;
00722 }
00723 
00724 
00725 //: Implicit conversion from a vnl_bignum to a short.
00726 vnl_bignum::operator short() const
00727 {
00728   int j = this->operator int();
00729   return (short)j;
00730 }
00731 
00732 
00733 //: Implicit conversion from a vnl_bignum to an int.
00734 vnl_bignum::operator int() const
00735 {
00736   int j = 0;
00737   for (Counter i = this->count; i > 0; )
00738     j = int(j*0x10000 + this->data[--i]);
00739   return (this->sign < 0) ? -j : j;
00740 }
00741 
00742 
00743 //: Implicit conversion from a vnl_bignum to a long.
00744 vnl_bignum::operator long() const
00745 {
00746   long l = 0;
00747   for (Counter i = this->count; i > 0; )
00748     l = l*0x10000L + this->data[--i];
00749   return (this->sign < 0) ? -l : l;
00750 }
00751 
00752 
00753 //: Implicit conversion from a vnl_bignum to a float.
00754 vnl_bignum::operator float() const
00755 {
00756   float f = 0.0f;
00757   for (Counter i = this->count; i > 0; )
00758     f = f*0x10000 + this->data[--i];
00759   if (this->is_infinity()) f = vcl_numeric_limits<float>::infinity();
00760   return (this->sign < 0) ? -f : f;
00761 }
00762 
00763 
00764 //: Implicit conversion from a vnl_bignum to a double.
00765 
00766 vnl_bignum::operator double() const
00767 {
00768   double d = 0.0;
00769   for (Counter i = this->count; i > 0; )
00770     d = d*0x10000 + this->data[--i];
00771   if (this->is_infinity()) d = vcl_numeric_limits<double>::infinity();
00772   return (this->sign < 0) ? -d : d;
00773 }
00774 
00775 //: Implicit conversion from a vnl_bignum to a long double.
00776 
00777 vnl_bignum::operator long double() const
00778 {
00779   long double d = 0.0;
00780   for (Counter i = this->count; i > 0; )
00781     d = d*0x10000 + this->data[--i];
00782   if (this->is_infinity()) d = vcl_numeric_limits<long double>::infinity();
00783   return (this->sign < 0) ? -d : d;
00784 }
00785 
00786 //: dump the contents of a vnl_bignum to a stream, default cout.
00787 
00788 void vnl_bignum::dump(vcl_ostream& os) const
00789 {
00790   os << "{count=" << this->count       // output count field
00791      << ", sign=" << this->sign        // output sign field
00792      << ", data=" << this->data        // output data pointer
00793      << ", value=" << *this
00794      << ", {";
00795   // format string == "%04X%s" or "%02X%s", etc.
00796   //  static char format_str[10] =
00797   //    {'%','0',char(2*2 + '0'),'X','%','s'};
00798   //  format_str[2] = char(2*2 + '0');
00799   if (this->count > 0) { // output data array
00800     os << vcl_hex << (this->data[this->count-1]);
00801     for (Counter i = this->count - 1; i > 0; --i) {
00802       os  << ',';
00803       if (this->data[i-1] < 0x10) os << '0';
00804       if (this->data[i-1] < 0x100) os << '0';
00805       if (this->data[i-1] < 0x1000) os << '0';
00806       os  << this->data[i-1];
00807     }
00808     os << vcl_dec;
00809   }
00810   os << "}}\n";                         // close brackets
00811 }
00812 
00813 
00814 //: Converts decimal string to a vnl_bignum.
00815 
00816 int vnl_bignum::dtoBigNum(const char *s)
00817 {
00818   this->resize(0); this->sign = 1;      // Reset number to 0.
00819   Counter len = 0;                      // No chars converted yet
00820   vnl_bignum sum;
00821   while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s; // skip whitespace
00822   if (*s == '-' || *s == '+') ++len;    // Skip over leading +,-
00823   while (s[len]>='0' && s[len]<='9') {  // While current char is digit
00824     *this *= vnl_bignum(10L),           // Shift vnl_bignum left a decimal
00825     add(*this,vnl_bignum(long(s[len++]-'0')),sum), // digit and add new digit
00826     *this = sum;
00827   }
00828   if (*s == '-') this->sign = -1;       // If s had leading -, note it
00829   return len;                           // Return # of chars processed
00830 }
00831 
00832 //: convert exponential string to a vnl_bignum
00833 
00834 void vnl_bignum::exptoBigNum(const char *s)
00835 {
00836   while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s; // skip whitespace
00837   Counter pos = Counter(this->dtoBigNum(s) + 1); // Convert the base, skip [eE]
00838   long pow = vcl_atol(s + pos);         // Convert the exponent to long
00839   while (pow-- > 0)                     // Raise vnl_bignum to the given
00840     *this = (*this) * 10L;              // power
00841 }
00842 
00843 
00844 //: convert hex character to integer hex value (ASCII or EBCDIC)
00845 // - Inputs:  character representation of a hex number
00846 // - Outputs: integer value of the hex number
00847 
00848 unsigned int ctox(int c)
00849 {
00850   if ('0' <= c && c <= '9')
00851     return c - '0';
00852   if ('a' <= c && c <= 'f')
00853     return c - 'a' + 10;
00854   return c - 'A' + 10;
00855 }
00856 
00857 //: convert hex string to vnl_bignum
00858 
00859 void vnl_bignum::xtoBigNum(const char *s)
00860 {
00861   this->resize(0); sign = 1;            // Reset number to 0.
00862   while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s; // skip whitespace
00863   Counter size = Counter(vcl_strlen(s));
00864   Counter len = 2;                      // skip leading "0x"
00865   while (len < size) {                  // While there are more chars
00866     (*this) = ((*this) * 16L) +         // Shift vnl_bignum left one hex
00867       vnl_bignum(long(ctox(s[len++]))); //   digit and add next digit
00868   }
00869 }
00870 
00871 
00872 //: convert octal string to vnl_bignum
00873 
00874 void vnl_bignum::otoBigNum(const char *s)
00875 {
00876   this->resize(0); sign = 1;           // Reset number to 0.
00877   while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s; // skip whitespace
00878   Counter size = Counter(vcl_strlen(s));
00879   Counter len = 0;                      // No chars converted yet
00880   while (len < size) {                  // While there are more chars
00881     (*this) = ((*this) * 8L) +          // Shift vnl_bignum left 1 oct dig.
00882       vnl_bignum(long(s[len++] - '0')); // Add next character value
00883   }
00884 }
00885 
00886 //: change the data allotment for a vnl_bignum
00887 
00888 void vnl_bignum::resize(short new_count)
00889 {
00890   assert(new_count >= 0);
00891   if (new_count == this->count) return;
00892   Data *new_data = (new_count > 0 ? new Data[new_count] : 0); // Allocate data if necessary
00893 
00894   if (this->count <= new_count) {       // Copy old data into new
00895     short i = 0;
00896     for (; i < this->count; i++)
00897       new_data[i] = this->data[i];
00898     for (; i < new_count; i++)
00899       new_data[i] = 0;
00900   }
00901   else {
00902     for (short i = 0; i < new_count; i++)
00903       new_data[i] = this->data[i];
00904   }
00905 
00906   delete [] this->data;                 // Get rid of old data
00907   this->data = new_data;                // Point to new data
00908   this->count = new_count;              // Save new count
00909 }
00910 
00911 
00912 //: trim non-infinite vnl_bignum of excess data allotment
00913 
00914 vnl_bignum& vnl_bignum::trim()
00915 {
00916   Counter i = this->count;
00917   for (; i > 0; i--)                    // Skip over high-order words
00918     if (this->data[i - 1] != 0) break;  //   that are zero
00919   if (i < this->count) {                // If there are some such words
00920     this->count = i;                    // Update the count
00921     Data *new_data = (i > 0 ? new Data[i] : 0); // Allocate data if necessary
00922     for (; i > 0; i--)                  // Copy old data into new
00923       new_data[i - 1] = this->data[i - 1];
00924     delete [] this->data;               // Delete old data
00925     this->data = new_data;              // Point to new data
00926   }
00927   return *this;                         // return reference to vnl_bignum
00928 }
00929 
00930 
00931 //: add two non-infinite vnl_bignum values and save their sum
00932 
00933 void add(const vnl_bignum& b1, const vnl_bignum& b2, vnl_bignum& sum)
00934 {
00935   const vnl_bignum *bmax, *bmin;        // Determine which of the two
00936   if (b1.count >= b2.count) {           //   addends has the most
00937     bmax = &b1;                         //   data.
00938     bmin = &b2;
00939   }
00940   else {
00941     bmax = &b2;
00942     bmin = &b1;
00943   }
00944   sum.resize(bmax->count);              // Allocate data for their sum
00945   unsigned long temp, carry = 0;
00946   Counter i = 0;
00947   while (i < bmin->count) {             // Add, element by element.
00948     // Add both elements and carry
00949     temp = (unsigned long)b1.data[i] + (unsigned long)b2.data[i] + carry;
00950     carry = temp/0x10000L;              // keep track of the carry
00951     sum.data[i] = Data(temp);           // store sum
00952     i++;                                // go to next element
00953   }
00954   while (i < bmax->count) {             // bmin has no more elements
00955     temp = bmax->data[i] + carry;       // propagate the carry through
00956     carry = temp/0x10000L;              // the rest of bmax's elements
00957     sum.data[i] = Data(temp);           // store sum
00958     i++;
00959   }
00960   if (carry) {                          // if carry left over
00961     sum.resize(bmax->count + 1);        //   allocate another word
00962     sum.data[bmax->count] = 1;          //   save the carry in it
00963   }
00964 }
00965 
00966 //: Add 1 to bnum (unsigned, non-infinite)
00967 void increment(vnl_bignum& bnum)
00968 {
00969   Counter i = 0;
00970   unsigned long carry = 1;
00971   while (i < bnum.count && carry) {             // increment, element by element.
00972     unsigned long temp = (unsigned long)bnum.data[i] + carry;
00973     carry = temp/0x10000L;
00974     bnum.data[i] = (Data)temp;
00975     ++i;
00976   }
00977   if (carry)
00978   {
00979     bnum.resize(bnum.count+1);
00980     bnum.data[bnum.count-1] = 1;
00981   }
00982 }
00983 
00984 
00985 //: subtract bmin from bmax (unsigned, non-infinite), result in diff
00986 
00987 void subtract(const vnl_bignum& bmax, const vnl_bignum& bmin, vnl_bignum& diff)
00988 {
00989   diff.resize(bmax.count);                      // Allocate data for difference
00990   unsigned long temp;
00991   int borrow = 0;
00992   Counter i = 0;
00993   for (; i < bmin.count; i++) {                 // Subtract word by word.
00994     temp = (unsigned long)bmax.data[i] + 0x10000L - borrow; // Add radix to bmax's data
00995     temp -= (unsigned long)bmin.data[i];        // Subtract off bmin's data
00996     borrow = (temp/0x10000L == 0);              // Did we have to borrow?
00997     diff.data[i] = (Data) temp;                 // Reduce modulo radix and save
00998   }
00999   for (; i < bmax.count; i++) {                 // No more data for bmin
01000     temp = (unsigned long)bmax.data[i] + 0x10000L - borrow; // Propagate the borrow through
01001     borrow = (temp/0x10000L == 0);              //   rest of bmax's data
01002     diff.data[i] = (Data) temp;
01003   }
01004   diff.trim();                                  // Done. Now trim excess data
01005 }
01006 
01007 
01008 //: Subtract 1 from bnum (unsigned, non-infinite, non-zero)
01009 void decrement(vnl_bignum& bnum)
01010 {
01011   Counter i = 0;
01012   unsigned long borrow = 1;
01013   while (i < bnum.count && borrow) {            // decrement, element by element.
01014     unsigned long temp = (unsigned long)bnum.data[i] + 0x10000L - borrow;
01015     borrow = (temp/0x10000L == 0);              // Did we have to borrow?
01016     bnum.data[i] = (Data)temp;                  // Reduce modulo radix and save
01017     ++i;
01018   }
01019   bnum.trim();                                  // Done. Now trim excess data
01020   if (bnum.count==0) bnum.sign=+1;              // Re-establish sign invariant
01021 }
01022 
01023 
01024 //: compare absolute values of two vnl_bignums
01025 // Outputs:  result of comparison:  -1 if abs(b1) < abs(b2)
01026 //                                   0 if abs(b1) == abs(b2)
01027 //                                  +1 if abs(b1) > abs(b2)
01028 
01029 int magnitude_cmp(const vnl_bignum& b1, const vnl_bignum& b2)
01030 {
01031   if (b1.is_infinity()) return b2.is_infinity() ? 0 : 1;
01032   if (b2.is_infinity()) return -1;
01033   if (b1.count > b2.count) return 1;    // If one has more data than
01034   if (b2.count > b1.count) return -1;   //   the other, it wins
01035   Counter i = b1.count;                 // Else same number of elements
01036   while (i > 0) {                       // Do lexicographic comparison
01037     if (b1.data[i - 1] > b2.data[i - 1])
01038       return 1;
01039     else if (b1.data[i - 1] < b2.data[i - 1])
01040       return -1;
01041     i--;
01042   }                                     // No data, or all elements same
01043   return 0;                             //  so must be equal
01044 }
01045 
01046 
01047 //: multiply a non-infinite vnl_bignum by a "single digit"
01048 // - Inputs:  vnl_bignum reference, single word multiplier, reference to the product,
01049 //            and index of starting storage location to use in product
01050 
01051 void multiply_aux(const vnl_bignum& b, Data d, vnl_bignum& prod, Counter i)
01052 {
01053   // this function works just like normal multiplication by hand, in that the
01054   // top number is multiplied by the first digit of the bottom number, then the
01055   // second digit, and so on.  The only difference is that instead of doing all
01056   // of the multiplication before adding the rows, addition is done
01057   // concurrently.
01058   if (i == 0) {                         // if index is zero
01059     Counter j = 0;                      //   then zero out all of
01060     while (j < prod.count)              //   prod's data elements
01061       prod.data[j++] = 0;
01062   }
01063   if (d != 0) {                         // if d == 0, nothing to do
01064     unsigned long temp;
01065     Data carry = 0;
01066 
01067     Counter j = 0;
01068     for (; j < b.count; j++) {
01069       // for each of b's data elements, multiply times d and add running product
01070       temp = (unsigned long)b.data[j] * (unsigned long)d
01071            + (unsigned long)prod.data[i + j] + carry;
01072       prod.data[i + j] = Data(temp % 0x10000L); //   store result in product
01073       carry = Data(temp/0x10000L);              //   keep track of carry
01074     }
01075     if (i+j < prod.count)
01076       prod.data[i + j] = carry;                 // Done.  Store the final carry
01077   }
01078 }
01079 
01080 //: normalize two vnl_bignums
01081 // (Refer to Knuth, V.2, Section 4.3.1, Algorithm D for details.
01082 //  A digit here is one data element in the radix 2**2.)
01083 // - Inputs:  references to two vnl_bignums b1, b2
01084 // - Outputs: their normalized counterparts u and v,
01085 //            and the integral normalization factor used
01086 
01087 Data normalize(const vnl_bignum& b1, const vnl_bignum& b2, vnl_bignum& u, vnl_bignum& v)
01088 {
01089   Data d = Data(0x10000L/((unsigned long)(b2.data[b2.count - 1]) + 1L)); // Calculate normalization factor.
01090   u.resize(b1.count + 1);                       // Get data for u (plus extra)
01091   v.resize(b2.count);                           // Get data for v
01092   u.data[b1.count] = 0;                         // Set u's leading digit to 0
01093   multiply_aux(b1,d,u,0);                       // u = b1 * d
01094   multiply_aux(b2,d,v,0);                       // v = b2 * d
01095   return d;                                     // return normalization factor
01096 }
01097 
01098 
01099 //: divide a vnl_bignum by a "single digit"
01100 // (Refer to Knuth, V.2, Section 4.3.2, exercise 16 for details.
01101 //  A digit here is one data element in the radix 2**2.)
01102 // - Inputs:  reference to vnl_bignum dividend, single digit divisor d, vnl_bignum
01103 //            quotient, and single digit remainder r
01104 
01105 void divide_aux(const vnl_bignum& b1, Data d, vnl_bignum& q, Data& r)
01106 {
01107   r = 0;                                        // init remainder to zero
01108   unsigned long temp;
01109   for (Counter j = b1.count; j > 0; j--) {
01110     temp = (unsigned long)r*0x10000L + (unsigned long)b1.data[j - 1]; // get remainder, append next
01111     if (j < 1 + q.count)
01112       q.data[j - 1] = Data(temp/d);             //   digit, then divide
01113     r = Data(temp % d);                         // calculate new remainder
01114   }
01115 }
01116 
01117 
01118 //: estimate next dividend
01119 // (Refer to Knuth, V.2, Section 4.3.1, Algorithm D for details.
01120 //  This function estimates how many times v goes into u, starting at u's
01121 //  jth digit.  A digit here is actually a data element, thought of as
01122 //  being in the radix 2**2.)
01123 // - Inputs:  reference to vnl_bignum dividend and divisor and starting digit j
01124 // - Outputs: estimated number of times v goes into u
01125 
01126 Data estimate_q_hat(const vnl_bignum& u, const vnl_bignum& v, Counter j)
01127 {
01128   Data q_hat,
01129        v1 = v.data[v.count - 1],                // localize frequent data
01130        v2 = v.data[v.count - 2],
01131        u0 = u.data[u.count - 1 - j],
01132        u1 = u.data[u.count - 2 - j],
01133        u2 = u.data[u.count - 3 - j];
01134 
01135   // Initial Knuth estimate, usually correct
01136   q_hat = (u0 == v1 ? Data(0xffff) : Data(((unsigned long)u0 * 0x10000L + (unsigned long)u1) / (unsigned long)v1));
01137 
01138   // high speed test to determine most of the cases in which initial
01139   // estimate is too large.  Eliminates most cases in which q_hat is one too
01140   // large.  Eliminates all cases in which q_hat is two too large.  The test
01141   // looks hairy because we have to watch out for overflow.  In the book, this
01142   // test is the simple inequality:
01143   //     v2*q_hat > (u0*0x10000L + u1 - q_hat*v1)*0x10000L + u2.
01144   // If the inequality is true, decrease q_hat by 1.  If inequality is still
01145   // true, decrease q_hat again.
01146   unsigned long lhs, rhs;               // lhs, rhs of Knuth inequality
01147   for (Counter i = 0; i < 2; i++) {     // loop at most twice
01148     lhs = (unsigned long)v2 * (unsigned long)q_hat;       // Calculate left-hand side of ineq.
01149     rhs = (unsigned long)u0 * 0x10000L +(unsigned long)u1;// Calculate part of right-hand side
01150     rhs -= ((unsigned long)q_hat * (unsigned long)v1);    // Now subtract off part
01151 
01152     if ( rhs >= 0x10000L )              // if multiplication with 0x10000L causes overflow
01153       break;                            //   then rhs > lhs, so test fails
01154     rhs *= 0x10000L;                    // No overflow:  ok to multiply
01155 
01156     if (rhs > rhs + (unsigned long)u2)  // if addition yields overflow
01157       break;                            //   then rhs > lhs, so test fails
01158     rhs += u2;                          // No overflow: ok to add.
01159     if (lhs <= rhs)                     // if lhs <= rhs
01160       break;                            //   test fails
01161     q_hat--;                            // Test passes:  decrement q_hat
01162   }                                     // Loop again
01163   return q_hat;                         // Return estimate
01164 }
01165 
01166 
01167 //: calculate u - v*q_hat
01168 // (Refer to Knuth, V. 2, Section 4.3.1, Algorithm D for details.
01169 //  A digit here is a data element, thought of as being in the radix 2**2.)
01170 // - Inputs:  reference to vnl_bignum dividend, divisor, estimated result, and index
01171 //            into jth digit of dividend
01172 // - Outputs: true number of times v goes into u
01173 
01174 Data multiply_subtract(vnl_bignum& u, const vnl_bignum& v, Data q_hat, Counter j)
01175 {
01176   // At this point it has been estimated that v goes into the jth and higher
01177   // digits of u about q_hat times, and in fact that q_hat is either the
01178   // correct number of times or one too large.
01179 
01180   if (q_hat == 0) return q_hat;         // if q_hat 0, nothing to do
01181   vnl_bignum rslt;                      // create a temporary vnl_bignum
01182   Counter tmpcnt;
01183   rslt.resize(v.count + 1u);            // allocate data for it
01184 
01185   // simultaneous computation of u - v*q_hat
01186   unsigned long prod, diff;
01187   Data carry = 0, borrow = 0;
01188   Counter i = 0;
01189   for (; i < v.count; ++i) {
01190     // for each digit of v, multiply it by q_hat and subtract the result
01191     prod = (unsigned long)v.data[i] * (unsigned long)q_hat + carry;
01192     diff = (unsigned long)u.data[u.count - v.count - 1 - j + i] + (0x10000L - (unsigned long)borrow);
01193     diff -= (unsigned long)Data(prod);  //   form proper digit of u
01194     rslt.data[i] = Data(diff);          //   save the result
01195     borrow = (diff/0x10000L == 0) ? 1 : 0; //   keep track of any borrows
01196     carry = Data(prod/0x10000L);        //   keep track of carries
01197   }
01198   tmpcnt = Counter(u.count - v.count + i - j - 1);
01199   diff = (unsigned long)u.data[tmpcnt]  //  special case for the last digit
01200          + (0x10000L - (unsigned long)borrow);
01201   diff -= (unsigned long)carry;
01202   rslt.data[i] = Data(diff);
01203   borrow = (diff/0x10000L == 0) ? 1 : 0;
01204 
01205   // A leftover borrow indicates that u - v*q_hat is negative, i.e., that
01206   // q_hat was one too large.  So to get correct result, decrement q_hat and
01207   // add back one multiple of v
01208   if (borrow) {
01209     q_hat--;
01210     carry = 0;
01211     unsigned long sum;
01212     for (i = 0; i < v.count; ++i) {
01213       sum = (unsigned long)rslt.data[i] + (unsigned long)v.data[i] + carry;
01214       carry = Data(sum/0x10000L);
01215       u.data[u.count - v.count + i - 1 - j] = Data(sum);
01216     }
01217     u.data[u.count - v.count + i - 1 - j] = rslt.data[i] + carry;
01218   }
01219   else {                                // otherwise, result is ok
01220     for (i = 0; i < rslt.count; ++i)    // store result back into u
01221       u.data[u.count - v.count + i - 1 - j] = rslt.data[i];
01222   }
01223   return q_hat;                         // return corrected q_hat
01224 }
01225 
01226 
01227 //: divide b2 into b1, getting quotient q and remainder r.
01228 // (Refer to Knuth, V.2, Section 4.3.1, Algorithm D for details.
01229 //  This function implements Algorithm D.)
01230 // - Inputs: references to a vnl_bignum dividend b1, divisor b2, quotient q, and
01231 //           remainder r.
01232 
01233 void divide(const vnl_bignum& b1, const vnl_bignum& b2, vnl_bignum& q, vnl_bignum& r)
01234 {
01235   // Note that q or r may *not* be identical to either b1 or b2 !
01236   assert(&b1 != &q && &b2 != &q && &b1 != &r && &b2 != &r);
01237   q = r = 0L;
01238   if (b1 == 0L)                      // If divisor is zero
01239     return;                          //   return zero quotient and remainder
01240   int mag = magnitude_cmp(b1,b2);    // Compare magnitudes
01241   if (mag < 0)                       // if abs(b1) < abs(b2)
01242     r = b1;                          //   return zero quotient, b1 remainder
01243   else if (mag == 0)                 // if abs(b1) == abs(b2)
01244     q = 1L;                          //   quotient is 1, remainder is 0
01245   else {                             // otherwise abs(b1) > abs(b2), so divide
01246     q.resize(b1.count + 1u - b2.count); // Allocate quotient
01247     r.resize(b2.count);              // Allocate remainder
01248     if (b2.count == 1) {                        // Single digit divisor?
01249       divide_aux(b1,b2.data[0],q,r.data[0]);    // Do single digit divide
01250     }
01251     else {                                      // Else full-blown divide
01252       vnl_bignum u,v;
01253 #ifdef DEBUG
01254       vcl_cerr << "\nvnl_bignum::divide: b1 ="; b1.dump(vcl_cerr);
01255       vcl_cerr << "vnl_bignum::divide: b2 ="; b2.dump(vcl_cerr);
01256 #endif
01257       Data d = normalize(b1,b2,u,v);            // Set u = b1*d, v = b2*d
01258 #ifdef DEBUG
01259       vcl_cerr << "vnl_bignum::divide: d = " << d << vcl_hex << " (0x" << d << ")\n" << vcl_dec;
01260       vcl_cerr << "vnl_bignum::divide: u = "; u.dump(vcl_cerr);
01261       vcl_cerr << "vnl_bignum::divide: v = "; v.dump(vcl_cerr);
01262 #endif
01263       Counter j = 0;
01264       while (j <= b1.count - b2.count) {        // Main division loop
01265         Data q_hat = estimate_q_hat(u,v,j);     // Estimate # times v divides u
01266         q.data[q.count - 1 - j] =               // Do division, get true answ.
01267           multiply_subtract(u,v,q_hat,j);
01268         j++;
01269 #ifdef DEBUG
01270         vcl_cerr << "vnl_bignum::divide: q_hat = " << q_hat << vcl_hex << " (0x" << q_hat << ")\n" << vcl_dec;
01271         vcl_cerr << "vnl_bignum::divide: u = "; u.dump(vcl_cerr);
01272 #endif
01273       }
01274       static Data dufus;                // dummy variable
01275       divide_aux(u,d,r,dufus);          // Unnormalize u for remainder
01276 
01277 #ifdef DEBUG
01278       vcl_cerr << "vnl_bignum::divide: q = "; q.dump(vcl_cerr);
01279       vcl_cerr << "vnl_bignum::divide: r = "; r.dump(vcl_cerr);
01280 #endif
01281     }
01282     q.trim();                           // Trim leading zeros of quot.
01283     r.trim();                           // Trim leading zeros of rem.
01284   }
01285   q.sign = r.sign = b1.sign * b2.sign;  // Calculate signs
01286 }
01287 
01288 
01289 //: left shift (arithmetic) non-infinite vnl_bignum by positive number.
01290 // - Inputs:  reference to vnl_bignum, positive shift value
01291 // - Outputs: vnl_bignum, multiplied by the corresponding power of two
01292 
01293 vnl_bignum left_shift(const vnl_bignum& b1, int l)
01294 {
01295   // to carry out this arithmetic left shift, we cheat.  Instead of physically
01296   // shifting the data array l bits to the left, we shift just enough to get
01297   // the correct word alignment, and then pad the array on the right with as
01298   // many zeros as we need.
01299   vnl_bignum rslt;                      // result of shift
01300   rslt.sign = b1.sign;                  // result follows sign of input
01301   Counter growth = Counter(l / 16);     // # of words rslt will grow by
01302   Data shift = Data(l % 16);            // amount to actually shift
01303   Data rshift = Data(16 - shift);       // amount to shift next word by
01304   Data carry = Data(                    // value that will be shifted
01305     b1.data[b1.count - 1] >> (16 - shift));// out end of current array
01306   rslt.resize(b1.count + growth + (carry ? 1u : 0u)); // allocate new data array
01307   Counter i = 0;
01308   while (i < growth)                            // zero out padded elements
01309     rslt.data[i++] = 0;
01310   rslt.data[i++] = b1.data[0] << shift;         // shift first non-zero element
01311   while (i < rslt.count - 1) {                  // for remaining data words
01312     rslt.data[i] = (b1.data[i - growth] << shift) + // shift current data word
01313       (b1.data[i - 1 - growth] >> rshift);          // propagate adjacent
01314     i++;                                        // carry into current word
01315   }
01316   if (i < rslt.count) {
01317     if (carry)                                  // if last word had overflow
01318       rslt.data[i] = carry;                     //   store it new data
01319     else                                        // otherwise,
01320       rslt.data[i] = (b1.data[i - growth] << shift) // do like the rest
01321                    + (b1.data[i - 1 - growth] >> rshift);
01322   }
01323   vnl_bignum& result = *((vnl_bignum*) &rslt);// same physical object
01324   return result;                              // shallow swap on return
01325 }
01326 
01327 
01328 //: right shift (arithmetic) non-infinite vnl_bignum by positive number.
01329 // - Inputs:  reference to vnl_bignum, positive shift value
01330 // - Outputs: vnl_bignum, divided by the corresponding power of two
01331 
01332 vnl_bignum right_shift(const vnl_bignum& b1, int l)
01333 {
01334   vnl_bignum rslt;                              // result of shift
01335   Counter shrinkage = Counter(l / 16);          // # of words rslt will shrink
01336   Data shift = Data(l % 16);                    // amount to actually shift
01337   Data dregs = Data(b1.data[b1.count-1] >> shift);// high end data to save
01338   if (shrinkage + (dregs == 0) < b1.count) {    // if not all data shifted out
01339     rslt.sign = b1.sign;                        // rslt follows sign of input
01340     rslt.resize(b1.count - shrinkage - (dregs == 0 ? 1 : 0)); // allocate new data
01341     Data lshift = Data(16 - shift);             // amount to shift high word
01342     Counter i = 0;
01343     while (i < rslt.count - 1) {                // shift current word
01344       rslt.data[i] = (b1.data[i + shrinkage] >> shift) + // propagate adjacent
01345         (b1.data[i + shrinkage + 1u] << lshift); // word into current word
01346       i++;
01347     }
01348     if (dregs)                                  // don't lose dregs
01349       rslt.data[i] = dregs;
01350     else {
01351       rslt.data[i] = (b1.data[i + shrinkage] >> shift) +
01352         (b1.data[i + shrinkage + 1u] << lshift);
01353     }
01354   }
01355   vnl_bignum& result = *((vnl_bignum*) &rslt);  // same physical object
01356   return result;                                // shallow swap on return
01357 }

Generated on Mon Mar 8 05:06:44 2010 for core/vnl by  doxygen 1.5.1