# pragma GCC target("avx2")
# pragma GCC optimize("O3")
# pragma GCC optimize("unroll-loops")


#ifndef CPPLIB_SRC_ALGORITHM_MATH_INTEGER_FAST_BIG_INTEGER_HPP_INCLUDED
#define CPPLIB_SRC_ALGORITHM_MATH_INTEGER_FAST_BIG_INTEGER_HPP_INCLUDED

#include<bits/stdc++.h>
namespace fbi{ template<class Integer> inline constexpr bool ni = std::is_integral_v<std::remove_cv_t<Integer>> || std::same_as<std::remove_cv_t<Integer>, __int128_t> || std::same_as<std::remove_cv_t<Integer>, __uint128_t>; template<class Integer> concept NI = ni<Integer>; template<class Integer> struct MU{ using type = std::make_unsigned_t<Integer>; }; template<> struct MU<__int128_t>{using type = __uint128_t;}; template<> struct MU<__uint128_t>{using type = __uint128_t;}; template<class Integer> using UT = typename MU<std::remove_cv_t<Integer>>::type; constexpr uint64_t pw( uint64_t value, uint64_t exponent, uint64_t modulus ){ uint64_t result = 1; while(exponent != 0){ if((exponent & 1U) != 0) result = result * value % modulus; value = value * value % modulus; exponent >>= 1; } return result; }
template<uint32_t M> class MI{ static_assert((M & 1U) != 0); static_assert(M < (uint32_t{1} << 30)); static consteval uint32_t mi(){ uint32_t inverse = M; for(int it = 0; it < 5; ++it){ inverse *= 2U - M * inverse; } return 0U - inverse; }
static constexpr uint32_t iv = mi(); static constexpr uint32_t r2 = static_cast<uint32_t>( (uint64_t{0} - M) % M ); static constexpr uint32_t m2 = M * 2U; uint32_t value_ = 0; static constexpr uint32_t reduce(uint64_t value) noexcept{ return static_cast<uint32_t>( (value + static_cast<uint64_t>( static_cast<uint32_t>(value) * iv ) * M) >> 32 ); } public: constexpr MI() = default; constexpr MI(uint64_t value) noexcept: value_(reduce(value % M * r2)){} constexpr MI& operator+=(const MI& o) noexcept{ value_ += o.value_; if(value_ >= m2) value_ -= m2; return *this; } constexpr MI& operator-=(const MI& o) noexcept{ if(value_ < o.value_) value_ += m2; value_ -= o.value_; return *this; } constexpr MI& operator*=(const MI& o) noexcept{ value_ = reduce(static_cast<uint64_t>(value_) * o.value_); return *this; }
friend constexpr MI operator+( MI x, const MI& y ) noexcept{return x += y;}
friend constexpr MI operator-( MI x, const MI& y ) noexcept{return x -= y;}
friend constexpr MI operator*( MI x, const MI& y ) noexcept{return x *= y;} constexpr MI pow(uint64_t exponent) const noexcept{ MI result(1), factor = *this; while(exponent != 0){ if((exponent & 1U) != 0) result *= factor; factor *= factor; exponent >>= 1; } return result; } constexpr uint32_t value() const noexcept{ uint32_t result = reduce(value_); if(result >= M) result -= M; return result; } }; template<uint32_t M, uint32_t G> void ntt( std::vector<MI<M>>& values, bool inverse ){ using Z = MI<M>; const size_t size = values.size(); if(size <= 1) return; for(size_t i = 1, reversed = 0; i < size; ++i){ size_t bit = size >> 1; while((reversed & bit) != 0){ reversed ^= bit; bit >>= 1; } reversed ^= bit; if(i < reversed) std::swap(values[i], values[reversed]); } const unsigned levels = std::countr_zero(size); size_t fs = 2; if((levels & 1U) != 0){ for(size_t block = 0; block < size; block += 2){ const Z second = values[block + 1]; values[block + 1] = values[block] - second; values[block] += second; } fs = 4; } for(; fs < size; fs <<= 2){ const size_t block_size = fs << 1; const size_t quarter = fs >> 1; Z root = Z(G).pow( (M - 1) / static_cast<uint64_t>(block_size) ); if(inverse) root = root.pow(M - 2); const Z imaginary = root.pow(quarter); for(size_t block = 0; block < size; block += block_size){ Z factor(1); for(size_t offset = 0; offset < quarter; ++offset){ const Z sqd_factor = factor * factor; const Z first = values[block + offset]; const Z second = values[block + quarter + offset] * sqd_factor; const Z third = values[block + fs + offset]; const Z fourth = values[block + fs + quarter + offset] * sqd_factor; const Z low_sum = first + second; const Z low_difference = first - second; const Z high_sum = (third + fourth) * factor; const Z high_difference = (third - fourth) * imaginary * factor; values[block + offset] = low_sum + high_sum; values[block + fs + offset] = low_sum - high_sum; values[block + quarter + offset] = low_difference + high_difference; values[block + fs + quarter + offset] = low_difference - high_difference; factor *= root; } } } if(inverse){ const Z inverse_size = Z(size).pow(M - 2); for(Z& value: values) value *= inverse_size; } }
template<uint32_t M, uint32_t G> std::vector<uint32_t> cv( std::span<const uint32_t> x, std::span<const uint32_t> y, size_t ts, bool sq ){ using Z = MI<M>; const size_t rs = x.size() + y.size() - 1; std::vector<Z> lv(ts); for(size_t i = 0; i < x.size(); ++i){ lv[i] = Z(x[i]); } ntt<M, G>(lv, false); if(sq){ for(Z& value: lv) value *= value; }else{ std::vector<Z> rv(ts); for(size_t i = 0; i < y.size(); ++i){ rv[i] = Z(y[i]); } ntt<M, G>(rv, false); for(size_t i = 0; i < ts; ++i){ lv[i] *= rv[i]; } } ntt<M, G>(lv, true); std::vector<uint32_t> result(rs); for(size_t i = 0; i < rs; ++i){ result[i] = lv[i].value(); } return result; } } class BigInteger{ using L = uint32_t; using W = __uint128_t; static constexpr L B = 1'000'000'000U; static constexpr size_t D = 9; static constexpr size_t I = 4; static constexpr size_t MT = 128; static constexpr size_t VT = 64; static constexpr size_t NT = size_t{1} << 24; class LS{ union Values{ std::array<L, I> a; std::vector<L> v; Values(){} ~Values(){} } x_; static constexpr std::uint8_t heap_mode = 255; std::uint8_t n_ = 0; bool heap() const noexcept{return n_ == heap_mode;} public: bool s = false; LS(){new(&x_.a) decltype(x_.a){};} ~LS(){if(heap()) x_.v.~vector();} LS(const LS& x): n_(x.n_), s(x.s){ if(heap()) new(&x_.v) decltype(x_.v)(x.x_.v); else new(&x_.a) decltype(x_.a)(x.x_.a); } LS(LS&& x) noexcept: n_(x.n_), s(x.s){ if(heap()) new(&x_.v) decltype(x_.v)(std::move(x.x_.v)); else new(&x_.a) decltype(x_.a)(x.x_.a); } LS& operator=(const LS& x){ if(this != &x){LS y(x); *this = std::move(y);} return *this; } LS& operator=(LS&& x) noexcept{ if(this != &x){ this->~LS(); new(this) LS(std::move(x)); } return *this; } size_t size() const noexcept{return heap() ? x_.v.size() : n_;} bool empty() const noexcept{return size() == 0;} L* data() noexcept{return heap() ? x_.v.data() : x_.a.data();} const L* data() const noexcept{ return heap() ? x_.v.data() : x_.a.data(); } L& operator[](size_t i) noexcept{return data()[i];} const L& operator[](size_t i) const noexcept{return data()[i];} L& back() noexcept{return data()[size() - 1];} const L& back() const noexcept{return data()[size() - 1];} L* begin() noexcept{return data();} L* end() noexcept{return data() + size();} const L* begin() const noexcept{return data();} const L* end() const noexcept{return data() + size();} void clear() noexcept{ if(heap()) x_.v.clear(); else n_ = 0; } void resize(size_t n, L z = 0){ if(heap()){ x_.v.resize(n, z); return; } if(n <= I){ if(n > n_) std::fill(x_.a.begin() + n_, x_.a.begin() + n, z); n_ = static_cast<std::uint8_t>(n); return; } std::vector<L> v(n, z); std::copy_n(x_.a.begin(), n_, v.begin()); new(&x_.v) decltype(x_.v)(std::move(v)); n_ = heap_mode; } void push_back(L z){ const size_t n = size(); resize(n + 1); (*this)[n] = z; } void pop_back() noexcept{ if(heap()) x_.v.pop_back(); else --n_; } void assign(std::vector<L>&& v){ if(v.size() <= I){ if(heap()){ x_.v.~vector(); new(&x_.a) decltype(x_.a){}; } n_ = static_cast<std::uint8_t>(v.size()); std::copy(v.begin(), v.end(), x_.a.begin()); }else if(heap()) x_.v = std::move(v); else{ new(&x_.v) decltype(x_.v)(std::move(v)); n_ = heap_mode; } } std::span<const L> span() const noexcept{return {data(), size()};} }; LS a_; static void shr(std::vector<L>& value) noexcept{ while(!value.empty() && value.back() == 0) value.pop_back(); } void norm() noexcept{ while(!a_.empty() && a_.back() == 0) a_.pop_back(); if(a_.empty()) a_.s = false; }
static int cmp( std::span<const L> x, std::span<const L> y ) noexcept{ while(!x.empty() && x.back() == 0){ x = x.first(x.size() - 1); } while(!y.empty() && y.back() == 0){ y = y.first(y.size() - 1); } if(x.size() != y.size()) return x.size() < y.size() ? -1 : 1; for(size_t i = x.size(); i-- > 0;){ if(x[i] != y[i]) return x[i] < y[i] ? -1 : 1; } return 0; }
static std::vector<L> add( std::span<const L> x, std::span<const L> y ){ const size_t size = std::max(x.size(), y.size()); std::vector<L> result(size + 1, 0); uint64_t cy = 0; for(size_t i = 0; i < size; ++i){ const uint64_t sum = cy + (i < x.size() ? x[i] : 0) + (i < y.size() ? y[i] : 0); if(sum >= B){ result[i] = static_cast<L>(sum - B); cy = 1; }else{ result[i] = static_cast<L>(sum); cy = 0; } } result[size] = static_cast<L>(cy); shr(result); return result; }
static std::vector<L> sub( std::span<const L> larger, std::span<const L> smaller ){ std::vector<L> result(larger.size(), 0); uint64_t bw = 0; for(size_t i = 0; i < larger.size(); ++i){ const uint64_t sb = (i < smaller.size() ? smaller[i] : 0) + bw; if(larger[i] < sb){ result[i] = static_cast<L>( static_cast<uint64_t>(larger[i]) + B - sb ); bw = 1; }else{ result[i] = static_cast<L>(larger[i] - sb); bw = 0; } } shr(result); return result; }
static std::vector<L> mul1( std::span<const L> value, L k ){ if(value.empty() || k == 0) return {}; if(k == 1) return {value.begin(), value.end()}; std::vector<L> result(value.size() + 1, 0); uint64_t cy = 0; for(size_t i = 0; i < value.size(); ++i){ const uint64_t p = static_cast<uint64_t>(value[i]) * k + cy; result[i] = static_cast<L>(p % B); cy = p / B; } result[value.size()] = static_cast<L>(cy); shr(result); return result; }
static uint64_t recip(L d) noexcept{ const uint64_t q = std::numeric_limits<uint64_t>::max() / d; return q + ( std::numeric_limits<uint64_t>::max() - q * d == d - 1 ); }
static L step( uint64_t x, L d, uint64_t rp, L& o ) noexcept{ uint64_t q = static_cast<uint64_t>( static_cast<W>(x) * rp >> 64 ); uint64_t r = x - q * d; if(r >= d){ ++q; r -= d; } o = static_cast<L>(r); return static_cast<L>(q); }
static std::pair<std::vector<L>, L> div1( std::span<const L> value, L y ){ if(y == 1) return {{value.begin(), value.end()}, 0}; std::vector<L> q(value.size(), 0); const uint64_t rp = recip(y); L r = 0; for(size_t i = value.size(); i-- > 0;){ const uint64_t c = static_cast<uint64_t>(r) * B + value[i]; q[i] = step( c, y, rp, r ); } shr(q); return {std::move(q), r}; } L mod1(L y) const noexcept{ if(y == 1) return 0; const uint64_t rp = recip(y); L r = 0; for(size_t i = a_.size(); i-- > 0;){ (void)step( static_cast<uint64_t>(r) * B + a_[i], y, rp, r ); } return r; } L divi(L y) noexcept{ if(y == 1) return 0; const uint64_t rp = recip(y); L r = 0; for(size_t i = a_.size(); i-- > 0;){ const uint64_t c = static_cast<uint64_t>(r) * B + a_[i]; a_[i] = step( c, y, rp, r ); } norm(); return r; } void muli(L k){ if(k == 0 || is_zero()){ a_.clear(); a_.s = false; return; } if(k == 1) return; const size_t size = a_.size(); a_.resize(size + 1, 0); uint64_t cy = 0; for(size_t i = 0; i < size; ++i){ const uint64_t p = static_cast<uint64_t>(a_[i]) * k + cy; a_[i] = static_cast<L>(p % B); cy = p / B; } a_[size] = static_cast<L>(cy); norm(); } void addi(const BigInteger& o){ const size_t ls = a_.size(); const size_t size = std::max(ls, o.a_.size()); a_.resize(size + 1, 0); uint64_t cy = 0; for(size_t i = 0; i < size; ++i){ const uint64_t sum = cy + (i < ls ? a_[i] : 0) + (i < o.a_.size() ? o.a_[i] : 0); if(sum >= B){ a_[i] = static_cast<L>(sum - B); cy = 1; }else{ a_[i] = static_cast<L>(sum); cy = 0; } } a_[size] = static_cast<L>(cy); norm(); } void subi(const BigInteger& o){ uint64_t bw = 0; for(size_t i = 0; i < a_.size(); ++i){ const uint64_t sb = bw + (i < o.a_.size() ? o.a_[i] : 0); if(a_[i] < sb){ a_[i] = static_cast<L>( static_cast<uint64_t>(a_[i]) + B - sb ); bw = 1; }else{ a_[i] = static_cast<L>(a_[i] - sb); bw = 0; } } norm(); }
static std::vector<L> smul( std::span<const L> x, std::span<const L> y ){ if(x.size() < y.size()) std::swap(x, y); std::vector<L> result(x.size() + y.size(), 0); for(size_t j = 0; j < y.size(); ++j){ uint64_t cy = 0; for(size_t i = 0; i < x.size(); ++i){ const size_t k = i + j; const uint64_t c = static_cast<uint64_t>(x[i]) * y[j] + result[k] + cy; result[k] = static_cast<L>(c % B); cy = c / B; } result[j + x.size()] = static_cast<L>(cy); } shr(result); return result; }
static std::vector<L> nmul( std::span<const L> x, std::span<const L> y, bool sq ){ constexpr uint32_t modulus0 = 167'772'161U; constexpr uint32_t modulus1 = 469'762'049U; constexpr uint32_t modulus2 = 754'974'721U; if(x.size() > (std::numeric_limits<size_t>::max)() - y.size() + 1){ throw std::length_error("BigInteger multiplication is too large"); } const size_t rs = x.size() + y.size() - 1; if(rs > NT){ throw std::length_error(std::string{}); } const size_t ts = std::bit_ceil(rs); auto residue0 = fbi::cv< modulus0, 3U >(x, y, ts, sq); auto residue1 = fbi::cv< modulus1, 3U >(x, y, ts, sq); auto residue2 = fbi::cv< modulus2, 11U >(x, y, ts, sq); constexpr uint64_t modulus01 = static_cast<uint64_t>(modulus0) * modulus1; constexpr uint64_t inverse0_mod1 = fbi::pw(modulus0, modulus1 - 2, modulus1); constexpr uint64_t inverse01_mod2 = fbi::pw( modulus01 % modulus2, modulus2 - 2, modulus2 ); std::vector<L> result; result.reserve(rs + 3); W cy = 0; for(size_t i = 0; i < rs; ++i){ const uint64_t first = residue0[i]; const uint64_t second_delta = (residue1[i] + modulus1 - first % modulus1) % modulus1; const uint64_t second = second_delta * inverse0_mod1 % modulus1; const uint64_t first_two = first + modulus0 * second; const uint64_t third_delta = (residue2[i] + modulus2 - first_two % modulus2) % modulus2; const uint64_t third = third_delta * inverse01_mod2 % modulus2; const W coefficient = static_cast<W>(first_two) + static_cast<W>(modulus01) * third; const W c = coefficient + cy; result.push_back(static_cast<L>(c % B)); cy = c / B; } while(cy != 0){ result.push_back(static_cast<L>(cy % B)); cy /= B; } shr(result); return result; }
static std::vector<L> mul( std::span<const L> x, std::span<const L> y, bool sq = false ){ if(x.empty() || y.empty()) return {}; if(x.size() == 1) return mul1(y, x.front()); if(y.size() == 1) return mul1(x, y.front()); if(std::min(x.size(), y.size()) <= MT){ return smul(x, y); } return nmul(x, y, sq); }
static void dec(std::vector<L>& value){ size_t i = 0; while(value[i] == 0){ value[i] = B - 1; ++i; } --value[i]; shr(value); }
static void inc(std::vector<L>& value){ size_t i = 0; while(i < value.size() && value[i] == B - 1){ value[i] = 0; ++i; } if(i == value.size()) value.push_back(1); else ++value[i]; }
static std::pair<std::vector<L>, std::vector<L>> ldiv( std::span<const L> x, std::span<const L> y, bool wq = true, bool wr = true ){ const int ord = cmp(x, y); if(ord < 0) return {{}, wr ? std::vector<L>(x.begin(), x.end()) : std::vector<L>{}}; if(ord == 0) return {wq ? std::vector<L>{1} : std::vector<L>{}, {}}; if(y.size() == 1){ auto [q, r] = div1(x, y[0]); return { wq ? std::move(q) : std::vector<L>{}, !wr || r == 0 ? std::vector<L>{} : std::vector<L>{r} }; } const L nm = static_cast<L>( B / (static_cast<uint64_t>(y.back()) + 1) ); std::vector<L> normd_x = mul1(x, nm); std::vector<L> normd_y = mul1(y, nm); normd_x.resize(x.size() + 1, 0); const size_t ds = normd_y.size(); const size_t qs = x.size() - y.size() + 1; std::vector<L> q( wq ? qs : 0, 0 ); const uint64_t hd = normd_y.back(); const uint64_t nd = normd_y[ds - 2]; for(size_t pos = qs; pos-- > 0;){ const uint64_t num = static_cast<uint64_t>( normd_x[pos + ds] ) * B + normd_x[pos + ds - 1]; uint64_t e = num / hd; uint64_t r = num % hd; if(e >= B){ e = B - 1; r = num - e * hd; } while(r < B && e * nd > r * B + normd_x[pos + ds - 2]){ --e; r += hd; } uint64_t cy = 0; uint64_t bw = 0; for(size_t i = 0; i < ds; ++i){ const uint64_t p = e * normd_y[i] + cy; cy = p / B; const uint64_t sb = p % B + bw; L& c = normd_x[pos + i]; if(c < sb){ c = static_cast<L>( static_cast<uint64_t>(c) + B - sb ); bw = 1; }else{ c = static_cast<L>(c - sb); bw = 0; } } L& high = normd_x[pos + ds]; const uint64_t hs = cy + bw; const bool overed = high < hs; if(overed){ high = static_cast<L>( static_cast<uint64_t>(high) + B - hs ); --e; uint64_t ac = 0; for(size_t i = 0; i < ds; ++i){ const uint64_t sum = static_cast<uint64_t>( normd_x[pos + i] ) + normd_y[i] + ac; if(sum >= B){ normd_x[pos + i] = static_cast<L>(sum - B); ac = 1; }else{ normd_x[pos + i] = static_cast<L>(sum); ac = 0; } } const uint64_t ch = high + ac; high = static_cast<L>( ch >= B ? ch - B : ch ); }else{ high = static_cast<L>(high - hs); } if(wq) q[pos] = static_cast<L>(e); } shr(q); if(!wr) return {std::move(q), {}}; std::vector<L> r( normd_x.begin(), normd_x.begin() + static_cast<ptrdiff_t>(ds) ); shr(r); if(nm != 1){ auto division = div1(r, nm); r = std::move(division.first); } return {std::move(q), std::move(r)}; }
static std::vector<L> inv( const std::vector<L>& value, size_t pr ){ size_t cp = pr; while(cp > VT){ cp = (cp + 1) / 2; } std::vector<L> num( value.size() + cp + 1, 0 ); num.back() = 1; std::vector<L> ap = ldiv(num, value, true, false).first; while(cp < pr){ std::vector<L> sq = mul( ap, ap, true ); sq.insert(sq.begin(), 0); const size_t ws = cp * 2 + 1; std::vector<L> hv(ws, 0); const size_t copied = std::min(value.size(), ws); std::copy( value.end() - static_cast<ptrdiff_t>(copied), value.end(), hv.end() - static_cast<ptrdiff_t>(copied) ); std::vector<L> co = mul( sq, hv ); if(co.size() <= ws){ co.clear(); }else{ co.erase( co.begin(), co.begin() + static_cast<ptrdiff_t>(ws) ); } std::vector<L> db = add( ap, ap ); std::vector<L> sc(cp + 1, 0); sc.insert(sc.end(), db.begin(), db.end()); ap = sub(sc, co); if(!ap.empty()) ap.erase(ap.begin()); shr(ap); cp *= 2; } if(cp > pr){ const size_t excess = cp - pr; ap.erase( ap.begin(), ap.begin() + static_cast<ptrdiff_t>(excess) ); } shr(ap); return ap; }
static std::pair<std::vector<L>, std::vector<L>> ndiv( std::span<const L> x, std::span<const L> y, bool wq = true, bool wr = true ){ const L nm = static_cast<L>( B / (static_cast<uint64_t>(y.back()) + 1) ); std::vector<L> normd_x = mul1(x, nm); std::vector<L> normd_y = mul1(y, nm); const size_t pr = normd_x.size() - normd_y.size() + 2; std::vector<L> rp = inv( normd_y, pr ); std::vector<L> q = mul( normd_x, rp ); const size_t discarded = normd_y.size() + pr; if(q.size() <= discarded){ q.clear(); }else{ q.erase( q.begin(), q.begin() + static_cast<ptrdiff_t>(discarded) ); } shr(q); std::vector<L> p = mul( normd_y, q ); while(cmp(normd_x, p) < 0){ dec(q); p = sub(p, normd_y); } std::vector<L> r = sub( normd_x, p ); while(cmp(r, normd_y) >= 0){ inc(q); r = sub(r, normd_y); } if(nm != 1 && wr){ auto division = div1(r, nm); r = std::move(division.first); } if(!wq) q.clear(); if(!wr) r.clear(); shr(q); shr(r); return {std::move(q), std::move(r)}; }
static std::pair<std::vector<L>, std::vector<L>> div( std::span<const L> x, std::span<const L> y, bool wq = true, bool wr = true ){ const int ord = cmp(x, y); if(ord < 0) return {{}, wr ? std::vector<L>(x.begin(), x.end()) : std::vector<L>{}}; if(ord == 0) return {wq ? std::vector<L>{1} : std::vector<L>{}, {}}; if(y.size() == 1 || y.size() <= VT || x.size() - y.size() <= VT){ return ldiv( x, y, wq, wr ); } return ndiv( x, y, wq, wr ); }
template<fbi::NI Integer> void set(Integer value){ using Value = std::remove_cv_t<Integer>; a_.clear(); a_.s = false; if constexpr(std::same_as<Value, bool>){ if(value) a_.push_back(1); }else{ using U = fbi::UT<Value>; U m = static_cast<U>(value); if constexpr(std::numeric_limits<Value>::is_signed){ if(value < 0){ a_.s = true; m = U{0} - m; } } std::vector<L> limbs; while(m != 0){ limbs.push_back(static_cast<L>(m % B)); m /= B; } a_.assign(std::move(limbs)); } }
static BigInteger make( std::vector<L>&& m, bool negative = false ){ shr(m); BigInteger result; result.a_.s = negative && !m.empty(); result.a_.assign(std::move(m)); return result; } public: BigInteger() = default; template<fbi::NI Integer> BigInteger(Integer value){set(value);} explicit BigInteger(std::string_view decimal){assign(decimal);}
template<fbi::NI Integer> BigInteger& operator=(Integer value){ set(value); return *this; }
BigInteger& assign(std::string_view decimal){ if(decimal.empty()){ throw std::invalid_argument("empty BigInteger literal"); } bool negative = false; size_t begin = 0; if(decimal.front() == '+' || decimal.front() == '-'){ negative = decimal.front() == '-'; begin = 1; } if(begin == decimal.size()){ throw std::invalid_argument("BigInteger literal has no digits"); } for(size_t i = begin; i < decimal.size(); ++i){ if(decimal[i] < '0' || decimal[i] > '9'){ throw std::invalid_argument("invalid BigInteger decimal digit"); } } while(begin < decimal.size() && decimal[begin] == '0') ++begin; if(begin == decimal.size()){ a_.clear(); a_.s = false; return *this; } std::vector<L> m; m.reserve( (decimal.size() - begin + D - 1) / D ); size_t end = decimal.size(); while(end > begin){ const size_t block_begin = end - begin > D ? end - D : begin; L block = 0; for(size_t i = block_begin; i < end; ++i){ block = static_cast<L>( block * 10U + static_cast<unsigned>(decimal[i] - '0') ); } m.push_back(block); end = block_begin; } a_.assign(std::move(m)); a_.s = negative; return *this; } bool is_zero() const noexcept{return a_.empty();} bool is_negative() const noexcept{return a_.s;}
BigInteger absolute() const{ BigInteger result = *this; result.a_.s = false; return result; } std::string to_string() const{ if(is_zero()) return "0"; char hb[16]; const auto [highest_end, error] = std::to_chars( hb, hb + sizeof(hb), a_.back() ); if(error != std::errc{}){ throw std::runtime_error("BigInteger decimal conversion failed"); } const size_t highest_size = static_cast<size_t>( highest_end - hb ); std::string result; result.reserve( static_cast<size_t>(a_.s) + highest_size + (a_.size() - 1) * D ); if(a_.s) result.push_back('-'); result.append(hb, highest_end); for(size_t i = a_.size() - 1; i-- > 0;){ char block[D]; L value = a_[i]; for(size_t pos = D; pos-- > 0;){ block[pos] = static_cast<char>('0' + value % 10U); value /= 10U; } result.append(block, block + D); } return result; }
template<fbi::NI Integer> Integer checked_to() const{ using Value = std::remove_cv_t<Integer>; if constexpr(std::same_as<Value, bool>){ if(*this == 0) return false; if(*this == 1) return true; throw std::overflow_error( "BigInteger does not fit target integer type" ); }else{ using U = fbi::UT<Value>; U limit; if constexpr(std::numeric_limits<Value>::is_signed){ const U positive_limit = static_cast<U>( (std::numeric_limits<Value>::max)() ); limit = a_.s ? positive_limit + U{1} : positive_limit; }else{ if(a_.s){ throw std::overflow_error( "BigInteger does not fit target integer type" ); } limit = (std::numeric_limits<Value>::max)(); } U m = 0; for(size_t i = a_.size(); i-- > 0;){ const U digit = static_cast<U>(a_[i]); if(digit > limit || m > (limit - digit) / B){ throw std::overflow_error( "BigInteger does not fit target integer type" ); } m = m * B + digit; } if constexpr(!std::numeric_limits<Value>::is_signed){ return static_cast<Value>(m); }else{ if(!a_.s) return static_cast<Value>(m); const U minimum_m = static_cast<U>((std::numeric_limits<Value>::max)()) + U{1}; if(m == minimum_m){ return (std::numeric_limits<Value>::min)(); } return static_cast<Value>(-static_cast<Value>(m)); } } }
static std::pair<BigInteger, BigInteger> divmod( const BigInteger& x, const BigInteger& y ){ if(y.is_zero()){ throw std::domain_error("BigInteger division by zero"); } auto [q, r] = div( x.a_.span(), y.a_.span() ); return { make( std::move(q), x.a_.s != y.a_.s ), make(std::move(r), x.a_.s) }; }
BigInteger operator-() const{ BigInteger result = *this; if(!result.is_zero()) result.a_.s = !result.a_.s; return result; }
BigInteger& operator+=(const BigInteger& o){ if(o.is_zero()) return *this; if(is_zero()){ *this = o; return *this; } if(a_.s == o.a_.s){ addi(o); return *this; } const int ord = cmp( a_.span(), o.a_.span() ); if(ord == 0){ a_.clear(); a_.s = false; }else if(ord > 0){ subi(o); }else{ a_.assign(sub( o.a_.span(), a_.span() )); a_.s = o.a_.s; } return *this; }
BigInteger& operator-=(const BigInteger& o){ if(this == &o){ a_.clear(); a_.s = false; return *this; } if(o.is_zero()) return *this; if(is_zero()){ *this = -o; return *this; } if(a_.s != o.a_.s){ addi(o); return *this; } const int ord = cmp( a_.span(), o.a_.span() ); if(ord == 0){ a_.clear(); a_.s = false; }else if(ord > 0){ subi(o); }else{ a_.assign(sub( o.a_.span(), a_.span() )); a_.s = !o.a_.s; } return *this; }
BigInteger& operator*=(const BigInteger& o){ const bool sn = a_.s != o.a_.s; if(o.a_.size() == 1){ const L k = o.a_[0]; muli(k); a_.s = sn && !a_.empty(); return *this; } const bool sq = this == &o; std::vector<L> p = mul( a_.span(), o.a_.span(), sq ); a_.assign(std::move(p)); a_.s = sn && !a_.empty(); return *this; }
BigInteger& operator/=(const BigInteger& o){ if(o.is_zero()){ throw std::domain_error("BigInteger division by zero"); } if(is_zero()) return *this; const bool sn = a_.s != o.a_.s; if(o.a_.size() == 1){ divi(o.a_[0]); a_.s = sn && !a_.empty(); return *this; } auto q = div( a_.span(), o.a_.span(), true, false ).first; a_.assign(std::move(q)); a_.s = sn && !a_.empty(); return *this; }
BigInteger& operator%=(const BigInteger& o){ if(o.is_zero()){ throw std::domain_error("BigInteger division by zero"); } if(is_zero()) return *this; if(o.a_.size() == 1){ const L r = mod1(o.a_[0]); if(r == 0){ a_.clear(); a_.s = false; }else{ a_.resize(1); a_[0] = r; } return *this; } auto r = div( a_.span(), o.a_.span(), false, true ).second; a_.assign(std::move(r)); a_.s = a_.s && !a_.empty(); return *this; }
BigInteger& operator++(){ *this += 1; return *this; }
BigInteger operator++(int){ BigInteger result = *this; ++*this; return result; }
BigInteger& operator--(){ *this -= 1; return *this; }
BigInteger operator--(int){ BigInteger result = *this; --*this; return result; }
friend BigInteger abs(const BigInteger& value){return value.absolute();}
friend bool operator==( const BigInteger& x, const BigInteger& y ){ if(x.a_.s != y.a_.s || x.a_.size() != y.a_.size()){ return false; } return std::equal( x.a_.begin(), x.a_.end(), y.a_.begin() ); }
friend std::strong_ordering operator<=> ( const BigInteger& x, const BigInteger& y ){ if(x.a_.s != y.a_.s){ return x.a_.s ? std::strong_ordering::less : std::strong_ordering::greater; } const int ord = cmp( x.a_.span(), y.a_.span() ); if(ord == 0) return std::strong_ordering::equal; const bool less = x.a_.s ? ord > 0 : ord < 0; return less ? std::strong_ordering::less : std::strong_ordering::greater; }
friend BigInteger operator+( BigInteger x, const BigInteger& y ){return x += y;}
friend BigInteger operator-( BigInteger x, const BigInteger& y ){return x -= y;}
friend BigInteger operator*( const BigInteger& x, const BigInteger& y ){ return make( mul( x.a_.span(), y.a_.span(), &x == &y ), x.a_.s != y.a_.s ); }
friend BigInteger operator/( BigInteger x, const BigInteger& y ){return x /= y;}
friend BigInteger operator%( const BigInteger& x, const BigInteger& y ){ if(y.is_zero()){ throw std::domain_error("BigInteger division by zero"); } if(x.is_zero()) return {}; if(y.a_.size() == 1){ const L r = x.mod1(y.a_[0]); BigInteger result; if(r != 0){ result.a_.push_back(r); result.a_.s = x.a_.s; } return result; } return divmod(x, y).second; }
friend std::ostream& operator<<( std::ostream& stream, const BigInteger& value ){return stream << value.to_string();}
friend std::istream& operator>>( std::istream& stream, BigInteger& value ){ std::string text; if(!(stream >> text)) return stream; try{ BigInteger parsed(text); value = std::move(parsed); }catch(const std::invalid_argument&){ stream.setstate(std::ios::failbit); } return stream; } };
#endif  // CPPLIB_SRC_ALGORITHM_MATH_INTEGER_FAST_BIG_INTEGER_HPP_INCLUDED

