big_int.cpp

Go to the documentation of this file.
00001 /*****************************************************************************
00002     
00003     big_int.cpp -- see the big_int.hpp file for description.
00004 
00005     This file is a part of the Arageli library.
00006     
00007     WARNIG. This file has no complate implementation.
00008 
00009     Copyright (C) Nikolai Yu. Zolotykh, 1999 -- 2005
00010     University of Nizhni Novgorod, Russia
00011 
00012 *****************************************************************************/
00013 
00014 #include "config.hpp"
00015 
00016 #if !defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE) || \
00017     defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_BIG_INT)
00018 
00019 
00020 #include <cstddef>
00021 #include <limits>
00022 #include <cmath>
00023 
00024 #include "big_int.hpp"
00025 #include "_utility.hpp"
00026 
00027 
00028 namespace Arageli
00029 {
00030 
00031 
00032 namespace _Internal
00033 {
00034 
00035 template <typename T> struct Unsigned { typedef T Type; };
00036 
00037 #define _ARAGELI_IMPL_UNSIGNED_SIGNED(TYPE) \
00038     template <> struct Unsigned<TYPE> { typedef unsigned TYPE Type; }
00039 
00040 _ARAGELI_IMPL_UNSIGNED_SIGNED(char);
00041 _ARAGELI_IMPL_UNSIGNED_SIGNED(short);
00042 _ARAGELI_IMPL_UNSIGNED_SIGNED(int);
00043 _ARAGELI_IMPL_UNSIGNED_SIGNED(long);
00044 
00045 #ifdef ARAGELI_INT64_SUPPORT
00046     _ARAGELI_IMPL_UNSIGNED_SIGNED(__int64);
00047 #endif
00048 
00049 template <> struct Unsigned<signed char> { typedef unsigned char Type; };
00050 
00051 }
00052 
00053 
00054 #ifdef ARAGELI_DISABLE_PARTICULAR_COMPILER_WARNINGS
00055     #pragma warning (push)
00056     #pragma warning (disable : 4146)
00057     #pragma warning (disable : 4244)
00058     #pragma warning (disable : 4804)
00059     #pragma warning (disable : 4800)
00060     #pragma warning (disable : 4806)
00061     #pragma warning (disable : 4305)
00062 #endif
00063 
00064 
00065 template <typename T>
00066 void big_int::from_native_int (const T& x)
00067 {
00068     typedef std::numeric_limits<T> Nl;
00069     ARAGELI_ASSERT_1(Nl::is_specialized);
00070     ARAGELI_ASSERT_1(Nl::is_integer);
00071     alloc_zero();
00072 
00073     if(Arageli::is_null(x))return;
00074     else if(is_negative(x))
00075     {
00076         // WARNING. The following expression is not portable.
00077         from_native_int(static_cast<typename _Internal::Unsigned<T>::Type>(-x));
00078         number->sign = -1;
00079     }
00080     else if(Nl::digits <= _Internal::bits_per_digit)
00081     {
00082         _Internal::digit* mem = big_int::get_mem_for_data(1);
00083         
00084         try
00085         {
00086             mem[0] = _Internal::digit(x);
00087             free_mem_and_alloc_number(1, mem, 1);
00088         }
00089         catch(...)
00090         {
00091             free_data(mem);
00092             throw;
00093         }
00094     }
00095     else
00096     {
00097         std::size_t n =
00098             Nl::digits / _Internal::bits_per_digit
00099             + bool(Nl::digits % _Internal::bits_per_digit);
00100 
00101         _Internal::digit* mem = get_mem_for_data(n);
00102 
00103         try
00104         {
00105             T xx = x;
00106             
00107             for(std::size_t i = 0; i < n; ++i, xx >>= _Internal::bits_per_digit)
00108                 mem[i] = xx & _Internal::max_digit;
00109 
00110             ARAGELI_ASSERT_1(Arageli::is_null(xx));
00111 
00112             std::size_t newlen;
00113             mem = optimize(newlen, mem, n);
00114             ARAGELI_ASSERT_1(newlen);
00115             free_mem_and_alloc_number(1, mem, newlen);
00116         }
00117         catch(...)
00118         {
00119             free_data(mem);
00120             throw;
00121         }
00122 
00123         number->sign = +1;
00124     }
00125 }
00126 
00127 
00128 template <typename T>
00129 void big_int::from_native_float (const T& x)
00130 {
00131     typedef std::numeric_limits<T> Nl;
00132     
00133     ARAGELI_ASSERT_1(Nl::is_specialized);
00134     ARAGELI_ASSERT_1(!Nl::is_integer);
00135     ARAGELI_ASSERT_0(!(Nl::has_infinity && x == Nl::infinity()));
00136     ARAGELI_ASSERT_0(!(Nl::has_quiet_NaN && x == Nl::quiet_NaN()));
00137     ARAGELI_ASSERT_0(!(Nl::has_signaling_NaN && x == Nl::signaling_NaN()));
00138     ARAGELI_ASSERT_0(!(Nl::has_denorm && x == Nl::denorm_min()));
00139 
00140     alloc_zero();
00141 
00142     if(x > opposite_unit(x) && x < unit(x))return;
00143     else if(is_negative(x))
00144     {
00145         from_native_float(-x);
00146         number->sign = -1;
00147     }
00148     else
00149     {
00150         T xx = std::floor(x);
00151         int expon;
00152         xx = std::frexp(xx, &expon);
00153         
00154         ARAGELI_ASSERT_1(expon > 0);
00155         
00156         xx *= std::pow(2 * unit(x), Nl::digits + 1);
00157         expon -= Nl::digits + 1;
00158 
00159         ARAGELI_ASSERT_1(std::floor(xx) == xx);
00160 
00161         std::size_t n =
00162             (Nl::digits + 1) / _Internal::bits_per_digit
00163             + bool((Nl::digits + 1) % _Internal::bits_per_digit);
00164 
00165         _Internal::digit* mem = get_mem_for_data(n);
00166 
00167         try
00168         {
00169             T module = _Internal::max_digit;
00170             module += unit(module);
00171 
00172             for(std::size_t i = 0; i < n; ++i)
00173             {
00174                 T remaind = std::fmod(xx, module);
00175 
00176                 ARAGELI_ASSERT_1(std::floor(remaind) == remaind);
00177                 ARAGELI_ASSERT_1(_Internal::digit(remaind) == remaind);
00178 
00179                 mem[i] = _Internal::digit(remaind);
00180                 xx = (xx - remaind) / module;
00181 
00182                 ARAGELI_ASSERT_1(std::floor(xx) == xx);
00183             }
00184 
00185             ARAGELI_ASSERT_1(Arageli::is_null(xx));
00186 
00187             std::size_t newlen;
00188             mem = big_int::optimize(newlen, mem, n);
00189             ARAGELI_ASSERT_1(newlen);
00190             free_mem_and_alloc_number(1, mem, newlen);
00191         }
00192         catch(...)
00193         {
00194             free_data(mem);
00195             throw;
00196         }
00197 
00198         number->sign = +1;
00199         
00200         if(expon > 0)*this <<= expon;
00201         else *this >>= -expon;
00202     }
00203 }
00204 
00205 
00206 template <typename T>
00207 T big_int::to_native_int () const
00208 {
00209     typedef std::numeric_limits<T> Nl;
00210     ARAGELI_ASSERT_1(Nl::is_specialized);
00211     ARAGELI_ASSERT_1(Nl::is_integer);
00212 
00213     ARAGELI_ASSERT_0(big_int(Nl::min()) <= *this && *this <= big_int(Nl::max()));
00214 
00215     if(is_null())return factory<T>::null();
00216     else if(sign() < 0)
00217     {
00218         // WARNING. The following expression is not portable.
00219         return
00220             -static_cast<T>
00221             (to_native_int_without_sign<typename _Internal::Unsigned<T>::Type>());
00222     }
00223     else return to_native_int_without_sign<T>();
00224 }
00225 
00226 
00227 template <typename T>
00228 T big_int::to_native_int_without_sign () const
00229 {
00230     typedef std::numeric_limits<T> Nl;
00231     ARAGELI_ASSERT_1(Nl::is_specialized);
00232     ARAGELI_ASSERT_1(Nl::is_integer);
00233     ARAGELI_ASSERT_1(*this <= big_int(Nl::max()));
00234     ARAGELI_ASSERT_1(!is_null());
00235 
00236     T res = number->data[0];
00237 
00238     for(std::size_t i = 1; i < number->len; ++i)
00239         res |= number->data[i] << (i * _Internal::bits_per_digit);
00240 
00241     return res;
00242 }
00243 
00244 
00245 template <typename T>
00246 T big_int::to_native_float () const
00247 {
00248     typedef std::numeric_limits<T> Nl;
00249     
00250     ARAGELI_ASSERT_1(Nl::is_specialized);
00251     ARAGELI_ASSERT_1(!Nl::is_integer);
00252     ARAGELI_ASSERT_0(Nl::has_infinity);
00253     
00254     if(is_null())return factory<T>::null();
00255     if(*this < big_int(-Nl::max()))return -Nl::infinity();
00256     else if(*this > big_int(Nl::max()))return Nl::infinity();
00257     else
00258     {
00259         big_int t = *this;
00260         int expon = 0;
00261         std::size_t blen = t.length();
00262 
00263         if(blen - 1 > Nl::digits)
00264         {
00265             expon = blen - 1 - Nl::digits;
00266             t >>= expon;
00267         }
00268 
00269         ARAGELI_ASSERT_1(t.length() - 1 <= Nl::digits);
00270 
00271         T res = factory<T>::null();
00272 
00273         T module = _Internal::max_digit;
00274         module += unit(module);
00275         T curscale = unit(module);
00276 
00277         for(std::size_t i = 0; i < t.number->len; ++i)
00278         {
00279             res += t.number->data[i]*curscale;
00280             curscale *= module;
00281         }
00282 
00283         res *= std::pow(2*unit(res), expon);
00284         if(sign() < 0)opposite(&res);
00285 
00286         ARAGELI_DEBUG_EXEC_1(big_int backres = res);
00287         ARAGELI_ASSERT_1(backres.length() == length());
00288         
00289         ARAGELI_ASSERT_1
00290         (
00291             (length() <= Nl::digits && backres - *this == 0)  ||
00292             (
00293                 length() > Nl::digits &&
00294                 (backres - *this).length() <= length() - Nl::digits
00295             )
00296         );
00297 
00298         return res;
00299     }
00300 }
00301 
00302 
00303 #ifdef ARAGELI_DISABLE_PARTICULAR_COMPILER_WARNINGS
00304     #pragma warning (pop)
00305 #endif
00306 
00307 
00308 }
00309 
00310 
00311 #else
00312 
00313 
00314 #include <cstdlib>
00315 #include <malloc.h>
00316 #include <sstream>
00317 #include <limits>
00318 #include <cctype>
00319 
00320 #include "big_int.hpp"
00321 #include "rational.hpp"
00322 
00323 
00324 namespace
00325 {
00326     typedef Arageli::_Internal::digit digit;
00327 }
00328 
00329 
00330 namespace Arageli
00331 {
00332 
00333 void big_arith_error(const char *s)
00334 { throw big_int::exception(std::string("Big arith error: ") + s); }
00335 
00336 
00337 /***************************/
00338 /*                         */
00339 /*    low level memory     */
00340 /*   managment routines    */
00341 /*                         */
00342 /***************************/
00343 
00344 
00345 digit* big_int::get_mem_for_data (std::size_t nitems)
00346 {
00347     digit* p = reinterpret_cast<digit*>(std::malloc(nitems * sizeof(digit)));
00348     if(!p)big_arith_error("the heap overflow");
00349     return p;
00350 }
00351 
00352 
00353 void big_int::free_data (digit *p) { free(p); }
00354 
00355 
00356 digit* big_int::realloc_data (digit* p, std::size_t newnitems)
00357 {
00358     if(newnitems)
00359     {
00360         p = reinterpret_cast<digit*>(realloc(p, newnitems * sizeof(digit)));
00361         if(!p)big_arith_error("the heap overflow");
00362     }
00363     else
00364     {
00365         free_data(p);
00366         p = 0;
00367     }
00368 
00369     return p;
00370 }
00371 
00372 
00373 void big_int::copy_data
00374 (digit* dest, const digit* source, std::size_t newnitems)
00375 { memmove(dest, source, newnitems * sizeof(digit)); }
00376 
00377 /***************************/
00378 /*                         */
00379 /*    number allocation    */
00380 /*         routines        */
00381 /*                         */
00382 /***************************/
00383 
00384 void big_int::alloc_number
00385 (int new_sign, digit * new_data, std::size_t new_len)
00386 {
00387     number = new big_struct;
00388 
00389     number->sign = new_sign;
00390     number->data = new_data;
00391     number->len = new_len;
00392     number->refs = 1;
00393 }
00394 
00395 void big_int::free_number()
00396 {
00397     if(!number)return;
00398     number->refs--;
00399     if(number->refs)return;
00400     free_data(number->data);
00401     delete number;
00402     number = 0;
00403 }
00404 
00405 void big_int::free_mem_and_alloc_number
00406 (int new_sign, digit* new_data, std::size_t new_len)
00407 {
00408     free_number();
00409     alloc_number(new_sign, new_data, new_len);
00410 }
00411 
00412 
00413 void big_int::free_mem_and_alloc_zero ()
00414 {
00415     free_number();
00416     alloc_zero();
00417 }
00418 
00419 /***************************/
00420 /*                         */
00421 /*      calc_bdn_radix     */
00422 /*                         */
00423 /***************************/
00424 
00425 void calc_bdn_radix (digit radix, digit& bdn_radix, std::size_t& chars_per_digit)
00426 {
00427     //bdn_radix = maximal power of radix fit in digit;
00428 
00429     bdn_radix = 1;
00430     chars_per_digit = 0;
00431     digit t = _Internal::max_digit/radix;
00432     
00433     while(t)
00434     {
00435         t /= radix;
00436         bdn_radix *= radix;
00437         chars_per_digit++;
00438     }
00439 }
00440 
00441 
00442 /***************************/
00443 /*                         */
00444 /*         optimize        */
00445 /*                         */
00446 /***************************/
00447 
00448 digit* big_int::optimize (std::size_t& new_len, digit * p, std::size_t len)
00449 {
00450     // deleting leading zeros
00451   
00452     new_len = _Internal::do_optimize(p, len);
00453 
00454     if(new_len != len)
00455         if(new_len)
00456         {
00457             digit* new_p = get_mem_for_data(new_len);
00458             copy_data(new_p, p, new_len);
00459             free_data(p);
00460             p = new_p;
00461         }
00462         else
00463         {
00464             free_data(p);
00465             p = 0;
00466         }
00467 
00468     return p;
00469 }
00470 
00471 /***************************/
00472 /*                         */
00473 /*      constructors       */
00474 /*    and destructors      */
00475 /*                         */
00476 /***************************/
00477 
00478 big_int::big_int (const char* str)
00479 {
00480     std::istringstream s(str);
00481     big_int b;
00482     s >> b;
00483     if(!s && !s.eof())
00484         throw incorrect_string(str);
00485     alloc_zero();
00486     
00487     try { *this = b; }
00488     catch(...)
00489     { free_number(); throw; }
00490 }
00491 
00492 
00493 /***************************/
00494 /*                         */
00495 /*        operator =       */
00496 /*                         */
00497 /***************************/
00498 
00499 
00500 big_int& big_int::operator= (const big_int & b)
00501 {
00502     // make a copy of a number, just increments the reference count
00503     if(number == b.number)return *this;
00504     free_number();
00505     b.number->refs++;
00506     number = b.number;
00507     return *this;
00508 }
00509 
00510 
00511 //big_int::big_int (const rational<big_int>& x)
00512 //{
00513 //  alloc_zero();
00514 //  *this = x.numerator()/x.denominator();
00515 //}
00516 
00517 /***************************/
00518 /*                         */
00519 /*      Unary operators    */
00520 /*                         */
00521 /***************************/
00522 
00523 
00524 big_int big_int::operator- () const // unary minus
00525 {
00526     big_int a;
00527     digit *result;
00528     std::size_t blen;
00529 
00530     blen = number->len;
00531 
00532     if(number->sign == 0)
00533         a.free_mem_and_alloc_zero();
00534     else
00535     {
00536         result = get_mem_for_data(blen);
00537         copy_data(result, number->data, blen);      // copy data
00538         a.free_mem_and_alloc_number(-number->sign, result, blen);
00539     }
00540 
00541     return a;
00542 }
00543 
00544 /***************************/
00545 /*                         */
00546 /*        operator +       */
00547 /*                         */
00548 /***************************/
00549 
00550 big_int operator+ (const big_int& b, const big_int& c)
00551 // Thanks to Max Alekseyev for a lot of suggestions
00552 {
00553     digit *sum;
00554     std::size_t sumlen;
00555     big_int a;
00556     int bsign = b.number->sign;
00557     int csign = c.number->sign;
00558     int sumsign;
00559     std::size_t blen = b.number->len;
00560     std::size_t clen = c.number->len;
00561     digit *bdata = b.number->data;
00562     digit *cdata = c.number->data;
00563     big_int::big_struct *u, *v;
00564 
00565     if(bsign == 0)a = c;
00566     else if(csign == 0)a = b;
00567     else
00568     {
00569         if(bsign == csign)
00570         {
00571             sumsign = bsign;
00572             if (blen >= clen)
00573             {
00574                 u = b.number;
00575                 v = c.number;
00576             }
00577             else
00578             {
00579                 u = c.number;
00580                 v = b.number;
00581             }
00582             sum = big_int::get_mem_for_data(u->len + 1);
00583             big_int::copy_data(sum, u->data, u->len);
00584             sumlen = _Internal::do_add(sum, v->data, u->len, v->len);
00585         }
00586         else // bsign != csign
00587         {
00588             if(blen >= clen)
00589             {
00590                 sum = big_int::get_mem_for_data(blen);
00591                 sumsign = bsign;
00592                 big_int::copy_data(sum, bdata, blen);
00593                 if(_Internal::do_sub(sum, cdata, blen, clen))     // был заем
00594                     goto c_ge_b;                // yes, it is
00595                 sum = big_int::optimize(sumlen, sum, blen);
00596             }
00597             else
00598             {
00599                 sum = big_int::get_mem_for_data(clen);
00600             c_ge_b:
00601                 sumsign = csign;
00602                 big_int::copy_data(sum, cdata, clen);
00603                 _Internal::do_sub(sum, bdata, clen, blen);
00604                 sum = big_int::optimize(sumlen, sum, clen);
00605             }
00606         }
00607 
00608         if(sumlen)
00609             a.free_mem_and_alloc_number(sumsign, sum, sumlen);
00610         else
00611             a.free_mem_and_alloc_zero();
00612     }
00613 
00614     return a;
00615 }
00616 
00617 
00618 /***************************/
00619 /*                         */
00620 /*        operator -       */
00621 /*                         */
00622 /***************************/
00623 
00624 
00625 big_int operator- (const big_int& b, const big_int& c)
00626 { return b + (-c); }
00627 
00628 
00629 /***************************/
00630 /*                         */
00631 /*        operator *       */
00632 /*                         */
00633 /***************************/
00634 
00635 
00636 big_int operator* (const big_int& b, const big_int& c)
00637 {
00638     digit *result;
00639     std::size_t blen, clen, resultlen;
00640     int bsign, csign;
00641     big_int a;
00642 
00643     blen = b.number->len;
00644     clen = c.number->len;
00645     bsign = b.number->sign;
00646     csign = c.number->sign;
00647 
00648     if((bsign == 0) || (csign == 0))
00649         a.free_mem_and_alloc_zero();
00650     else
00651     {
00652         resultlen = blen + clen;
00653         result = big_int::get_mem_for_data(resultlen);
00654         if(_Internal::do_mult(b.number->data, c.number->data, result, blen, clen) == 0)
00655             resultlen--;
00656         a.free_mem_and_alloc_number(bsign * csign, result, resultlen);
00657     }
00658 
00659     return a;
00660 }
00661 
00662 
00663 /***************************/
00664 /*                         */
00665 /*         xdivide          */
00666 /*                         */
00667 /***************************/
00668 
00669 
00670 void _Internal::xdivide (big_int& a, const big_int& b, const big_int& c, big_int& res)
00671 {
00672     // a = b / c; res = b % c
00673     
00674     digit *u, *v, *q, *r;
00675     std::size_t alen, blen, clen, rlen;
00676     digit runint;
00677     int bsign, csign;
00678 
00679     blen = b.number->len;
00680     clen = c.number->len;
00681     bsign = b.number->sign;
00682     csign = c.number->sign;
00683 
00684     if(!csign)big_arith_error("divide by zero");
00685     else if(!bsign)
00686     {
00687         a.free_mem_and_alloc_zero();
00688         res.free_mem_and_alloc_zero();
00689     }
00690     else if (blen < clen)
00691     {
00692         a.free_mem_and_alloc_zero();
00693         res = b;
00694     }
00695     else if (clen == 1)
00696     {
00697         q = big_int::get_mem_for_data(blen);
00698         runint = _Internal::do_divide_by_digit
00699             (b.number->data, q, blen, c.number->data[0]);
00700         
00701         if(q[blen - 1])
00702             alen = blen;
00703         else
00704             alen = blen - 1;
00705         
00706         if(alen)
00707             a.free_mem_and_alloc_number(bsign * csign, q, alen);
00708         else
00709         {
00710             a.free_mem_and_alloc_zero();
00711             big_int::free_data(q);
00712         }
00713 
00714         if(runint == 0)
00715             res.free_mem_and_alloc_zero();
00716         else
00717         {
00718             u = big_int::get_mem_for_data(1);
00719             u[0] = runint;
00720             res.free_mem_and_alloc_number(bsign, u, 1);
00721         }
00722 
00723     }
00724     else
00725     {
00726         u = big_int::get_mem_for_data(blen + 1);
00727         v = big_int::get_mem_for_data(clen);
00728         q = big_int::get_mem_for_data(blen - clen + 1);
00729         big_int::copy_data(u, b.number->data, blen);
00730         big_int::copy_data(v, c.number->data, clen);
00731 
00732         rlen = _Internal::do_divide(u, v, q, blen, clen);
00733         if(q[blen - clen])
00734             alen = blen - clen + 1;
00735         else
00736             alen = blen - clen;
00737 
00738         big_int::free_data(v);
00739         if(rlen == 0)
00740             res.free_mem_and_alloc_zero();
00741         else
00742         {
00743             r = big_int::get_mem_for_data(rlen);
00744             big_int::copy_data(r, u, rlen);
00745             res.free_mem_and_alloc_number(bsign, r, rlen);
00746         }
00747 
00748         big_int::free_data(u);
00749         if(alen != 0)
00750             a.free_mem_and_alloc_number(bsign * csign, q, alen);
00751         else
00752         {
00753             a.free_mem_and_alloc_zero();
00754             big_int::free_data(q);
00755         }
00756     }
00757 
00758     ARAGELI_ASSERT_1(b == c*a + res);
00759 }
00760 
00761 
00762 big_int operator% (const big_int& b, const big_int& c)
00763 {
00764     big_int a, r;
00765     _Internal::xdivide(a, b, c, r);
00766     return r;
00767 }
00768 
00769 
00770 big_int operator/ (const big_int& b, const big_int& c)
00771 {
00772     big_int a, r;
00773     _Internal::xdivide(a, b, c, r);
00774     return a;
00775 }
00776 
00777 
00778 /***************************/
00779 /*                         */
00780 /*    Compare functions    */
00781 /*                         */
00782 /***************************/
00783 
00784 
00785 int cmp (const big_int & a, const big_int & b)
00786 {
00787     //sign(a-b)
00788     
00789     int result;
00790 
00791     int asign, bsign;
00792     std::size_t alen, blen;
00793     big_int c;
00794 
00795     asign = a.number->sign;
00796     bsign = b.number->sign;
00797     alen = a.number->len;
00798     blen = b.number->len;
00799 
00800     if(asign < bsign)result = -1;
00801     else if(asign > bsign)result = 1;
00802     else if(asign == 0)result = 0;          // asign == bsign == 0
00803     else if(alen < blen)result = -asign;    // asign == bsign != 0
00804     else if(alen > blen)result = asign;
00805     else
00806     { // asign == bsign != 0, alen == blen
00807         c = a - b;
00808         result = c.number->sign;
00809     }
00810 
00811     return result;
00812 }
00813 
00814 
00815 /***************************/
00816 /*                         */
00817 /*      Random numbers     */
00818 /*                         */
00819 /***************************/
00820 
00821 #ifdef ARAGELI_DISABLE_PARTICULAR_COMPILER_WARNINGS
00822     #pragma warning (push)
00823     #pragma warning (disable : 4800)
00824 #endif
00825 
00826 digit random_digit ()
00827 {
00828     if(_Internal::max_digit <= RAND_MAX)return rand();
00829 
00830     static const std::size_t rand_max_len = big_int(RAND_MAX).length(); // WARNING! It is not efficient.
00831     const std::size_t n =
00832         _Internal::bits_per_digit / rand_max_len
00833         + bool(_Internal::bits_per_digit % rand_max_len);
00834     
00835     digit res = 0;
00836     for(std::size_t i = 0; i < n; ++i)
00837         (res <<= rand_max_len) |= rand();   // WARNING! Lower bits from rand is
00838                                             // placed to higher bits of the result.
00839     return res;
00840 }
00841 
00842 #ifdef ARAGELI_DISABLE_PARTICULAR_COMPILER_WARNINGS
00843     #pragma warning (pop)
00844 #endif
00845 
00846 
00847 //big_int random_number(std::size_t length)
00848 //{
00849 //  big_int a;
00850 //  digit *p = big_int::get_mem_for_data(length);
00851 //  
00852 //  try
00853 //  {
00854 //      for(std::size_t i = 0; i < length; i++)
00855 //          p[i] = random_digit();
00856 //      
00857 //      std::size_t new_len;
00858 //      p = big_int::optimize(new_len, p, length);
00859 //      
00860 //      if(new_len)
00861 //          a.free_mem_and_alloc_number(1, p, new_len);
00862 //  }
00863 //  catch(...)
00864 //  { big_int::free_data(p); throw; }
00865 //  
00866 //  return a;
00867 //}
00868 
00869 
00870 big_int big_int::random_in_range (const big_int& max)
00871 {
00872     if(max.is_null())return max;
00873     
00874     std::size_t len = max.length();
00875 
00876     for(;;)
00877     {
00878         big_int res = random_with_length_or_less(len);
00879         if(res <= max)return res;
00880     }
00881 }
00882 
00883 
00884 big_int big_int::random_with_length_or_less (std::size_t length) //bits
00885 {
00886     int bits_in_highest = length % _Internal::bits_per_digit;
00887     std::size_t len = length/_Internal::bits_per_digit + 1;
00888     digit *p = big_int::get_mem_for_data(len);
00889 
00890     try
00891     {
00892         for(std::size_t i = 0; i < len - 1; i++)
00893             p[i] = random_digit();
00894 
00895         if(bits_in_highest)
00896             p[len - 1] = random_digit() >>
00897                 (_Internal::bits_per_digit - bits_in_highest);
00898         else
00899             p[len - 1] = 0;
00900 
00901         std::size_t new_len;
00902 
00903         p = big_int::optimize(new_len, p, len);
00904         big_int a;
00905 
00906         if(new_len)
00907             a.free_mem_and_alloc_number(1, p, new_len);
00908 
00909         return a;
00910     }
00911     catch(...)
00912     { big_int::free_data(p); throw; }
00913 }
00914 
00915 /***************************/
00916 /*                         */
00917 /*     Stream operators    */
00918 /*                         */
00919 /***************************/
00920 
00921 digit stream_radix(std::ios & s)
00922 {
00923     if(s.flags() & std::ios::dec)return 10;
00924     if(s.flags() & std::ios::hex)return 16;
00925     if(s.flags() & std::ios::oct)return 8;
00926     else return 10; // it is very strange
00927 }
00928 
00929 
00930 std::ostream& operator<< (std::ostream& s, const big_int& x)
00931 {
00932     if(s.flags() & std::ios_base::showpos && x.sign() > 0)
00933         s << '+';
00934     
00935     if(!x.number->sign)return s << "0";
00936 
00937     digit bdn_radix;
00938     std::size_t chars_per_block;
00939     calc_bdn_radix(stream_radix(s), bdn_radix, chars_per_block);
00940 
00941     std::size_t bdnlen = x.number->len;
00942     digit *bdn = big_int::get_mem_for_data(2 * bdnlen); //bdnlen + bdnlen/chars_per_block
00943     
00944     try
00945     {
00946         digit *numberdata = big_int::get_mem_for_data(x.number->len);
00947         
00948         try
00949         {
00950             big_int::copy_data(numberdata, x.number->data, x.number->len);
00951             bdnlen = _Internal::do_big_int_to_bdn(numberdata, bdn, x.number->len, bdn_radix);
00952         }
00953         catch(...)
00954         { big_int::free_data(numberdata); throw; }
00955         
00956         big_int::free_data(numberdata);
00957 
00958         _Internal::Auto_stream_state _ass(s, s.flags() & ~std::iostream::showpos);
00959 
00960         if(x.number->sign == -1)s << '-';
00961         s << bdn[bdnlen - 1];
00962 
00963         for (std::size_t i = bdnlen - 1; i > 0; i--)
00964         {
00965             s.width(chars_per_block);
00966             s.fill('0');
00967             s << bdn[i - 1];
00968         }
00969     }
00970     catch(...)
00971     { big_int::free_data(bdn); throw; }
00972 
00973     big_int::free_data(bdn);
00974     return s;
00975 }
00976 
00977 
00978 inline void set_stream_radix (std::ios& s, digit radix)
00979 {
00980     switch(radix)
00981     {
00982         case 10: s.setf(std::ios::dec, std::ios::basefield); break;
00983         case 16: s.setf(std::ios::hex, std::ios::basefield); break;
00984         case 8:  s.setf(std::ios::oct, std::ios::basefield); break;
00985         default: throw big_int::exception("Value for radix is invalid");
00986     }
00987 }
00988 
00989 
00990 std::istream& operator>> (std::istream& s, big_int& x)
00991 {
00992 #if 0   // old version
00993 
00994     // WARNING!!! This function is incorrect!
00995 
00996     char ch;
00997     int sign = 1;
00998     digit radix = 10;
00999     bool brackets = false;
01000 
01001     while(s.get(ch) && isspace(ch));    // pass leading spaces
01002 
01003     if(!s)big_arith_error("empty string for reading as big int");
01004 
01005     if(ch == '-')
01006     {
01007         sign = -1;
01008         s.get(ch);
01009     }
01010     else if(ch == '+')s.get(ch);
01011     else if(ch == '(')
01012     {
01013         s.get(ch);
01014         brackets = true;
01015     }
01016 
01017     if(!s)big_arith_error("invalid string format for reading as big int");
01018 
01019     // now in ch the first digit of number notation
01020 
01021     int zero = 0;
01022     if (ch == '0')
01023     {
01024         zero = 1;
01025         s.get(ch);
01026         switch (ch)
01027         {
01028             case 'o': radix = 8; break;
01029             case 'x': radix = 16; break;
01030             default: s.putback(ch);
01031         }
01032     }
01033     else s.putback(ch);
01034 
01035     do s.get(ch); while (ch == '0');    // pass non significant first zeros
01036     s.putback(ch);
01037 
01038     std::ostringstream buffer;
01039     std::size_t tellp = 0;
01040 
01041     while(s.get(ch))
01042     {
01043         if
01044         (
01045             radix == 10 && isdigit(ch)
01046             || radix == 16 && isxdigit(ch)
01047             || radix == 8  && isdigit(ch) && ch < '8'
01048         )
01049         {
01050             buffer << ch;
01051             ++tellp;
01052         }
01053         else
01054         {
01055             s.putback(ch);
01056             break;
01057         }
01058     }
01059 
01060     if(tellp)
01061     {
01062         if(!zero)
01063             big_arith_error("empty string for reading as big int");
01064         x = 0;
01065         return s;
01066     }
01067 
01068     std::string buffer_cpp_str = buffer.str();
01069     const char* buffer_str = buffer_cpp_str.c_str();
01070 
01071     digit bdn_radix;
01072     std::size_t chars_per_block;
01073     calc_bdn_radix(radix, bdn_radix, chars_per_block);
01074 
01075     std::size_t length_of_str = buffer_cpp_str.length();
01076     std::size_t number_of_blocks = (length_of_str - 1)/chars_per_block + 1;
01077     digit* bdn = big_int::get_mem_for_data(number_of_blocks);
01078     digit* data = 0;
01079     std::size_t len = 0;
01080 
01081     try
01082     {
01083         std::size_t length_of_block = length_of_str % chars_per_block;
01084                                 // for the first block
01085         if(!length_of_block)length_of_block = chars_per_block;
01086 
01087         const char * buffer_str_cur = buffer_str;
01088         for(std::size_t j = number_of_blocks; j > 0; j--)
01089         {
01090             std::istringstream stream_block(std::string(buffer_str_cur).substr(0, length_of_block));
01091             set_stream_radix(stream_block, radix);
01092             stream_block >> bdn[j-1];
01093             buffer_str_cur += length_of_block;
01094             length_of_block = chars_per_block;  // for the last blocks
01095         }
01096 
01097         data = big_int::get_mem_for_data(number_of_blocks);
01098         len = _Internal::do_bdn_to_big_int(bdn, data, number_of_blocks, bdn_radix);
01099     }
01100     catch(...)
01101     { big_int::free_data(bdn); big_int::free_data(data); throw; }
01102         
01103     big_int::free_data(bdn);
01104 
01105     try
01106     {
01107         data = big_int::realloc_data(data, len);
01108 
01109         if(data)x.free_mem_and_alloc_number(sign, data, len);
01110         else x.free_mem_and_alloc_zero();
01111     }
01112     catch(...)
01113     { big_int::free_data(data); throw; }
01114 
01115     if(brackets)
01116     {
01117         s >> ch;
01118         if(ch != ')')
01119             big_arith_error("invalid string format for reading as big int");
01120     }
01121 
01122     s.clear();
01123     return s;
01124 
01125 #else
01126 
01127     char ch = 0;
01128     int sign = 1;
01129     digit radix = 10;
01130     bool brackets = false;
01131 
01132     std::size_t len;
01133 
01134     //do s.get(ch); while (isspace(ch));
01135     s >> ch;    // pass leading spaces (if this feature is turn on for s)
01136                 // and read first non space character
01137     if(!s)
01138     {
01139         s.clear(std::ios_base::badbit);
01140         return s;
01141     }
01142 
01143     if(ch == '(')
01144     {
01145         brackets = true;
01146         s >> ch;
01147     }
01148 
01149     if(ch == '-')sign = -1;
01150     else if(ch != '+')s.putback(ch);
01151 
01152     s.get(ch);
01153 
01154     int zero = 0;
01155     if(ch == '0')
01156     {
01157         zero = 1;
01158         s.get(ch);
01159 
01160         switch (ch)
01161         {
01162             case 'o': radix = 8; break;
01163             case 'x': radix = 16; break;
01164             default: s.putback(ch);
01165         }
01166     }
01167     else s.putback(ch);
01168 
01169     s.clear();
01170     do s.get(ch); while (ch == '0' && s);
01171     
01172     if(s)s.putback(ch);
01173     else
01174     {
01175         x = big_int();
01176         return s;
01177     }
01178 
01179     std::string buffer;
01180 
01181     while(s.get(ch))
01182     {
01183         if
01184         (
01185             radix == 10 && isdigit(ch) ||
01186             radix == 16 && isxdigit(ch) ||
01187             radix == 8  && isdigit(ch) && ch < '8'
01188         )
01189             buffer += ch;
01190         else
01191         {
01192             s.putback(ch);
01193             break;
01194         }
01195     }
01196 
01197     if(buffer.empty())
01198     {
01199         if(!zero)throw big_int::incorrect_string(std::string(1, s.peek()));
01200         
01201         if(brackets && !_Internal::read_literal(s, ")"))
01202             throw big_int::incorrect_string(std::string(1, s.peek()));
01203         
01204         x = big_int();
01205         return s;
01206     }
01207 
01208     digit bdn_radix;
01209     std::size_t chars_per_block;
01210     calc_bdn_radix(radix, bdn_radix, chars_per_block);
01211 
01212     std::size_t length_of_str = buffer.length();
01213     std::size_t number_of_blocks = (length_of_str - 1)/chars_per_block + 1;
01214     _Internal::auto_array<digit> bdn(big_int::get_mem_for_data(number_of_blocks));
01215     std::size_t length_of_block = length_of_str % chars_per_block;
01216                             // for the first block
01217     if(!length_of_block)length_of_block = chars_per_block;
01218 
01219     std::size_t buffer_cur = 0;
01220     for(std::size_t j = number_of_blocks; j > 0; j--)
01221     {
01222         std::string block = buffer.substr(buffer_cur, length_of_block);
01223         std::istringstream stream_block(block); 
01224         set_stream_radix(stream_block, radix);
01225         stream_block >> bdn.get()[j-1];
01226         buffer_cur += length_of_block;
01227         length_of_block = chars_per_block;  // for the last blocks
01228     }
01229 
01230     _Internal::auto_array<digit> data(big_int::get_mem_for_data(number_of_blocks));
01231     
01232     len = _Internal::do_bdn_to_big_int
01233     (
01234         bdn.get(), data.get(),
01235         number_of_blocks, bdn_radix
01236     );
01237     
01238     big_int::free_data(bdn.release());
01239     digit* d = big_int::realloc_data(data.release(), len);
01240 
01241     if(d)x.free_mem_and_alloc_number(sign, d, len);
01242     else x.free_mem_and_alloc_zero();
01243 
01244     if(brackets && !_Internal::read_literal(s, ")"))
01245         throw big_int::incorrect_string(std::string(1, s.peek()));
01246 
01247     return s;
01248 
01249 #endif
01250 }
01251 
01252 
01253 std::size_t big_int::length () const
01254 {
01255     if(!number->len)return 0;
01256     std::size_t l = (number->len - 1) * _Internal::bits_per_digit;
01257     digit highest = number->data[number->len - 1];
01258     while(highest >>= 1)l++;
01259     return l + 1;
01260 }
01261 
01262 
01263 #ifdef ARAGELI_DISABLE_PARTICULAR_COMPILER_WARNINGS
01264     #pragma warning (push)
01265     #pragma warning (disable : 4800)
01266 #endif
01267 
01268 bool big_int::operator[] (std::size_t k) const
01269 {
01270     ARAGELI_ASSERT_0(k < length());
01271     return (number->data[k / _Internal::bits_per_digit] >>
01272         (k % _Internal::bits_per_digit)) % 2;
01273 }
01274 
01275 #ifdef ARAGELI_DISABLE_PARTICULAR_COMPILER_WARNINGS
01276     #pragma warning (pop)
01277 #endif
01278 
01279 
01280 
01281 big_int operator<< (const big_int& a, std::size_t n)
01282 {
01283     std::size_t a_len = a.number->len;
01284     if(!n || !a_len) return a;
01285     std::size_t l = a.length() + n;
01286     std::size_t res_len = l/_Internal::bits_per_digit;
01287     if (l%_Internal::bits_per_digit) 
01288         res_len++;
01289     std::size_t m = n/_Internal::bits_per_digit;
01290     n %= _Internal::bits_per_digit;
01291     digit t = 0;
01292     std::size_t k = _Internal::bits_per_digit - n;
01293     digit* res_data = big_int::get_mem_for_data(res_len);
01294     const digit* data = a.number->data;
01295 
01296     //  memset(res_data, 0, m * sizeof(digit));
01297     std::fill_n(res_data, m, digit(0));
01298 
01299     if(n)
01300     {
01301         for(std::size_t i = 0; i < a_len; i++) 
01302         { 
01303             res_data[m + i] = (data[i] << n) | t; 
01304             t = data[i] >> k;
01305         }
01306         if(t)res_data[m + a_len] = t;
01307     }
01308     else
01309     {
01310         std::copy(data, data + a_len, res_data + m);
01311     }
01312     
01313     big_int res;
01314     res.free_mem_and_alloc_number(a.number->sign, res_data, res_len);
01315     return res;
01316 }
01317 
01318 
01319 big_int operator>> (const big_int& a, std::size_t n)
01320 {
01321     std::size_t a_len = a.number->len;
01322     if (!n || !a_len) return a;
01323     std::size_t l = a.length(); 
01324     if (l <= n) return big_int();
01325     l -= n;
01326     std::size_t res_len = l/_Internal::bits_per_digit;
01327     if (l % _Internal::bits_per_digit) res_len++;
01328     std::size_t m = n / _Internal::bits_per_digit;
01329     n %= _Internal::bits_per_digit;
01330     std::size_t k = _Internal::bits_per_digit - n; 
01331     digit t = 0;
01332     digit* res_data = big_int::get_mem_for_data(res_len);
01333     const digit* data = a.number->data;
01334 
01335     if(n)
01336     {
01337         if(a_len > res_len + m)t = data[a_len - 1] << k;
01338         for(std::size_t i = res_len; i > 0; i--) 
01339         { 
01340             res_data[i - 1] = (data[i + m - 1] >> n) | t; 
01341             t = data[i + m - 1] << k;
01342         }
01343     }
01344     else
01345     {
01346         std::copy(data + m, data + m + res_len, res_data);
01347     }
01348 
01349     big_int res;
01350     res.free_mem_and_alloc_number(a.number->sign, res_data, res_len);
01351     return res;
01352 }
01353 
01354 
01355 big_int operator& (const big_int& a, const big_int& b)
01356 {
01357     std::size_t reslen = std::min(a.number->len, b.number->len);
01358     if(reslen == 0)return big_int();
01359     digit *res = big_int::get_mem_for_data(reslen);
01360 
01361     for(std::size_t i = 0; i < reslen; ++i)
01362         res[i] = a.number->data[i] & b.number->data[i];
01363 
01364     big_int resbi;
01365     resbi.free_mem_and_alloc_number(1, res, reslen);
01366     return resbi;
01367 }
01368 
01369 
01370 //big_int operator| (const big_int& a, const big_int& b)
01371 //{
01372 //  std::size_t reslen = std::min(a.number->len, b.number->len);
01373 //  if(reslen == 0)return big_int();
01374 //  digit *res = big_int::get_mem_for_data(reslen);
01375 //
01376 //  for(std::size_t i = 0; i < reslen; ++i)
01377 //      res[i] = a.number->data[i] | b.number->data[i];
01378 //
01379 //  big_int resbi;
01380 //  resbi.free_mem_and_alloc_number(1, res, reslen);
01381 //  return resbi;
01382 //}
01383 //
01384 //
01385 //big_int operator^ (const big_int& a, const big_int& b)
01386 //{
01387 //  std::size_t reslen = std::min(a.number->len, b.number->len);
01388 //  if(reslen == 0)return big_int();
01389 //  digit *res = big_int::get_mem_for_data(reslen);
01390 //
01391 //  for(std::size_t i = 0; i < reslen; ++i)
01392 //      res[i] = a.number->data[i] ^ b.number->data[i];
01393 //
01394 //  big_int resbi;
01395 //  resbi.free_mem_and_alloc_number(1, res, reslen);
01396 //  return resbi;
01397 //}
01398 
01399 
01400 
01401 
01402 
01403 //bool big_int::is_even() const
01404 //{
01405 //  if(!number->len)return true;
01406 //  return !(operator[](0));
01407 //}
01408 
01409 //#ifdef ARAGELI_DISABLE_PARTICULAR_COMPILER_WARNINGS
01410 //  #pragma warning (push)
01411 //  #pragma warning (disable : 4800)
01412 //#endif
01413 //
01414 //
01415 //bool big_int::is_odd () const
01416 //{
01417 //  if(!number->len)return false;
01418 //  return (operator[](0));
01419 //}
01420 //
01421 //
01422 //#ifdef ARAGELI_DISABLE_PARTICULAR_COMPILER_WARNINGS
01423 //  #pragma warning (pop)
01424 //#endif
01425 
01426 
01427 }
01428 
01429 
01430 #endif

Generated on Thu Aug 31 17:38:04 2006 for Arageli by  doxygen 1.4.7