#include #include #include #include #include #include #include #include using namespace std; // base and base_digits must be consistent const int base = 1000000000; const int base_digits = 9; struct bigint { vector a; int sign; bigint() : sign(1) { } bigint(long long v) { *this = v; } bigint(const string &s) { read(s); } void operator=(const bigint &v) { sign = v.sign; a = v.a; } void operator=(long long v) { sign = 1; if (v < 0) sign = -1, v = -v; a.clear(); for (; v > 0; v = v / base) a.push_back(v % base); } bigint operator+(const bigint &v) const { if (sign == v.sign) { bigint res = v; for (int i = 0, carry = 0; i < (int) max(a.size(), v.a.size()) || carry; ++i) { if (i == (int) res.a.size()) res.a.push_back(0); res.a[i] += carry + (i < (int) a.size() ? a[i] : 0); carry = res.a[i] >= base; if (carry) res.a[i] -= base; } return res; } return *this - (-v); } bigint operator-(const bigint &v) const { if (sign == v.sign) { if (abs() >= v.abs()) { bigint res = *this; for (int i = 0, carry = 0; i < (int) v.a.size() || carry; ++i) { res.a[i] -= carry + (i < (int) v.a.size() ? v.a[i] : 0); carry = res.a[i] < 0; if (carry) res.a[i] += base; } res.trim(); return res; } return -(v - *this); } return *this + (-v); } void operator*=(int v) { if (v < 0) sign = -sign, v = -v; for (int i = 0, carry = 0; i < (int) a.size() || carry; ++i) { if (i == (int) a.size()) a.push_back(0); long long cur = a[i] * (long long) v + carry; carry = (int) (cur / base); a[i] = (int) (cur % base); //asm("divl %%ecx" : "=a"(carry), "=d"(a[i]) : "A"(cur), "c"(base)); /* int val; __asm { lea esi, cur mov eax, [esi] mov edx, [esi+4] mov ecx, base div ecx mov carry, eax mov val, edx; } a[i] = val; */ } trim(); } bigint operator*(int v) const { bigint res = *this; res *= v; return res; } friend pair divmod(const bigint &a1, const bigint &b1) { int norm = base / (b1.a.back() + 1); bigint a = a1.abs() * norm; bigint b = b1.abs() * norm; bigint q, r; q.a.resize(a.a.size()); for (int i = a.a.size() - 1; i >= 0; i--) { r *= base; r += a.a[i]; int s1 = r.a.size() <= b.a.size() ? 0 : r.a[b.a.size()]; int s2 = r.a.size() <= b.a.size() - 1 ? 0 : r.a[b.a.size() - 1]; int d = ((long long) base * s1 + s2) / b.a.back(); r -= b * d; while (r < 0) r += b, --d; q.a[i] = d; } q.sign = a1.sign * b1.sign; r.sign = a1.sign; q.trim(); r.trim(); return make_pair(q, r / norm); } bigint operator/(const bigint &v) const { return divmod(*this, v).first; } bigint operator%(const bigint &v) const { return divmod(*this, v).second; } void operator/=(int v) { if (v < 0) sign = -sign, v = -v; for (int i = (int) a.size() - 1, rem = 0; i >= 0; --i) { long long cur = a[i] + rem * (long long) base; a[i] = (int) (cur / v); rem = (int) (cur % v); } trim(); } bigint operator/(int v) const { bigint res = *this; res /= v; return res; } int operator%(int v) const { if (v < 0) v = -v; int m = 0; for (int i = a.size() - 1; i >= 0; --i) m = (a[i] + m * (long long) base) % v; return m * sign; } void operator+=(const bigint &v) { *this = *this + v; } void operator-=(const bigint &v) { *this = *this - v; } void operator*=(const bigint &v) { *this = *this * v; } void operator/=(const bigint &v) { *this = *this / v; } bool operator<(const bigint &v) const { if (sign != v.sign) return sign < v.sign; if (a.size() != v.a.size()) return a.size() * sign < v.a.size() * v.sign; for (int i = a.size() - 1; i >= 0; i--) if (a[i] != v.a[i]) return a[i] * sign < v.a[i] * sign; return false; } bool operator>(const bigint &v) const { return v < *this; } bool operator<=(const bigint &v) const { return !(v < *this); } bool operator>=(const bigint &v) const { return !(*this < v); } bool operator==(const bigint &v) const { return !(*this < v) && !(v < *this); } bool operator!=(const bigint &v) const { return *this < v || v < *this; } void trim() { while (!a.empty() && a.back() == 0) a.pop_back(); if (a.empty()) sign = 1; } bool isZero() const { return a.empty() || (a.size() == 1 && !a[0]); } bigint operator-() const { bigint res = *this; res.sign = -sign; return res; } bigint abs() const { bigint res = *this; res.sign *= res.sign; return res; } long long longValue() const { long long res = 0; for (int i = a.size() - 1; i >= 0; i--) res = res * base + a[i]; return res * sign; } friend bigint gcd(const bigint &a, const bigint &b) { return b.isZero() ? a : gcd(b, a % b); } friend bigint lcm(const bigint &a, const bigint &b) { return a / gcd(a, b) * b; } void read(const string &s) { sign = 1; a.clear(); int pos = 0; while (pos < (int) s.size() && (s[pos] == '-' || s[pos] == '+')) { if (s[pos] == '-') sign = -sign; ++pos; } for (int i = s.size() - 1; i >= pos; i -= base_digits) { int x = 0; for (int j = max(pos, i - base_digits + 1); j <= i; j++) x = x * 10 + s[j] - '0'; a.push_back(x); } trim(); } friend istream& operator>>(istream &stream, bigint &v) { string s; stream >> s; v.read(s); return stream; } friend ostream& operator<<(ostream &stream, const bigint &v) { if (v.sign == -1) stream << '-'; stream << (v.a.empty() ? 0 : v.a.back()); for (int i = (int) v.a.size() - 2; i >= 0; --i) stream << setw(base_digits) << setfill('0') << v.a[i]; return stream; } static vector convert_base(const vector &a, int old_digits, int new_digits) { vector p(max(old_digits, new_digits) + 1); p[0] = 1; for (int i = 1; i < (int) p.size(); i++) p[i] = p[i - 1] * 10; vector res; long long cur = 0; int cur_digits = 0; for (int i = 0; i < (int) a.size(); i++) { cur += a[i] * p[cur_digits]; cur_digits += old_digits; while (cur_digits >= new_digits) { res.push_back(int(cur % p[new_digits])); cur /= p[new_digits]; cur_digits -= new_digits; } } res.push_back((int) cur); while (!res.empty() && res.back() == 0) res.pop_back(); return res; } void fft(vector > & a, bool invert) const { int n = (int) a.size(); for (int i = 1, j = 0; i < n; ++i) { int bit = n >> 1; for (; j >= bit; bit >>= 1) j -= bit; j += bit; if (i < j) swap(a[i], a[j]); } for (int len = 2; len <= n; len <<= 1) { double ang = 2 * 3.14159265358979323846 / len * (invert ? -1 : 1); complex wlen(cos(ang), sin(ang)); for (int i = 0; i < n; i += len) { complex w(1); for (int j = 0; j < len / 2; ++j) { complex u = a[i + j]; complex v = a[i + j + len / 2] * w; a[i + j] = u + v; a[i + j + len / 2] = u - v; w *= wlen; } } } if (invert) for (int i = 0; i < n; ++i) a[i] /= n; } void multiply_fft(const vector &a, const vector &b, vector &res) const { vector > fa(a.begin(), a.end()); vector > fb(b.begin(), b.end()); int n = 1; while (n < (int) max(a.size(), b.size())) n <<= 1; n <<= 1; fa.resize(n); fb.resize(n); fft(fa, false); fft(fb, false); for (int i = 0; i < n; ++i) fa[i] *= fb[i]; fft(fa, true); res.resize(n); for (int i = 0, carry = 0; i < n; ++i) { res[i] = int(fa[i].real() + 0.5) + carry; carry = res[i] / 1000; res[i] %= 1000; } } bigint operator*(const bigint &v) const { bigint res; res.sign = sign * v.sign; multiply_fft(convert_base(a, base_digits, 3), convert_base(v.a, base_digits, 3), res.a); res.a = convert_base(res.a, 3, base_digits); res.trim(); return res; } bigint mul_simple(const bigint &v) const { bigint res; res.sign = sign * v.sign; res.a.resize(a.size() + v.a.size()); for (int i = 0; i < (int) a.size(); ++i) if (a[i]) for (int j = 0, carry = 0; j < (int) v.a.size() || carry; ++j) { long long cur = res.a[i + j] + (long long) a[i] * (j < (int) v.a.size() ? v.a[j] : 0) + carry; carry = (int) (cur / base); res.a[i + j] = (int) (cur % base); } res.trim(); return res; } typedef vector vll; static vll karatsubaMultiply(const vll &a, const vll &b) { int n = a.size(); vll res(n + n); if (n <= 32) { for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) res[i + j] += a[i] * b[j]; return res; } int k = n >> 1; vll a1(a.begin(), a.begin() + k); vll a2(a.begin() + k, a.end()); vll b1(b.begin(), b.begin() + k); vll b2(b.begin() + k, b.end()); vll a1b1 = karatsubaMultiply(a1, b1); vll a2b2 = karatsubaMultiply(a2, b2); for (int i = 0; i < k; i++) a2[i] += a1[i]; for (int i = 0; i < k; i++) b2[i] += b1[i]; vll r = karatsubaMultiply(a2, b2); for (int i = 0; i < (int) a1b1.size(); i++) r[i] -= a1b1[i]; for (int i = 0; i < (int) a2b2.size(); i++) r[i] -= a2b2[i]; for (int i = 0; i < (int) r.size(); i++) res[i + k] += r[i]; for (int i = 0; i < (int) a1b1.size(); i++) res[i] += a1b1[i]; for (int i = 0; i < (int) a2b2.size(); i++) res[i + n] += a2b2[i]; return res; } bigint mul_karatsuba(const bigint &v) const { vector a6 = convert_base(this->a, base_digits, 6); vector b6 = convert_base(v.a, base_digits, 6); vll a(a6.begin(), a6.end()); vll b(b6.begin(), b6.end()); while (a.size() < b.size()) a.push_back(0); while (b.size() < a.size()) b.push_back(0); while (a.size() & (a.size() - 1)) a.push_back(0), b.push_back(0); vll c = karatsubaMultiply(a, b); bigint res; res.sign = sign * v.sign; for (int i = 0, carry = 0; i < (int) c.size(); i++) { long long cur = c[i] + carry; res.a.push_back((int) (cur % 1000000)); carry = (int) (cur / 1000000); } res.a = convert_base(res.a, 6, base_digits); res.trim(); return res; } }; bigint getRandomBigint(int len) { string s; for (int i = 0; i < len; i++) s += rand() % 10 + '0'; return bigint(s); } int main() { srand(1); bigint a("99999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999"); bigint b("19999999999999999999999999999999999999999999999999999999999999999999999999999999999999998"); cout << a * b << endl; cout << a.mul_karatsuba(b) << endl; cout << a / b << endl; for (int step = 0; step < 100; step++) { bigint a = getRandomBigint(1000); bigint b = getRandomBigint(1000); bigint x = a * b; bigint y = a.mul_simple(b); bigint z = a.mul_karatsuba(b); if (x != y || x != z) { cout << a << " " << b << " " << x << " " << y << " " << z << endl; } } int steps = 1; vector x(steps), y(steps); for (int i = 0; i < steps; i++) { x[i] = getRandomBigint(60000); y[i] = getRandomBigint(60000); } clock_t start = clock(); for (int i = 0; i < steps; i++) bigint z = x[i] * y[i]; fprintf(stderr, "time=%.3lfsec\n", 0.001 * (clock() - start)); start = clock(); for (int i = 0; i < steps; i++) bigint z = x[i].mul_karatsuba(y[i]); fprintf(stderr, "time=%.3lfsec\n", 0.001 * (clock() - start)); a = getRandomBigint(10000); b = getRandomBigint(2000); start = clock(); bigint c = a / b; fprintf(stderr, "time=%.3lfsec\n", 0.001 * (clock() - start)); bigint z = 5; z = 6; cout << z << endl; }