#ifndef CPPLIB_SRC_STRUCTURE_IO_FASTIO_HPP_INCLUDED
#define CPPLIB_SRC_STRUCTURE_IO_FASTIO_HPP_INCLUDED

#include <limits>
#include <stdio.h>

#define read_size 1000000
#define write_size 1000000

struct fastio {
private:
  char read_data[read_size];
  int read_pos = 0;
  int read_len = 0;

  char write_data[write_size];
  int write_pos = 0;

  int getch() {
    if (read_pos == read_len) {
      read_len = fread(read_data, 1, read_size, stdin);
      read_pos = 0;
      if (read_len == 0) return EOF;
    }
    return read_data[read_pos++];
  }

  void ungetch() {
    if (read_pos > 0) read_pos--;
  }

  void readspeoln() {
    int c;
    while (true) {
      c = getch();
      if (c == EOF) return;
      if (c != ' ' && c != '\n' && c != '\r' && c != '\t') {
        ungetch();
        return;
      }
    }
  }

  void flush() {
    if (write_pos != 0) {
      fwrite(write_data, 1, write_pos, stdout);
      write_pos = 0;
    }
  }

public:
  fastio() {}

  ~fastio() {
    flush();
  }

  void readint(int &x) {
    readspeoln();

    int c = getch();
    bool negative = false;
    unsigned int value = 0;

    if (c == '-') {
      negative = true;
      c = getch();
    }

    while ('0' <= c && c <= '9') {
      value = value * 10U + static_cast<unsigned int>(c & 15);
      c = getch();
    }

    if (c != EOF) ungetch();

    if (negative && value == static_cast<unsigned int>(
        std::numeric_limits<int>::max()) + 1U) {
      x = std::numeric_limits<int>::min();
    } else {
      x = negative ? -static_cast<int>(value) : static_cast<int>(value);
    }
  }

