#5 14.09.11 23:01
Re: библиотека больших чисел
хм, тоже скоро понадобится такая библиотека, ссылка на гугл конечно очень помогла)
Offline
#6 14.09.11 23:31
Re: библиотека больших чисел
ну да, и мне тоже. я нашел вот это:
large.h
#pragma once
#include "stdafx.h"
#include <iostream>
#include <vector>
#include <iterator>
using namespace std;
typedef unsigned int uint;
typedef vector<uint> vec;
class large
{
public:
large (const char *str);
large(int i);
large(uint i=0);
large(long i);
large operator-()const;
large &operator+=(const large &y);
large &operator-=(const large &y);
large &operator*=(int y);
large &operator*=(uint y);
large &operator*=(large y);
large &operator/=(const large &divisor);
large &operator%=(const large &divisor);
large &operator <<= (uint k);
large &operator >>= (uint k);
void divide(large denom, large ", large &rem, bool RemDesired)const;
// Функция num2char преобразует объект х класса large
//в его символьное представление s в обратном порядке:
void num2char(vector<char> &s) const;
int compare(const large &y)const;
private:
vec P;
bool neg;
void SetLen(int n);
void reduce();
void DDproduct(uint A, uint B, uint &Hi, uint &Lo)const;
uint DDquotient(uint A, uint B, uint d)const;
void subtractmul(uint *a, uint *b, int n, uint &q)const;
bool normalize(large &denom, large &num,int &x)const;
void unnormalize(large &rem, int x, bool SecondDone)const;
};
large operator+(large x, const large &y);
large operator-(large x, const large &y);
large operator*(large x, const large &y);
large operator/(large x, const large &y);
large operator%(large x, const large &y);
large operator<<(large u, uint k);
large operator>>(large u, uint k);
ostream &operator<<(ostream &os, const large &x);
istream &operator>>(istream &os, large &x) ;
bool operator==(const large &x, const large &y);
bool operator<(const large &x, const large &y);
bool operator!=(const large &x, const large &y);
bool operator>(const large &x, const large &y);
large abs(large x);
large sqrt(const large &a);
large power(large x, uint n) ;
#if (defined(_MSC_VER) && !defined(_SGI_MSVC))
namespace flinn
{
template <class T>
inline const T& min(const T& a, const T& b)
{
return a < b ? a : b;
}
template <class T>
inline const T& max(const T& a, const T& b)
{
return a < b ? b : a;
}
}
#endif
large.cpp
#include "StdAfx.h"
#include "large.h"
// large.cpp: Многоразрядная целочисленная арифметика,
// значения 10, 20, 30 и 40 использующая STL.
#include <iostream>
#include <strstream>
#include <iomanip>
#include <stdlib.h>
#include <limits.h>
#include <string>
#include <ctype.h>
#include "large.h"
const uint uintmax = UINT_MAX;
const int wLen = sizeof(uint) * 8; // Количество бит
const int hLen = wLen/2;
const uint rMask = (1 << hLen) - 1;
const uint lMask = uintmax - rMask;
const uint lBit = uintmax - (uintmax >> 1);
// Добавить нули или удалить элементы в конце
//то есть в позициях старших разрядов.
void large::SetLen(int LenNew)
{
int LenOld = P.size();
if (LenNew > LenOld)
P.insert(P.end(), LenNew - LenOld, 0);
else if (LenNew < LenOld)
P.erase(P.begin()+LenNew, P.end());
};
// Удалить все старшие нули (в конце вектора Р):
void large::reduce()
{
int L = P.size();
while (L > 0 && P[L-1] == 0) L--;
P. erase (P. begin () + L, P.end());
if (L == 0) neg = 0;
};
// Разные конструкторы:
large::large(const char *str)
{
bool Neg = *str == '=';
int i = (Neg ? 1 : 0);
*this = 0;
while (str[i])
*this = *this * 10 + (str[i++] - '0');
neg = Neg;
reduce();
};
large::large(int i)
{
if (i > 0) P.push_back(i); else
if (i < 0) P.push_back(-i);
neg = i < 0;
};
large::large(long L)
{
neg = L < 0;
unsigned long UL = neg ? -L : L;
if (UL != 0)
{
#if LONG_MAX == INT_MAX
// sizeof(long) == sizeof(int):
P.push_back(uint(UL));
#else
// sizeof(long) > sizeof(int):
while (UL != 0)
{
P.push_back(uint(UL & UINT_MAX));
UL >>= wLen; }
#endif
};
};
large::large(uint i)
{
neg = 0;
if (i) P.push_back(i);
};
// Operator + определяется через += и т. п.:
large operator+(large x, const large &y){return x+=y;}
large operator-(large x, const large &y){return x-=y;}
large operator*(large x, const large &y){return x*=y;}
large operator/(large x, const large &y){return x/=y;}
large operator%(large x, const large &y){return x%=y;}
large operator<<(large u, uint k){return u <<= k;}
large operator>>(large u, uint k) {return u >>= k;}
large large::operator-()const
{
large v = *this;
if (v.P.size()) v.neg = !v.neg;
return v;
};
large &large::operator+=(const large &y)
{
if (neg != y.neg) return *this -= -y;
int i, Ly = y.P.size();
uint d, carry = 0;
// Длина результата = max. длин операндов + 1:
SetLen(flinn::max(Ly, (int)P.size()) + 1);
int L = P.size();
for (i=0; i<L; i++)
{
if (i >= Ly && carry == 0) break;
d = P[i] + carry;
// Сбросить признак переноса 'carry' в 0, если
// последнее сложение не привело
//к дополнительному переносу:
carry = d < carry;
if (i < Ly)
{
P[i] = d + y.P[i];
// Последнее сложение может привести к переносу,
// только если текущий перенос 'carry' равен 0
// (и d > 0):
if (P[i] < d) carry = 1;
} else P[i] = d;
}
reduce();
return *this;
}
large &large::operator-=(const large &y)
{
if (neg != y.neg) return *this += -y;
if (!neg && y > *this || neg && y < *this)
return *this = -(y - *this);
int i, borrow = 0, Ly = y.P.size(), L = P.size();
uint d;
for (i=0; i<L; i++)
{
if (i >= Ly && borrow == 0) break;
d = P[i] - borrow;
borrow = d > P[i];
if (i < Ly)
{
P[i] = d - y.P[i];
if (P[i] > d) borrow = 1;
} else P[i] = d;
}
reduce();
return *this;
}
large &large::operator*=(int y)
{
bool Neg = y < 0;
if (Neg) y = -y;
*this *= uint(y);
if (Neg)
neg = !neg;
return *this;
};
large &large::operator*= (uint y)
{
int len0 = P.size(), i;
uint Hi, Lo, dig = P[0], nextdig = 0;
SetLen(len0 + 1) ;
for (i=0; i<len0; i++)
{ // Вычислим произведение двух цифр dig * у;
// результатом является (Hi, Lo):
DDproduct(dig, y, Hi, Lo);
P[i] = Lo + nextdig;
dig = P[i+1];
nextdig = Hi + (P[i] < Lo) ;
}
P[i] = nextdig;
reduce();
return *this;
};
large &large::operator*=(large y)
{
int L = P.size(), Ly = y.P.size();
if (L == 0 || Ly == 0) return *this = 0;
bool DifSigns = neg != y.neg;
if (L + Ly == 2)
// L = Ly = 1: Произведение длиной в одну или две цифры:
{
uint a = P[0], b = y.P[0];
P[0] = a * b; // Предположим произведение длиной
//в одну цифру:
if (P[0] / a != b)
{
P.push_back(0);
DDproduct(a, b, P[1], P[0]);
reduce();
}
neg = DifSigns;;
return *this;
}
if (L == 1) // && Ly > 1
{
uint digit = P[0];
*this = y; // Поменять местами операнды
*this *= digit; // и вызвать operator*=(uint у)
} else
if (Ly == 1) // && L > 1
{ *this *= y.P[0]; // Вызвать operator*=(uint y)
} else
// Длины обоих операндов L и Ly больше 1:
{
int lenProd = L + Ly, i, jA, jB;
uint sumHi = 0, sumLo, Hi, Lo, sumLoOld, sumHiOld, carry=0;
large x = *this;
SetLen(lenProd); // Установить длину *this
// равной lenProd
for (i=0; i<lenProd; i++)
{
sumLo = sumHi; sumHi = carry; carry = 0;
int max_jA = flinn::min(i, L-1);
// jA <= i гарантирует jB >= 0
// jA < L, поскольку в *this имеется лишь L цифр
for (jA = flinn::max(0, i+1-Ly); jA<=max_jA; jA++)
// jA > i - Ly гарантирует jB < Ly
{
jB = i - jA;
// Другое произведение цифр, влияющее на
// позицию i (= jA + jB) произведения
// всего числа:
DDproduct(x.P[jA], y.P[jB], Hi, Lo);
sumLoOld = sumLo; sumHiOld = sumHi;
sumLo += Lo;
if (sumLo < sumLoOld)
sumHi++;
sumHi += Hi;
carry += (sumHi < sumHiOld);
}
P[i] = sumLo;
}
}
reduce();
neg = DifSigns;
return *this;
};
large &large::operator/=(const large &divisor)
{
large r;
// Разделить *this на делитель, сохраняя результат
// в *this; 0 означает, что остаток в г
//не обязан быть правильным:
divide(divisor, *this, r, 0);
return *this;
};
large &large::operator %= (const large &divisor)
{
large q;
// // Разделить *this на делитель, сохраняя результат
// // в *this; 1 означает, что остаток в г
// // должен быть правильным:
divide(divisor, q, *this, 1);
return *this;
};
large &large::operator<<= (uint k)
{
int q = k / wLen; // Количество полных слов,
if (q)
{
int i; // Увеличить длину на q:
SetLen(P.size() + q);
// Сдвинуть *this на q слов влево:
for (i=P.size()-1; i>=0; i--)
P[i] = (i < q ? 0 : P[i - q] );
k %= wLen; // Теперь к содержит оставшиеся
} // позиции сдвига.
if (k) // 0 < k < wLen:
{
int kl = wLen - k;
uint mask = (1 << k) - 1;
// маска: 00..Oil..1 (k единичных бит)
SetLen(P.size() + 1);
// Каждый из P[i] сдвигается на к позиций влево,
// а затем комбинируется по "или" с к самыми левыми
// битами правого соседа P[i-1]:
for (int i=P.size()-1; i>=0; i--)
{
P[i] <<= k;
if (i > 0)
P[i] |= (P[i-1] >> kl) & mask;
}
}
reduce();
return *this;
}
large &large::operator>>=(uint k) // Аналогично «=
{
int q = k / wLen, L = P. size();
if (q >= L){*this = 0; return *this;}
if (q)
{
for (int i=q; i<L; i++) P[i-q] = P[i];
SetLen(L - q);
k %= wLen;
if (k == 0) {reduce(); return *this;}
}
int n = P.size() - 1, kl = wLen - k;
uint mask = (1 << k) - 1;
for (int i=0; i<=n; i++)
{
P[i] >>= k;
if (i < n) P[i] |= ((P[i+1] & mask) << kl);
}
reduce();
return *this;
};
// compare возвращает: отрицательное число, если *this < у,
// ноль, если *this == у, и положительное, если *this > у.
int large::compare(const large &y)const
{
if (neg != y.neg) return y.neg - neg;
int code = 0, L = P.size(), Ly = y.P.size();
if (L == 0 || Ly == 0) code = L - Ly; else
if (L < Ly) code = -1; else
if (L > Ly) code = +1; else
for (int i = L - 1; i >= 0; i--)
{
if (P[i] > y.P[i]) {code = 1; break;} else
if (P[i] < y.P[i]) {code = -1; break;}
}
return neg ? -code : code;
};
// Двуцифровое произведение (Hi, Lo) = A * B:
void large::DDproduct(uint A, uint B, uint &Hi, uint &Lo)const
{
uint hiA = A >> hLen, loA = A & rMask,
hiB = B >> hLen, loB = B & rMask,
midl, mid2, old;
Lo = loA * loB; Hi = hiA * hiB;
midl = loA * hiB; mid2 = hiA * loB;
old = Lo;
Lo += midl << hLen;
Hi += (Lo < old) + (midl >> hLen);
old = Lo;
Lo += mid2 << hLen;
Hi += (Lo < old) + (mid2 >> hLen);
}
// Двуцифровое значение (А, В) делится на d:
uint large::DDquotient(uint A, uint B, uint d)const
{
uint left, middle, right, qHi, qLo, x, dLol,
dHi = d >> hLen, dLo = d & rMask;
qHi = A/(dHi + 1);
// Это начальное приближение к qHi может оказаться
// слишком малым,
middle = qHi * dLo;
left = qHi * dHi;
x = B - (middle << hLen);
A -= (middle >> hLen) + left + (x > B); B = x;
dLol = dLo << hLen;
// Увеличить qHi при необходимости:
while (A > dHi || (A == dHi && B >= dLol))
{
x = B - dLol;
A -= dHi + (x > B);
B = x;
qHi++;
}
qLo = ((A << hLen) | (B >> hLen))/(dHi + 1);
// Это начальное приближение к qLo может оказаться
// слишком малым.
right = qLo * dLo; middle = qLo * dHi;
x = B - right;
A -= (x > B);
B = x;
x = B - (middle << hLen);
A -= (middle >> hLen) + (x > B);
B = x;
// Увеличить qLo при необходимости:
while (A || B >= d)
{
x = B - d;
A -= (x > B);
B = x;
qLo++;
}
uint result = (qHi << hLen) + qLo;
return result == 0 && qHi > 0 ? uintmax : result;
}
// Вычесть произведение q * b из а, где а и b -
// значения длиной n цифр. Остаток а - q * b
// будет меньше, чем b, и должен быть неотрицательным.
// Последнее условие может потребовать
// уменьшения q на 1:
void large::subtractmul(uint *a, uint *b, int n, uint &q)const
{
uint Hi, Lo, d, carry = 0;
int i;
for (i=0; i<n; i++)
{
DDproduct(b[i], q, Hi, Lo);
d = a[i];
a[i] -= Lo;
if (a[i] > d) carry++;
d = a[i + 1];
a[i + 1] -= Hi + carry;
carry = a[i + 1] > d;
}
if (carry) // q was too large
{
q--; carry = 0;
for (i=0; i<n; i++)
{
d = a[i] + carry;
carry = d < carry;
a[i] = d + b[i];
if (a[i] < d)
carry = 1;
}
a[n] = 0;
}
}
// Нормализация путем сдвига denom и num влево,
// так что самая левая позиция denom станет 1;
// оба операнда выражения num/denom сдвигаются
// на х битовых позиций:
bool large::normalize(large &denom, large &num, int &x)const
{
int r = denom.P.size() - 1;
uint y = denom.P[r]; x = 0;
while ( (y & lBit) == 0) {y <<= 1; x++;}
denom <<= x; num <<= x;
// Возможно второе действие согласно К. Дж. Мифсуду
// (см. библиографию):
if (r > 0 && denom.P[r] < denom.P[r-1])
{
denom *= uintmax; num *= uintmax;
return 1;
}
return 0;
}
// Откатить нормализацию (включая поправку Мифсуда,
// если SecondDone ==1), чтобы получить
// правильный остаток:
void large::unnormalize(large &rem, int x, bool SecondDone)const
{
if (SecondDone) rem /= uintmax;
if (x > 0) rem >>= x; else rem.reduce();
}
// Разделить *this на denom, получив quot = num / denom,
// и, если RemDesired == 1, rem = num % denom:
void large::divide(large denom, large ", large &rem, bool RemDesired)const
{
int L = P.size(), Ld = denom.P.size();
if (Ld == 0) {cout << "Zero divide.\n"; return;}
bool QuotNeg = neg ^ denom.neg;
int i, r, x = 0, n;
bool RemNeg = neg, SecondDone = false;
uint q, d;
large num = *this;
num.neg = denom.neg = 0;
if (num < denom)
{
quot = 0; rem = num; rem.neg = RemNeg; return;
}
if (Ld == 1 && L == 1)
{
quot = uint(num.P[0]/denom.P[0]);
rem = uint(num.P[0]%denom.P[0]);
quot.neg = QuotNeg; rem.neg = RemNeg;
return;
} else
if (Ld == 1 && (denom.P[0] & lMask) == 0)
{ // Делитель умещается в половину слова
uint divisor = denom.P[0], dHi = 0,
ql, r, q2, dividend;
quot.SetLen(L);
for (int i=L-1; i>=0; i--)
{
dividend = (dHi << hLen) | (P[i] >> hLen);
ql = dividend/divisor; r = dividend % divisor;
dividend = (r << hLen) | (P[i] & rMask);
q2 = dividend/divisor;
dHi = dividend % divisor;
quot.P[i] = (ql << hLen) | q2;
}
quot.reduce(); rem = dHi;
quot.neg = QuotNeg; rem.neg = RemNeg;
return;
}
large numO = num, denomO = denom;
SecondDone = normalize(denom, num, x);
r = denom.P.size() - 1; n = num.P.size() - 1;
quot.SetLen(n - r);
int Lq = quot.P.size();
for (i=Lq-1; i>=0; i--) quot.P[i] = 0;
rem = num;
if (rem.P[n] >= denom.P[r])
{
rem.SetLen(rem.P.size() + 1);
n++;
quot.SetLen(Lq + 1);
}
d = denom.P[r];
for (int k=n; k>r; k--)
{
q = DDquotient(rem.P[k], rem.P[k-1], d);
subtractmul(&rem.P[k - r - 1], &denom.P[0],
r + 1, q) ;
quot.P[k - r - 1] = q;
}
quot.reduce(); quot.neg = QuotNeg;
if (RemDesired)
{
unnormalize(rem, x, SecondDone);
rem.neg = RemNeg;
}
};
bool operator==(const large &x, const large &y)
{
return x.compare(y) == 0;
};
bool operator<(const large &x, const large &y)
{
return x.compare(y) < 0;
};
bool operator!=(const large& x, const large& y)
{
return x.compare(y) != 0;
};
bool operator>(const large& x, const large& y)
{
return x.compare(y) > 0;
};
// Функция num2char преобразует объект х класса large
//в его символьное представление s в обратном порядке:
void large::num2char(vector<char> &s) const
{
large x = *this;
static uint pl0 = 1, ipl0 = 0;
if (x.P.size() == 0) s.push_back('0'); else
{
uint r;
if (pl0 == 1)
{
while (pl0 <= UINT_MAX/10)
{
pl0 *= 10;
ipl0++;
}
} // plO является максимальной степенью 10,
// представляемой с помощью uint
// LP10 = plO = pow(10, iplO)
large R, LP10 = pl0;
bool neg = x.neg;
do
{
x.divide(LP10, x, R, 1) ;
r = (R.P.size() ? R.P[0] : 0);
for (uint j=0; j<ipl0; j++)
{
s.push_back(char(r % 10 + '0'));
r /= 10;
if (r + x.P.size() == 0) break;
}
} while (x.P.size());
if (neg) s.push_back('-');
// s содержит строку в обратном порядке
}
}
ostream &operator<< (ostream &os, const large &x)
{
vector<char> s;
x.num2char(s);
vector<char>::reverse_iterator i;
for (i=s.rbegin(); i != s.rend(); ++i) os << *i;
return os;
};
istream &operator>>(istream &is, large &x)
{
char ch; x = 0;
bool neg = 0; is >> ch;
if (ch == '-' )
{
neg = 1;
is.get(ch);
}
while (isdigit(ch))
{
x = x * 10 + (ch - '0');
is.get(ch);
}
if (neg) x = -x;
is.putback(ch);
return is;
};
large abs(large a)
{
if (a < 0) a = -a;
return a;
};
large sqrt(const large &a)
{
large x = a, b = a, q;
b <<= 1;
while (b >>= 2, b > large(0)) x >>= 1;
while (x > (q = a/x) + 1 || x < q - 1)
{
x += q; x >>= 1;
}
return x < q ? x : q;
};
large power(large x, uint n)
{
large y = 1;
while (n)
{
if (n & 1) y *= x;
x *= x; n >>= 1;
}
return y;
};
Offline
