hex_big_integer_division.hpp¶
对非负十六进制大整数求商和余数；使用 Karatsuba 与 Burnikel–Ziegler 分治除法。
\[
\displaystyle A = BQ + R,\; 0\le R < B
\]
Complexity: Time: O(n^(log_2 3)) for n hexadecimal digits. Space: O(n log n).
AC 记录：division_of_hex_big_integers。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(n^(log_2 3)) for n hexadecimal digits.
/// Space: O(n log n).
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <string_view>
#include <tuple>
#include <type_traits>
#include <utility>
#include <vector>
namespace noya {
namespace hex_big_integer_division_internal {
using namespace std;
struct BigInteger {
using M = BigInteger;
bool neg;
vector<uint32_t> bit;
static constexpr int log = 8;
BigInteger() : neg(false), bit() {}
BigInteger(bool n, const vector<uint32_t> &d) : neg(n), bit(d) {}
BigInteger(uint32_t x) : neg(false) { bit = _integer_to_vector(x); }
BigInteger(int32_t x) : neg(false) {
if (x < 0)
neg = true, x = -x;
bit = _integer_to_vector((uint32_t)x);
}
template <typename I, enable_if_t<is_integral_v<I> ||
is_same_v<I, __int128_t>> * = nullptr>
BigInteger(I x) : neg(false) {
if constexpr (is_signed_v<I> || is_same_v<I, __int128_t>) {
if (x < 0)
neg = true, x = -x;
}
bit = _integer_to_vector(x);
}
BigInteger(const string &S) : neg(false) {
assert(!S.empty());
if (S.size() == 1u && S[0] == '0')
return;
int l = 0;
if (S[0] == '-')
++l, neg = true;
for (int ie = S.size(); l < ie; ie -= log) {
int is = max(l, ie - log);
uint32_t x = 0;
for (int i = is; i < ie; i++) {
x <<= 4;
if ('0' <= S[i] and S[i] <= '9') {
x |= S[i] - '0';
} else if ('A' <= S[i] and S[i] <= 'F') {
x |= S[i] - 'A' + 10;
} else if ('a' <= S[i] and S[i] <= 'f') {
x |= S[i] - 'a' + 10;
} else {
assert(false);
}
}
bit.push_back(x);
}
}
friend M operator+(const M &lhs, const M &rhs) {
if (lhs.neg == rhs.neg)
return {lhs.neg, _add(lhs.bit, rhs.bit)};
if (_leq(lhs.bit, rhs.bit)) {
// |l| <= |r|
auto c = _sub(rhs.bit, lhs.bit);
bool n = _is_zero(c) ? false : rhs.neg;
return {n, c};
}
auto c = _sub(lhs.bit, rhs.bit);
bool n = _is_zero(c) ? false : lhs.neg;
return {n, c};
}
friend M operator-(const M &lhs, const M &rhs) { return lhs + (-rhs); }
friend M operator*(const M &lhs, const M &rhs) {
auto c = _mul(lhs.bit, rhs.bit);
bool n = _is_zero(c) ? false : (lhs.neg ^ rhs.neg);
return {n, c};
}
friend pair<M, M> divmod(const M &lhs, const M &rhs) {
auto dm = _divmod(lhs.bit, rhs.bit);
bool dn = _is_zero(dm.first) ? false : lhs.neg != rhs.neg;
bool mn = _is_zero(dm.second) ? false : lhs.neg;
return {M{dn, dm.first}, M{mn, dm.second}};
}
friend M operator/(const M &lhs, const M &rhs) {
return divmod(lhs, rhs).first;
}
friend M operator%(const M &lhs, const M &rhs) {
return divmod(lhs, rhs).second;
}
M &operator+=(const M &rhs) { return (*this) = (*this) + rhs; }
M &operator-=(const M &rhs) { return (*this) = (*this) - rhs; }
M &operator*=(const M &rhs) { return (*this) = (*this) * rhs; }
M &operator/=(const M &rhs) { return (*this) = (*this) / rhs; }
M &operator%=(const M &rhs) { return (*this) = (*this) % rhs; }
M operator-() const {
if (is_zero())
return *this;
return {!neg, bit};
}
M operator+() const { return *this; }
friend M abs(const M &m) { return {false, m.bit}; }
bool is_zero() const { return _is_zero(bit); }
friend bool operator==(const M &lhs, const M &rhs) {
return lhs.neg == rhs.neg && lhs.bit == rhs.bit;
}
friend bool operator!=(const M &lhs, const M &rhs) {
return lhs.neg != rhs.neg || lhs.bit != rhs.bit;
}
friend bool operator<(const M &lhs, const M &rhs) {
if (lhs == rhs)
return false;
return _neq_lt(lhs, rhs);
}
friend bool operator<=(const M &lhs, const M &rhs) {
if (lhs == rhs)
return true;
return _neq_lt(lhs, rhs);
}
friend bool operator>(const M &lhs, const M &rhs) {
if (lhs == rhs)
return false;
return _neq_lt(rhs, lhs);
}
friend bool operator>=(const M &lhs, const M &rhs) {
if (lhs == rhs)
return true;
return _neq_lt(rhs, lhs);
}
string to_string() const {
if (is_zero())
return "0";
string res;
if (neg)
res.push_back('-');
for (int i = _size() - 1; i >= 0; i--) {
res += _itos(bit[i], i != _size() - 1);
}
return res;
}
friend istream &operator>>(istream &is, M &m) {
string s;
is >> s;
m = M{s};
return is;
}
friend ostream &operator<<(ostream &os, const M &m) {
return os << m.to_string();
}
private:
// size
int _size() const { return bit.size(); }
// a == b
static bool _eq(const vector<uint32_t> &a, const vector<uint32_t> &b) {
return a == b;
}
// a < b
static bool _lt(const vector<uint32_t> &a, const vector<uint32_t> &b) {
return _lt(a.cbegin(), a.cend(), b.cbegin(), b.cend());
}
// a <= b
static bool _leq(const vector<uint32_t> &a, const vector<uint32_t> &b) {
return _eq(a, b) || _lt(a, b);
}
// a < b (s.t. a != b)
static bool _neq_lt(const M &lhs, const M &rhs) {
assert(lhs != rhs);
if (lhs.neg != rhs.neg)
return lhs.neg;
bool f = _lt(lhs.bit, rhs.bit);
if (f)
return !lhs.neg;
return lhs.neg;
}
// a == 0
static bool _is_zero(const vector<uint32_t> &a) { return a.empty(); }
// a == 1
static bool _is_one(const vector<uint32_t> &a) {
return (int)a.size() == 1 && a[0] == 1;
}
// 末尾 0 を削除
static void _shrink(vector<uint32_t> &a) {
while (a.size() && a.back() == 0)
a.pop_back();
}
// 末尾 0 を削除
void _shrink() {
while (_size() && bit.back() == 0)
bit.pop_back();
}
// a + b
static vector<uint32_t> _add(const vector<uint32_t> &a,
const vector<uint32_t> &b) {
vector<uint32_t> c(max<int>(a.size(), b.size()) + 1);
_add(a.cbegin(), a.cend(), b.cbegin(), b.cend(), c.begin(), c.end());
_shrink(c);
return c;
}
// a - b
static vector<uint32_t> _sub(const vector<uint32_t> &a,
const vector<uint32_t> &b) {
assert(_leq(b, a));
vector<uint32_t> c{a};
_sub(c.begin(), c.end(), b.cbegin(), b.cend());
_shrink(c);
return c;
}
// a * b
static vector<uint32_t> _mul(const vector<uint32_t> &a,
const vector<uint32_t> &b) {
if (_is_zero(a) || _is_zero(b))
return {};
if (_is_one(a))
return b;
if (_is_one(b))
return a;
vector<uint32_t> c(a.size() + b.size());
_mul(a.cbegin(), a.cend(), b.cbegin(), b.cend(), c.begin(), c.end());
_shrink(c);
return c;
}
// a / b
static pair<vector<uint32_t>, vector<uint32_t>>
_divmod(const vector<uint32_t> &a, const vector<uint32_t> &b) {
if (_is_zero(b)) {
cerr << "Divide by Zero Exception" << endl;
exit(1);
}
if (a.size() < b.size()) {
return {{}, a};
}
vector<uint32_t> q(a.size() - b.size() + 1);
vector<uint32_t> r(b.size());
_divmod(a.cbegin(), a.cend(), b.cbegin(), b.cend(), q.begin(), q.end(),
r.begin(), r.end());
_shrink(q);
_shrink(r);
return {q, r};
}
using iter = typename vector<uint32_t>::iterator;
using citer = typename vector<uint32_t>::const_iterator;
// a + b
static void _add(citer a, citer ae, citer b, citer be, iter c, iter ce) {
if (ae - a < be - b) {
swap(a, b);
swap(ae, be);
}
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t csi = ce - c;
assert(asi <= csi);
uint32_t car = 0;
for (size_t i = 0; i < bsi; i++) {
uint64_t v = (uint64_t)a[i] + b[i] + car;
c[i] = (uint32_t)v;
car = v >> 32;
}
for (size_t i = bsi; i < asi; i++) {
uint64_t v = (uint64_t)a[i] + car;
c[i] = (uint32_t)v;
car = v >> 32;
}
if (car != 0) {
assert(c + asi < ce);
c[asi] = car;
}
}
// a += b
static void _add(iter a, iter ae, citer b, citer be) {
const size_t asi = ae - a;
const size_t bsi = be - b;
uint32_t car = 0;
for (size_t i = 0; i < bsi; i++) {
uint64_t v = (uint64_t)a[i] + b[i] + car;
a[i] = (uint32_t)v;
car = v >> 32;
}
for (size_t i = bsi; car != 0 && i < asi; i++) {
uint64_t v = (uint64_t)a[i] + car;
a[i] = (uint32_t)v;
car = v >> 32;
}
}
// a -= b
static void _sub(iter a, iter ae, citer b, citer be) {
const size_t asi = ae - a;
const size_t bsi = be - b;
int32_t car = 0;
for (size_t i = 0; i < bsi; i++) {
int64_t v = (int64_t)a[i] - b[i] + car;
a[i] = v;
car = v >> 32;
}
for (size_t i = bsi; car != 0 && i < asi; i++) {
int64_t v = (int64_t)a[i] + car;
a[i] = v;
car = v >> 32;
}
assert(car == 0);
}
// a * b
static void _mul(citer a, citer ae, citer b, citer be, iter c, iter ce) {
if (ae - a < be - b) {
swap(a, b);
swap(ae, be);
}
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t csi = ce - c;
assert(csi == asi + bsi);
if (bsi <= 128) {
_mul_naive(a, ae, b, be, c, ce);
return;
}
// Karatsuba reduces four half-size products to three.
// | A_hi | A_lo |
// | B_hi | B_lo |
// z0 := A_lo * B_lo
// z2 := A_hi * B_hi
// z1 := A_hi * B_lo + A_lo * B_hi
// z1 = (A_hi+A_lo) * (B_hi+B_lo) - z0 - z2
// A * B = (z2<<(2*shf)) + (z1<<shf) + z0
const size_t n = (asi + 1) >> 1;
if (bsi <= n) {
// | A_hi | A_lo |
// | 0 | B |
// A_lo * B
_mul(a, a + n, b, be, c, c + n + bsi);
vector<uint32_t> car(c + n, c + n + bsi);
fill(c + n, c + n + bsi, 0);
// A_hi * B
_mul(a + n, ae, b, be, c + n, ce);
_add(c + n, ce, car.cbegin(), car.cend());
} else {
// A_lo * B_lo
_mul(a, a + n, b, b + n, c, c + n + n);
// A_hi * B_hi
_mul(a + n, ae, b + n, be, c + n + n, ce);
vector<uint32_t> a1(n + 1);
vector<uint32_t> b1(n + 1);
vector<uint32_t> z1(2 * n + 2);
_add(a, a + n, a + n, ae, a1.begin(), a1.end());
_add(b, b + n, b + n, be, b1.begin(), b1.end());
_mul(a1.cbegin(), a1.cend(), b1.cbegin(), b1.cend(), z1.begin(),
z1.end());
_sub(z1.begin(), z1.end(), c, c + n + n);
_sub(z1.begin(), z1.end(), c + n + n, ce);
_shrink(z1);
_add(c + n, ce, z1.begin(), z1.end());
}
}
// a * b (naive)
static void _mul_naive(citer a, citer ae, citer b, citer be, iter c,
iter ce) {
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t csi = ce - c;
if (asi == 0 || bsi == 0)
return;
assert(csi == asi + bsi);
for (size_t i = 0; i < asi; i++) {
uint32_t car = 0;
for (size_t j = 0; j < bsi; j++) {
uint64_t p = 1LL * a[i] * b[j] + car + c[i + j];
c[i + j] = p;
car = p >> 32;
}
c[i + bsi] = car;
}
}
static const int DNT = 64;
// a / b
static void _divmod(citer a, citer ae, citer b, citer be, iter quo, iter qe,
iter rem, iter re) {
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t qsi = qe - quo;
const size_t rsi = re - rem;
assert(bsi > 0);
assert(qsi == asi - bsi + 1);
assert(rsi == bsi);
if (min(asi - bsi, bsi) <= DNT) {
_divmod_naive(a, ae, b, be, quo, qe, rem, re);
return;
}
// Round the divisor length to balanced blocks before recursive division.
size_t n;
{
size_t m = (bsi + DNT - 1) / DNT;
if (m > 1)
m = 1 << (32 - __builtin_clz(m - 1));
size_t j = (bsi + m - 1) / m;
n = j * m;
}
const int sd = n - bsi;
const int ss = __builtin_clz(*(be - 1));
vector<uint32_t> x(asi + sd + (__builtin_clz(*(ae - 1)) <= ss ? 1 : 0));
vector<uint32_t> y(n);
vector<uint32_t> r(n + 1);
vector<uint32_t> z(2 * n);
copy(a, ae, x.begin() + sd);
copy(b, be, y.begin() + sd);
_left_shift(x.begin() + sd, x.end(), ss);
_left_shift(y.begin() + sd, y.end(), ss);
size_t t = max<size_t>(2, (x.size() + n - 1) / n);
copy(x.cbegin() + (t - 2) * n, x.cend(), z.begin());
const size_t ql = qe - (quo + (t - 2) * n);
if (ql < n) {
vector<uint32_t> qq(n);
_divmod_d2n1n(z.begin(), z.end(), y.begin(), y.end(), qq.begin(),
qq.end(), r.begin(), r.end());
copy(qq.cbegin(), qq.cbegin() + ql, quo + (t - 2) * n);
} else {
_divmod_d2n1n(z.begin(), z.end(), y.begin(), y.end(), quo + (t - 2) * n,
quo + (t - 1) * n, r.begin(), r.end());
}
for (int i = t - 3; i >= 0; --i) {
copy(x.begin() + i * n, x.begin() + (i + 1) * n, z.begin());
copy(r.begin(), r.begin() + n, z.begin() + n);
fill(r.begin(), r.end(), 0);
_divmod_d2n1n(z.begin(), z.end(), y.begin(), y.end(), quo + i * n,
quo + (i + 1) * n, r.begin(), r.end());
}
_shrink(r);
copy(r.cbegin() + sd, r.cbegin() + sd + rsi, rem);
_right_shift(rem, re, ss);
}
// a(2n bits) / b(n bits)
static void _divmod_d2n1n(citer a, citer ae, citer b, citer be, iter quo,
iter qe, iter rem, iter re) {
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t qsi = qe - quo;
const size_t rsi = re - rem;
const size_t n = bsi;
assert(asi == 2 * n);
assert(qsi == n);
assert(rsi == n + 1);
assert(_lt(a + n, ae, b, be));
if (n % 2 != 0 || n <= DNT) {
_divmod_naive(a, ae, b, be, quo, qe, rem, re);
return;
}
const size_t hal = n >> 1;
vector<uint32_t> r1(n + hal + 1);
copy(a, a + hal, r1.begin());
_divmod_d3n2n(a + hal, ae, b, be, quo + hal, qe, r1.begin() + hal,
r1.end());
_divmod_d3n2n(r1.cbegin(), r1.cend() - 1, b, be, quo, quo + hal, rem, re);
}
// a(3n bits) / b(2n bits)
static void _divmod_d3n2n(citer a, citer ae, citer b, citer be, iter quo,
iter qe, iter rem, iter re) {
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t qsi = qe - quo;
const size_t rsi = re - rem;
const size_t n = bsi >> 1;
assert(asi == 3 * n);
assert(bsi == 2 * n);
assert(qsi == n);
assert(rsi == 2 * n + 1);
assert(_lt(a + n, ae, b, be));
vector<uint32_t> d(2 * n);
if (_lt(a + n + n, ae, b + n, be)) {
_divmod_d2n1n(a + n, ae, b + n, be, quo, qe, rem + n, re);
_mul(quo, qe, b, b + n, d.begin(), d.end());
} else {
fill(quo, qe, UINT32_MAX);
_add(a + n, a + n + n, b + n, be, rem + n, re);
copy(b, b + n, d.begin() + n);
_sub(d.begin(), d.end(), b, b + n);
}
copy(a, a + n, rem);
while (_lt(rem, re, d.cbegin(), d.cend())) {
_add(rem, re, b, be);
vector<uint32_t> one{1};
_sub(quo, qe, one.cbegin(), one.cend());
}
_sub(rem, re, d.cbegin(), d.cend());
}
// a / b (naive)
static void _divmod_naive(citer a, citer ae, citer b, citer be, iter quo,
iter qe, iter rem, iter re) {
const size_t asi = ae - a;
const size_t bsi = be - b;
if (bsi == 0) {
cerr << "Divide by Zero Exception" << endl;
exit(1);
}
if (bsi == 1) {
uint64_t car = 0;
for (size_t i = asi - 1; 0 <= (int)i; --i) {
uint64_t v = (car << 32) | a[i];
quo[i] = v / (*b);
car = v - quo[i] * (*b);
}
*rem = (uint32_t)car;
return;
}
if (max(asi, bsi) <= 2) {
uint64_t a64 = _iter_to_integer<uint64_t>(a, ae),
b64 = _iter_to_integer<uint64_t>(b, be);
_integer_to_iter(a64 / b64, quo, qe);
_integer_to_iter(a64 % b64, rem, re);
return;
}
if (_lt(a, ae, b, be)) {
if (asi > bsi)
ae = a + bsi;
copy(a, ae, rem);
return;
}
// B > 0xFFFFFFFF, A >= B
// 割る数のビットを大きくする
const int shf = __builtin_clz(*(be - 1));
vector<uint32_t> x(asi + (__builtin_clz(*(ae - 1)) < shf ? 1 : 0));
copy(a, ae, x.begin());
vector<uint32_t> y(b, be);
_left_shift(x.begin(), x.end(), shf);
_left_shift(y.begin(), y.end(), shf);
uint32_t yb = y.back();
vector<uint32_t> qv(x.size() - y.size() + 1);
vector<uint32_t> rv(x.end() - y.size(), x.end());
for (int i = qv.size() - 1; i >= 0; i--) {
if (rv.size() < y.size()) {
// do nothing
} else if (rv.size() == y.size()) {
if (_leq(y, rv)) {
qv[i] = 1, rv = _sub(rv, y);
}
} else {
assert(y.size() + 1 == rv.size());
uint64_t rb = rv[rv.size() - 1];
rb <<= 32;
rb |= rv[rv.size() - 2];
uint64_t q = rb / yb;
if (q > UINT32_MAX)
q = UINT32_MAX;
vector<uint32_t> yq = _mul(y, {(uint32_t)q});
// 真の商は q-2 以上 q+1 以下だが自信が無いので念のため while を回す
while (_lt(rv, yq))
q--, yq = _sub(yq, y);
rv = _sub(rv, yq);
while (_leq(y, rv))
q++, rv = _sub(rv, y);
qv[i] = q;
}
if (i)
rv.insert(begin(rv), x[i - 1]);
_shrink(rv);
}
_shrink(qv);
_right_shift(rv.begin(), rv.end(), shf);
_shrink(rv);
copy(qv.cbegin(), qv.cend(), quo);
copy(rv.cbegin(), rv.cend(), rem);
}
// int -> string
// 先頭かどうかに応じて zero padding するかを決める
static string _itos(uint32_t x, bool pad) {
string res;
for (int i = 0; i < log; i++) {
uint32_t d = x % 16;
if (d < 10)
res.push_back('0' + d);
else if (d < 16)
res.push_back('A' + d - 10);
else
assert(false);
x /= 16;
}
if (!pad) {
while (res.size() && res.back() == '0')
res.pop_back();
assert(!res.empty());
}
reverse(begin(res), end(res));
return res;
}
template <typename I, enable_if_t<is_integral_v<I> ||
is_same_v<I, __int128_t>> * = nullptr>
static I _iter_to_integer(citer it, const citer end) {
I res = 0;
int shf = 0;
while (it != end) {
res |= (I)*it << shf;
shf += 32;
++it;
}
return res;
}
template <typename I, enable_if_t<is_integral_v<I> ||
is_same_v<I, __int128_t>> * = nullptr>
static void _integer_to_iter(I x, iter it, const iter end) {
if constexpr (is_signed_v<I> || is_same_v<I, __int128_t>) {
assert(x >= 0);
}
while (x) {
assert(it < end);
*it = (uint32_t)x;
++it;
x >>= 32;
}
}
static vector<uint32_t> _integer_to_vector(uint32_t x) {
return vector<uint32_t>{x};
}
template <typename I, enable_if_t<is_integral_v<I> ||
is_same_v<I, __int128_t>> * = nullptr>
static vector<uint32_t> _integer_to_vector(I x) {
if constexpr (is_signed_v<I> || is_same_v<I, __int128_t>) {
assert(x >= 0);
}
vector<uint32_t> rt;
while (x) {
rt.push_back((uint32_t)x);
x >>= 32;
}
return rt;
}
// a << shf
static void _left_shift(const iter a, const iter ae, int shf) {
assert(shf < 32);
if (shf == 0)
return;
const int bac = 32 - shf;
uint32_t car = 0;
for (iter it = a; it < ae; ++it) {
uint32_t v = (*it << shf) | car;
car = *it >> bac;
*it = v;
}
assert(car == 0);
}
// a >> shf
static void _right_shift(const iter a, const iter ae, int shf) {
assert(shf < 32);
if (shf == 0)
return;
const int bac = 32 - shf;
uint32_t car = 0;
for (iter it = ae - 1; a <= it; --it) {
uint32_t v = (*it >> shf) | car;
car = *it << bac;
*it = v;
}
assert(car == 0);
}
// a < b
static bool _lt(citer a, citer ae, citer b, citer be) {
size_t asi = ae - a;
size_t bsi = be - b;
if (asi != bsi) {
if (asi < bsi) {
auto nz = std::find_if(b + asi, be, [](uint32_t x) { return x != 0; });
if (nz != be)
return true;
be = b + asi;
bsi = asi;
} else {
auto nz = std::find_if(a + bsi, ae, [](uint32_t x) { return x != 0; });
if (nz != ae)
return false;
ae = a + bsi;
asi = bsi;
}
}
for (int i = asi - 1; i >= 0; i--) {
if (a[i] != b[i])
return a[i] < b[i];
}
return false;
}
static void _dump(const vector<uint32_t> &a, string s = "") {
_dump(a.cbegin(), a.cend(), s);
}
static void _dump(citer a, citer ae, string s = "") {
if (!s.empty())
cerr << s << " : ";
cerr << "{ ";
for (citer it = a; it < ae; ++it) {
cerr << *it << ", ";
}
cerr << "}" << endl;
}
};
using bigint = BigInteger;
} // namespace hex_big_integer_division_internal
/// @brief Divide nonnegative hexadecimal integers stored in 32-bit limbs.
/// Karatsuba supplies subquadratic multiplication, while Burnikel-Ziegler
/// splits a normalized 2n-by-n division into balanced 3n/2-by-n and
/// n-by-n/2 subproblems; short blocks fall back to Knuth-style long division.
inline std::pair<std::string, std::string>
divide_hex_big_integers(std::string_view div, std::string_view dvs) {
using hex_big_integer_division_internal::bigint;
bigint fir{std::string(div)};
bigint sec{std::string(dvs)};
auto [qtn, rmn] = divmod(fir, sec);
return {qtn.to_string(), rmn.to_string()};
}
} // namespace noya
#ifndef NOYA_HEX_BIG_INTEGER_DIVISION_HPP
#define NOYA_HEX_BIG_INTEGER_DIVISION_HPP 1
/// @complexity Time: O(n^(log_2 3)) for n hexadecimal digits.
/// Space: O(n log n).
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <string_view>
#include <tuple>
#include <type_traits>
#include <utility>
#include <vector>
namespace noya {
namespace hex_big_integer_division_internal {
using namespace std;
struct BigInteger {
using M = BigInteger;
bool neg;
vector<uint32_t> bit;
static constexpr int log = 8;
BigInteger() : neg(false), bit() {}
BigInteger(bool n, const vector<uint32_t> &d) : neg(n), bit(d) {}
BigInteger(uint32_t x) : neg(false) { bit = _integer_to_vector(x); }
BigInteger(int32_t x) : neg(false) {
if (x < 0)
neg = true, x = -x;
bit = _integer_to_vector((uint32_t)x);
}
template <typename I, enable_if_t<is_integral_v<I> ||
is_same_v<I, __int128_t>> * = nullptr>
BigInteger(I x) : neg(false) {
if constexpr (is_signed_v<I> || is_same_v<I, __int128_t>) {
if (x < 0)
neg = true, x = -x;
}
bit = _integer_to_vector(x);
}
BigInteger(const string &S) : neg(false) {
assert(!S.empty());
if (S.size() == 1u && S[0] == '0')
return;
int l = 0;
if (S[0] == '-')
++l, neg = true;
for (int ie = S.size(); l < ie; ie -= log) {
int is = max(l, ie - log);
uint32_t x = 0;
for (int i = is; i < ie; i++) {
x <<= 4;
if ('0' <= S[i] and S[i] <= '9') {
x |= S[i] - '0';
} else if ('A' <= S[i] and S[i] <= 'F') {
x |= S[i] - 'A' + 10;
} else if ('a' <= S[i] and S[i] <= 'f') {
x |= S[i] - 'a' + 10;
} else {
assert(false);
}
}
bit.push_back(x);
}
}
friend M operator+(const M &lhs, const M &rhs) {
if (lhs.neg == rhs.neg)
return {lhs.neg, _add(lhs.bit, rhs.bit)};
if (_leq(lhs.bit, rhs.bit)) {
// |l| <= |r|
auto c = _sub(rhs.bit, lhs.bit);
bool n = _is_zero(c) ? false : rhs.neg;
return {n, c};
}
auto c = _sub(lhs.bit, rhs.bit);
bool n = _is_zero(c) ? false : lhs.neg;
return {n, c};
}
friend M operator-(const M &lhs, const M &rhs) { return lhs + (-rhs); }
friend M operator*(const M &lhs, const M &rhs) {
auto c = _mul(lhs.bit, rhs.bit);
bool n = _is_zero(c) ? false : (lhs.neg ^ rhs.neg);
return {n, c};
}
friend pair<M, M> divmod(const M &lhs, const M &rhs) {
auto dm = _divmod(lhs.bit, rhs.bit);
bool dn = _is_zero(dm.first) ? false : lhs.neg != rhs.neg;
bool mn = _is_zero(dm.second) ? false : lhs.neg;
return {M{dn, dm.first}, M{mn, dm.second}};
}
friend M operator/(const M &lhs, const M &rhs) {
return divmod(lhs, rhs).first;
}
friend M operator%(const M &lhs, const M &rhs) {
return divmod(lhs, rhs).second;
}
M &operator+=(const M &rhs) { return (*this) = (*this) + rhs; }
M &operator-=(const M &rhs) { return (*this) = (*this) - rhs; }
M &operator*=(const M &rhs) { return (*this) = (*this) * rhs; }
M &operator/=(const M &rhs) { return (*this) = (*this) / rhs; }
M &operator%=(const M &rhs) { return (*this) = (*this) % rhs; }
M operator-() const {
if (is_zero())
return *this;
return {!neg, bit};
}
M operator+() const { return *this; }
friend M abs(const M &m) { return {false, m.bit}; }
bool is_zero() const { return _is_zero(bit); }
friend bool operator==(const M &lhs, const M &rhs) {
return lhs.neg == rhs.neg && lhs.bit == rhs.bit;
}
friend bool operator!=(const M &lhs, const M &rhs) {
return lhs.neg != rhs.neg || lhs.bit != rhs.bit;
}
friend bool operator<(const M &lhs, const M &rhs) {
if (lhs == rhs)
return false;
return _neq_lt(lhs, rhs);
}
friend bool operator<=(const M &lhs, const M &rhs) {
if (lhs == rhs)
return true;
return _neq_lt(lhs, rhs);
}
friend bool operator>(const M &lhs, const M &rhs) {
if (lhs == rhs)
return false;
return _neq_lt(rhs, lhs);
}
friend bool operator>=(const M &lhs, const M &rhs) {
if (lhs == rhs)
return true;
return _neq_lt(rhs, lhs);
}
string to_string() const {
if (is_zero())
return "0";
string res;
if (neg)
res.push_back('-');
for (int i = _size() - 1; i >= 0; i--) {
res += _itos(bit[i], i != _size() - 1);
}
return res;
}
friend istream &operator>>(istream &is, M &m) {
string s;
is >> s;
m = M{s};
return is;
}
friend ostream &operator<<(ostream &os, const M &m) {
return os << m.to_string();
}
private:
// size
int _size() const { return bit.size(); }
// a == b
static bool _eq(const vector<uint32_t> &a, const vector<uint32_t> &b) {
return a == b;
}
// a < b
static bool _lt(const vector<uint32_t> &a, const vector<uint32_t> &b) {
return _lt(a.cbegin(), a.cend(), b.cbegin(), b.cend());
}
// a <= b
static bool _leq(const vector<uint32_t> &a, const vector<uint32_t> &b) {
return _eq(a, b) || _lt(a, b);
}
// a < b (s.t. a != b)
static bool _neq_lt(const M &lhs, const M &rhs) {
assert(lhs != rhs);
if (lhs.neg != rhs.neg)
return lhs.neg;
bool f = _lt(lhs.bit, rhs.bit);
if (f)
return !lhs.neg;
return lhs.neg;
}
// a == 0
static bool _is_zero(const vector<uint32_t> &a) { return a.empty(); }
// a == 1
static bool _is_one(const vector<uint32_t> &a) {
return (int)a.size() == 1 && a[0] == 1;
}
// 末尾 0 を削除
static void _shrink(vector<uint32_t> &a) {
while (a.size() && a.back() == 0)
a.pop_back();
}
// 末尾 0 を削除
void _shrink() {
while (_size() && bit.back() == 0)
bit.pop_back();
}
// a + b
static vector<uint32_t> _add(const vector<uint32_t> &a,
const vector<uint32_t> &b) {
vector<uint32_t> c(max<int>(a.size(), b.size()) + 1);
_add(a.cbegin(), a.cend(), b.cbegin(), b.cend(), c.begin(), c.end());
_shrink(c);
return c;
}
// a - b
static vector<uint32_t> _sub(const vector<uint32_t> &a,
const vector<uint32_t> &b) {
assert(_leq(b, a));
vector<uint32_t> c{a};
_sub(c.begin(), c.end(), b.cbegin(), b.cend());
_shrink(c);
return c;
}
// a * b
static vector<uint32_t> _mul(const vector<uint32_t> &a,
const vector<uint32_t> &b) {
if (_is_zero(a) || _is_zero(b))
return {};
if (_is_one(a))
return b;
if (_is_one(b))
return a;
vector<uint32_t> c(a.size() + b.size());
_mul(a.cbegin(), a.cend(), b.cbegin(), b.cend(), c.begin(), c.end());
_shrink(c);
return c;
}
// a / b
static pair<vector<uint32_t>, vector<uint32_t>>
_divmod(const vector<uint32_t> &a, const vector<uint32_t> &b) {
if (_is_zero(b)) {
cerr << "Divide by Zero Exception" << endl;
exit(1);
}
if (a.size() < b.size()) {
return {{}, a};
}
vector<uint32_t> q(a.size() - b.size() + 1);
vector<uint32_t> r(b.size());
_divmod(a.cbegin(), a.cend(), b.cbegin(), b.cend(), q.begin(), q.end(),
r.begin(), r.end());
_shrink(q);
_shrink(r);
return {q, r};
}
using iter = typename vector<uint32_t>::iterator;
using citer = typename vector<uint32_t>::const_iterator;
// a + b
static void _add(citer a, citer ae, citer b, citer be, iter c, iter ce) {
if (ae - a < be - b) {
swap(a, b);
swap(ae, be);
}
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t csi = ce - c;
assert(asi <= csi);
uint32_t car = 0;
for (size_t i = 0; i < bsi; i++) {
uint64_t v = (uint64_t)a[i] + b[i] + car;
c[i] = (uint32_t)v;
car = v >> 32;
}
for (size_t i = bsi; i < asi; i++) {
uint64_t v = (uint64_t)a[i] + car;
c[i] = (uint32_t)v;
car = v >> 32;
}
if (car != 0) {
assert(c + asi < ce);
c[asi] = car;
}
}
// a += b
static void _add(iter a, iter ae, citer b, citer be) {
const size_t asi = ae - a;
const size_t bsi = be - b;
uint32_t car = 0;
for (size_t i = 0; i < bsi; i++) {
uint64_t v = (uint64_t)a[i] + b[i] + car;
a[i] = (uint32_t)v;
car = v >> 32;
}
for (size_t i = bsi; car != 0 && i < asi; i++) {
uint64_t v = (uint64_t)a[i] + car;
a[i] = (uint32_t)v;
car = v >> 32;
}
}
// a -= b
static void _sub(iter a, iter ae, citer b, citer be) {
const size_t asi = ae - a;
const size_t bsi = be - b;
int32_t car = 0;
for (size_t i = 0; i < bsi; i++) {
int64_t v = (int64_t)a[i] - b[i] + car;
a[i] = v;
car = v >> 32;
}
for (size_t i = bsi; car != 0 && i < asi; i++) {
int64_t v = (int64_t)a[i] + car;
a[i] = v;
car = v >> 32;
}
assert(car == 0);
}
// a * b
static void _mul(citer a, citer ae, citer b, citer be, iter c, iter ce) {
if (ae - a < be - b) {
swap(a, b);
swap(ae, be);
}
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t csi = ce - c;
assert(csi == asi + bsi);
if (bsi <= 128) {
_mul_naive(a, ae, b, be, c, ce);
return;
}
// Karatsuba reduces four half-size products to three.
// | A_hi | A_lo |
// | B_hi | B_lo |
// z0 := A_lo * B_lo
// z2 := A_hi * B_hi
// z1 := A_hi * B_lo + A_lo * B_hi
// z1 = (A_hi+A_lo) * (B_hi+B_lo) - z0 - z2
// A * B = (z2<<(2*shf)) + (z1<<shf) + z0
const size_t n = (asi + 1) >> 1;
if (bsi <= n) {
// | A_hi | A_lo |
// | 0 | B |
// A_lo * B
_mul(a, a + n, b, be, c, c + n + bsi);
vector<uint32_t> car(c + n, c + n + bsi);
fill(c + n, c + n + bsi, 0);
// A_hi * B
_mul(a + n, ae, b, be, c + n, ce);
_add(c + n, ce, car.cbegin(), car.cend());
} else {
// A_lo * B_lo
_mul(a, a + n, b, b + n, c, c + n + n);
// A_hi * B_hi
_mul(a + n, ae, b + n, be, c + n + n, ce);
vector<uint32_t> a1(n + 1);
vector<uint32_t> b1(n + 1);
vector<uint32_t> z1(2 * n + 2);
_add(a, a + n, a + n, ae, a1.begin(), a1.end());
_add(b, b + n, b + n, be, b1.begin(), b1.end());
_mul(a1.cbegin(), a1.cend(), b1.cbegin(), b1.cend(), z1.begin(),
z1.end());
_sub(z1.begin(), z1.end(), c, c + n + n);
_sub(z1.begin(), z1.end(), c + n + n, ce);
_shrink(z1);
_add(c + n, ce, z1.begin(), z1.end());
}
}
// a * b (naive)
static void _mul_naive(citer a, citer ae, citer b, citer be, iter c,
iter ce) {
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t csi = ce - c;
if (asi == 0 || bsi == 0)
return;
assert(csi == asi + bsi);
for (size_t i = 0; i < asi; i++) {
uint32_t car = 0;
for (size_t j = 0; j < bsi; j++) {
uint64_t p = 1LL * a[i] * b[j] + car + c[i + j];
c[i + j] = p;
car = p >> 32;
}
c[i + bsi] = car;
}
}
static const int DNT = 64;
// a / b
static void _divmod(citer a, citer ae, citer b, citer be, iter quo, iter qe,
iter rem, iter re) {
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t qsi = qe - quo;
const size_t rsi = re - rem;
assert(bsi > 0);
assert(qsi == asi - bsi + 1);
assert(rsi == bsi);
if (min(asi - bsi, bsi) <= DNT) {
_divmod_naive(a, ae, b, be, quo, qe, rem, re);
return;
}
// Round the divisor length to balanced blocks before recursive division.
size_t n;
{
size_t m = (bsi + DNT - 1) / DNT;
if (m > 1)
m = 1 << (32 - __builtin_clz(m - 1));
size_t j = (bsi + m - 1) / m;
n = j * m;
}
const int sd = n - bsi;
const int ss = __builtin_clz(*(be - 1));
vector<uint32_t> x(asi + sd + (__builtin_clz(*(ae - 1)) <= ss ? 1 : 0));
vector<uint32_t> y(n);
vector<uint32_t> r(n + 1);
vector<uint32_t> z(2 * n);
copy(a, ae, x.begin() + sd);
copy(b, be, y.begin() + sd);
_left_shift(x.begin() + sd, x.end(), ss);
_left_shift(y.begin() + sd, y.end(), ss);
size_t t = max<size_t>(2, (x.size() + n - 1) / n);
copy(x.cbegin() + (t - 2) * n, x.cend(), z.begin());
const size_t ql = qe - (quo + (t - 2) * n);
if (ql < n) {
vector<uint32_t> qq(n);
_divmod_d2n1n(z.begin(), z.end(), y.begin(), y.end(), qq.begin(),
qq.end(), r.begin(), r.end());
copy(qq.cbegin(), qq.cbegin() + ql, quo + (t - 2) * n);
} else {
_divmod_d2n1n(z.begin(), z.end(), y.begin(), y.end(), quo + (t - 2) * n,
quo + (t - 1) * n, r.begin(), r.end());
}
for (int i = t - 3; i >= 0; --i) {
copy(x.begin() + i * n, x.begin() + (i + 1) * n, z.begin());
copy(r.begin(), r.begin() + n, z.begin() + n);
fill(r.begin(), r.end(), 0);
_divmod_d2n1n(z.begin(), z.end(), y.begin(), y.end(), quo + i * n,
quo + (i + 1) * n, r.begin(), r.end());
}
_shrink(r);
copy(r.cbegin() + sd, r.cbegin() + sd + rsi, rem);
_right_shift(rem, re, ss);
}
// a(2n bits) / b(n bits)
static void _divmod_d2n1n(citer a, citer ae, citer b, citer be, iter quo,
iter qe, iter rem, iter re) {
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t qsi = qe - quo;
const size_t rsi = re - rem;
const size_t n = bsi;
assert(asi == 2 * n);
assert(qsi == n);
assert(rsi == n + 1);
assert(_lt(a + n, ae, b, be));
if (n % 2 != 0 || n <= DNT) {
_divmod_naive(a, ae, b, be, quo, qe, rem, re);
return;
}
const size_t hal = n >> 1;
vector<uint32_t> r1(n + hal + 1);
copy(a, a + hal, r1.begin());
_divmod_d3n2n(a + hal, ae, b, be, quo + hal, qe, r1.begin() + hal,
r1.end());
_divmod_d3n2n(r1.cbegin(), r1.cend() - 1, b, be, quo, quo + hal, rem, re);
}
// a(3n bits) / b(2n bits)
static void _divmod_d3n2n(citer a, citer ae, citer b, citer be, iter quo,
iter qe, iter rem, iter re) {
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t qsi = qe - quo;
const size_t rsi = re - rem;
const size_t n = bsi >> 1;
assert(asi == 3 * n);
assert(bsi == 2 * n);
assert(qsi == n);
assert(rsi == 2 * n + 1);
assert(_lt(a + n, ae, b, be));
vector<uint32_t> d(2 * n);
if (_lt(a + n + n, ae, b + n, be)) {
_divmod_d2n1n(a + n, ae, b + n, be, quo, qe, rem + n, re);
_mul(quo, qe, b, b + n, d.begin(), d.end());
} else {
fill(quo, qe, UINT32_MAX);
_add(a + n, a + n + n, b + n, be, rem + n, re);
copy(b, b + n, d.begin() + n);
_sub(d.begin(), d.end(), b, b + n);
}
copy(a, a + n, rem);
while (_lt(rem, re, d.cbegin(), d.cend())) {
_add(rem, re, b, be);
vector<uint32_t> one{1};
_sub(quo, qe, one.cbegin(), one.cend());
}
_sub(rem, re, d.cbegin(), d.cend());
}
// a / b (naive)
static void _divmod_naive(citer a, citer ae, citer b, citer be, iter quo,
iter qe, iter rem, iter re) {
const size_t asi = ae - a;
const size_t bsi = be - b;
if (bsi == 0) {
cerr << "Divide by Zero Exception" << endl;
exit(1);
}
if (bsi == 1) {
uint64_t car = 0;
for (size_t i = asi - 1; 0 <= (int)i; --i) {
uint64_t v = (car << 32) | a[i];
quo[i] = v / (*b);
car = v - quo[i] * (*b);
}
*rem = (uint32_t)car;
return;
}
if (max(asi, bsi) <= 2) {
uint64_t a64 = _iter_to_integer<uint64_t>(a, ae),
b64 = _iter_to_integer<uint64_t>(b, be);
_integer_to_iter(a64 / b64, quo, qe);
_integer_to_iter(a64 % b64, rem, re);
return;
}
if (_lt(a, ae, b, be)) {
if (asi > bsi)
ae = a + bsi;
copy(a, ae, rem);
return;
}
// B > 0xFFFFFFFF, A >= B
// 割る数のビットを大きくする
const int shf = __builtin_clz(*(be - 1));
vector<uint32_t> x(asi + (__builtin_clz(*(ae - 1)) < shf ? 1 : 0));
copy(a, ae, x.begin());
vector<uint32_t> y(b, be);
_left_shift(x.begin(), x.end(), shf);
_left_shift(y.begin(), y.end(), shf);
uint32_t yb = y.back();
vector<uint32_t> qv(x.size() - y.size() + 1);
vector<uint32_t> rv(x.end() - y.size(), x.end());
for (int i = qv.size() - 1; i >= 0; i--) {
if (rv.size() < y.size()) {
// do nothing
} else if (rv.size() == y.size()) {
if (_leq(y, rv)) {
qv[i] = 1, rv = _sub(rv, y);
}
} else {
assert(y.size() + 1 == rv.size());
uint64_t rb = rv[rv.size() - 1];
rb <<= 32;
rb |= rv[rv.size() - 2];
uint64_t q = rb / yb;
if (q > UINT32_MAX)
q = UINT32_MAX;
vector<uint32_t> yq = _mul(y, {(uint32_t)q});
// 真の商は q-2 以上 q+1 以下だが自信が無いので念のため while を回す
while (_lt(rv, yq))
q--, yq = _sub(yq, y);
rv = _sub(rv, yq);
while (_leq(y, rv))
q++, rv = _sub(rv, y);
qv[i] = q;
}
if (i)
rv.insert(begin(rv), x[i - 1]);
_shrink(rv);
}
_shrink(qv);
_right_shift(rv.begin(), rv.end(), shf);
_shrink(rv);
copy(qv.cbegin(), qv.cend(), quo);
copy(rv.cbegin(), rv.cend(), rem);
}
// int -> string
// 先頭かどうかに応じて zero padding するかを決める
static string _itos(uint32_t x, bool pad) {
string res;
for (int i = 0; i < log; i++) {
uint32_t d = x % 16;
if (d < 10)
res.push_back('0' + d);
else if (d < 16)
res.push_back('A' + d - 10);
else
assert(false);
x /= 16;
}
if (!pad) {
while (res.size() && res.back() == '0')
res.pop_back();
assert(!res.empty());
}
reverse(begin(res), end(res));
return res;
}
template <typename I, enable_if_t<is_integral_v<I> ||
is_same_v<I, __int128_t>> * = nullptr>
static I _iter_to_integer(citer it, const citer end) {
I res = 0;
int shf = 0;
while (it != end) {
res |= (I)*it << shf;
shf += 32;
++it;
}
return res;
}
template <typename I, enable_if_t<is_integral_v<I> ||
is_same_v<I, __int128_t>> * = nullptr>
static void _integer_to_iter(I x, iter it, const iter end) {
if constexpr (is_signed_v<I> || is_same_v<I, __int128_t>) {
assert(x >= 0);
}
while (x) {
assert(it < end);
*it = (uint32_t)x;
++it;
x >>= 32;
}
}
static vector<uint32_t> _integer_to_vector(uint32_t x) {
return vector<uint32_t>{x};
}
template <typename I, enable_if_t<is_integral_v<I> ||
is_same_v<I, __int128_t>> * = nullptr>
static vector<uint32_t> _integer_to_vector(I x) {
if constexpr (is_signed_v<I> || is_same_v<I, __int128_t>) {
assert(x >= 0);
}
vector<uint32_t> rt;
while (x) {
rt.push_back((uint32_t)x);
x >>= 32;
}
return rt;
}
// a << shf
static void _left_shift(const iter a, const iter ae, int shf) {
assert(shf < 32);
if (shf == 0)
return;
const int bac = 32 - shf;
uint32_t car = 0;
for (iter it = a; it < ae; ++it) {
uint32_t v = (*it << shf) | car;
car = *it >> bac;
*it = v;
}
assert(car == 0);
}
// a >> shf
static void _right_shift(const iter a, const iter ae, int shf) {
assert(shf < 32);
if (shf == 0)
return;
const int bac = 32 - shf;
uint32_t car = 0;
for (iter it = ae - 1; a <= it; --it) {
uint32_t v = (*it >> shf) | car;
car = *it << bac;
*it = v;
}
assert(car == 0);
}
// a < b
static bool _lt(citer a, citer ae, citer b, citer be) {
size_t asi = ae - a;
size_t bsi = be - b;
if (asi != bsi) {
if (asi < bsi) {
auto nz = std::find_if(b + asi, be, [](uint32_t x) { return x != 0; });
if (nz != be)
return true;
be = b + asi;
bsi = asi;
} else {
auto nz = std::find_if(a + bsi, ae, [](uint32_t x) { return x != 0; });
if (nz != ae)
return false;
ae = a + bsi;
asi = bsi;
}
}
for (int i = asi - 1; i >= 0; i--) {
if (a[i] != b[i])
return a[i] < b[i];
}
return false;
}
static void _dump(const vector<uint32_t> &a, string s = "") {
_dump(a.cbegin(), a.cend(), s);
}
static void _dump(citer a, citer ae, string s = "") {
if (!s.empty())
cerr << s << " : ";
cerr << "{ ";
for (citer it = a; it < ae; ++it) {
cerr << *it << ", ";
}
cerr << "}" << endl;
}
};
using bigint = BigInteger;
} // namespace hex_big_integer_division_internal
/// @brief Divide nonnegative hexadecimal integers stored in 32-bit limbs.
/// Karatsuba supplies subquadratic multiplication, while Burnikel-Ziegler
/// splits a normalized 2n-by-n division into balanced 3n/2-by-n and
/// n-by-n/2 subproblems; short blocks fall back to Knuth-style long division.
inline std::pair<std::string, std::string>
divide_hex_big_integers(std::string_view div, std::string_view dvs) {
using hex_big_integer_division_internal::bigint;
bigint fir{std::string(div)};
bigint sec{std::string(dvs)};
auto [qtn, rmn] = divmod(fir, sec);
return {qtn.to_string(), rmn.to_string()};
}
} // namespace noya
#endif // NOYA_HEX_BIG_INTEGER_DIVISION_HPP
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <string_view>
#include <tuple>
#include <type_traits>
#include <utility>
#include <vector>
/// @complexity Time: O(n^(log_2 3)) for n hexadecimal digits.
/// Space: O(n log n).
namespace noya {
namespace hex_big_integer_division_internal {
using namespace std;
struct BigInteger {
using M = BigInteger;
bool neg;
vector<uint32_t> bit;
static constexpr int log = 8;
BigInteger() : neg(false), bit() {}
BigInteger(bool n, const vector<uint32_t> &d) : neg(n), bit(d) {}
BigInteger(uint32_t x) : neg(false) { bit = _integer_to_vector(x); }
BigInteger(int32_t x) : neg(false) {
if (x < 0)
neg = true, x = -x;
bit = _integer_to_vector((uint32_t)x);
}
template <typename I, enable_if_t<is_integral_v<I> ||
is_same_v<I, __int128_t>> * = nullptr>
BigInteger(I x) : neg(false) {
if constexpr (is_signed_v<I> || is_same_v<I, __int128_t>) {
if (x < 0)
neg = true, x = -x;
}
bit = _integer_to_vector(x);
}
BigInteger(const string &S) : neg(false) {
assert(!S.empty());
if (S.size() == 1u && S[0] == '0')
return;
int l = 0;
if (S[0] == '-')
++l, neg = true;
for (int ie = S.size(); l < ie; ie -= log) {
int is = max(l, ie - log);
uint32_t x = 0;
for (int i = is; i < ie; i++) {
x <<= 4;
if ('0' <= S[i] and S[i] <= '9') {
x |= S[i] - '0';
} else if ('A' <= S[i] and S[i] <= 'F') {
x |= S[i] - 'A' + 10;
} else if ('a' <= S[i] and S[i] <= 'f') {
x |= S[i] - 'a' + 10;
} else {
assert(false);
}
}
bit.push_back(x);
}
}
friend M operator+(const M &lhs, const M &rhs) {
if (lhs.neg == rhs.neg)
return {lhs.neg, _add(lhs.bit, rhs.bit)};
if (_leq(lhs.bit, rhs.bit)) {
// |l| <= |r|
auto c = _sub(rhs.bit, lhs.bit);
bool n = _is_zero(c) ? false : rhs.neg;
return {n, c};
}
auto c = _sub(lhs.bit, rhs.bit);
bool n = _is_zero(c) ? false : lhs.neg;
return {n, c};
}
friend M operator-(const M &lhs, const M &rhs) { return lhs + (-rhs); }
friend M operator*(const M &lhs, const M &rhs) {
auto c = _mul(lhs.bit, rhs.bit);
bool n = _is_zero(c) ? false : (lhs.neg ^ rhs.neg);
return {n, c};
}
friend pair<M, M> divmod(const M &lhs, const M &rhs) {
auto dm = _divmod(lhs.bit, rhs.bit);
bool dn = _is_zero(dm.first) ? false : lhs.neg != rhs.neg;
bool mn = _is_zero(dm.second) ? false : lhs.neg;
return {M{dn, dm.first}, M{mn, dm.second}};
}
friend M operator/(const M &lhs, const M &rhs) {
return divmod(lhs, rhs).first;
}
friend M operator%(const M &lhs, const M &rhs) {
return divmod(lhs, rhs).second;
}
M &operator+=(const M &rhs) { return (*this) = (*this) + rhs; }
M &operator-=(const M &rhs) { return (*this) = (*this) - rhs; }
M &operator*=(const M &rhs) { return (*this) = (*this) * rhs; }
M &operator/=(const M &rhs) { return (*this) = (*this) / rhs; }
M &operator%=(const M &rhs) { return (*this) = (*this) % rhs; }
M operator-() const {
if (is_zero())
return *this;
return {!neg, bit};
}
M operator+() const { return *this; }
friend M abs(const M &m) { return {false, m.bit}; }
bool is_zero() const { return _is_zero(bit); }
friend bool operator==(const M &lhs, const M &rhs) {
return lhs.neg == rhs.neg && lhs.bit == rhs.bit;
}
friend bool operator!=(const M &lhs, const M &rhs) {
return lhs.neg != rhs.neg || lhs.bit != rhs.bit;
}
friend bool operator<(const M &lhs, const M &rhs) {
if (lhs == rhs)
return false;
return _neq_lt(lhs, rhs);
}
friend bool operator<=(const M &lhs, const M &rhs) {
if (lhs == rhs)
return true;
return _neq_lt(lhs, rhs);
}
friend bool operator>(const M &lhs, const M &rhs) {
if (lhs == rhs)
return false;
return _neq_lt(rhs, lhs);
}
friend bool operator>=(const M &lhs, const M &rhs) {
if (lhs == rhs)
return true;
return _neq_lt(rhs, lhs);
}
string to_string() const {
if (is_zero())
return "0";
string res;
if (neg)
res.push_back('-');
for (int i = _size() - 1; i >= 0; i--) {
res += _itos(bit[i], i != _size() - 1);
}
return res;
}
friend istream &operator>>(istream &is, M &m) {
string s;
is >> s;
m = M{s};
return is;
}
friend ostream &operator<<(ostream &os, const M &m) {
return os << m.to_string();
}
private:
// size
int _size() const { return bit.size(); }
// a == b
static bool _eq(const vector<uint32_t> &a, const vector<uint32_t> &b) {
return a == b;
}
// a < b
static bool _lt(const vector<uint32_t> &a, const vector<uint32_t> &b) {
return _lt(a.cbegin(), a.cend(), b.cbegin(), b.cend());
}
// a <= b
static bool _leq(const vector<uint32_t> &a, const vector<uint32_t> &b) {
return _eq(a, b) || _lt(a, b);
}
// a < b (s.t. a != b)
static bool _neq_lt(const M &lhs, const M &rhs) {
assert(lhs != rhs);
if (lhs.neg != rhs.neg)
return lhs.neg;
bool f = _lt(lhs.bit, rhs.bit);
if (f)
return !lhs.neg;
return lhs.neg;
}
// a == 0
static bool _is_zero(const vector<uint32_t> &a) { return a.empty(); }
// a == 1
static bool _is_one(const vector<uint32_t> &a) {
return (int)a.size() == 1 && a[0] == 1;
}
// 末尾 0 を削除
static void _shrink(vector<uint32_t> &a) {
while (a.size() && a.back() == 0)
a.pop_back();
}
// 末尾 0 を削除
void _shrink() {
while (_size() && bit.back() == 0)
bit.pop_back();
}
// a + b
static vector<uint32_t> _add(const vector<uint32_t> &a,
const vector<uint32_t> &b) {
vector<uint32_t> c(max<int>(a.size(), b.size()) + 1);
_add(a.cbegin(), a.cend(), b.cbegin(), b.cend(), c.begin(), c.end());
_shrink(c);
return c;
}
// a - b
static vector<uint32_t> _sub(const vector<uint32_t> &a,
const vector<uint32_t> &b) {
assert(_leq(b, a));
vector<uint32_t> c{a};
_sub(c.begin(), c.end(), b.cbegin(), b.cend());
_shrink(c);
return c;
}
// a * b
static vector<uint32_t> _mul(const vector<uint32_t> &a,
const vector<uint32_t> &b) {
if (_is_zero(a) || _is_zero(b))
return {};
if (_is_one(a))
return b;
if (_is_one(b))
return a;
vector<uint32_t> c(a.size() + b.size());
_mul(a.cbegin(), a.cend(), b.cbegin(), b.cend(), c.begin(), c.end());
_shrink(c);
return c;
}
// a / b
static pair<vector<uint32_t>, vector<uint32_t>>
_divmod(const vector<uint32_t> &a, const vector<uint32_t> &b) {
if (_is_zero(b)) {
cerr << "Divide by Zero Exception" << endl;
exit(1);
}
if (a.size() < b.size()) {
return {{}, a};
}
vector<uint32_t> q(a.size() - b.size() + 1);
vector<uint32_t> r(b.size());
_divmod(a.cbegin(), a.cend(), b.cbegin(), b.cend(), q.begin(), q.end(),
r.begin(), r.end());
_shrink(q);
_shrink(r);
return {q, r};
}
using iter = typename vector<uint32_t>::iterator;
using citer = typename vector<uint32_t>::const_iterator;
// a + b
static void _add(citer a, citer ae, citer b, citer be, iter c, iter ce) {
if (ae - a < be - b) {
swap(a, b);
swap(ae, be);
}
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t csi = ce - c;
assert(asi <= csi);
uint32_t car = 0;
for (size_t i = 0; i < bsi; i++) {
uint64_t v = (uint64_t)a[i] + b[i] + car;
c[i] = (uint32_t)v;
car = v >> 32;
}
for (size_t i = bsi; i < asi; i++) {
uint64_t v = (uint64_t)a[i] + car;
c[i] = (uint32_t)v;
car = v >> 32;
}
if (car != 0) {
assert(c + asi < ce);
c[asi] = car;
}
}
// a += b
static void _add(iter a, iter ae, citer b, citer be) {
const size_t asi = ae - a;
const size_t bsi = be - b;
uint32_t car = 0;
for (size_t i = 0; i < bsi; i++) {
uint64_t v = (uint64_t)a[i] + b[i] + car;
a[i] = (uint32_t)v;
car = v >> 32;
}
for (size_t i = bsi; car != 0 && i < asi; i++) {
uint64_t v = (uint64_t)a[i] + car;
a[i] = (uint32_t)v;
car = v >> 32;
}
}
// a -= b
static void _sub(iter a, iter ae, citer b, citer be) {
const size_t asi = ae - a;
const size_t bsi = be - b;
int32_t car = 0;
for (size_t i = 0; i < bsi; i++) {
int64_t v = (int64_t)a[i] - b[i] + car;
a[i] = v;
car = v >> 32;
}
for (size_t i = bsi; car != 0 && i < asi; i++) {
int64_t v = (int64_t)a[i] + car;
a[i] = v;
car = v >> 32;
}
assert(car == 0);
}
// a * b
static void _mul(citer a, citer ae, citer b, citer be, iter c, iter ce) {
if (ae - a < be - b) {
swap(a, b);
swap(ae, be);
}
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t csi = ce - c;
assert(csi == asi + bsi);
if (bsi <= 128) {
_mul_naive(a, ae, b, be, c, ce);
return;
}
// Karatsuba reduces four half-size products to three.
// | A_hi | A_lo |
// | B_hi | B_lo |
// z0 := A_lo * B_lo
// z2 := A_hi * B_hi
// z1 := A_hi * B_lo + A_lo * B_hi
// z1 = (A_hi+A_lo) * (B_hi+B_lo) - z0 - z2
// A * B = (z2<<(2*shf)) + (z1<<shf) + z0
const size_t n = (asi + 1) >> 1;
if (bsi <= n) {
// | A_hi | A_lo |
// | 0 | B |
// A_lo * B
_mul(a, a + n, b, be, c, c + n + bsi);
vector<uint32_t> car(c + n, c + n + bsi);
fill(c + n, c + n + bsi, 0);
// A_hi * B
_mul(a + n, ae, b, be, c + n, ce);
_add(c + n, ce, car.cbegin(), car.cend());
} else {
// A_lo * B_lo
_mul(a, a + n, b, b + n, c, c + n + n);
// A_hi * B_hi
_mul(a + n, ae, b + n, be, c + n + n, ce);
vector<uint32_t> a1(n + 1);
vector<uint32_t> b1(n + 1);
vector<uint32_t> z1(2 * n + 2);
_add(a, a + n, a + n, ae, a1.begin(), a1.end());
_add(b, b + n, b + n, be, b1.begin(), b1.end());
_mul(a1.cbegin(), a1.cend(), b1.cbegin(), b1.cend(), z1.begin(),
z1.end());
_sub(z1.begin(), z1.end(), c, c + n + n);
_sub(z1.begin(), z1.end(), c + n + n, ce);
_shrink(z1);
_add(c + n, ce, z1.begin(), z1.end());
}
}
// a * b (naive)
static void _mul_naive(citer a, citer ae, citer b, citer be, iter c,
iter ce) {
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t csi = ce - c;
if (asi == 0 || bsi == 0)
return;
assert(csi == asi + bsi);
for (size_t i = 0; i < asi; i++) {
uint32_t car = 0;
for (size_t j = 0; j < bsi; j++) {
uint64_t p = 1LL * a[i] * b[j] + car + c[i + j];
c[i + j] = p;
car = p >> 32;
}
c[i + bsi] = car;
}
}
static const int DNT = 64;
// a / b
static void _divmod(citer a, citer ae, citer b, citer be, iter quo, iter qe,
iter rem, iter re) {
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t qsi = qe - quo;
const size_t rsi = re - rem;
assert(bsi > 0);
assert(qsi == asi - bsi + 1);
assert(rsi == bsi);
if (min(asi - bsi, bsi) <= DNT) {
_divmod_naive(a, ae, b, be, quo, qe, rem, re);
return;
}
// Round the divisor length to balanced blocks before recursive division.
size_t n;
{
size_t m = (bsi + DNT - 1) / DNT;
if (m > 1)
m = 1 << (32 - __builtin_clz(m - 1));
size_t j = (bsi + m - 1) / m;
n = j * m;
}
const int sd = n - bsi;
const int ss = __builtin_clz(*(be - 1));
vector<uint32_t> x(asi + sd + (__builtin_clz(*(ae - 1)) <= ss ? 1 : 0));
vector<uint32_t> y(n);
vector<uint32_t> r(n + 1);
vector<uint32_t> z(2 * n);
copy(a, ae, x.begin() + sd);
copy(b, be, y.begin() + sd);
_left_shift(x.begin() + sd, x.end(), ss);
_left_shift(y.begin() + sd, y.end(), ss);
size_t t = max<size_t>(2, (x.size() + n - 1) / n);
copy(x.cbegin() + (t - 2) * n, x.cend(), z.begin());
const size_t ql = qe - (quo + (t - 2) * n);
if (ql < n) {
vector<uint32_t> qq(n);
_divmod_d2n1n(z.begin(), z.end(), y.begin(), y.end(), qq.begin(),
qq.end(), r.begin(), r.end());
copy(qq.cbegin(), qq.cbegin() + ql, quo + (t - 2) * n);
} else {
_divmod_d2n1n(z.begin(), z.end(), y.begin(), y.end(), quo + (t - 2) * n,
quo + (t - 1) * n, r.begin(), r.end());
}
for (int i = t - 3; i >= 0; --i) {
copy(x.begin() + i * n, x.begin() + (i + 1) * n, z.begin());
copy(r.begin(), r.begin() + n, z.begin() + n);
fill(r.begin(), r.end(), 0);
_divmod_d2n1n(z.begin(), z.end(), y.begin(), y.end(), quo + i * n,
quo + (i + 1) * n, r.begin(), r.end());
}
_shrink(r);
copy(r.cbegin() + sd, r.cbegin() + sd + rsi, rem);
_right_shift(rem, re, ss);
}
// a(2n bits) / b(n bits)
static void _divmod_d2n1n(citer a, citer ae, citer b, citer be, iter quo,
iter qe, iter rem, iter re) {
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t qsi = qe - quo;
const size_t rsi = re - rem;
const size_t n = bsi;
assert(asi == 2 * n);
assert(qsi == n);
assert(rsi == n + 1);
assert(_lt(a + n, ae, b, be));
if (n % 2 != 0 || n <= DNT) {
_divmod_naive(a, ae, b, be, quo, qe, rem, re);
return;
}
const size_t hal = n >> 1;
vector<uint32_t> r1(n + hal + 1);
copy(a, a + hal, r1.begin());
_divmod_d3n2n(a + hal, ae, b, be, quo + hal, qe, r1.begin() + hal,
r1.end());
_divmod_d3n2n(r1.cbegin(), r1.cend() - 1, b, be, quo, quo + hal, rem, re);
}
// a(3n bits) / b(2n bits)
static void _divmod_d3n2n(citer a, citer ae, citer b, citer be, iter quo,
iter qe, iter rem, iter re) {
const size_t asi = ae - a;
const size_t bsi = be - b;
const size_t qsi = qe - quo;
const size_t rsi = re - rem;
const size_t n = bsi >> 1;
assert(asi == 3 * n);
assert(bsi == 2 * n);
assert(qsi == n);
assert(rsi == 2 * n + 1);
assert(_lt(a + n, ae, b, be));
vector<uint32_t> d(2 * n);
if (_lt(a + n + n, ae, b + n, be)) {
_divmod_d2n1n(a + n, ae, b + n, be, quo, qe, rem + n, re);
_mul(quo, qe, b, b + n, d.begin(), d.end());
} else {
fill(quo, qe, UINT32_MAX);
_add(a + n, a + n + n, b + n, be, rem + n, re);
copy(b, b + n, d.begin() + n);
_sub(d.begin(), d.end(), b, b + n);
}
copy(a, a + n, rem);
while (_lt(rem, re, d.cbegin(), d.cend())) {
_add(rem, re, b, be);
vector<uint32_t> one{1};
_sub(quo, qe, one.cbegin(), one.cend());
}
_sub(rem, re, d.cbegin(), d.cend());
}
// a / b (naive)
static void _divmod_naive(citer a, citer ae, citer b, citer be, iter quo,
iter qe, iter rem, iter re) {
const size_t asi = ae - a;
const size_t bsi = be - b;
if (bsi == 0) {
cerr << "Divide by Zero Exception" << endl;
exit(1);
}
if (bsi == 1) {
uint64_t car = 0;
for (size_t i = asi - 1; 0 <= (int)i; --i) {
uint64_t v = (car << 32) | a[i];
quo[i] = v / (*b);
car = v - quo[i] * (*b);
}
*rem = (uint32_t)car;
return;
}
if (max(asi, bsi) <= 2) {
uint64_t a64 = _iter_to_integer<uint64_t>(a, ae),
b64 = _iter_to_integer<uint64_t>(b, be);
_integer_to_iter(a64 / b64, quo, qe);
_integer_to_iter(a64 % b64, rem, re);
return;
}
if (_lt(a, ae, b, be)) {
if (asi > bsi)
ae = a + bsi;
copy(a, ae, rem);
return;
}
// B > 0xFFFFFFFF, A >= B
// 割る数のビットを大きくする
const int shf = __builtin_clz(*(be - 1));
vector<uint32_t> x(asi + (__builtin_clz(*(ae - 1)) < shf ? 1 : 0));
copy(a, ae, x.begin());
vector<uint32_t> y(b, be);
_left_shift(x.begin(), x.end(), shf);
_left_shift(y.begin(), y.end(), shf);
uint32_t yb = y.back();
vector<uint32_t> qv(x.size() - y.size() + 1);
vector<uint32_t> rv(x.end() - y.size(), x.end());
for (int i = qv.size() - 1; i >= 0; i--) {
if (rv.size() < y.size()) {
// do nothing
} else if (rv.size() == y.size()) {
if (_leq(y, rv)) {
qv[i] = 1, rv = _sub(rv, y);
}
} else {
assert(y.size() + 1 == rv.size());
uint64_t rb = rv[rv.size() - 1];
rb <<= 32;
rb |= rv[rv.size() - 2];
uint64_t q = rb / yb;
if (q > UINT32_MAX)
q = UINT32_MAX;
vector<uint32_t> yq = _mul(y, {(uint32_t)q});
// 真の商は q-2 以上 q+1 以下だが自信が無いので念のため while を回す
while (_lt(rv, yq))
q--, yq = _sub(yq, y);
rv = _sub(rv, yq);
while (_leq(y, rv))
q++, rv = _sub(rv, y);
qv[i] = q;
}
if (i)
rv.insert(begin(rv), x[i - 1]);
_shrink(rv);
}
_shrink(qv);
_right_shift(rv.begin(), rv.end(), shf);
_shrink(rv);
copy(qv.cbegin(), qv.cend(), quo);
copy(rv.cbegin(), rv.cend(), rem);
}
// int -> string
// 先頭かどうかに応じて zero padding するかを決める
static string _itos(uint32_t x, bool pad) {
string res;
for (int i = 0; i < log; i++) {
uint32_t d = x % 16;
if (d < 10)
res.push_back('0' + d);
else if (d < 16)
res.push_back('A' + d - 10);
else
assert(false);
x /= 16;
}
if (!pad) {
while (res.size() && res.back() == '0')
res.pop_back();
assert(!res.empty());
}
reverse(begin(res), end(res));
return res;
}
template <typename I, enable_if_t<is_integral_v<I> ||
is_same_v<I, __int128_t>> * = nullptr>
static I _iter_to_integer(citer it, const citer end) {
I res = 0;
int shf = 0;
while (it != end) {
res |= (I)*it << shf;
shf += 32;
++it;
}
return res;
}
template <typename I, enable_if_t<is_integral_v<I> ||
is_same_v<I, __int128_t>> * = nullptr>
static void _integer_to_iter(I x, iter it, const iter end) {
if constexpr (is_signed_v<I> || is_same_v<I, __int128_t>) {
assert(x >= 0);
}
while (x) {
assert(it < end);
*it = (uint32_t)x;
++it;
x >>= 32;
}
}
static vector<uint32_t> _integer_to_vector(uint32_t x) {
return vector<uint32_t>{x};
}
template <typename I, enable_if_t<is_integral_v<I> ||
is_same_v<I, __int128_t>> * = nullptr>
static vector<uint32_t> _integer_to_vector(I x) {
if constexpr (is_signed_v<I> || is_same_v<I, __int128_t>) {
assert(x >= 0);
}
vector<uint32_t> rt;
while (x) {
rt.push_back((uint32_t)x);
x >>= 32;
}
return rt;
}
// a << shf
static void _left_shift(const iter a, const iter ae, int shf) {
assert(shf < 32);
if (shf == 0)
return;
const int bac = 32 - shf;
uint32_t car = 0;
for (iter it = a; it < ae; ++it) {
uint32_t v = (*it << shf) | car;
car = *it >> bac;
*it = v;
}
assert(car == 0);
}
// a >> shf
static void _right_shift(const iter a, const iter ae, int shf) {
assert(shf < 32);
if (shf == 0)
return;
const int bac = 32 - shf;
uint32_t car = 0;
for (iter it = ae - 1; a <= it; --it) {
uint32_t v = (*it >> shf) | car;
car = *it << bac;
*it = v;
}
assert(car == 0);
}
// a < b
static bool _lt(citer a, citer ae, citer b, citer be) {
size_t asi = ae - a;
size_t bsi = be - b;
if (asi != bsi) {
if (asi < bsi) {
auto nz = std::find_if(b + asi, be, [](uint32_t x) { return x != 0; });
if (nz != be)
return true;
be = b + asi;
bsi = asi;
} else {
auto nz = std::find_if(a + bsi, ae, [](uint32_t x) { return x != 0; });
if (nz != ae)
return false;
ae = a + bsi;
asi = bsi;
}
}
for (int i = asi - 1; i >= 0; i--) {
if (a[i] != b[i])
return a[i] < b[i];
}
return false;
}
static void _dump(const vector<uint32_t> &a, string s = "") {
_dump(a.cbegin(), a.cend(), s);
}
static void _dump(citer a, citer ae, string s = "") {
if (!s.empty())
cerr << s << " : ";
cerr << "{ ";
for (citer it = a; it < ae; ++it) {
cerr << *it << ", ";
}
cerr << "}" << endl;
}
};
using bigint = BigInteger;
} // namespace hex_big_integer_division_internal
/// @brief Divide nonnegative hexadecimal integers stored in 32-bit limbs.
/// Karatsuba supplies subquadratic multiplication, while Burnikel-Ziegler
/// splits a normalized 2n-by-n division into balanced 3n/2-by-n and
/// n-by-n/2 subproblems; short blocks fall back to Knuth-style long division.
inline std::pair<std::string, std::string>
divide_hex_big_integers(std::string_view div, std::string_view dvs) {
using hex_big_integer_division_internal::bigint;
bigint fir{std::string(div)};
bigint sec{std::string(dvs)};
auto [qtn, rmn] = divmod(fir, sec);
return {qtn.to_string(), rmn.to_string()};
}
} // namespace noya