  void readll(long long &x) {
    readspeoln();

    int c = getch();
    bool negative = false;
    unsigned long long value = 0;

    if (c == '-') {
      negative = true;
      c = getch();
    }

    while ('0' <= c && c <= '9') {
      value = value * 10ULL + static_cast<unsigned long long>(c & 15);
      c = getch();
    }

    if (c != EOF) ungetch();

    if (negative && value == static_cast<unsigned long long>(
        std::numeric_limits<long long>::max()) + 1ULL) {
      x = std::numeric_limits<long long>::min();
    } else {
      x = negative ? -static_cast<long long>(value) : static_cast<long long>(value);
    }
  }

  // [a, z]
  void readstr(char *s) {
    readspeoln();

    int c = getch();
    
    while (c != EOF && c != ' ' && c != '\n' && c != '\r' && c != '\t') {
      *s++ = (char)c;
      c = getch();
    }

    if (c != EOF) ungetch();

    *s = '\0';
  }

  void readstr(std::string &s) {
    readspeoln();

    s.clear();

    int c = getch();

    while (
      c != EOF &&
      c != ' ' &&
      c != '\n' &&
      c != '\r' &&
      c != '\t'
    ) {
      s.push_back(static_cast<char>(c));
      c = getch();
    }

    if (c != EOF) ungetch();
  }

