00001
00002 #include "vnl_bignum.h"
00003
00004
00005
00006 #include <vcl_cstdlib.h>
00007 #include <vcl_cstring.h>
00008 #include <vcl_cmath.h>
00009 #include <vcl_algorithm.h>
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>
00015
00016 typedef unsigned short Counter;
00017 typedef unsigned short Data;
00018
00019
00020
00021 vnl_bignum::vnl_bignum()
00022 : count(0), sign(1), data(0)
00023 {
00024 }
00025
00026
00027
00028 vnl_bignum::vnl_bignum(long l)
00029 : count(0), sign(1), data(0)
00030 {
00031 if (l < 0) {
00032 l = -l;
00033 this->sign = -1;
00034 }
00035 Data buf[sizeof(l)];
00036 Counter i = 0;
00037 while (l) {
00038 assert(i < sizeof(l));
00039 buf[i] = Data(l);
00040 l >>= 16;
00041 ++i;
00042 }
00043 if (i > 0)
00044 this->data = new Data[this->count=i];
00045
00046 while (i--)
00047 this->data[i] = buf[i];
00048 }
00049
00050
00051
00052 vnl_bignum::vnl_bignum(int l)
00053 : count(0), sign(1), data(0)
00054 {
00055 if (l < 0) {
00056 l = -l;
00057 this->sign = -1;
00058 }
00059 Data buf[sizeof(l)];
00060 Counter i = 0;
00061 while (l) {
00062 assert(i < sizeof(l));
00063 buf[i] = Data(l);
00064 l >>= 16;
00065 i++;
00066 }
00067 if (i > 0)
00068 this->data = new Data[this->count=i];
00069
00070 while (i--)
00071 this->data[i] = buf[i];
00072 }
00073
00074
00075
00076 vnl_bignum::vnl_bignum(unsigned long l)
00077 : count(0), sign(1), data(0)
00078 {
00079 Data buf[sizeof(l)];
00080 Counter i = 0;
00081 while (l) {
00082 assert(i < sizeof(l));
00083 buf[i] = Data(l);
00084 l >>= 16;
00085 i++;
00086 }
00087 if (i > 0)
00088 this->data = new Data[this->count=i];
00089
00090 while (i--)
00091 this->data[i] = buf[i];
00092 }
00093
00094
00095
00096 vnl_bignum::vnl_bignum(unsigned int l)
00097 : count(0), sign(1), data(0)
00098 {
00099 Data buf[sizeof(l)];
00100 Counter i = 0;
00101 while (l) {
00102 assert(i < sizeof(l));
00103 buf[i] = Data(l);
00104 l >>= 16;
00105 i++;
00106 }
00107 if (i > 0)
00108 this->data = new Data[this->count=i];
00109
00110 while (i--)
00111 this->data[i] = buf[i];
00112 }
00113
00114
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) {
00121 d = -d;
00122 this->sign = -1;
00123 }
00124
00125 if (!vnl_math_isfinite(d)) {
00126
00127
00128 this->count = 1;
00129 this->data = new Data[1];
00130 this->data[0] = 0;
00131 }
00132 else if (d >= 1.0) {
00133
00134 vcl_vector<Data> buf;
00135 while (d >= 1.0) {
00136 buf.push_back( Data(vcl_fmod(d,0x10000L)) );
00137 d /= 0x10000L;
00138 }
00139
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
00147
00148 vnl_bignum::vnl_bignum(double d)
00149 : count(0), sign(1), data(0)
00150 {
00151 if (d < 0.0) {
00152 d = -d;
00153 this->sign = -1;
00154 }
00155
00156 if (!vnl_math_isfinite(d)) {
00157
00158
00159 this->count = 1;
00160 this->data = new Data[1];
00161 this->data[0] = 0;
00162 }
00163 else if (d >= 1.0) {
00164
00165 vcl_vector<Data> buf;
00166 while (d >= 1.0) {
00167 buf.push_back( Data(vcl_fmod(d,0x10000L)) );
00168 d /= 0x10000L;
00169 }
00170
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
00178
00179 vnl_bignum::vnl_bignum(long double d)
00180 : count(0), sign(1), data(0)
00181 {
00182 if (d < 0.0) {
00183 d = -d;
00184 this->sign = -1;
00185 }
00186
00187 if (!vnl_math_isfinite(d)) {
00188
00189
00190 this->count = 1;
00191 this->data = new Data[1];
00192 this->data[0] = 0;
00193 }
00194 else if (d >= 1.0) {
00195
00196 vcl_vector<Data> buf;
00197 while (d >= 1.0) {
00198 buf.push_back( Data(vcl_fmod(d,0x10000L)) );
00199 d /= 0x10000L;
00200 }
00201
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
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]);
00268 if (*s) ++s;
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);
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
00375
00376 vnl_bignum::vnl_bignum(const char *s)
00377 : count(0), sign(1), data(0)
00378 {
00379
00380
00381
00382
00383
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))
00390 this->dtoBigNum(s);
00391 else if (is_exponential(s))
00392 this->exptoBigNum(s);
00393 else if (is_hexadecimal(s))
00394 this->xtoBigNum(s);
00395 else if (is_octal(s))
00396 this->otoBigNum(s);
00397 else {
00398 vcl_cerr << "Cannot convert string " << s << " to vnl_bignum\n";
00399 }
00400 }
00401
00402
00403
00404 vcl_istream& operator>>(vcl_istream& is, vnl_bignum& x)
00405 {
00406
00407
00408
00409
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))
00419 x.exptoBigNum(rt);
00420 else if (is_decimal(rt,&isp))
00421 x.dtoBigNum(rt);
00422 else if (is_hexadecimal(rt,&isp))
00423 x.xtoBigNum(rt);
00424 else if (is_octal(rt,&isp))
00425 x.otoBigNum(rt);
00426 else {
00427 vcl_cerr << "Cannot convert string " << rt << " to vnl_bignum\n";
00428 }
00429 return is;
00430 }
00431
00432
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];
00439 for (Counter i = 0; i < this->count; ++i)
00440 this->data[i] = b.data[i];
00441 }
00442 else {
00443 this->data = 0;
00444 }
00445 }
00446
00447
00448
00449 vnl_bignum::~vnl_bignum()
00450 {
00451 delete [] this->data; this->count = 0;
00452 }
00453
00454
00455
00456 vnl_bignum& vnl_bignum::operator=(const vnl_bignum& rhs)
00457 {
00458 if (this != &rhs) {
00459 delete[] this->data;
00460 this->count = rhs.count;
00461 this->data = rhs.data ? new Data[rhs.count] : 0;
00462 for (Counter i = 0; i < rhs.count; ++i)
00463 this->data[i] = rhs.data[i];
00464 this->sign = rhs.sign;
00465 }
00466 return *this;
00467 }
00468
00469
00470
00471 vnl_bignum vnl_bignum::operator-() const
00472 {
00473 vnl_bignum neg(*this);
00474 if (neg.count)
00475 neg.sign = -neg.sign;
00476 return neg;
00477 }
00478
00479
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
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
00518
00519 vnl_bignum vnl_bignum::operator+(const vnl_bignum& b) const
00520 {
00521
00522 assert (! b.is_minus_infinity() || ! this->is_plus_infinity() );
00523 assert (! b.is_plus_infinity() || ! this->is_minus_infinity() );
00524 if (b.is_infinity()) { return b; }
00525 if (this->is_infinity()) { return *this; }
00526
00527 vnl_bignum sum;
00528 if (this->sign == b.sign) {
00529 add(*this,b,sum);
00530 sum.sign = this->sign;
00531 }
00532 else {
00533 int mag = magnitude_cmp(*this,b);
00534 if (mag > 0) {
00535 subtract(*this,b,sum);
00536 sum.sign = this->sign;
00537 }
00538 else if (mag < 0) {
00539 subtract(b,*this,sum);
00540 sum.sign = b.sign;
00541 }
00542 }
00543 return sum;
00544 }
00545
00546
00547
00548 vnl_bignum& vnl_bignum::operator*=(const vnl_bignum& b)
00549 {
00550
00551 assert (! b.is_infinity() || this->count != 0 );
00552 assert (! this->is_infinity() || b.count != 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);
00560 for (Counter i = 0; i < b.count; i++)
00561 multiply_aux(*this, b.data[i], prod, i);
00562 prod.sign = this->sign * b.sign;
00563 prod.trim();
00564 return (*this)=prod;
00565 }
00566
00567
00568
00569 vnl_bignum& vnl_bignum::operator/=(const vnl_bignum& b)
00570 {
00571
00572 assert (! b.is_infinity() || ! this->is_infinity() );
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);
00576 if (b.count == 0)
00577 return (*this) = (this->sign < 0 ? vnl_bignum("-Inf") : vnl_bignum("+Inf"));
00578
00579 vnl_bignum quot, r;
00580 divide(*this,b,quot,r);
00581 return (*this) = quot;
00582 }
00583
00584
00585
00586 vnl_bignum& vnl_bignum::operator%=(const vnl_bignum& b)
00587 {
00588
00589 assert (! b.is_infinity() || ! this->is_infinity() );
00590 if (b.is_infinity()) return *this;
00591 if (this->is_infinity()) return (*this) = 0L;
00592 assert (b.count!=0 || this->count != 0);
00593 if (b.count == 0) return (*this) = 0L;
00594
00595 vnl_bignum remain, q;
00596 divide(*this,b,q,remain);
00597 return (*this) = remain;
00598 }
00599
00600
00601
00602 vnl_bignum vnl_bignum::operator<<(int l) const
00603 {
00604
00605 if (this->is_infinity()) return *this;
00606
00607 if (l == 0 || *this == 0L)
00608 return *this;
00609 if (l < 0)
00610 return right_shift(*this,-l);
00611 else
00612 return left_shift(*this,l);
00613 }
00614
00615
00616
00617 vnl_bignum vnl_bignum::operator>>(int l) const
00618 {
00619
00620 if (this->is_infinity()) return *this;
00621
00622 if (l == 0 || *this == 0L)
00623 return *this;
00624 if (l < 0)
00625 return left_shift(*this,-l);
00626 else
00627 return right_shift(*this,l);
00628 }
00629
00630
00631
00632 bool vnl_bignum::operator==(const vnl_bignum& rhs) const
00633 {
00634 if (this != &rhs) {
00635 if (this->sign != rhs.sign) return false;
00636 if (this->count != rhs.count) return false;
00637 for (Counter i = 0; i < this->count; i++)
00638 if (this->data[i] != rhs.data[i]) return false;
00639 }
00640 return true;
00641 }
00642
00643
00644
00645 bool vnl_bignum::operator<(const vnl_bignum& rhs) const
00646 {
00647 if (this->sign < rhs.sign) return true;
00648 if (this->sign > rhs.sign) return false;
00649 if (this->sign == 1)
00650 return magnitude_cmp(*this,rhs) < 0;
00651 else
00652 return magnitude_cmp(*this,rhs) > 0;
00653 }
00654
00655
00656
00657 vcl_ostream& operator<<(vcl_ostream& os, const vnl_bignum& b)
00658 {
00659 vnl_bignum d = b;
00660 if (d.sign == -1) {
00661 os << '-';
00662 d.sign = 1;
00663 }
00664 if (d.is_infinity()) return os<<"Inf";
00665 vnl_bignum q,r;
00666 char *cbuf = new char[5 * (b.count+1)];
00667 Counter i = 0;
00668 do {
00669 divide(d,10L,q,r);
00670 cbuf[i++] = char(long(r) + '0');
00671 d = q;
00672 q = r = 0L;
00673 } while (d != 0L);
00674 do {
00675 os << cbuf[--i];
00676 } while (i);
00677 delete [] cbuf;
00678 return os;
00679 }
00680
00681
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;
00687
00688 vnl_bignum d = b;
00689 if (d.sign == -1) {
00690 s.insert(insert_point,"-");
00691 d.sign = 1;
00692 ++insert_point;
00693 }
00694 if (d.is_infinity()) return s+="Inf";
00695 vnl_bignum q,r;
00696 do {
00697 divide(d,10L,q,r);
00698 s.insert(insert_point, 1, char('0'+long(r)));
00699 d = q;
00700 q = r = 0L;
00701 } while (d != 0L);
00702 return s;
00703 }
00704
00705
00706
00707 vnl_bignum& vnl_bignum_from_string(vnl_bignum& b, const vcl_string& s)
00708 {
00709
00710
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());
00718 return b;
00719 }
00720
00721
00722
00723 vnl_bignum::operator short() const
00724 {
00725 int j = this->operator int();
00726 return (short)j;
00727 }
00728
00729
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
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
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
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
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
00783
00784 void vnl_bignum::dump(vcl_ostream& os) const
00785 {
00786 os << "{count=" << this->count
00787 << ", sign=" << this->sign
00788 << ", data=" << this->data
00789 << ", value=" << *this
00790 << ", {";
00791
00792
00793
00794
00795 if (this->count > 0) {
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";
00807 }
00808
00809
00810
00811 int vnl_bignum::dtoBigNum(const char *s)
00812 {
00813 this->resize(0); this->sign = 1;
00814 Counter len = 0;
00815 vnl_bignum sum;
00816 while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s;
00817 if (*s == '-' || *s == '+') ++len;
00818 while (s[len]>='0' && s[len]<='9') {
00819 *this *= vnl_bignum(10L),
00820 add(*this,vnl_bignum(long(s[len++]-'0')),sum),
00821 *this = sum;
00822 }
00823 if (*s == '-') this->sign = -1;
00824 return len;
00825 }
00826
00827
00828
00829 void vnl_bignum::exptoBigNum(const char *s)
00830 {
00831 while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s;
00832 Counter pos = Counter(this->dtoBigNum(s) + 1);
00833 long pow = vcl_atol(s + pos);
00834 while (pow-- > 0)
00835 *this = (*this) * 10L;
00836 }
00837
00838
00839
00840
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
00852
00853 void vnl_bignum::xtoBigNum(const char *s)
00854 {
00855 this->resize(0); sign = 1;
00856 while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s;
00857 Counter size = Counter(vcl_strlen(s));
00858 Counter len = 2;
00859 while (len < size) {
00860 (*this) = ((*this) * 16L) +
00861 vnl_bignum(long(ctox(s[len++])));
00862 }
00863 }
00864
00865
00866
00867 void vnl_bignum::otoBigNum(const char *s)
00868 {
00869 this->resize(0); sign = 1;
00870 while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') ++s;
00871 Counter size = Counter(vcl_strlen(s));
00872 Counter len = 0;
00873 while (len < size) {
00874 (*this) = ((*this) * 8L) +
00875 vnl_bignum(long(s[len++] - '0'));
00876 }
00877 }
00878
00879
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);
00886
00887 if (this->count <= new_count) {
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;
00900 this->data = new_data;
00901 this->count = new_count;
00902 }
00903
00904
00905
00906 vnl_bignum& vnl_bignum::trim()
00907 {
00908 Counter i = this->count;
00909 for (; i > 0; i--)
00910 if (this->data[i - 1] != 0) break;
00911 if (i < this->count) {
00912 this->count = i;
00913 Data *new_data = (i > 0 ? new Data[i] : 0);
00914 for (; i > 0; i--)
00915 new_data[i - 1] = this->data[i - 1];
00916 delete [] this->data;
00917 this->data = new_data;
00918 }
00919 return *this;
00920 }
00921
00922
00923
00924 void add(const vnl_bignum& b1, const vnl_bignum& b2, vnl_bignum& sum)
00925 {
00926 const vnl_bignum *bmax, *bmin;
00927 if (b1.count >= b2.count) {
00928 bmax = &b1;
00929 bmin = &b2;
00930 }
00931 else {
00932 bmax = &b2;
00933 bmin = &b1;
00934 }
00935 sum.resize(bmax->count);
00936 unsigned long temp, carry = 0;
00937 Counter i = 0;
00938 while (i < bmin->count) {
00939
00940 temp = (unsigned long)b1.data[i] + (unsigned long)b2.data[i] + carry;
00941 carry = temp/0x10000L;
00942 sum.data[i] = Data(temp);
00943 i++;
00944 }
00945 while (i < bmax->count) {
00946 temp = bmax->data[i] + carry;
00947 carry = temp/0x10000L;
00948 sum.data[i] = Data(temp);
00949 i++;
00950 }
00951 if (carry) {
00952 sum.resize(bmax->count + 1);
00953 sum.data[bmax->count] = 1;
00954 }
00955 }
00956
00957
00958
00959 void increment(vnl_bignum& bnum)
00960 {
00961 Counter i = 0;
00962 unsigned long carry = 1;
00963 while (i < bnum.count && carry) {
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
00977
00978 void subtract(const vnl_bignum& bmax, const vnl_bignum& bmin, vnl_bignum& diff)
00979 {
00980 diff.resize(bmax.count);
00981 unsigned long temp;
00982 int borrow = 0;
00983 Counter i = 0;
00984 for (; i < bmin.count; i++) {
00985 temp = (unsigned long)bmax.data[i] + 0x10000L - borrow;
00986 temp -= (unsigned long)bmin.data[i];
00987 borrow = (temp/0x10000L == 0);
00988 diff.data[i] = (Data) temp;
00989 }
00990 for (; i < bmax.count; i++) {
00991 temp = (unsigned long)bmax.data[i] + 0x10000L - borrow;
00992 borrow = (temp/0x10000L == 0);
00993 diff.data[i] = (Data) temp;
00994 }
00995 diff.trim();
00996 }
00997
00998
00999
01000 void decrement(vnl_bignum& bnum)
01001 {
01002 Counter i = 0;
01003 unsigned long borrow = 1;
01004 while (i < bnum.count && borrow) {
01005 unsigned long temp = (unsigned long)bnum.data[i] + 0x10000L - borrow;
01006 borrow = (temp/0x10000L == 0);
01007 bnum.data[i] = (Data)temp;
01008 ++i;
01009 }
01010 bnum.trim();
01011 if (bnum.count==0) bnum.sign=+1;
01012 }
01013
01014
01015
01016
01017
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;
01024 if (b2.count > b1.count) return -1;
01025 Counter i = b1.count;
01026 while (i > 0) {
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 }
01033 return 0;
01034 }
01035
01036
01037
01038
01039
01040 void multiply_aux(const vnl_bignum& b, Data d, vnl_bignum& prod, Counter i)
01041 {
01042
01043
01044
01045
01046
01047 if (i == 0) {
01048 Counter j = 0;
01049 while (j < prod.count)
01050 prod.data[j++] = 0;
01051 }
01052 if (d != 0) {
01053 unsigned long temp;
01054 Data carry = 0;
01055
01056 Counter j = 0;
01057 for (; j < b.count; j++) {
01058
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);
01062 carry = Data(temp/0x10000L);
01063 }
01064 if (i+j < prod.count)
01065 prod.data[i + j] = carry;
01066 }
01067 }
01068
01069
01070
01071
01072
01073
01074
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));
01079 u.resize(b1.count + 1);
01080 v.resize(b2.count);
01081 u.data[b1.count] = 0;
01082 multiply_aux(b1,d,u,0);
01083 multiply_aux(b2,d,v,0);
01084 return d;
01085 }
01086
01087
01088
01089
01090
01091
01092
01093 void divide_aux(const vnl_bignum& b1, Data d, vnl_bignum& q, Data& r)
01094 {
01095 r = 0;
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];
01099 if (j < 1 + q.count)
01100 q.data[j - 1] = Data(temp/d);
01101 r = Data(temp % d);
01102 }
01103 }
01104
01105
01106
01107
01108
01109
01110
01111
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],
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
01123 q_hat = (u0 == v1 ? Data(0xffff) : Data(((unsigned long)u0 * 0x10000L + (unsigned long)u1) / (unsigned long)v1));
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133 unsigned long lhs, rhs;
01134 for (Counter i = 0; i < 2; i++) {
01135 lhs = (unsigned long)v2 * (unsigned long)q_hat;
01136 rhs = (unsigned long)u0 * 0x10000L +(unsigned long)u1;
01137 rhs -= ((unsigned long)q_hat * (unsigned long)v1);
01138
01139 if ( rhs >= 0x10000L )
01140 break;
01141 rhs *= 0x10000L;
01142
01143 if (rhs > rhs + (unsigned long)u2)
01144 break;
01145 rhs += u2;
01146 if (lhs <= rhs)
01147 break;
01148 q_hat--;
01149 }
01150 return q_hat;
01151 }
01152
01153
01154
01155
01156
01157
01158
01159
01160 Data multiply_subtract(vnl_bignum& u, const vnl_bignum& v, Data q_hat, Counter j)
01161 {
01162
01163
01164
01165
01166 if (q_hat == 0) return q_hat;
01167 vnl_bignum rslt;
01168 Counter tmpcnt;
01169 rslt.resize(v.count + 1u);
01170
01171
01172 unsigned long prod, diff;
01173 Data carry = 0, borrow = 0;
01174 Counter i = 0;
01175 for (; i < v.count; ++i) {
01176
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);
01180 rslt.data[i] = Data(diff);
01181 borrow = (diff/0x10000L == 0) ? 1 : 0;
01182 carry = Data(prod/0x10000L);
01183 }
01184 tmpcnt = Counter(u.count - v.count + i - j - 1);
01185 diff = (unsigned long)u.data[tmpcnt]
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
01192
01193
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 {
01206 for (i = 0; i < rslt.count; ++i)
01207 u.data[u.count - v.count + i - 1 - j] = rslt.data[i];
01208 }
01209 return q_hat;
01210 }
01211
01212
01213
01214
01215
01216
01217
01218 void divide(const vnl_bignum& b1, const vnl_bignum& b2, vnl_bignum& q, vnl_bignum& r)
01219 {
01220
01221 assert(&b1 != &q && &b2 != &q && &b1 != &r && &b2 != &r);
01222 q = r = 0L;
01223 if (b1 == 0L)
01224 return;
01225 int mag = magnitude_cmp(b1,b2);
01226 if (mag < 0)
01227 r = b1;
01228 else if (mag == 0)
01229 q = 1L;
01230 else {
01231 q.resize(b1.count + 1u - b2.count);
01232 r.resize(b2.count);
01233 if (b2.count == 1) {
01234 divide_aux(b1,b2.data[0],q,r.data[0]);
01235 }
01236 else {
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);
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) {
01250 Data q_hat = estimate_q_hat(u,v,j);
01251 q.data[q.count - 1 - j] =
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;
01260 divide_aux(u,d,r,dufus);
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();
01268 r.trim();
01269 }
01270 q.sign = r.sign = b1.sign * b2.sign;
01271 }
01272
01273
01274
01275
01276
01277 vnl_bignum left_shift(const vnl_bignum& b1, int l)
01278 {
01279
01280
01281
01282
01283 vnl_bignum rslt;
01284 rslt.sign = b1.sign;
01285 Counter growth = Counter(l / 16);
01286 Data shift = Data(l % 16);
01287 Data rshift = Data(16 - shift);
01288 Data carry = Data(
01289 b1.data[b1.count - 1] >> (16 - shift));
01290 rslt.resize(b1.count + growth + (carry ? 1u : 0u));
01291 Counter i = 0;
01292 while (i < growth)
01293 rslt.data[i++] = 0;
01294 rslt.data[i++] = b1.data[0] << shift;
01295 while (i < rslt.count - 1) {
01296 rslt.data[i] = (b1.data[i - growth] << shift) +
01297 (b1.data[i - 1 - growth] >> rshift);
01298 i++;
01299 }
01300 if (i < rslt.count) {
01301 if (carry)
01302 rslt.data[i] = carry;
01303 else
01304 rslt.data[i] = (b1.data[i - growth] << shift)
01305 + (b1.data[i - 1 - growth] >> rshift);
01306 }
01307 vnl_bignum& result = *((vnl_bignum*) &rslt);
01308 return result;
01309 }
01310
01311
01312
01313
01314
01315 vnl_bignum right_shift(const vnl_bignum& b1, int l)
01316 {
01317 vnl_bignum rslt;
01318 Counter shrinkage = Counter(l / 16);
01319 Data shift = Data(l % 16);
01320 Data dregs = Data(b1.data[b1.count-1] >> shift);
01321 if (shrinkage + (dregs == 0) < b1.count) {
01322 rslt.sign = b1.sign;
01323 rslt.resize(b1.count - shrinkage - (dregs == 0 ? 1 : 0));
01324 Data lshift = Data(16 - shift);
01325 Counter i = 0;
01326 while (i < rslt.count - 1) {
01327 rslt.data[i] = (b1.data[i + shrinkage] >> shift) +
01328 (b1.data[i + shrinkage + 1u] << lshift);
01329 i++;
01330 }
01331 if (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);
01339 return result;
01340 }