CP-Algorithms Library
This documentation is automatically generated by competitive-verifier/competitive-verifier
cp-algo/math/decimal.hpp
- View this file on GitHub
- Last update: 2025-12-17 02:23:31+01:00
- Include:
#include "cp-algo/math/decimal.hpp"
Depends on
- cp-algo/math/bigint.hpp
- cp-algo/math/common.hpp
- cp-algo/math/cvector.hpp
- cp-algo/math/fft_simple.hpp
- cp-algo/random/rng.hpp
- cp-algo/util/big_alloc.hpp
- cp-algo/util/checkpoint.hpp
- cp-algo/util/complex.hpp
- cp-algo/util/simd.hpp
Verified with
- Division of Big Integers (verify/bigint/divmod.test.cpp)
- Division of Hex Big Integers (verify/bigint/hex_divmod.test.cpp)
Code
#ifndef CP_ALGO_MATH_DECIMAL_HPP
#define CP_ALGO_MATH_DECIMAL_HPP
#include "bigint.hpp"
#include <utility>
namespace cp_algo::math {
template<base_v base = x10>
struct decimal {
bigint<base> value;
int64_t scale; // value * base^scale
decimal(int64_t v=0, int64_t s=0): value(bigint<base>(v)), scale(s) {}
decimal(bigint<base> v, int64_t s=0): value(v), scale(s) {}
decimal& operator *= (const decimal &other) {
value *= other.value;
scale += other.scale;
return *this;
}
decimal& operator += (decimal const& other) {
if (scale < other.scale) {
value += other.value.pad(other.scale - scale);
} else {
value.pad_inplace(scale - other.scale);
value += other.value;
scale = other.scale;
}
return *this;
}
decimal& operator -= (decimal const& other) {
if (scale < other.scale) {
value -= other.value.pad(other.scale - scale);
} else {
value.pad_inplace(scale - other.scale);
value -= other.value;
scale = other.scale;
}
return *this;
}
decimal operator * (const decimal &other) const {
return decimal(*this) *= other;
}
decimal operator + (const decimal &other) const {
return decimal(*this) += other;
}
decimal operator - (const decimal &other) const {
return decimal(*this) -= other;
}
auto split() const {
auto int_part = scale >= -ssize(value.digits) ? value.top(ssize(value.digits) + scale) : bigint<base>(0);
auto frac_part = *this - decimal(int_part);
return std::pair{int_part, frac_part};
}
void print() {
auto [int_part, frac_part] = split();
print_bigint(std::cout, int_part);
if (frac_part.value != bigint<base>(0)) {
std::cout << '.';
std::cout << std::string(bigint<base>::digit_length * (-frac_part.magnitude()), '0');
frac_part.value.negative = false;
print_bigint<true>(std::cout, frac_part.value);
}
std::cout << std::endl;
}
bigint<base> trunc() const {
if (scale >= 0) {
return value.pad(scale);
} else if (-scale >= ssize(value.digits)) {
return 0;
} else {
return value.top(ssize(value.digits) + scale);
}
}
bigint<base> round() const {
if (scale >= 0) {
return value.pad(scale);
} else if (-scale > ssize(value.digits)) {
return 0;
} else {
auto res = value.top(ssize(value.digits) + scale);
if (value.digits[-scale - 1] * 2 >= bigint<base>::Base) {
res += 1;
}
return res;
}
}
decimal trunc(size_t digits) const {
digits = std::min(digits, size(value.digits));
return decimal(
value.top(digits),
scale + ssize(value.digits) - digits
);
}
auto magnitude() const {
static constexpr int64_t inf = 1e18;
if (value.digits.empty()) return -inf;
return ssize(value.digits) + scale;
}
decimal inv(int64_t precision) {
assert(precision >= 0);
int64_t lead = llround((double)bigint<base>::Base / (double)value.digits.back());
decimal d(bigint<base>(lead), -ssize(value.digits));
size_t cur = 2;
decimal amend = decimal(1) - trunc(cur) * d;
while(-amend.magnitude() < precision) {
d += d * amend;
cur = 2 * (1 - amend.magnitude());
d = d.trunc(cur);
amend = decimal(1) - trunc(cur) * d;
}
return d;
}
};
template<base_v base>
auto divmod_fast(bigint<base> const& a, int64_t b) {
// Optimized divmod for small divisors that fit in int64_t
if (b == 0) {
assert(false && "Division by zero");
}
bool neg_a = a.negative;
bool neg_b = b < 0;
b = std::abs(b);
bigint<base> quotient;
uint64_t remainder = 0;
auto n = ssize(a.digits);
for (auto i = n - 1; i >= 0; i--) {
__uint128_t val = (__uint128_t)remainder * bigint<base>::Base + a.digits[i];
uint64_t q = uint64_t(val / b);
remainder = uint64_t(val % b);
quotient.digits.push_back(q);
}
std::ranges::reverse(quotient.digits);
quotient.negative = (neg_a ^ neg_b);
quotient.normalize();
bigint<base> rem{int64_t(remainder)};
rem.negative = neg_a;
return std::pair{quotient, rem};
}
template<base_v base>
auto divmod(bigint<base> const& a, bigint<base> const& b) {
if (a < b) {
return std::pair{bigint<base>(0), a};
}
// Use fast path if b fits in int64_t
if (size(b.digits) == 1) {
int64_t b_val = b.negative ? -int64_t(b.digits[0]) : int64_t(b.digits[0]);
return divmod_fast(a, b_val);
}
// General case using decimal arithmetic
auto A = decimal<base>(a);
auto B = decimal<base>(b);
auto d = (A * B.inv(A.magnitude() - B.magnitude() + 1)).round();
auto r = a - d * b;
if (r >= b) {
d += 1;
r -= b;
}
if (r < bigint<base>(0)) {
d -= 1;
r += b;
}
return std::pair{d, r};
}
}
#endif // CP_ALGO_MATH_DECIMAL_HPP
#line 1 "cp-algo/math/decimal.hpp"
#line 1 "cp-algo/math/bigint.hpp"
#line 1 "cp-algo/util/big_alloc.hpp"
#include <set>
#include <map>
#include <deque>
#include <stack>
#include <queue>
#include <vector>
#include <string>
#include <cstddef>
#include <iostream>
#include <generator>
#include <forward_list>
// Single macro to detect POSIX platforms (Linux, Unix, macOS)
#if defined(__linux__) || defined(__unix__) || (defined(__APPLE__) && defined(__MACH__))
# define CP_ALGO_USE_MMAP 1
# include <sys/mman.h>
#else
# define CP_ALGO_USE_MMAP 0
#endif
namespace cp_algo {
template <typename T, size_t Align = 32>
class big_alloc {
static_assert( Align >= alignof(void*), "Align must be at least pointer-size");
static_assert(std::popcount(Align) == 1, "Align must be a power of two");
public:
using value_type = T;
template <class U> struct rebind { using other = big_alloc<U, Align>; };
constexpr bool operator==(const big_alloc&) const = default;
constexpr bool operator!=(const big_alloc&) const = default;
big_alloc() noexcept = default;
template <typename U, std::size_t A>
big_alloc(const big_alloc<U, A>&) noexcept {}
[[nodiscard]] T* allocate(std::size_t n) {
std::size_t padded = round_up(n * sizeof(T));
std::size_t align = std::max<std::size_t>(alignof(T), Align);
#if CP_ALGO_USE_MMAP
if (padded >= MEGABYTE) {
void* raw = mmap(nullptr, padded,
PROT_READ | PROT_WRITE,
MAP_PRIVATE | MAP_ANONYMOUS, -1, 0);
madvise(raw, padded, MADV_HUGEPAGE);
madvise(raw, padded, MADV_POPULATE_WRITE);
return static_cast<T*>(raw);
}
#endif
return static_cast<T*>(::operator new(padded, std::align_val_t(align)));
}
void deallocate(T* p, std::size_t n) noexcept {
if (!p) return;
std::size_t padded = round_up(n * sizeof(T));
std::size_t align = std::max<std::size_t>(alignof(T), Align);
#if CP_ALGO_USE_MMAP
if (padded >= MEGABYTE) { munmap(p, padded); return; }
#endif
::operator delete(p, padded, std::align_val_t(align));
}
private:
static constexpr std::size_t MEGABYTE = 1 << 20;
static constexpr std::size_t round_up(std::size_t x) noexcept {
return (x + Align - 1) / Align * Align;
}
};
template<typename T> using big_vector = std::vector<T, big_alloc<T>>;
template<typename T> using big_basic_string = std::basic_string<T, std::char_traits<T>, big_alloc<T>>;
template<typename T> using big_deque = std::deque<T, big_alloc<T>>;
template<typename T> using big_stack = std::stack<T, big_deque<T>>;
template<typename T> using big_queue = std::queue<T, big_deque<T>>;
template<typename T> using big_priority_queue = std::priority_queue<T, big_vector<T>>;
template<typename T> using big_forward_list = std::forward_list<T, big_alloc<T>>;
using big_string = big_basic_string<char>;
template<typename Key, typename Value, typename Compare = std::less<Key>>
using big_map = std::map<Key, Value, Compare, big_alloc<std::pair<const Key, Value>>>;
template<typename T, typename Compare = std::less<T>>
using big_multiset = std::multiset<T, Compare, big_alloc<T>>;
template<typename T, typename Compare = std::less<T>>
using big_set = std::set<T, Compare, big_alloc<T>>;
template<typename Ref, typename V = void>
using big_generator = std::generator<Ref, V, big_alloc<std::byte>>;
}
// Deduction guide to make elements_of with big_generator default to big_alloc
namespace std::ranges {
template<typename Ref, typename V>
elements_of(cp_algo::big_generator<Ref, V>&&) -> elements_of<cp_algo::big_generator<Ref, V>&&, cp_algo::big_alloc<std::byte>>;
}
#line 1 "cp-algo/math/fft_simple.hpp"
#line 1 "cp-algo/random/rng.hpp"
#include <chrono>
#include <random>
namespace cp_algo::random {
std::mt19937_64 gen(
std::chrono::steady_clock::now().time_since_epoch().count()
);
uint64_t rng() {
return gen();
}
}
#line 1 "cp-algo/math/common.hpp"
#include <functional>
#include <cstdint>
#include <cassert>
namespace cp_algo::math {
#ifdef CP_ALGO_MAXN
const int maxn = CP_ALGO_MAXN;
#else
const int maxn = 1 << 19;
#endif
const int magic = 64; // threshold for sizes to run the naive algo
auto bpow(auto const& x, auto n, auto const& one, auto op) {
if(n == 0) {
return one;
} else {
auto t = bpow(x, n / 2, one, op);
t = op(t, t);
if(n % 2) {
t = op(t, x);
}
return t;
}
}
auto bpow(auto x, auto n, auto ans) {
return bpow(x, n, ans, std::multiplies{});
}
template<typename T>
T bpow(T const& x, auto n) {
return bpow(x, n, T(1));
}
inline constexpr auto inv2(auto x) {
assert(x % 2);
std::make_unsigned_t<decltype(x)> y = 1;
while(y * x != 1) {
y *= 2 - x * y;
}
return y;
}
}
#line 1 "cp-algo/math/cvector.hpp"
#line 1 "cp-algo/util/simd.hpp"
#include <experimental/simd>
#line 6 "cp-algo/util/simd.hpp"
#include <memory>
#if defined(__x86_64__) && !defined(CP_ALGO_DISABLE_AVX2)
#define CP_ALGO_SIMD_AVX2_TARGET _Pragma("GCC target(\"avx2\")")
#else
#define CP_ALGO_SIMD_AVX2_TARGET
#endif
#define CP_ALGO_SIMD_PRAGMA_PUSH \
_Pragma("GCC push_options") \
CP_ALGO_SIMD_AVX2_TARGET
CP_ALGO_SIMD_PRAGMA_PUSH
namespace cp_algo {
template<typename T, size_t len>
using simd [[gnu::vector_size(len * sizeof(T))]] = T;
using u64x8 = simd<uint64_t, 8>;
using u32x16 = simd<uint32_t, 16>;
using i64x4 = simd<int64_t, 4>;
using u64x4 = simd<uint64_t, 4>;
using u32x8 = simd<uint32_t, 8>;
using u16x16 = simd<uint16_t, 16>;
using i32x4 = simd<int32_t, 4>;
using u32x4 = simd<uint32_t, 4>;
using u16x8 = simd<uint16_t, 8>;
using u16x4 = simd<uint16_t, 4>;
using i16x4 = simd<int16_t, 4>;
using u8x32 = simd<uint8_t, 32>;
using u8x8 = simd<uint8_t, 8>;
using u8x4 = simd<uint8_t, 4>;
using dx4 = simd<double, 4>;
inline dx4 abs(dx4 a) {
return dx4{
std::abs(a[0]),
std::abs(a[1]),
std::abs(a[2]),
std::abs(a[3])
};
}
// https://stackoverflow.com/a/77376595
// works for ints in (-2^51, 2^51)
static constexpr dx4 magic = dx4() + (3ULL << 51);
inline i64x4 lround(dx4 x) {
return i64x4(x + magic) - i64x4(magic);
}
inline dx4 to_double(i64x4 x) {
return dx4(x + i64x4(magic)) - magic;
}
inline dx4 round(dx4 a) {
return dx4{
std::nearbyint(a[0]),
std::nearbyint(a[1]),
std::nearbyint(a[2]),
std::nearbyint(a[3])
};
}
inline u64x4 low32(u64x4 x) {
return x & uint32_t(-1);
}
inline auto swap_bytes(auto x) {
return decltype(x)(__builtin_shufflevector(u32x8(x), u32x8(x), 1, 0, 3, 2, 5, 4, 7, 6));
}
inline u64x4 montgomery_reduce(u64x4 x, uint32_t mod, uint32_t imod) {
#ifdef __AVX2__
auto x_ninv = u64x4(_mm256_mul_epu32(__m256i(x), __m256i() + imod));
x += u64x4(_mm256_mul_epu32(__m256i(x_ninv), __m256i() + mod));
#else
auto x_ninv = u64x4(u32x8(low32(x)) * imod);
x += x_ninv * uint64_t(mod);
#endif
return swap_bytes(x);
}
inline u64x4 montgomery_mul(u64x4 x, u64x4 y, uint32_t mod, uint32_t imod) {
#ifdef __AVX2__
return montgomery_reduce(u64x4(_mm256_mul_epu32(__m256i(x), __m256i(y))), mod, imod);
#else
return montgomery_reduce(x * y, mod, imod);
#endif
}
inline u32x8 montgomery_mul(u32x8 x, u32x8 y, uint32_t mod, uint32_t imod) {
return u32x8(montgomery_mul(u64x4(x), u64x4(y), mod, imod)) |
u32x8(swap_bytes(montgomery_mul(u64x4(swap_bytes(x)), u64x4(swap_bytes(y)), mod, imod)));
}
inline dx4 rotate_right(dx4 x) {
static constexpr u64x4 shuffler = {3, 0, 1, 2};
return __builtin_shuffle(x, shuffler);
}
template<std::size_t Align = 32>
inline bool is_aligned(const auto* p) noexcept {
return (reinterpret_cast<std::uintptr_t>(p) % Align) == 0;
}
template<class Target>
inline Target& vector_cast(auto &&p) {
return *reinterpret_cast<Target*>(std::assume_aligned<alignof(Target)>(&p));
}
}
#pragma GCC pop_options
#line 1 "cp-algo/util/complex.hpp"
#line 4 "cp-algo/util/complex.hpp"
#include <cmath>
#include <type_traits>
#line 7 "cp-algo/util/complex.hpp"
CP_ALGO_SIMD_PRAGMA_PUSH
namespace cp_algo {
// Custom implementation, since std::complex is UB on non-floating types
template<typename T>
struct complex {
using value_type = T;
T x, y;
inline constexpr complex(): x(), y() {}
inline constexpr complex(T const& x): x(x), y() {}
inline constexpr complex(T const& x, T const& y): x(x), y(y) {}
inline complex& operator *= (T const& t) {x *= t; y *= t; return *this;}
inline complex& operator /= (T const& t) {x /= t; y /= t; return *this;}
inline complex operator * (T const& t) const {return complex(*this) *= t;}
inline complex operator / (T const& t) const {return complex(*this) /= t;}
inline complex& operator += (complex const& t) {x += t.x; y += t.y; return *this;}
inline complex& operator -= (complex const& t) {x -= t.x; y -= t.y; return *this;}
inline complex operator * (complex const& t) const {return {x * t.x - y * t.y, x * t.y + y * t.x};}
inline complex operator / (complex const& t) const {return *this * t.conj() / t.norm();}
inline complex operator + (complex const& t) const {return complex(*this) += t;}
inline complex operator - (complex const& t) const {return complex(*this) -= t;}
inline complex& operator *= (complex const& t) {return *this = *this * t;}
inline complex& operator /= (complex const& t) {return *this = *this / t;}
inline complex operator - () const {return {-x, -y};}
inline complex conj() const {return {x, -y};}
inline T norm() const {return x * x + y * y;}
inline T abs() const {return std::sqrt(norm());}
inline T const real() const {return x;}
inline T const imag() const {return y;}
inline T& real() {return x;}
inline T& imag() {return y;}
inline static constexpr complex polar(T r, T theta) {return {T(r * cos(theta)), T(r * sin(theta))};}
inline auto operator <=> (complex const& t) const = default;
};
template<typename T> inline complex<T> conj(complex<T> const& x) {return x.conj();}
template<typename T> inline T norm(complex<T> const& x) {return x.norm();}
template<typename T> inline T abs(complex<T> const& x) {return x.abs();}
template<typename T> inline T& real(complex<T> &x) {return x.real();}
template<typename T> inline T& imag(complex<T> &x) {return x.imag();}
template<typename T> inline T const real(complex<T> const& x) {return x.real();}
template<typename T> inline T const imag(complex<T> const& x) {return x.imag();}
template<typename T>
inline constexpr complex<T> polar(T r, T theta) {
return complex<T>::polar(r, theta);
}
template<typename T>
inline std::ostream& operator << (std::ostream &out, complex<T> const& x) {
return out << x.real() << ' ' << x.imag();
}
}
#pragma GCC pop_options
#line 1 "cp-algo/util/checkpoint.hpp"
#line 8 "cp-algo/util/checkpoint.hpp"
namespace cp_algo {
#ifdef CP_ALGO_CHECKPOINT
big_map<big_string, double> checkpoints;
double last;
#endif
template<bool final = false>
void checkpoint([[maybe_unused]] auto const& _msg) {
#ifdef CP_ALGO_CHECKPOINT
big_string msg = _msg;
double now = (double)clock() / CLOCKS_PER_SEC;
double delta = now - last;
last = now;
if(msg.size() && !final) {
checkpoints[msg] += delta;
}
if(final) {
for(auto const& [key, value] : checkpoints) {
std::cerr << key << ": " << value * 1000 << " ms\n";
}
std::cerr << "Total: " << now * 1000 << " ms\n";
}
#endif
}
template<bool final = false>
void checkpoint() {
checkpoint<final>("");
}
}
#line 7 "cp-algo/math/cvector.hpp"
#include <ranges>
#include <bit>
CP_ALGO_SIMD_PRAGMA_PUSH
namespace stdx = std::experimental;
namespace cp_algo::math::fft {
static constexpr size_t flen = 4;
using ftype = double;
using vftype = dx4;
using point = complex<ftype>;
using vpoint = complex<vftype>;
static constexpr vftype vz = {};
vpoint vi(vpoint const& r) {
return {-imag(r), real(r)};
}
struct cvector {
big_vector<vpoint> r;
cvector(size_t n) {
n = std::max(flen, std::bit_ceil(n));
r.resize(n / flen);
checkpoint("cvector create");
}
vpoint& at(size_t k) {return r[k / flen];}
vpoint at(size_t k) const {return r[k / flen];}
template<class pt = point>
inline void set(size_t k, pt const& t) {
if constexpr(std::is_same_v<pt, point>) {
real(r[k / flen])[k % flen] = real(t);
imag(r[k / flen])[k % flen] = imag(t);
} else {
at(k) = t;
}
}
template<class pt = point>
inline pt get(size_t k) const {
if constexpr(std::is_same_v<pt, point>) {
return {real(r[k / flen])[k % flen], imag(r[k / flen])[k % flen]};
} else {
return at(k);
}
}
size_t size() const {
return flen * r.size();
}
static constexpr size_t eval_arg(size_t n) {
if(n < pre_evals) {
return eval_args[n];
} else {
return eval_arg(n / 2) | (n & 1) << (std::bit_width(n) - 1);
}
}
static constexpr point eval_point(size_t n) {
if(n % 2) {
return -eval_point(n - 1);
} else if(n % 4) {
return eval_point(n - 2) * point(0, 1);
} else if(n / 4 < pre_evals) {
return evalp[n / 4];
} else {
return polar<ftype>(1., std::numbers::pi / (ftype)std::bit_floor(n) * (ftype)eval_arg(n));
}
}
static constexpr std::array<point, 32> roots = []() {
std::array<point, 32> res;
for(size_t i = 2; i < 32; i++) {
res[i] = polar<ftype>(1., std::numbers::pi / (1ull << (i - 2)));
}
return res;
}();
static constexpr point root(size_t n) {
return roots[std::bit_width(n)];
}
template<int step>
static void exec_on_eval(size_t n, size_t k, auto &&callback) {
callback(k, root(4 * step * n) * eval_point(step * k));
}
template<int step>
static void exec_on_evals(size_t n, auto &&callback) {
point factor = root(4 * step * n);
for(size_t i = 0; i < n; i++) {
callback(i, factor * eval_point(step * i));
}
}
static void do_dot_iter(point rt, vpoint& Bv, vpoint const& Av, vpoint& res) {
res += Av * Bv;
real(Bv) = rotate_right(real(Bv));
imag(Bv) = rotate_right(imag(Bv));
auto x = real(Bv)[0], y = imag(Bv)[0];
real(Bv)[0] = x * real(rt) - y * imag(rt);
imag(Bv)[0] = x * imag(rt) + y * real(rt);
}
void dot(cvector const& t) {
size_t n = this->size();
exec_on_evals<1>(n / flen, [&](size_t k, point rt) __attribute__((always_inline)) {
k *= flen;
auto [Ax, Ay] = at(k);
auto Bv = t.at(k);
vpoint res = vz;
for (size_t i = 0; i < flen; i++) {
vpoint Av = vpoint(vz + Ax[i], vz + Ay[i]);
do_dot_iter(rt, Bv, Av, res);
}
set(k, res);
});
checkpoint("dot");
}
template<bool partial = true>
void ifft() {
size_t n = size();
if constexpr (!partial) {
point pi(0, 1);
exec_on_evals<4>(n / 4, [&](size_t k, point rt) __attribute__((always_inline)) {
k *= 4;
point v1 = conj(rt);
point v2 = v1 * v1;
point v3 = v1 * v2;
auto A = get(k);
auto B = get(k + 1);
auto C = get(k + 2);
auto D = get(k + 3);
set(k, (A + B) + (C + D));
set(k + 2, ((A + B) - (C + D)) * v2);
set(k + 1, ((A - B) - pi * (C - D)) * v1);
set(k + 3, ((A - B) + pi * (C - D)) * v3);
});
}
bool parity = std::countr_zero(n) % 2;
if(parity) {
exec_on_evals<2>(n / (2 * flen), [&](size_t k, point rt) __attribute__((always_inline)) {
k *= 2 * flen;
vpoint cvrt = {vz + real(rt), vz - imag(rt)};
auto B = at(k) - at(k + flen);
at(k) += at(k + flen);
at(k + flen) = B * cvrt;
});
}
for(size_t leaf = 3 * flen; leaf < n; leaf += 4 * flen) {
size_t level = std::countr_one(leaf + 3);
for(size_t lvl = 4 + parity; lvl <= level; lvl += 2) {
size_t i = (1 << lvl) / 4;
exec_on_eval<4>(n >> lvl, leaf >> lvl, [&](size_t k, point rt) __attribute__((always_inline)) {
k <<= lvl;
vpoint v1 = {vz + real(rt), vz - imag(rt)};
vpoint v2 = v1 * v1;
vpoint v3 = v1 * v2;
for(size_t j = k; j < k + i; j += flen) {
auto A = at(j);
auto B = at(j + i);
auto C = at(j + 2 * i);
auto D = at(j + 3 * i);
at(j) = ((A + B) + (C + D));
at(j + 2 * i) = ((A + B) - (C + D)) * v2;
at(j + i) = ((A - B) - vi(C - D)) * v1;
at(j + 3 * i) = ((A - B) + vi(C - D)) * v3;
}
});
}
}
checkpoint("ifft");
for(size_t k = 0; k < n; k += flen) {
if constexpr (partial) {
set(k, get<vpoint>(k) /= vz + ftype(n / flen));
} else {
set(k, get<vpoint>(k) /= vz + ftype(n));
}
}
}
template<bool partial = true>
void fft() {
size_t n = size();
bool parity = std::countr_zero(n) % 2;
for(size_t leaf = 0; leaf < n; leaf += 4 * flen) {
size_t level = std::countr_zero(n + leaf);
level -= level % 2 != parity;
for(size_t lvl = level; lvl >= 4; lvl -= 2) {
size_t i = (1 << lvl) / 4;
exec_on_eval<4>(n >> lvl, leaf >> lvl, [&](size_t k, point rt) __attribute__((always_inline)) {
k <<= lvl;
vpoint v1 = {vz + real(rt), vz + imag(rt)};
vpoint v2 = v1 * v1;
vpoint v3 = v1 * v2;
for(size_t j = k; j < k + i; j += flen) {
auto A = at(j);
auto B = at(j + i) * v1;
auto C = at(j + 2 * i) * v2;
auto D = at(j + 3 * i) * v3;
at(j) = (A + C) + (B + D);
at(j + i) = (A + C) - (B + D);
at(j + 2 * i) = (A - C) + vi(B - D);
at(j + 3 * i) = (A - C) - vi(B - D);
}
});
}
}
if(parity) {
exec_on_evals<2>(n / (2 * flen), [&](size_t k, point rt) __attribute__((always_inline)) {
k *= 2 * flen;
vpoint vrt = {vz + real(rt), vz + imag(rt)};
auto t = at(k + flen) * vrt;
at(k + flen) = at(k) - t;
at(k) += t;
});
}
if constexpr (!partial) {
point pi(0, 1);
exec_on_evals<4>(n / 4, [&](size_t k, point rt) __attribute__((always_inline)) {
k *= 4;
point v1 = rt;
point v2 = v1 * v1;
point v3 = v1 * v2;
auto A = get(k);
auto B = get(k + 1) * v1;
auto C = get(k + 2) * v2;
auto D = get(k + 3) * v3;
set(k, (A + C) + (B + D));
set(k + 1, (A + C) - (B + D));
set(k + 2, (A - C) + pi * (B - D));
set(k + 3, (A - C) - pi * (B - D));
});
}
checkpoint("fft");
}
static constexpr size_t pre_evals = 1 << 16;
static const std::array<size_t, pre_evals> eval_args;
static const std::array<point, pre_evals> evalp;
};
const std::array<size_t, cvector::pre_evals> cvector::eval_args = []() {
std::array<size_t, pre_evals> res = {};
for(size_t i = 1; i < pre_evals; i++) {
res[i] = res[i >> 1] | (i & 1) << (std::bit_width(i) - 1);
}
return res;
}();
const std::array<point, cvector::pre_evals> cvector::evalp = []() {
std::array<point, pre_evals> res = {};
res[0] = 1;
for(size_t n = 1; n < pre_evals; n++) {
res[n] = polar<ftype>(1., std::numbers::pi * ftype(eval_args[n]) / ftype(4 * std::bit_floor(n)));
}
return res;
}();
}
#pragma GCC pop_options
#line 6 "cp-algo/math/fft_simple.hpp"
CP_ALGO_SIMD_PRAGMA_PUSH
namespace cp_algo::math::fft {
struct dft_simple {
cp_algo::math::fft::cvector cv;
dft_simple(auto const& a, size_t n): cv(n) {
for(size_t i = 0; i < std::min(std::size(a), n); i++) {
real(cv.at(i))[i % 4] = ftype(a[i]);
imag(cv.at(i))[i % 4] = ftype(i + n < std::size(a) ? a[i + n] : 0);
}
checkpoint("dft64 init");
cv.fft();
}
void dot(dft_simple const& t) {
cv.dot(t.cv);
}
void recover_mod(auto &res, size_t k) {
cv.ifft();
size_t n = cv.size();
for(size_t i = 0; i < std::min(k, n); i++) {
res[i] = llround(real(cv.get(i)));
}
for(size_t i = n; i < k; i++) {
res[i] = llround(imag(cv.get(i - n)));
}
cp_algo::checkpoint("recover mod");
}
};
// Multiplies a and b, assuming perfect precision and no overflow
void conv_simple(auto& a, auto const& b) {
if (empty(a) || empty(b)) {
a.clear();
return;
}
size_t n = a.size(), m = b.size();
size_t N = std::max(flen, std::bit_ceil(n + m - 1) / 2);
dft_simple A(a, N), B(b, N);
A.dot(B);
a.resize(n + m - 1);
A.recover_mod(a, n + m - 1);
}
}
#pragma GCC pop_options
#line 5 "cp-algo/math/bigint.hpp"
#include <bits/stdc++.h>
namespace cp_algo::math {
enum base_v {
x10 = uint64_t(1e16),
x16 = uint64_t(1ull << 60)
};
template<base_v base = x10>
struct bigint {
static constexpr uint64_t Base = uint64_t(base);
static constexpr uint16_t digit_length = base == x10 ? 16 : 15;
static constexpr uint16_t sub_base = base == x10 ? 10 : 16;
static constexpr uint32_t meta_base = base == x10 ? uint32_t(1e4) : uint32_t(1 << 15);
big_basic_string<uint64_t> digits;
bool negative;
auto operator <=> (bigint const& other) const {
// Handle zero cases
if (digits.empty() && other.digits.empty()) {
return std::strong_ordering::equal;
}
if (digits.empty()) {
return other.negative ? std::strong_ordering::greater : std::strong_ordering::less;
}
if (other.digits.empty()) {
return negative ? std::strong_ordering::less : std::strong_ordering::greater;
}
// Handle sign differences
if (negative != other.negative) {
return negative ? std::strong_ordering::less : std::strong_ordering::greater;
}
// Both have the same sign - compare magnitudes
if (digits.size() != other.digits.size()) {
auto size_cmp = digits.size() <=> other.digits.size();
// If both negative, reverse the comparison
return negative ? 0 <=> size_cmp : size_cmp;
}
// Same size, compare digits from most significant to least
for (auto i = ssize(digits) - 1; i >= 0; i--) {
auto digit_cmp = digits[i] <=> other.digits[i];
if (digit_cmp != std::strong_ordering::equal) {
return negative ? 0 <=> digit_cmp : digit_cmp;
}
}
return std::strong_ordering::equal;
}
bigint() {}
bigint(big_basic_string<uint64_t> d, bool neg): digits(std::move(d)), negative(neg) {
normalize();
}
bigint& pad_inplace(size_t to_add) {
digits.insert(0, to_add, 0);
return normalize();
}
bigint& drop_inplace(size_t to_drop) {
digits.erase(0, std::min(to_drop, size(digits)));
return normalize();
}
bigint& take_inplace(size_t to_keep) {
digits.erase(std::min(to_keep, size(digits)), std::string::npos);
return normalize();
}
bigint& top_inplace(size_t to_keep) {
if (to_keep >= size(digits)) {
return pad_inplace(to_keep - size(digits));
} else {
return drop_inplace(size(digits) - to_keep);
}
}
bigint pad(size_t to_add) const {
return bigint{big_basic_string<uint64_t>(to_add, 0) + digits, negative}.normalize();
}
bigint drop(size_t to_drop) const {
return bigint{digits.substr(std::min(to_drop, size(digits))), negative}.normalize();
}
bigint take(size_t to_keep) const {
return bigint{digits.substr(0, std::min(to_keep, size(digits))), negative}.normalize();
}
bigint top(size_t to_keep) const {
if (to_keep >= size(digits)) {
return pad(to_keep - size(digits));
} else {
return drop(size(digits) - to_keep);
}
}
bigint& normalize() {
while (!empty(digits) && digits.back() == 0) {
digits.pop_back();
}
if (digits.empty()) {
negative = false;
}
return *this;
}
bigint& negate() {
negative ^= 1;
return *this;
}
bigint operator -() {
return bigint(*this).negate();
}
bigint& operator -= (const bigint& other) {
if (negative != other.negative) {
return (negate() += other).negate().normalize();
}
digits.resize(std::max(size(digits), size(other.digits)));
bool carry = false;
auto d_ptr = std::assume_aligned<32>(digits.data());
auto o_ptr = std::assume_aligned<32>(other.digits.data());
size_t N = size(other.digits);
size_t i = 0;
for (; i < N; i++) {
d_ptr[i] -= o_ptr[i] + carry;
carry = d_ptr[i] >= base;
d_ptr[i] += carry ? uint64_t(base) : 0;
}
if (carry) {
N = size(digits);
for (; i < N && d_ptr[i] == 0; i++) {
d_ptr[i] = base - 1;
}
if (i < N) {
d_ptr[i]--;
} else {
// Two's complement: flip all digits then add 1
for (i = 0; i < N; i++) {
d_ptr[i] = base - d_ptr[i] - 1;
}
bool carry = true;
for (i = 0; i < N && carry; i++) {
d_ptr[i]++;
carry = d_ptr[i] >= base;
d_ptr[i] -= carry * base;
}
negate();
}
}
return normalize();
}
bigint& operator += (const bigint& other) {
if (negative != other.negative) {
return (negate() -= other).negate().normalize();
}
digits.resize(std::max(size(digits), size(other.digits)));
bool carry = false;
auto d_ptr = std::assume_aligned<32>(digits.data());
auto o_ptr = std::assume_aligned<32>(other.digits.data());
size_t N = size(other.digits);
size_t i = 0;
for (; i < N; i++) {
d_ptr[i] += o_ptr[i] + carry;
carry = d_ptr[i] >= base;
d_ptr[i] -= carry ? uint64_t(base) : 0;
}
if (carry) {
N = size(digits);
for (; i < N && d_ptr[i] == uint64_t(base) - 1; i++) {
d_ptr[i] = 0;
}
if (i < N) {
d_ptr[i]++;
} else {
digits.push_back(1);
}
}
return *this;
}
bigint(int64_t x) {
negative = x < 0;
x = negative ? -x : x;
digits = x ? big_basic_string<uint64_t>{uint64_t(x)} : big_basic_string<uint64_t>{};
}
bigint(std::span<char> s): negative(false) {
if (size(s) < digit_length) {
int64_t val = 0;
std::from_chars(s.data(), s.data() + size(s), val, sub_base);
*this = bigint(val);
return;
}
if (!empty(s) && s[0] == '-') {
negative = true;
s = s.subspan(1);
}
size_t len = size(s);
assert(len > 0);
size_t num_digits = (len + digit_length - 1) / digit_length;
digits.resize(num_digits);
size_t i = len;
for (size_t j = 0; j < num_digits - 1; j++) {
std::from_chars(s.data() + i - digit_length, s.data() + i, digits[j], sub_base);
i -= digit_length;
}
std::from_chars(s.data(), s.data() + i, digits.back(), sub_base);
normalize();
}
bigint operator + (const bigint& other) const {
return bigint(*this) += other;
}
bigint operator - (const bigint& other) const {
return bigint(*this) -= other;
}
void to_metabase() {
auto N = ssize(digits);
digits.resize(4 * N);
for (auto i = N - 1; i >= 0; i--) {
uint64_t val = digits[i];
digits[4 * i] = val % meta_base;
val /= meta_base;
digits[4 * i + 1] = val % meta_base;
val /= meta_base;
digits[4 * i + 2] = val % meta_base;
val /= meta_base;
digits[4 * i + 3] = val;
}
}
void from_metabase() {
auto N = (ssize(digits) + 3) / 4;
digits.resize(4 * N);
uint64_t carry = 0;
for (int i = 0; i < N; i++) {
__uint128_t val = digits[4 * i + 3];
val = val * meta_base + digits[4 * i + 2];
val = val * meta_base + digits[4 * i + 1];
val = val * meta_base + digits[4 * i];
val += carry;
digits[i] = uint64_t(val % base);
carry = uint64_t(val / base);
}
digits.resize(N);
while (carry) {
digits.push_back(carry % base);
carry /= base;
}
}
bigint& operator *= (int64_t other) {
if (other < 0) {
negative ^= 1;
other = -other;
}
if (other == 0) {
return *this = bigint(0);
} else if (other == 1) {
return *this;
}
uint64_t carry = 0;
for (auto &d: digits) {
__uint128_t val = __uint128_t(d) * other + carry;
d = uint64_t(val % base);
carry = uint64_t(val / base);
}
if (carry) {
digits.push_back(carry % base);
carry /= base;
}
return *this;
}
bigint operator * (int64_t other) const {
return bigint(*this) *= other;
}
friend bigint operator * (int64_t lhs, const bigint& rhs) {
return bigint(rhs) *= lhs;
}
bigint& mul_inplace(auto &&other) {
negative ^= other.negative;
auto n = size(digits), m = size(other.digits);
if (n < m) {
std::swap(n, m);
std::swap(digits, other.digits);
}
if (m <= 1) {
return *this *= int64_t(m == 0 ? 0 : other.digits[0]);
}
// Small m: use schoolbook long multiplication in base `Base`
// Threshold chosen empirically to avoid FFT overhead on small sizes
constexpr size_t SMALL_M_THRESHOLD = 32;
if (m <= SMALL_M_THRESHOLD) {
big_basic_string<uint64_t> res;
res.assign(n + m, 0);
for (size_t i = 0; i < n; i++) {
if (digits[i] == 0) continue;
uint64_t carry = 0;
for (size_t j = 0; j < m; j++) {
__uint128_t cur = res[i + j]
+ (__uint128_t)digits[i] * other.digits[j]
+ carry;
res[i + j] = uint64_t(cur % Base);
carry = uint64_t(cur / Base);
}
size_t k = i + m;
if (carry) {
uint64_t cur = res[k] + carry;
res[k] = cur % Base;
carry = cur / Base;
k++;
}
}
digits = std::move(res);
return normalize();
}
to_metabase();
other.to_metabase();
fft::conv_simple(digits, other.digits);
from_metabase();
return normalize();
}
bigint& operator *= (bigint const& other) {
return mul_inplace(bigint(other));
}
bigint operator * (const bigint& other) const {
return bigint(*this).mul_inplace(bigint(other));
}
};
template<base_v base>
decltype(std::cin)& operator >> (decltype(std::cin) &in, cp_algo::math::bigint<base> &x) {
std::string s;
in >> s;
x = {s};
return in;
}
template<base_v base, bool fill = true>
auto& print_digit(auto &out, uint64_t d) {
char buf[16];
auto [ptr, ec] = std::to_chars(buf, buf + sizeof(buf), d, bigint<base>::sub_base);
if constexpr (base == x16) {
std::ranges::transform(buf, buf, toupper);
}
auto len = ptr - buf;
if constexpr (fill) {
out << std::string(bigint<base>::digit_length - len, '0');
}
return out << std::string_view(buf, len);
}
template<bool fill_all = false, base_v base>
auto& print_bigint(auto &out, cp_algo::math::bigint<base> const& x) {
if (x.negative) {
out << '-';
}
if (empty(x.digits)) {
return print_digit<base, fill_all>(out, 0);
}
print_digit<base, fill_all>(out, x.digits.back());
for (auto d: x.digits | std::views::reverse | std::views::drop(1)) {
print_digit<base, true>(out, d);
}
return out;
}
template<base_v base>
decltype(std::cout)& operator << (decltype(std::cout) &out, cp_algo::math::bigint<base> const& x) {
return print_bigint(out, x);
}
}
#line 5 "cp-algo/math/decimal.hpp"
namespace cp_algo::math {
template<base_v base = x10>
struct decimal {
bigint<base> value;
int64_t scale; // value * base^scale
decimal(int64_t v=0, int64_t s=0): value(bigint<base>(v)), scale(s) {}
decimal(bigint<base> v, int64_t s=0): value(v), scale(s) {}
decimal& operator *= (const decimal &other) {
value *= other.value;
scale += other.scale;
return *this;
}
decimal& operator += (decimal const& other) {
if (scale < other.scale) {
value += other.value.pad(other.scale - scale);
} else {
value.pad_inplace(scale - other.scale);
value += other.value;
scale = other.scale;
}
return *this;
}
decimal& operator -= (decimal const& other) {
if (scale < other.scale) {
value -= other.value.pad(other.scale - scale);
} else {
value.pad_inplace(scale - other.scale);
value -= other.value;
scale = other.scale;
}
return *this;
}
decimal operator * (const decimal &other) const {
return decimal(*this) *= other;
}
decimal operator + (const decimal &other) const {
return decimal(*this) += other;
}
decimal operator - (const decimal &other) const {
return decimal(*this) -= other;
}
auto split() const {
auto int_part = scale >= -ssize(value.digits) ? value.top(ssize(value.digits) + scale) : bigint<base>(0);
auto frac_part = *this - decimal(int_part);
return std::pair{int_part, frac_part};
}
void print() {
auto [int_part, frac_part] = split();
print_bigint(std::cout, int_part);
if (frac_part.value != bigint<base>(0)) {
std::cout << '.';
std::cout << std::string(bigint<base>::digit_length * (-frac_part.magnitude()), '0');
frac_part.value.negative = false;
print_bigint<true>(std::cout, frac_part.value);
}
std::cout << std::endl;
}
bigint<base> trunc() const {
if (scale >= 0) {
return value.pad(scale);
} else if (-scale >= ssize(value.digits)) {
return 0;
} else {
return value.top(ssize(value.digits) + scale);
}
}
bigint<base> round() const {
if (scale >= 0) {
return value.pad(scale);
} else if (-scale > ssize(value.digits)) {
return 0;
} else {
auto res = value.top(ssize(value.digits) + scale);
if (value.digits[-scale - 1] * 2 >= bigint<base>::Base) {
res += 1;
}
return res;
}
}
decimal trunc(size_t digits) const {
digits = std::min(digits, size(value.digits));
return decimal(
value.top(digits),
scale + ssize(value.digits) - digits
);
}
auto magnitude() const {
static constexpr int64_t inf = 1e18;
if (value.digits.empty()) return -inf;
return ssize(value.digits) + scale;
}
decimal inv(int64_t precision) {
assert(precision >= 0);
int64_t lead = llround((double)bigint<base>::Base / (double)value.digits.back());
decimal d(bigint<base>(lead), -ssize(value.digits));
size_t cur = 2;
decimal amend = decimal(1) - trunc(cur) * d;
while(-amend.magnitude() < precision) {
d += d * amend;
cur = 2 * (1 - amend.magnitude());
d = d.trunc(cur);
amend = decimal(1) - trunc(cur) * d;
}
return d;
}
};
template<base_v base>
auto divmod_fast(bigint<base> const& a, int64_t b) {
// Optimized divmod for small divisors that fit in int64_t
if (b == 0) {
assert(false && "Division by zero");
}
bool neg_a = a.negative;
bool neg_b = b < 0;
b = std::abs(b);
bigint<base> quotient;
uint64_t remainder = 0;
auto n = ssize(a.digits);
for (auto i = n - 1; i >= 0; i--) {
__uint128_t val = (__uint128_t)remainder * bigint<base>::Base + a.digits[i];
uint64_t q = uint64_t(val / b);
remainder = uint64_t(val % b);
quotient.digits.push_back(q);
}
std::ranges::reverse(quotient.digits);
quotient.negative = (neg_a ^ neg_b);
quotient.normalize();
bigint<base> rem{int64_t(remainder)};
rem.negative = neg_a;
return std::pair{quotient, rem};
}
template<base_v base>
auto divmod(bigint<base> const& a, bigint<base> const& b) {
if (a < b) {
return std::pair{bigint<base>(0), a};
}
// Use fast path if b fits in int64_t
if (size(b.digits) == 1) {
int64_t b_val = b.negative ? -int64_t(b.digits[0]) : int64_t(b.digits[0]);
return divmod_fast(a, b_val);
}
// General case using decimal arithmetic
auto A = decimal<base>(a);
auto B = decimal<base>(b);
auto d = (A * B.inv(A.magnitude() - B.magnitude() + 1)).round();
auto r = a - d * b;
if (r >= b) {
d += 1;
r -= b;
}
if (r < bigint<base>(0)) {
d -= 1;
r += b;
}
return std::pair{d, r};
}
}
#ifndef CP_ALGO_MATH_DECIMAL_HPP
#define CP_ALGO_MATH_DECIMAL_HPP
#include "bigint.hpp"
#include <utility>
namespace cp_algo::math{template<base_v base=x10>struct decimal{bigint<base>value;int64_t scale;decimal(int64_t v=0,int64_t s=0):value(bigint<base>(v)),scale(s){}decimal(bigint<base>v,int64_t s=0):value(v),scale(s){}decimal&operator*=(const decimal&other){value*=other.value;scale+=other.scale;return*this;}decimal&operator+=(decimal const&other){if(scale<other.scale){value+=other.value.pad(other.scale-scale);}else{value.pad_inplace(scale-other.scale);value+=other.value;scale=other.scale;}return*this;}decimal&operator-=(decimal const&other){if(scale<other.scale){value-=other.value.pad(other.scale-scale);}else{value.pad_inplace(scale-other.scale);value-=other.value;scale=other.scale;}return*this;}decimal operator*(const decimal&other)const{return decimal(*this)*=other;}decimal operator+(const decimal&other)const{return decimal(*this)+=other;}decimal operator-(const decimal&other)const{return decimal(*this)-=other;}auto split()const{auto int_part=scale>=-ssize(value.digits)?value.top(ssize(value.digits)+scale):bigint<base>(0);auto frac_part=*this-decimal(int_part);return std::pair{int_part,frac_part};}void print(){auto[int_part,frac_part]=split();print_bigint(std::cout,int_part);if(frac_part.value!=bigint<base>(0)){std::cout<<'.';std::cout<<std::string(bigint<base>::digit_length*(-frac_part.magnitude()),'0');frac_part.value.negative=false;print_bigint<true>(std::cout,frac_part.value);}std::cout<<std::endl;}bigint<base>trunc()const{if(scale>=0){return value.pad(scale);}else if(-scale>=ssize(value.digits)){return 0;}else{return value.top(ssize(value.digits)+scale);}}bigint<base>round()const{if(scale>=0){return value.pad(scale);}else if(-scale>ssize(value.digits)){return 0;}else{auto res=value.top(ssize(value.digits)+scale);if(value.digits[-scale-1]*2>=bigint<base>::Base){res+=1;}return res;}}decimal trunc(size_t digits)const{digits=std::min(digits,size(value.digits));return decimal(value.top(digits),scale+ssize(value.digits)-digits);}auto magnitude()const{static constexpr int64_t inf=1e18;if(value.digits.empty())return-inf;return ssize(value.digits)+scale;}decimal inv(int64_t precision){assert(precision>=0);int64_t lead=llround((double)bigint<base>::Base/(double)value.digits.back());decimal d(bigint<base>(lead),-ssize(value.digits));size_t cur=2;decimal amend=decimal(1)-trunc(cur)*d;while(-amend.magnitude()<precision){d+=d*amend;cur=2*(1-amend.magnitude());d=d.trunc(cur);amend=decimal(1)-trunc(cur)*d;}return d;}};template<base_v base>auto divmod_fast(bigint<base>const&a,int64_t b){if(b==0){assert(false&&"Division by zero");}bool neg_a=a.negative;bool neg_b=b<0;b=std::abs(b);bigint<base>quotient;uint64_t remainder=0;auto n=ssize(a.digits);for(auto i=n-1;i>=0;i--){__uint128_t val=(__uint128_t)remainder*bigint<base>::Base+a.digits[i];uint64_t q=uint64_t(val/b);remainder=uint64_t(val%b);quotient.digits.push_back(q);}std::ranges::reverse(quotient.digits);quotient.negative=(neg_a^neg_b);quotient.normalize();bigint<base>rem{int64_t(remainder)};rem.negative=neg_a;return std::pair{quotient,rem};}template<base_v base>auto divmod(bigint<base>const&a,bigint<base>const&b){if(a<b){return std::pair{bigint<base>(0),a};}if(size(b.digits)==1){int64_t b_val=b.negative?-int64_t(b.digits[0]):int64_t(b.digits[0]);return divmod_fast(a,b_val);}auto A=decimal<base>(a);auto B=decimal<base>(b);auto d=(A*B.inv(A.magnitude()-B.magnitude()+1)).round();auto r=a-d*b;if(r>=b){d+=1;r-=b;}if(r<bigint<base>(0)){d-=1;r+=b;}return std::pair{d,r};}}
#endif
#line 1 "cp-algo/math/decimal.hpp"
#line 1 "cp-algo/math/bigint.hpp"
#line 1 "cp-algo/util/big_alloc.hpp"
#include <set>
#include <map>
#include <deque>
#include <stack>
#include <queue>
#include <vector>
#include <string>
#include <cstddef>
#include <iostream>
#include <generator>
#include <forward_list>
#if defined(__linux__) || defined(__unix__) || (defined(__APPLE__) && defined(__MACH__))
# define CP_ALGO_USE_MMAP 1
# include <sys/mman.h>
#else
# define CP_ALGO_USE_MMAP 0
#endif
namespace cp_algo{template<typename T,size_t Align=32>class big_alloc{static_assert(Align>=alignof(void*),"Align must be at least pointer-size");static_assert(std::popcount(Align)==1,"Align must be a power of two");public:using value_type=T;template<class U>struct rebind{using other=big_alloc<U,Align>;};constexpr bool operator==(const big_alloc&)const=default;constexpr bool operator!=(const big_alloc&)const=default;big_alloc()noexcept=default;template<typename U,std::size_t A>big_alloc(const big_alloc<U,A>&)noexcept{}[[nodiscard]]T*allocate(std::size_t n){std::size_t padded=round_up(n*sizeof(T));std::size_t align=std::max<std::size_t>(alignof(T),Align);
#if CP_ALGO_USE_MMAP
if(padded>=MEGABYTE){void*raw=mmap(nullptr,padded,PROT_READ|PROT_WRITE,MAP_PRIVATE|MAP_ANONYMOUS,-1,0);madvise(raw,padded,MADV_HUGEPAGE);madvise(raw,padded,MADV_POPULATE_WRITE);return static_cast<T*>(raw);}
#endif
return static_cast<T*>(::operator new(padded,std::align_val_t(align)));}void deallocate(T*p,std::size_t n)noexcept{if(!p)return;std::size_t padded=round_up(n*sizeof(T));std::size_t align=std::max<std::size_t>(alignof(T),Align);
#if CP_ALGO_USE_MMAP
if(padded>=MEGABYTE){munmap(p,padded);return;}
#endif
::operator delete(p,padded,std::align_val_t(align));}private:static constexpr std::size_t MEGABYTE=1<<20;static constexpr std::size_t round_up(std::size_t x)noexcept{return(x+Align-1)/Align*Align;}};template<typename T>using big_vector=std::vector<T,big_alloc<T>>;template<typename T>using big_basic_string=std::basic_string<T,std::char_traits<T>,big_alloc<T>>;template<typename T>using big_deque=std::deque<T,big_alloc<T>>;template<typename T>using big_stack=std::stack<T,big_deque<T>>;template<typename T>using big_queue=std::queue<T,big_deque<T>>;template<typename T>using big_priority_queue=std::priority_queue<T,big_vector<T>>;template<typename T>using big_forward_list=std::forward_list<T,big_alloc<T>>;using big_string=big_basic_string<char>;template<typename Key,typename Value,typename Compare=std::less<Key>>using big_map=std::map<Key,Value,Compare,big_alloc<std::pair<const Key,Value>>>;template<typename T,typename Compare=std::less<T>>using big_multiset=std::multiset<T,Compare,big_alloc<T>>;template<typename T,typename Compare=std::less<T>>using big_set=std::set<T,Compare,big_alloc<T>>;template<typename Ref,typename V=void>using big_generator=std::generator<Ref,V,big_alloc<std::byte>>;}namespace std::ranges{template<typename Ref,typename V>elements_of(cp_algo::big_generator<Ref,V>&&)->elements_of<cp_algo::big_generator<Ref,V>&&,cp_algo::big_alloc<std::byte>>;}
#line 1 "cp-algo/math/fft_simple.hpp"
#line 1 "cp-algo/random/rng.hpp"
#include <chrono>
#include <random>
namespace cp_algo::random{std::mt19937_64 gen(std::chrono::steady_clock::now().time_since_epoch().count());uint64_t rng(){return gen();}}
#line 1 "cp-algo/math/common.hpp"
#include <functional>
#include <cstdint>
#include <cassert>
namespace cp_algo::math{
#ifdef CP_ALGO_MAXN
const int maxn=CP_ALGO_MAXN;
#else
const int maxn=1<<19;
#endif
const int magic=64;auto bpow(auto const&x,auto n,auto const&one,auto op){if(n==0){return one;}else{auto t=bpow(x,n/2,one,op);t=op(t,t);if(n%2){t=op(t,x);}return t;}}auto bpow(auto x,auto n,auto ans){return bpow(x,n,ans,std::multiplies{});}template<typename T>T bpow(T const&x,auto n){return bpow(x,n,T(1));}inline constexpr auto inv2(auto x){assert(x%2);std::make_unsigned_t<decltype(x)>y=1;while(y*x!=1){y*=2-x*y;}return y;}}
#line 1 "cp-algo/math/cvector.hpp"
#line 1 "cp-algo/util/simd.hpp"
#include <experimental/simd>
#line 6 "cp-algo/util/simd.hpp"
#include <memory>
#if defined(__x86_64__) && !defined(CP_ALGO_DISABLE_AVX2)
#define CP_ALGO_SIMD_AVX2_TARGET _Pragma("GCC target(\"avx2\")")
#else
#define CP_ALGO_SIMD_AVX2_TARGET
#endif
#define CP_ALGO_SIMD_PRAGMA_PUSH \
_Pragma("GCC push_options")\CP_ALGO_SIMD_AVX2_TARGETCP_ALGO_SIMD_PRAGMA_PUSHnamespace cp_algo{template<typename T,size_t len>using simd[[gnu::vector_size(len*sizeof(T))]]=T;using u64x8=simd<uint64_t,8>;using u32x16=simd<uint32_t,16>;using i64x4=simd<int64_t,4>;using u64x4=simd<uint64_t,4>;using u32x8=simd<uint32_t,8>;using u16x16=simd<uint16_t,16>;using i32x4=simd<int32_t,4>;using u32x4=simd<uint32_t,4>;using u16x8=simd<uint16_t,8>;using u16x4=simd<uint16_t,4>;using i16x4=simd<int16_t,4>;using u8x32=simd<uint8_t,32>;using u8x8=simd<uint8_t,8>;using u8x4=simd<uint8_t,4>;using dx4=simd<double,4>;inline dx4 abs(dx4 a){return dx4{std::abs(a[0]),std::abs(a[1]),std::abs(a[2]),std::abs(a[3])};}static constexpr dx4 magic=dx4()+(3ULL<<51);inline i64x4 lround(dx4 x){return i64x4(x+magic)-i64x4(magic);}inline dx4 to_double(i64x4 x){return dx4(x+i64x4(magic))-magic;}inline dx4 round(dx4 a){return dx4{std::nearbyint(a[0]),std::nearbyint(a[1]),std::nearbyint(a[2]),std::nearbyint(a[3])};}inline u64x4 low32(u64x4 x){return x&uint32_t(-1);}inline auto swap_bytes(auto x){return decltype(x)(__builtin_shufflevector(u32x8(x),u32x8(x),1,0,3,2,5,4,7,6));}inline u64x4 montgomery_reduce(u64x4 x,uint32_t mod,uint32_t imod){
#ifdef __AVX2__
auto x_ninv=u64x4(_mm256_mul_epu32(__m256i(x),__m256i()+imod));x+=u64x4(_mm256_mul_epu32(__m256i(x_ninv),__m256i()+mod));
#else
auto x_ninv=u64x4(u32x8(low32(x))*imod);x+=x_ninv*uint64_t(mod);
#endif
return swap_bytes(x);}inline u64x4 montgomery_mul(u64x4 x,u64x4 y,uint32_t mod,uint32_t imod){
#ifdef __AVX2__
return montgomery_reduce(u64x4(_mm256_mul_epu32(__m256i(x),__m256i(y))),mod,imod);
#else
return montgomery_reduce(x*y,mod,imod);
#endif
}inline u32x8 montgomery_mul(u32x8 x,u32x8 y,uint32_t mod,uint32_t imod){return u32x8(montgomery_mul(u64x4(x),u64x4(y),mod,imod))|u32x8(swap_bytes(montgomery_mul(u64x4(swap_bytes(x)),u64x4(swap_bytes(y)),mod,imod)));}inline dx4 rotate_right(dx4 x){static constexpr u64x4 shuffler={3,0,1,2};return __builtin_shuffle(x,shuffler);}template<std::size_t Align=32>inline bool is_aligned(const auto*p)noexcept{return(reinterpret_cast<std::uintptr_t>(p)%Align)==0;}template<class Target>inline Target&vector_cast(auto&&p){return*reinterpret_cast<Target*>(std::assume_aligned<alignof(Target)>(&p));}}
#pragma GCC pop_options
#line 1 "cp-algo/util/complex.hpp"
#line 4 "cp-algo/util/complex.hpp"
#include <cmath>
#include <type_traits>
#line 7 "cp-algo/util/complex.hpp"
CP_ALGO_SIMD_PRAGMA_PUSHnamespace cp_algo{template<typename T>struct complex{using value_type=T;T x,y;inline constexpr complex():x(),y(){}inline constexpr complex(T const&x):x(x),y(){}inline constexpr complex(T const&x,T const&y):x(x),y(y){}inline complex&operator*=(T const&t){x*=t;y*=t;return*this;}inline complex&operator/=(T const&t){x/=t;y/=t;return*this;}inline complex operator*(T const&t)const{return complex(*this)*=t;}inline complex operator/(T const&t)const{return complex(*this)/=t;}inline complex&operator+=(complex const&t){x+=t.x;y+=t.y;return*this;}inline complex&operator-=(complex const&t){x-=t.x;y-=t.y;return*this;}inline complex operator*(complex const&t)const{return{x*t.x-y*t.y,x*t.y+y*t.x};}inline complex operator/(complex const&t)const{return*this*t.conj()/t.norm();}inline complex operator+(complex const&t)const{return complex(*this)+=t;}inline complex operator-(complex const&t)const{return complex(*this)-=t;}inline complex&operator*=(complex const&t){return*this=*this*t;}inline complex&operator/=(complex const&t){return*this=*this/t;}inline complex operator-()const{return{-x,-y};}inline complex conj()const{return{x,-y};}inline T norm()const{return x*x+y*y;}inline T abs()const{return std::sqrt(norm());}inline T const real()const{return x;}inline T const imag()const{return y;}inline T&real(){return x;}inline T&imag(){return y;}inline static constexpr complex polar(T r,T theta){return{T(r*cos(theta)),T(r*sin(theta))};}inline auto operator<=>(complex const&t)const=default;};template<typename T>inline complex<T>conj(complex<T>const&x){return x.conj();}template<typename T>inline T norm(complex<T>const&x){return x.norm();}template<typename T>inline T abs(complex<T>const&x){return x.abs();}template<typename T>inline T&real(complex<T>&x){return x.real();}template<typename T>inline T&imag(complex<T>&x){return x.imag();}template<typename T>inline T const real(complex<T>const&x){return x.real();}template<typename T>inline T const imag(complex<T>const&x){return x.imag();}template<typename T>inline constexpr complex<T>polar(T r,T theta){return complex<T>::polar(r,theta);}template<typename T>inline std::ostream&operator<<(std::ostream&out,complex<T>const&x){return out<<x.real()<<' '<<x.imag();}}
#pragma GCC pop_options
#line 1 "cp-algo/util/checkpoint.hpp"
#line 8 "cp-algo/util/checkpoint.hpp"
namespace cp_algo{
#ifdef CP_ALGO_CHECKPOINT
big_map<big_string,double>checkpoints;double last;
#endif
template<bool final=false>void checkpoint([[maybe_unused]]auto const&_msg){
#ifdef CP_ALGO_CHECKPOINT
big_string msg=_msg;double now=(double)clock()/CLOCKS_PER_SEC;double delta=now-last;last=now;if(msg.size()&&!final){checkpoints[msg]+=delta;}if(final){for(auto const&[key,value]:checkpoints){std::cerr<<key<<": "<<value*1000<<" ms\n";}std::cerr<<"Total: "<<now*1000<<" ms\n";}
#endif
}template<bool final=false>void checkpoint(){checkpoint<final>("");}}
#line 7 "cp-algo/math/cvector.hpp"
#include <ranges>
#include <bit>
CP_ALGO_SIMD_PRAGMA_PUSHnamespace stdx=std::experimental;namespace cp_algo::math::fft{static constexpr size_t flen=4;using ftype=double;using vftype=dx4;using point=complex<ftype>;using vpoint=complex<vftype>;static constexpr vftype vz={};vpoint vi(vpoint const&r){return{-imag(r),real(r)};}struct cvector{big_vector<vpoint>r;cvector(size_t n){n=std::max(flen,std::bit_ceil(n));r.resize(n/flen);checkpoint("cvector create");}vpoint&at(size_t k){return r[k/flen];}vpoint at(size_t k)const{return r[k/flen];}template<class pt=point>inline void set(size_t k,pt const&t){if constexpr(std::is_same_v<pt,point>){real(r[k/flen])[k%flen]=real(t);imag(r[k/flen])[k%flen]=imag(t);}else{at(k)=t;}}template<class pt=point>inline pt get(size_t k)const{if constexpr(std::is_same_v<pt,point>){return{real(r[k/flen])[k%flen],imag(r[k/flen])[k%flen]};}else{return at(k);}}size_t size()const{return flen*r.size();}static constexpr size_t eval_arg(size_t n){if(n<pre_evals){return eval_args[n];}else{return eval_arg(n/2)|(n&1)<<(std::bit_width(n)-1);}}static constexpr point eval_point(size_t n){if(n%2){return-eval_point(n-1);}else if(n%4){return eval_point(n-2)*point(0,1);}else if(n/4<pre_evals){return evalp[n/4];}else{return polar<ftype>(1.,std::numbers::pi/(ftype)std::bit_floor(n)*(ftype)eval_arg(n));}}static constexpr std::array<point,32>roots=[](){std::array<point,32>res;for(size_t i=2;i<32;i++){res[i]=polar<ftype>(1.,std::numbers::pi/(1ull<<(i-2)));}return res;}();static constexpr point root(size_t n){return roots[std::bit_width(n)];}template<int step>static void exec_on_eval(size_t n,size_t k,auto&&callback){callback(k,root(4*step*n)*eval_point(step*k));}template<int step>static void exec_on_evals(size_t n,auto&&callback){point factor=root(4*step*n);for(size_t i=0;i<n;i++){callback(i,factor*eval_point(step*i));}}static void do_dot_iter(point rt,vpoint&Bv,vpoint const&Av,vpoint&res){res+=Av*Bv;real(Bv)=rotate_right(real(Bv));imag(Bv)=rotate_right(imag(Bv));auto x=real(Bv)[0],y=imag(Bv)[0];real(Bv)[0]=x*real(rt)-y*imag(rt);imag(Bv)[0]=x*imag(rt)+y*real(rt);}void dot(cvector const&t){size_t n=this->size();exec_on_evals<1>(n/flen,[&](size_t k,point rt)__attribute__((always_inline)){k*=flen;auto[Ax,Ay]=at(k);auto Bv=t.at(k);vpoint res=vz;for(size_t i=0;i<flen;i++){vpoint Av=vpoint(vz+Ax[i],vz+Ay[i]);do_dot_iter(rt,Bv,Av,res);}set(k,res);});checkpoint("dot");}template<bool partial=true>void ifft(){size_t n=size();if constexpr(!partial){point pi(0,1);exec_on_evals<4>(n/4,[&](size_t k,point rt)__attribute__((always_inline)){k*=4;point v1=conj(rt);point v2=v1*v1;point v3=v1*v2;auto A=get(k);auto B=get(k+1);auto C=get(k+2);auto D=get(k+3);set(k,(A+B)+(C+D));set(k+2,((A+B)-(C+D))*v2);set(k+1,((A-B)-pi*(C-D))*v1);set(k+3,((A-B)+pi*(C-D))*v3);});}bool parity=std::countr_zero(n)%2;if(parity){exec_on_evals<2>(n/(2*flen),[&](size_t k,point rt)__attribute__((always_inline)){k*=2*flen;vpoint cvrt={vz+real(rt),vz-imag(rt)};auto B=at(k)-at(k+flen);at(k)+=at(k+flen);at(k+flen)=B*cvrt;});}for(size_t leaf=3*flen;leaf<n;leaf+=4*flen){size_t level=std::countr_one(leaf+3);for(size_t lvl=4+parity;lvl<=level;lvl+=2){size_t i=(1<<lvl)/4;exec_on_eval<4>(n>>lvl,leaf>>lvl,[&](size_t k,point rt)__attribute__((always_inline)){k<<=lvl;vpoint v1={vz+real(rt),vz-imag(rt)};vpoint v2=v1*v1;vpoint v3=v1*v2;for(size_t j=k;j<k+i;j+=flen){auto A=at(j);auto B=at(j+i);auto C=at(j+2*i);auto D=at(j+3*i);at(j)=((A+B)+(C+D));at(j+2*i)=((A+B)-(C+D))*v2;at(j+i)=((A-B)-vi(C-D))*v1;at(j+3*i)=((A-B)+vi(C-D))*v3;}});}}checkpoint("ifft");for(size_t k=0;k<n;k+=flen){if constexpr(partial){set(k,get<vpoint>(k)/=vz+ftype(n/flen));}else{set(k,get<vpoint>(k)/=vz+ftype(n));}}}template<bool partial=true>void fft(){size_t n=size();bool parity=std::countr_zero(n)%2;for(size_t leaf=0;leaf<n;leaf+=4*flen){size_t level=std::countr_zero(n+leaf);level-=level%2!=parity;for(size_t lvl=level;lvl>=4;lvl-=2){size_t i=(1<<lvl)/4;exec_on_eval<4>(n>>lvl,leaf>>lvl,[&](size_t k,point rt)__attribute__((always_inline)){k<<=lvl;vpoint v1={vz+real(rt),vz+imag(rt)};vpoint v2=v1*v1;vpoint v3=v1*v2;for(size_t j=k;j<k+i;j+=flen){auto A=at(j);auto B=at(j+i)*v1;auto C=at(j+2*i)*v2;auto D=at(j+3*i)*v3;at(j)=(A+C)+(B+D);at(j+i)=(A+C)-(B+D);at(j+2*i)=(A-C)+vi(B-D);at(j+3*i)=(A-C)-vi(B-D);}});}}if(parity){exec_on_evals<2>(n/(2*flen),[&](size_t k,point rt)__attribute__((always_inline)){k*=2*flen;vpoint vrt={vz+real(rt),vz+imag(rt)};auto t=at(k+flen)*vrt;at(k+flen)=at(k)-t;at(k)+=t;});}if constexpr(!partial){point pi(0,1);exec_on_evals<4>(n/4,[&](size_t k,point rt)__attribute__((always_inline)){k*=4;point v1=rt;point v2=v1*v1;point v3=v1*v2;auto A=get(k);auto B=get(k+1)*v1;auto C=get(k+2)*v2;auto D=get(k+3)*v3;set(k,(A+C)+(B+D));set(k+1,(A+C)-(B+D));set(k+2,(A-C)+pi*(B-D));set(k+3,(A-C)-pi*(B-D));});}checkpoint("fft");}static constexpr size_t pre_evals=1<<16;static const std::array<size_t,pre_evals>eval_args;static const std::array<point,pre_evals>evalp;};const std::array<size_t,cvector::pre_evals>cvector::eval_args=[](){std::array<size_t,pre_evals>res={};for(size_t i=1;i<pre_evals;i++){res[i]=res[i>>1]|(i&1)<<(std::bit_width(i)-1);}return res;}();const std::array<point,cvector::pre_evals>cvector::evalp=[](){std::array<point,pre_evals>res={};res[0]=1;for(size_t n=1;n<pre_evals;n++){res[n]=polar<ftype>(1.,std::numbers::pi*ftype(eval_args[n])/ftype(4*std::bit_floor(n)));}return res;}();}
#pragma GCC pop_options
#line 6 "cp-algo/math/fft_simple.hpp"
CP_ALGO_SIMD_PRAGMA_PUSHnamespace cp_algo::math::fft{struct dft_simple{cp_algo::math::fft::cvector cv;dft_simple(auto const&a,size_t n):cv(n){for(size_t i=0;i<std::min(std::size(a),n);i++){real(cv.at(i))[i%4]=ftype(a[i]);imag(cv.at(i))[i%4]=ftype(i+n<std::size(a)?a[i+n]:0);}checkpoint("dft64 init");cv.fft();}void dot(dft_simple const&t){cv.dot(t.cv);}void recover_mod(auto&res,size_t k){cv.ifft();size_t n=cv.size();for(size_t i=0;i<std::min(k,n);i++){res[i]=llround(real(cv.get(i)));}for(size_t i=n;i<k;i++){res[i]=llround(imag(cv.get(i-n)));}cp_algo::checkpoint("recover mod");}};void conv_simple(auto&a,auto const&b){if(empty(a)||empty(b)){a.clear();return;}size_t n=a.size(),m=b.size();size_t N=std::max(flen,std::bit_ceil(n+m-1)/2);dft_simple A(a,N),B(b,N);A.dot(B);a.resize(n+m-1);A.recover_mod(a,n+m-1);}}
#pragma GCC pop_options
#line 5 "cp-algo/math/bigint.hpp"
#include <bits/stdc++.h>
namespace cp_algo::math{enum base_v{x10=uint64_t(1e16),x16=uint64_t(1ull<<60)};template<base_v base=x10>struct bigint{static constexpr uint64_t Base=uint64_t(base);static constexpr uint16_t digit_length=base==x10?16:15;static constexpr uint16_t sub_base=base==x10?10:16;static constexpr uint32_t meta_base=base==x10?uint32_t(1e4):uint32_t(1<<15);big_basic_string<uint64_t>digits;bool negative;auto operator<=>(bigint const&other)const{if(digits.empty()&&other.digits.empty()){return std::strong_ordering::equal;}if(digits.empty()){return other.negative?std::strong_ordering::greater:std::strong_ordering::less;}if(other.digits.empty()){return negative?std::strong_ordering::less:std::strong_ordering::greater;}if(negative!=other.negative){return negative?std::strong_ordering::less:std::strong_ordering::greater;}if(digits.size()!=other.digits.size()){auto size_cmp=digits.size()<=>other.digits.size();return negative?0<=>size_cmp:size_cmp;}for(auto i=ssize(digits)-1;i>=0;i--){auto digit_cmp=digits[i]<=>other.digits[i];if(digit_cmp!=std::strong_ordering::equal){return negative?0<=>digit_cmp:digit_cmp;}}return std::strong_ordering::equal;}bigint(){}bigint(big_basic_string<uint64_t>d,bool neg):digits(std::move(d)),negative(neg){normalize();}bigint&pad_inplace(size_t to_add){digits.insert(0,to_add,0);return normalize();}bigint&drop_inplace(size_t to_drop){digits.erase(0,std::min(to_drop,size(digits)));return normalize();}bigint&take_inplace(size_t to_keep){digits.erase(std::min(to_keep,size(digits)),std::string::npos);return normalize();}bigint&top_inplace(size_t to_keep){if(to_keep>=size(digits)){return pad_inplace(to_keep-size(digits));}else{return drop_inplace(size(digits)-to_keep);}}bigint pad(size_t to_add)const{return bigint{big_basic_string<uint64_t>(to_add,0)+digits,negative}.normalize();}bigint drop(size_t to_drop)const{return bigint{digits.substr(std::min(to_drop,size(digits))),negative}.normalize();}bigint take(size_t to_keep)const{return bigint{digits.substr(0,std::min(to_keep,size(digits))),negative}.normalize();}bigint top(size_t to_keep)const{if(to_keep>=size(digits)){return pad(to_keep-size(digits));}else{return drop(size(digits)-to_keep);}}bigint&normalize(){while(!empty(digits)&&digits.back()==0){digits.pop_back();}if(digits.empty()){negative=false;}return*this;}bigint&negate(){negative^=1;return*this;}bigint operator-(){return bigint(*this).negate();}bigint&operator-=(const bigint&other){if(negative!=other.negative){return(negate()+=other).negate().normalize();}digits.resize(std::max(size(digits),size(other.digits)));bool carry=false;auto d_ptr=std::assume_aligned<32>(digits.data());auto o_ptr=std::assume_aligned<32>(other.digits.data());size_t N=size(other.digits);size_t i=0;for(;i<N;i++){d_ptr[i]-=o_ptr[i]+carry;carry=d_ptr[i]>=base;d_ptr[i]+=carry?uint64_t(base):0;}if(carry){N=size(digits);for(;i<N&&d_ptr[i]==0;i++){d_ptr[i]=base-1;}if(i<N){d_ptr[i]--;}else{for(i=0;i<N;i++){d_ptr[i]=base-d_ptr[i]-1;}bool carry=true;for(i=0;i<N&&carry;i++){d_ptr[i]++;carry=d_ptr[i]>=base;d_ptr[i]-=carry*base;}negate();}}return normalize();}bigint&operator+=(const bigint&other){if(negative!=other.negative){return(negate()-=other).negate().normalize();}digits.resize(std::max(size(digits),size(other.digits)));bool carry=false;auto d_ptr=std::assume_aligned<32>(digits.data());auto o_ptr=std::assume_aligned<32>(other.digits.data());size_t N=size(other.digits);size_t i=0;for(;i<N;i++){d_ptr[i]+=o_ptr[i]+carry;carry=d_ptr[i]>=base;d_ptr[i]-=carry?uint64_t(base):0;}if(carry){N=size(digits);for(;i<N&&d_ptr[i]==uint64_t(base)-1;i++){d_ptr[i]=0;}if(i<N){d_ptr[i]++;}else{digits.push_back(1);}}return*this;}bigint(int64_t x){negative=x<0;x=negative?-x:x;digits=x?big_basic_string<uint64_t>{uint64_t(x)}:big_basic_string<uint64_t>{};}bigint(std::span<char>s):negative(false){if(size(s)<digit_length){int64_t val=0;std::from_chars(s.data(),s.data()+size(s),val,sub_base);*this=bigint(val);return;}if(!empty(s)&&s[0]=='-'){negative=true;s=s.subspan(1);}size_t len=size(s);assert(len>0);size_t num_digits=(len+digit_length-1)/digit_length;digits.resize(num_digits);size_t i=len;for(size_t j=0;j<num_digits-1;j++){std::from_chars(s.data()+i-digit_length,s.data()+i,digits[j],sub_base);i-=digit_length;}std::from_chars(s.data(),s.data()+i,digits.back(),sub_base);normalize();}bigint operator+(const bigint&other)const{return bigint(*this)+=other;}bigint operator-(const bigint&other)const{return bigint(*this)-=other;}void to_metabase(){auto N=ssize(digits);digits.resize(4*N);for(auto i=N-1;i>=0;i--){uint64_t val=digits[i];digits[4*i]=val%meta_base;val/=meta_base;digits[4*i+1]=val%meta_base;val/=meta_base;digits[4*i+2]=val%meta_base;val/=meta_base;digits[4*i+3]=val;}}void from_metabase(){auto N=(ssize(digits)+3)/4;digits.resize(4*N);uint64_t carry=0;for(int i=0;i<N;i++){__uint128_t val=digits[4*i+3];val=val*meta_base+digits[4*i+2];val=val*meta_base+digits[4*i+1];val=val*meta_base+digits[4*i];val+=carry;digits[i]=uint64_t(val%base);carry=uint64_t(val/base);}digits.resize(N);while(carry){digits.push_back(carry%base);carry/=base;}}bigint&operator*=(int64_t other){if(other<0){negative^=1;other=-other;}if(other==0){return*this=bigint(0);}else if(other==1){return*this;}uint64_t carry=0;for(auto&d:digits){__uint128_t val=__uint128_t(d)*other+carry;d=uint64_t(val%base);carry=uint64_t(val/base);}if(carry){digits.push_back(carry%base);carry/=base;}return*this;}bigint operator*(int64_t other)const{return bigint(*this)*=other;}friend bigint operator*(int64_t lhs,const bigint&rhs){return bigint(rhs)*=lhs;}bigint&mul_inplace(auto&&other){negative^=other.negative;auto n=size(digits),m=size(other.digits);if(n<m){std::swap(n,m);std::swap(digits,other.digits);}if(m<=1){return*this*=int64_t(m==0?0:other.digits[0]);}constexpr size_t SMALL_M_THRESHOLD=32;if(m<=SMALL_M_THRESHOLD){big_basic_string<uint64_t>res;res.assign(n+m,0);for(size_t i=0;i<n;i++){if(digits[i]==0)continue;uint64_t carry=0;for(size_t j=0;j<m;j++){__uint128_t cur=res[i+j]+(__uint128_t)digits[i]*other.digits[j]+carry;res[i+j]=uint64_t(cur%Base);carry=uint64_t(cur/Base);}size_t k=i+m;if(carry){uint64_t cur=res[k]+carry;res[k]=cur%Base;carry=cur/Base;k++;}}digits=std::move(res);return normalize();}to_metabase();other.to_metabase();fft::conv_simple(digits,other.digits);from_metabase();return normalize();}bigint&operator*=(bigint const&other){return mul_inplace(bigint(other));}bigint operator*(const bigint&other)const{return bigint(*this).mul_inplace(bigint(other));}};template<base_v base>decltype(std::cin)&operator>>(decltype(std::cin)&in,cp_algo::math::bigint<base>&x){std::string s;in>>s;x={s};return in;}template<base_v base,bool fill=true>auto&print_digit(auto&out,uint64_t d){char buf[16];auto[ptr,ec]=std::to_chars(buf,buf+sizeof(buf),d,bigint<base>::sub_base);if constexpr(base==x16){std::ranges::transform(buf,buf,toupper);}auto len=ptr-buf;if constexpr(fill){out<<std::string(bigint<base>::digit_length-len,'0');}return out<<std::string_view(buf,len);}template<bool fill_all=false,base_v base>auto&print_bigint(auto&out,cp_algo::math::bigint<base>const&x){if(x.negative){out<<'-';}if(empty(x.digits)){return print_digit<base,fill_all>(out,0);}print_digit<base,fill_all>(out,x.digits.back());for(auto d:x.digits|std::views::reverse|std::views::drop(1)){print_digit<base,true>(out,d);}return out;}template<base_v base>decltype(std::cout)&operator<<(decltype(std::cout)&out,cp_algo::math::bigint<base>const&x){return print_bigint(out,x);}}
#line 5 "cp-algo/math/decimal.hpp"
namespace cp_algo::math{template<base_v base=x10>struct decimal{bigint<base>value;int64_t scale;decimal(int64_t v=0,int64_t s=0):value(bigint<base>(v)),scale(s){}decimal(bigint<base>v,int64_t s=0):value(v),scale(s){}decimal&operator*=(const decimal&other){value*=other.value;scale+=other.scale;return*this;}decimal&operator+=(decimal const&other){if(scale<other.scale){value+=other.value.pad(other.scale-scale);}else{value.pad_inplace(scale-other.scale);value+=other.value;scale=other.scale;}return*this;}decimal&operator-=(decimal const&other){if(scale<other.scale){value-=other.value.pad(other.scale-scale);}else{value.pad_inplace(scale-other.scale);value-=other.value;scale=other.scale;}return*this;}decimal operator*(const decimal&other)const{return decimal(*this)*=other;}decimal operator+(const decimal&other)const{return decimal(*this)+=other;}decimal operator-(const decimal&other)const{return decimal(*this)-=other;}auto split()const{auto int_part=scale>=-ssize(value.digits)?value.top(ssize(value.digits)+scale):bigint<base>(0);auto frac_part=*this-decimal(int_part);return std::pair{int_part,frac_part};}void print(){auto[int_part,frac_part]=split();print_bigint(std::cout,int_part);if(frac_part.value!=bigint<base>(0)){std::cout<<'.';std::cout<<std::string(bigint<base>::digit_length*(-frac_part.magnitude()),'0');frac_part.value.negative=false;print_bigint<true>(std::cout,frac_part.value);}std::cout<<std::endl;}bigint<base>trunc()const{if(scale>=0){return value.pad(scale);}else if(-scale>=ssize(value.digits)){return 0;}else{return value.top(ssize(value.digits)+scale);}}bigint<base>round()const{if(scale>=0){return value.pad(scale);}else if(-scale>ssize(value.digits)){return 0;}else{auto res=value.top(ssize(value.digits)+scale);if(value.digits[-scale-1]*2>=bigint<base>::Base){res+=1;}return res;}}decimal trunc(size_t digits)const{digits=std::min(digits,size(value.digits));return decimal(value.top(digits),scale+ssize(value.digits)-digits);}auto magnitude()const{static constexpr int64_t inf=1e18;if(value.digits.empty())return-inf;return ssize(value.digits)+scale;}decimal inv(int64_t precision){assert(precision>=0);int64_t lead=llround((double)bigint<base>::Base/(double)value.digits.back());decimal d(bigint<base>(lead),-ssize(value.digits));size_t cur=2;decimal amend=decimal(1)-trunc(cur)*d;while(-amend.magnitude()<precision){d+=d*amend;cur=2*(1-amend.magnitude());d=d.trunc(cur);amend=decimal(1)-trunc(cur)*d;}return d;}};template<base_v base>auto divmod_fast(bigint<base>const&a,int64_t b){if(b==0){assert(false&&"Division by zero");}bool neg_a=a.negative;bool neg_b=b<0;b=std::abs(b);bigint<base>quotient;uint64_t remainder=0;auto n=ssize(a.digits);for(auto i=n-1;i>=0;i--){__uint128_t val=(__uint128_t)remainder*bigint<base>::Base+a.digits[i];uint64_t q=uint64_t(val/b);remainder=uint64_t(val%b);quotient.digits.push_back(q);}std::ranges::reverse(quotient.digits);quotient.negative=(neg_a^neg_b);quotient.normalize();bigint<base>rem{int64_t(remainder)};rem.negative=neg_a;return std::pair{quotient,rem};}template<base_v base>auto divmod(bigint<base>const&a,bigint<base>const&b){if(a<b){return std::pair{bigint<base>(0),a};}if(size(b.digits)==1){int64_t b_val=b.negative?-int64_t(b.digits[0]):int64_t(b.digits[0]);return divmod_fast(a,b_val);}auto A=decimal<base>(a);auto B=decimal<base>(b);auto d=(A*B.inv(A.magnitude()-B.magnitude()+1)).round();auto r=a-d*b;if(r>=b){d+=1;r-=b;}if(r<bigint<base>(0)){d-=1;r+=b;}return std::pair{d,r};}}