  void write(char c) {
    if (write_pos == write_size) flush();
    write_data[write_pos++] = c;
  }

  void write(const char *s) {
    while (*s) write(*s++);
  }

  void writeint(int x) {
    if (x == 0) {
      write('0');
      return;
    }

    unsigned int value;
    if (x < 0) {
      write('-');
      value = 0U - static_cast<unsigned int>(x);
    } else {
      value = static_cast<unsigned int>(x);
    }

    char s[20];
    int n = 0;

    while (value > 0) {
      s[n++] = char('0' + value % 10U);
      value /= 10U;
    }

    while (n--) write(s[n]);
  }

  void writell(long long x) {
    if (x == 0) {
      write('0');
      return;
    }

    unsigned long long value;
    if (x < 0) {
      write('-');
      value = 0ULL - static_cast<unsigned long long>(x);
    } else {
      value = static_cast<unsigned long long>(x);
    }

    char s[30];
    int n = 0;

    while (value > 0) {
      s[n++] = char('0' + value % 10ULL);
      value /= 10ULL;
    }

    while (n--) write(s[n]);
  }
};

#endif  // CPPLIB_SRC_STRUCTURE_IO_FASTIO_HPP_INCLUDED

#include <bits/stdc++.h>
using namespace std;

using cpp_int = BigInteger;


namespace {

constexpr int MAX_AB = 300;

struct RawEdge {
    int u;
    int v;
    int numerator;
    int denominator;
};

struct Edge {
    int to;
    cpp_int weight;
};

struct State {
    cpp_int distance;
    int vertex;
};

struct StateGreater {
    bool operator()(const State& lhs, const State& rhs) const {
        if (lhs.distance != rhs.distance) {
            return lhs.distance > rhs.distance;
        }
        return lhs.vertex > rhs.vertex;
    }
};

}  // namespace

int main() {
    fastio io;
    int n, m;
    io.readint(n);
    io.readint(m);

    std::vector<RawEdge> rawEdges;
    rawEdges.reserve(m);

    // maximumExponent[p] is the largest exponent of p occurring in any b_i.
    std::vector<int> maximumExponent(MAX_AB + 1, 0);

    for (int i = 0; i < m; ++i) {
        int u, v, a, b;
        io.readint(u);
        io.readint(v);
        io.readint(a);
        io.readint(b);
        --u;
        --v;
        rawEdges.push_back({u, v, a, b});

        int value = b;
        for (int prime = 2; prime * prime <= value; ++prime) {
            if (value % prime != 0) {
                continue;
            }

            int exponent = 0;
            while (value % prime == 0) {
                value /= prime;
                ++exponent;
            }
            maximumExponent[prime] =
                std::max(maximumExponent[prime], exponent);
        }
        if (value > 1) {
            maximumExponent[value] =
                std::max(maximumExponent[value], 1);
        }
    }

    // Every b_i divides this common denominator.
    cpp_int commonDenominator = 1;
    std::vector<std::pair<int, int>> primePowers;

    for (int prime = 2; prime <= MAX_AB; ++prime) {
        if (maximumExponent[prime] == 0) {
            continue;
        }

        primePowers.push_back({prime, maximumExponent[prime]});
        for (int exponent = 0;
             exponent < maximumExponent[prime];
             ++exponent) {
            commonDenominator *= prime;
        }
    }

    std::vector<std::vector<Edge>> graph(n);

    for (const RawEdge& raw : rawEdges) {
        cpp_int scaledWeight = commonDenominator;
        scaledWeight /= raw.denominator;
        scaledWeight *= raw.numerator;

        // One copy and one move are necessary because the graph is undirected.
        graph[raw.u].push_back({raw.v, scaledWeight});
        graph[raw.v].push_back({raw.u, std::move(scaledWeight)});
    }

    std::vector<cpp_int> distance(n);
    std::vector<char> reached(n, false);

    std::priority_queue<State, std::vector<State>, StateGreater> queue;

    reached[0] = true;
    distance[0] = 0;
    queue.push({cpp_int(0), 0});

    while (!queue.empty()) {
        State current = queue.top();
        queue.pop();

        if (!reached[current.vertex] ||
            current.distance != distance[current.vertex]) {
            continue;
        }

        for (const Edge& edge : graph[current.vertex]) {
            // Evaluate the arbitrary-precision addition exactly once.
            cpp_int nextDistance = current.distance + edge.weight;

            if (!reached[edge.to] || nextDistance < distance[edge.to]) {
                reached[edge.to] = true;
                distance[edge.to] = nextDistance;
                queue.push({std::move(nextDistance), edge.to});
            }
        }
    }

    // The official output limit is approximately 8 MiB, so buffering the
    // complete answer avoids tens of thousands of formatted stream writes.
    std::string output;
    output.reserve(8U * 1024U * 1024U);

    for (int vertex = 1; vertex < n; ++vertex) {
        cpp_int numerator = distance[vertex];
        cpp_int denominator = commonDenominator;

        // The complete prime factorization of the denominator is already
        // known. Reduce by small primes instead of running a general cpp_int
        // Euclidean gcd for every output vertex.
        for (const auto& [prime, exponent] : primePowers) {
            for (int count = 0; count < exponent; ++count) {
                if (numerator % prime != 0) {
                    break;
                }
                numerator /= prime;
                denominator /= prime;
            }
        }

        output += numerator.to_string();
        output.push_back(' ');
        output += denominator.to_string();
        output.push_back('\n');
    }

    for(const char& i:output){
        io.write(i);
    }
    return 0;
}