#line 1 "verify/math/rational.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"
#line 1 "math/rational.hpp"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <compare>
#include <concepts>
#include <iostream>
#include <limits>
#include <sstream>
#include <string>
#include <type_traits>
#include <utility>
namespace m1une {
namespace math {
namespace rational_detail {
template <class T>
concept IntegerLike =
std::signed_integral<T> ||
(!std::integral<T> && std::copyable<T> && requires(T first, T second) {
T(0);
T(1);
{ -first } -> std::same_as<T>;
{ first + second } -> std::same_as<T>;
{ first - second } -> std::same_as<T>;
{ first * second } -> std::same_as<T>;
{ first / second } -> std::same_as<T>;
{ first % second } -> std::same_as<T>;
{ first += second } -> std::same_as<T&>;
{ first -= second } -> std::same_as<T&>;
{ first /= second } -> std::same_as<T&>;
{ first == second } -> std::convertible_to<bool>;
{ first < second } -> std::convertible_to<bool>;
});
} // namespace rational_detail
template <rational_detail::IntegerLike T = long long>
struct Rational {
static_assert(!std::signed_integral<T> || sizeof(T) <= sizeof(long long));
private:
static constexpr bool BUILTIN_INTEGER = std::signed_integral<T>;
using Wide = std::conditional_t<BUILTIN_INTEGER, __int128_t, T>;
using Magnitude = std::conditional_t<BUILTIN_INTEGER, __uint128_t, T>;
T _numerator;
T _denominator;
static constexpr Magnitude magnitude(Wide value) {
if constexpr (BUILTIN_INTEGER) {
if (value < 0) {
return static_cast<Magnitude>(-(value + 1)) + 1;
}
return static_cast<Magnitude>(value);
} else {
return value < 0 ? -value : value;
}
}
static constexpr Magnitude gcd(Magnitude first, Magnitude second) {
while (second != 0) {
Magnitude remainder = first % second;
first = second;
second = remainder;
}
return first;
}
static constexpr T narrow(Wide value) {
if constexpr (BUILTIN_INTEGER) {
assert(Wide(std::numeric_limits<T>::min()) <= value);
assert(value <= Wide(std::numeric_limits<T>::max()));
return static_cast<T>(value);
} else {
return value;
}
}
constexpr void assign_normalized(Wide numerator, Wide denominator) {
assert(denominator != 0);
if (numerator == 0) {
_numerator = 0;
_denominator = 1;
return;
}
Magnitude divisor = gcd(magnitude(numerator), magnitude(denominator));
numerator /= static_cast<Wide>(divisor);
denominator /= static_cast<Wide>(divisor);
if (denominator < 0) {
numerator = -numerator;
denominator = -denominator;
}
_numerator = narrow(numerator);
_denominator = narrow(denominator);
}
static constexpr Rational from_wide(Wide numerator, Wide denominator) {
Rational result;
result.assign_normalized(numerator, denominator);
return result;
}
static std::pair<long double, long long> decimal_scientific(const T& value) {
std::ostringstream output;
output << value;
const std::string text = output.str();
std::size_t begin = 0;
int sign = 1;
if (!text.empty() && (text[0] == '-' || text[0] == '+')) {
if (text[0] == '-') sign = -1;
begin = 1;
}
while (begin < text.size() && text[begin] == '0') ++begin;
if (begin == text.size()) return std::make_pair(0.0L, 0LL);
constexpr int DIGITS = std::numeric_limits<long double>::digits10 + 1;
const std::size_t used = std::min<std::size_t>(DIGITS, text.size() - begin);
long double significand = 0;
for (std::size_t i = 0; i < used; ++i) {
assert('0' <= text[begin + i] && text[begin + i] <= '9');
significand = significand * 10 + (text[begin + i] - '0');
}
for (std::size_t i = 1; i < used; ++i) significand /= 10;
const long long exponent = static_cast<long long>(text.size() - begin - 1);
return std::make_pair(sign * significand, exponent);
}
public:
constexpr Rational() : _numerator(0), _denominator(1) {}
constexpr Rational(T integer) : _numerator(integer), _denominator(1) {}
template <std::integral U>
requires std::constructible_from<T, U> &&
(!std::same_as<std::remove_cv_t<U>, T>)
constexpr Rational(U integer) : Rational(T(integer)) {}
constexpr Rational(T numerator, T denominator) {
assign_normalized(Wide(numerator), Wide(denominator));
}
constexpr T numerator() const {
return _numerator;
}
constexpr T denominator() const {
return _denominator;
}
constexpr bool is_integer() const {
return _denominator == 1;
}
constexpr int sign() const {
return (_numerator > 0) - (_numerator < 0);
}
constexpr Rational reciprocal() const {
assert(_numerator != 0);
return from_wide(Wide(_denominator), Wide(_numerator));
}
constexpr Rational abs() const {
return _numerator < 0 ? -*this : *this;
}
constexpr long double to_long_double() const
requires requires(const T& value) { static_cast<long double>(value); }
{
return static_cast<long double>(_numerator) / static_cast<long double>(_denominator);
}
long double to_long_double() const
requires(!requires(const T& value) { static_cast<long double>(value); })
{
const auto [numerator, numerator_exponent] = decimal_scientific(_numerator);
const auto [denominator, denominator_exponent] = decimal_scientific(_denominator);
return numerator / denominator *
std::pow(10.0L, numerator_exponent - denominator_exponent);
}
template <std::floating_point F>
explicit constexpr operator F() const
requires requires(const T& value) { static_cast<long double>(value); }
{
return static_cast<F>(to_long_double());
}
template <std::floating_point F>
explicit operator F() const
requires(!requires(const T& value) { static_cast<long double>(value); })
{
return static_cast<F>(to_long_double());
}
constexpr T trunc() const {
return _numerator / _denominator;
}
constexpr T floor() const {
T quotient = _numerator / _denominator;
if (_numerator < 0 && _numerator % _denominator != 0) quotient -= T(1);
return quotient;
}
constexpr T ceil() const {
T quotient = _numerator / _denominator;
if (0 < _numerator && _numerator % _denominator != 0) quotient += T(1);
return quotient;
}
constexpr Rational operator+() const {
return *this;
}
constexpr Rational operator-() const {
return from_wide(-Wide(_numerator), Wide(_denominator));
}
constexpr Rational& operator+=(const Rational& other) {
Magnitude common =
gcd(static_cast<Magnitude>(_denominator), static_cast<Magnitude>(other._denominator));
Wide left_scale = Wide(other._denominator) / static_cast<Wide>(common);
Wide right_scale = Wide(_denominator) / static_cast<Wide>(common);
Wide numerator =
Wide(_numerator) * left_scale + Wide(other._numerator) * right_scale;
// With both operands already reduced, every factor shared by the new
// numerator and denominator must divide `common`. Restricting the
// second gcd to that value avoids a full-size gcd against the product
// of both denominators, which is especially important for BigInt.
Magnitude reduction = common == Magnitude(1)
? Magnitude(1)
: gcd(magnitude(numerator), common);
if (reduction != Magnitude(1)) {
numerator /= static_cast<Wide>(reduction);
}
Wide remaining_denominator = Wide(other._denominator);
if (reduction != Magnitude(1)) {
remaining_denominator /= static_cast<Wide>(reduction);
}
_numerator = narrow(numerator);
_denominator = narrow(right_scale * remaining_denominator);
return *this;
}
constexpr Rational& operator-=(const Rational& other) {
return *this += -other;
}
constexpr Rational& operator*=(const Rational& other) {
Magnitude first_gcd = gcd(magnitude(Wide(_numerator)), static_cast<Magnitude>(other._denominator));
Magnitude second_gcd = gcd(magnitude(Wide(other._numerator)), static_cast<Magnitude>(_denominator));
assign_normalized((Wide(_numerator) / static_cast<Wide>(first_gcd)) *
(Wide(other._numerator) / static_cast<Wide>(second_gcd)),
(Wide(_denominator) / static_cast<Wide>(second_gcd)) *
(Wide(other._denominator) / static_cast<Wide>(first_gcd)));
return *this;
}
constexpr Rational& operator/=(const Rational& other) {
return *this *= other.reciprocal();
}
friend constexpr Rational operator+(Rational left, const Rational& right) {
return left += right;
}
friend constexpr Rational operator-(Rational left, const Rational& right) {
return left -= right;
}
friend constexpr Rational operator*(Rational left, const Rational& right) {
return left *= right;
}
friend constexpr Rational operator/(Rational left, const Rational& right) {
return left /= right;
}
friend constexpr bool operator==(const Rational& left, const Rational& right) {
return left._numerator == right._numerator && left._denominator == right._denominator;
}
friend constexpr std::strong_ordering operator<=>(const Rational& left, const Rational& right) {
Wide first = Wide(left._numerator) * Wide(right._denominator);
Wide second = Wide(right._numerator) * Wide(left._denominator);
if (first < second) return std::strong_ordering::less;
if (second < first) return std::strong_ordering::greater;
return std::strong_ordering::equal;
}
friend std::ostream& operator<<(std::ostream& output, const Rational& value) {
output << value._numerator;
if (value._denominator != 1) {
output << '/' << value._denominator;
}
return output;
}
friend std::istream& operator>>(std::istream& input, Rational& value) {
std::string token;
if (!(input >> token)) return input;
std::size_t slash = token.find('/');
if (slash != std::string::npos && token.find('/', slash + 1) != std::string::npos) {
input.setstate(std::ios::failbit);
return input;
}
T numerator = 0;
T denominator = 1;
std::istringstream numerator_input(token.substr(0, slash));
if (!(numerator_input >> numerator) || numerator_input.peek() != std::char_traits<char>::eof()) {
input.setstate(std::ios::failbit);
return input;
}
if (slash != std::string::npos) {
std::istringstream denominator_input(token.substr(slash + 1));
if (!(denominator_input >> denominator) ||
denominator_input.peek() != std::char_traits<char>::eof()) {
input.setstate(std::ios::failbit);
return input;
}
}
value = Rational(numerator, denominator);
return input;
}
};
template <rational_detail::IntegerLike T>
constexpr Rational<T> abs(const Rational<T>& value) {
return value.abs();
}
} // namespace math
} // namespace m1une
namespace std {
// Integer/rational common types already follow from implicit integer
// construction. Mixing a floating scalar explicitly chooses approximation.
template <m1une::math::rational_detail::IntegerLike T, floating_point F>
struct common_type<m1une::math::Rational<T>, F> {
using type = long double;
};
template <floating_point F, m1une::math::rational_detail::IntegerLike T>
struct common_type<F, m1une::math::Rational<T>> {
using type = long double;
};
} // namespace std
#line 1 "utilities/bigint.hpp"
#line 5 "utilities/bigint.hpp"
#include <array>
#include <bit>
#line 8 "utilities/bigint.hpp"
#include <charconv>
#line 10 "utilities/bigint.hpp"
#include <cstdint>
#include <cstring>
#line 13 "utilities/bigint.hpp"
#include <memory>
#include <numbers>
#include <stdexcept>
#line 18 "utilities/bigint.hpp"
#include <vector>
#if (defined(__GNUC__) || defined(__clang__)) && \
(defined(__x86_64__) || defined(__i386__))
#include <immintrin.h>
#define M1UNE_BIGINT_HAS_X86_SIMD 1
#endif
#line 1 "math/fps/convolution.hpp"
#line 9 "math/fps/convolution.hpp"
#include <new>
#line 13 "math/fps/convolution.hpp"
#if defined(__GNUC__) && !defined(__clang__) && \
(defined(__x86_64__) || defined(__i386__)) && \
!defined(M1UNE_FPS_DISABLE_X86_SIMD)
#include <immintrin.h>
#define M1UNE_FPS_HAS_X86_SIMD 1
#pragma GCC push_options
#pragma GCC target("avx2,bmi")
#endif
#line 1 "math/fps/internal/ntt998_faster.hpp"
#ifdef M1UNE_FPS_HAS_X86_SIMD
#line 9 "math/fps/internal/ntt998_faster.hpp"
#include <immintrin.h>
namespace m1une {
namespace fps {
namespace internal {
namespace fast998_v2 {
// Fixed-modulus AVX2 transform with an in-register degree-8 residue product.
using u32=unsigned;
using u64=unsigned long long;
using idt=std::size_t;
using I256=__m256i;
inline void store256(void*p,I256 x){
_mm256_store_si256((I256*)p,x);
}
inline I256 load256(const void*p){
return _mm256_load_si256((const I256*)p);
}
constexpr u32 shrk(u32 x,u32 M){
return std::min(x,x-M);
}
constexpr u32 dilt(u32 x,u32 M){
return std::min(x,x+M);
}
constexpr u32 reduce(u64 x,u32 niv,u32 M){
return (x+u64(u32(x)*niv)*M)>>32;
}
constexpr u32 mul(u32 x,u32 y,u32 niv,u32 M){
return reduce(u64(x)*y,niv,M);
}
constexpr u32 mul_s(u32 x,u32 y,u32 niv,u32 M){
return shrk(reduce(u64(x)*y,niv,M),M);
}
constexpr u32 qpw(u32 a,u32 b,u32 niv,u32 M,u32 r){
for(;b;b>>=1,a=mul(a,a,niv,M)){
if(b&1){
r=mul(r,a,niv,M);
}
}
return r;
}
constexpr u32 qpw_s(u32 a,u32 b,u32 niv,u32 M,u32 r){
return shrk(qpw(a,b,niv,M,r),M);
}
inline I256 shrk32(I256 x,I256 M){
return _mm256_min_epu32(x,_mm256_sub_epi32(x,M));
}
inline I256 dilt32(I256 x,I256 M){
return _mm256_min_epu32(x,_mm256_add_epi32(x,M));
}
inline I256 Ladd32(I256 x,I256 y,I256){
return _mm256_add_epi32(x,y);
}
inline I256 Lsub32(I256 x,I256 y,I256 M){
return _mm256_add_epi32(_mm256_sub_epi32(x,y),M);
}
inline I256 add32(I256 x,I256 y,I256 M){
return shrk32(_mm256_add_epi32(x,y),M);
}
inline I256 sub32(I256 x,I256 y,I256 M){
return dilt32(_mm256_sub_epi32(x,y),M);
}
template<int msk>inline I256 neg32_m(I256 x,I256 M){
return _mm256_blend_epi32(x,_mm256_sub_epi32(M,x),msk);
}
inline I256 reduce(I256 a,I256 b,I256 niv,I256 M){
I256 c=_mm256_mul_epu32(a,niv),d=_mm256_mul_epu32(b,niv);
c=_mm256_mul_epu32(c,M),d=_mm256_mul_epu32(d,M);
return _mm256_blend_epi32(_mm256_srli_epi64(_mm256_add_epi64(a,c),32),_mm256_add_epi64(b,d),0xaa);
}
inline I256 mul(I256 a,I256 b,I256 niv,I256 M){
return reduce(_mm256_mul_epu32(a,b),_mm256_mul_epu32(_mm256_srli_epi64(a,32),_mm256_srli_epi64(b,32)),niv,M);
}
inline I256 mul_s(I256 a,I256 b,I256 niv,I256 M){
return shrk32(mul(a,b,niv,M),M);
}
inline I256 mul_bsm(I256 a,I256 b,I256 niv,I256 M){
return reduce(_mm256_mul_epu32(a,b),_mm256_mul_epu32(_mm256_srli_epi64(a,32),b),niv,M);
}
inline I256 mul_bsmfxd(I256 a,I256 b,I256 bniv,I256 M){
I256 cc=_mm256_mul_epu32(a,bniv),dd=_mm256_mul_epu32(_mm256_srli_epi64(a,32),bniv);
I256 c=_mm256_mul_epu32(a,b),d=_mm256_mul_epu32(_mm256_srli_epi64(a,32),b);
cc=_mm256_mul_epu32(cc,M),dd=_mm256_mul_epu32(dd,M);
return _mm256_blend_epi32(_mm256_srli_epi64(_mm256_add_epi64(c,cc),32),_mm256_add_epi64(d,dd),0xaa);
}
inline I256 mul_bfxd(I256 a,I256 b,I256 bniv,I256 M){
I256 cc=_mm256_mul_epu32(a,bniv),dd=_mm256_mul_epu32(_mm256_srli_epi64(a,32),_mm256_srli_epi64(bniv,32));
I256 c=_mm256_mul_epu32(a,b),d=_mm256_mul_epu32(_mm256_srli_epi64(a,32),_mm256_srli_epi64(b,32));
cc=_mm256_mul_epu32(cc,M),dd=_mm256_mul_epu32(dd,M);
return _mm256_blend_epi32(_mm256_srli_epi64(_mm256_add_epi64(c,cc),32),_mm256_add_epi64(d,dd),0xaa);
}
inline I256 mul_upd_rt(I256 a,I256 bu,I256 M){
I256 cc=_mm256_mul_epu32(a,bu),c=_mm256_mul_epu32(a,_mm256_srli_epi64(bu,32));
cc=_mm256_mul_epu32(cc,M);
return shrk32(_mm256_srli_epi64(_mm256_add_epi64(c,cc),32),M);
}
constexpr auto _mxlg=26,_lg_itth=6;
constexpr auto _itth=idt(1)<<_lg_itth;
static_assert(_lg_itth%2==0);
struct FNTT32_info{
u32 mod,mod2,niv,one,r2,r3,img,imgniv,RT1[_mxlg];
alignas(32) std::array<u32,8> rt3[_mxlg-2],rt3i[_mxlg-2],bwbr,bwb,bwbi,rt4[_mxlg-3],rt4niv[_mxlg-3],rt4i[_mxlg-3],rt4iniv[_mxlg-3],pr2,pr4,pr2niv,pr4niv,pr2i,pr2iniv,pr4i,pr4iniv;
constexpr FNTT32_info(const u32 m):mod(m),mod2(m*2),niv([&]{u32 n=2+m;for(int i=0;i<4;++i){n*=2+m*n;}return n;}()),one((-m)%m),r2((-u64(m))%m),r3(mul_s(r2,r2,niv,m)),img{},imgniv{},RT1{},rt3{},rt3i{},bwbr{},bwb{},bwbi{},rt4{},rt4niv{},rt4i{},rt4iniv{},pr2{},pr4{},pr2niv{},pr4niv{},pr2i{},pr2iniv{},pr4i{},pr4iniv{}{
const int k=__builtin_ctz(m-1);
u32 _g=mul(3,r2,niv,mod);
for(;;++_g){
if(qpw_s(_g,mod>>1,niv,mod,one)!=one){
break;
}
}
_g=qpw(_g,mod>>k,niv,mod,one);
u32 rt1[_mxlg-1],rt1i[_mxlg-1];
rt1[k-2]=_g,rt1i[k-2]=qpw(_g,mod-2,niv,mod,one);
for(int i=k-2;i>0;--i){
rt1[i-1]=mul(rt1[i],rt1[i],niv,mod);
rt1i[i-1]=mul(rt1i[i],rt1i[i],niv,mod);
}
RT1[k-1]=qpw_s(_g,3,niv,mod,one);
for(int i=k-1;i>0;--i){
RT1[i-1]=mul_s(RT1[i],RT1[i],niv,mod);
}
img=rt1[0],imgniv=img*niv;
bwbr={one,0,one,0,one};
bwb={rt1[1],0,rt1[0],0,mod-mul_s(rt1[0],rt1[1],niv,mod)};
bwbi={rt1i[1],0,rt1i[0],0,mul_s(rt1i[0],rt1i[1],niv,mod)};
u32 pr=one,pri=one;
for(int i=0;i<k-2;++i){
const u32 r=mul_s(pr,rt1[i+1],niv,mod),ri=mul_s(pri,rt1i[i+1],niv,mod);
const u32 r2=mul_s(r,r,niv,mod),r2i=mul_s(ri,ri,niv,mod);
const u32 r3=mul_s(r,r2,niv,mod),r3i=mul_s(ri,r2i,niv,mod);
rt3[i]={r*niv,r,r2*niv,r2,r3*niv,r3};
rt3i[i]={ri*niv,ri,r2i*niv,r2i,r3i*niv,r3i};
pr=mul(pr,rt1i[i+1],niv,mod),pri=mul(pri,rt1[i+1],niv,mod);
}
pr=one,pri=one;
for(int i=0;i<k-3;++i){
const u32 r=mul_s(pr,rt1[i+2],niv,mod),ri=mul_s(pri,rt1i[i+2],niv,mod);
rt4[i][0]=rt4i[i][0]=one;
for(int j=1;j<8;++j){
rt4[i][j]=mul_s(rt4[i][j-1],r,niv,mod);
rt4i[i][j]=mul_s(rt4i[i][j-1],ri,niv,mod);
}
for(int j=0;j<8;++j){
rt4niv[i][j]=rt4[i][j]*niv;
rt4iniv[i][j]=rt4i[i][j]*niv;
}
pr=mul(pr,rt1i[i+2],niv,mod),pri=mul(pri,rt1[i+2],niv,mod);
}
pr2={one,one,one,img,one,one,one,img};
pr4={one,one,one,one,one,rt1[1],img,mul_s(img,rt1[1],niv,mod)};
const u32 nr2=mod-r2,imgr2=mul_s(img,r2,niv,mod);
pr2i={nr2,nr2,nr2,imgr2,nr2,nr2,nr2,imgr2};
pr4i={one,one,one,one,one,rt1i[1],rt1i[0],mul_s(rt1i[0],rt1i[1],niv,mod)};
for(int j=0;j<8;++j){
pr2niv[j]=pr2[j]*niv,pr4niv[j]=pr4[j]*niv;
pr2iniv[j]=pr2i[j]*niv,pr4iniv[j]=pr4i[j]*niv;
}
}
};
inline void vector_dif(I256*const f,const idt n,const FNTT32_info*info){
alignas(32) std::array<u32,8> st_1[_mxlg>>1];
const I256 Mod=_mm256_set1_epi32(info->mod),Mod2=_mm256_set1_epi32(info->mod2),Niv=_mm256_set1_epi32(info->niv);
const I256 Img=_mm256_set1_epi32(info->img),ImgNiv=_mm256_set1_epi32(info->imgniv),id=_mm256_setr_epi32(0,2,0,4,0,2,0,4);
const int lgn=__builtin_ctzll(n);
std::fill(st_1,st_1+(lgn>>1),info->bwb);
const idt nn=n>>(lgn&1),m=std::min(n,_itth),mm=std::min(nn,_itth);
// I256 rr=_mm256_set1_epi32(info->one);
if(nn!=n){
for(idt i=0;i<nn;++i){
auto const p0=f+i,p1=f+nn+i;
const auto f0=load256(p0),f1=load256(p1);
const auto g0=add32(f0,f1,Mod2),g1=Lsub32(f0,f1,Mod2);
store256(p0,g0),store256(p1,g1);
}
}
for(idt L=nn>>2;L>0;L>>=2){
for(idt i=0;i<L;++i){
auto const p0=f+i,p1=p0+L,p2=p1+L,p3=p2+L;
const auto f1=load256(p1),f3=load256(p3),f2=load256(p2),f0=load256(p0);
const auto g3=mul_bsmfxd(Lsub32(f1,f3,Mod2),Img,ImgNiv,Mod),g1=add32(f1,f3,Mod2);
const auto g0=add32(f0,f2,Mod2),g2=sub32(f0,f2,Mod2);
const auto h0=add32(g0,g1,Mod2),h1=Lsub32(g0,g1,Mod2);
const auto h2=Ladd32(g2,g3,Mod2),h3=Lsub32(g2,g3,Mod2);
store256(p0,h0),store256(p1,h1),store256(p2,h2),store256(p3,h3);
}
}
for(idt j=0;j<n;j+=m){
int t=((j==0)?std::min(_lg_itth,lgn):__builtin_ctzll(j))&-2,p=(t-2)>>1;
for(idt L=(idt(1)<<t)>>2;L>=_itth;L>>=2,t-=2,--p){
auto rt=load256(st_1+p);
const auto r1=_mm256_permutevar8x32_epi32(rt,id);
const auto r1Niv=_mm256_permutevar8x32_epi32(_mm256_mul_epu32(rt,Niv),id);
rt=mul_upd_rt(rt,load256(info->rt3+__builtin_ctzll(~j>>t)),Mod);
const auto r2=_mm256_shuffle_epi32(r1,_MM_PERM_BBBB),nr3=_mm256_shuffle_epi32(r1,_MM_PERM_DDDD);
const auto r2Niv=_mm256_shuffle_epi32(r1Niv,_MM_PERM_BBBB),nr3Niv=_mm256_shuffle_epi32(r1Niv,_MM_PERM_DDDD);
store256(st_1+p,rt);
for(idt i=0;i<L;++i){
auto const p0=f+i+j,p1=p0+L,p2=p1+L,p3=p2+L;
const auto f1=load256(p1),f3=load256(p3),f2=load256(p2),f0=load256(p0);
const auto g1=mul_bsmfxd(f1,r1,r1Niv,Mod),ng3=mul_bsmfxd(f3,nr3,nr3Niv,Mod);
const auto g2=mul_bsmfxd(f2,r2,r2Niv,Mod),g0=shrk32(f0,Mod2);
const auto h3=mul_bsmfxd(Ladd32(g1,ng3,Mod2),Img,ImgNiv,Mod),h1=sub32(g1,ng3,Mod2);
const auto h0=add32(g0,g2,Mod2),h2=sub32(g0,g2,Mod2);
const auto u0=Ladd32(h0,h1,Mod2),u1=Lsub32(h0,h1,Mod2);
const auto u2=Ladd32(h2,h3,Mod2),u3=Lsub32(h2,h3,Mod2);
store256(p0,u0),store256(p1,u1),store256(p2,u2),store256(p3,u3);
}
}
I256*const g=f+j;
for(idt l=mm,L=mm>>2;L;l=L,L>>=2,t-=2,--p){
auto rt=load256(st_1+p);
for(idt i=(j==0?l:0),k=(j+i)>>t;i<m;i+=l,++k){
const auto r1=_mm256_permutevar8x32_epi32(rt,id);
const auto r2=_mm256_shuffle_epi32(r1,_MM_PERM_BBBB);
const auto nr3=_mm256_shuffle_epi32(r1,_MM_PERM_DDDD);
for(idt j=0;j<L;++j){
auto const p0=g+i+j,p1=p0+L,p2=p1+L,p3=p2+L;
const auto f1=load256(p1),f3=load256(p3),f2=load256(p2),f0=load256(p0);
const auto g1=mul_bsm(f1,r1,Niv,Mod),ng3=mul_bsm(f3,nr3,Niv,Mod);
const auto g2=mul_bsm(f2,r2,Niv,Mod),g0=shrk32(f0,Mod2);
const auto h3=mul_bsmfxd(Ladd32(g1,ng3,Mod2),Img,ImgNiv,Mod),h1=sub32(g1,ng3,Mod2);
const auto h0=add32(g0,g2,Mod2),h2=sub32(g0,g2,Mod2);
const auto u0=Ladd32(h0,h1,Mod2),u1=Lsub32(h0,h1,Mod2);
const auto u2=Ladd32(h2,h3,Mod2),u3=Lsub32(h2,h3,Mod2);
store256(p0,u0),store256(p1,u1),store256(p2,u2),store256(p3,u3);
}
rt=mul_upd_rt(rt,load256(info->rt3+__builtin_ctzll(~k)),Mod);
}
store256(st_1+p,rt);
}
// const auto pr2=load256(&info->pr2),pr4=load256(&info->pr4);
// const auto pr2Niv=load256(&info->pr2niv),pr4Niv=load256(&info->pr4niv);
// for(idt i=j;i<j+m;++i){
// auto fi=load256(f+i);
// fi=mul(fi,rr,Niv,Mod);
// rr=shrk32(mul_bfxd(rr,load256(info->rt4+__builtin_ctzll(~i)),load256(info->rt4niv+__builtin_ctzll(~i)),Mod),Mod);
// fi=mul_bfxd(Ladd32(neg32_m<0xf0>(fi,Mod2),_mm256_permute2x128_si256(fi,fi,1),Mod2),pr4,pr4Niv,Mod);
// fi=mul_bfxd(Ladd32(neg32_m<0xcc>(fi,Mod2),_mm256_shuffle_epi32(fi,0x4e),Mod2),pr2,pr2Niv,Mod);
// fi=sub32(_mm256_shuffle_epi32(fi,0xb1),neg32_m<0x55>(fi,Mod2),Mod2);
// store256(f+i,fi);
// }
}
}
template<bool shrk=false>inline void vector_dit(I256*const f,idt n,const FNTT32_info*const info){
alignas(32) std::array<u32,8> st_1[_mxlg>>1];
const I256 Mod=_mm256_set1_epi32(info->mod),Mod2=_mm256_set1_epi32(info->mod2),Niv=_mm256_set1_epi32(info->niv);
const I256 Img=_mm256_set1_epi32(info->img),ImgNiv=_mm256_set1_epi32(info->imgniv),id=_mm256_setr_epi32(0,2,0,4,0,2,0,4);
const int lgn=__builtin_ctzll(n);
std::fill(st_1,st_1+(_lg_itth>>1),info->bwbr);
std::fill(st_1+(_lg_itth>>1),st_1+(_mxlg>>1),info->bwbi);
const idt nn=n>>(lgn&1),mm=std::min(nn,_itth);
// I256 rr=_mm256_set1_epi32((info->mod-1)>>(lgn+3));
for(idt j=0;j<n;j+=mm){
// const auto pr2=load256(&info->pr2i),pr4=load256(&info->pr4i);
// const auto pr2Niv=load256(&info->pr2iniv),pr4Niv=load256(&info->pr4iniv);
// for(idt i=j;i<j+mm;++i){
// auto fi=load256(f+i);
// const auto rt=rr;
// rr=shrk32(mul_bfxd(rr,load256(info->rt4i+__builtin_ctzll(~i)),load256(info->rt4iniv+__builtin_ctzll(~i)),Mod),Mod);
// fi=mul_bfxd(Ladd32(neg32_m<0xaa>(fi,Mod2),_mm256_shuffle_epi32(fi,0xb1),Mod2),pr2,pr2Niv,Mod);
// fi=mul_bfxd(Ladd32(neg32_m<0xcc>(fi,Mod2),_mm256_shuffle_epi32(fi,0x4e),Mod2),pr4,pr4Niv,Mod);
// fi=mul(Ladd32(neg32_m<0xf0>(fi,Mod2),_mm256_permute2x128_si256(fi,fi,1),Mod2),rt,Niv,Mod);
// store256(f+i,fi);
// }
I256*const g=f+j;
int t=2,p=0;
for(idt l=4,L=1;l<=mm;L=l,l<<=2,t+=2,++p){
auto rt=load256(st_1+p);
for(idt i=0,k=j>>t;i<mm;i+=l,++k){
const auto r1=_mm256_permutevar8x32_epi32(rt,id);
const auto r2=_mm256_shuffle_epi32(r1,_MM_PERM_BBBB);
const auto r3=_mm256_shuffle_epi32(r1,_MM_PERM_DDDD);
for(idt j=0;j<L;++j){
auto const p0=g+i+j,p1=p0+L,p2=p1+L,p3=p2+L;
const auto f0=load256(p0),f1=load256(p1),f2=load256(p2),f3=load256(p3);
const auto g0=add32(f0,f1,Mod2),g1=sub32(f0,f1,Mod2);
const auto g2=add32(f2,f3,Mod2),g3=mul_bsmfxd(Lsub32(f3,f2,Mod2),Img,ImgNiv,Mod);
const auto h0=Ladd32(g0,g2,Mod2),h1=Ladd32(g1,g3,Mod2);
const auto h2=Lsub32(g0,g2,Mod2),h3=Lsub32(g1,g3,Mod2);
const auto u0=shrk32(h0,Mod2),u1=mul_bsm(h1,r1,Niv,Mod);
const auto u2=mul_bsm(h2,r2,Niv,Mod),u3=mul_bsm(h3,r3,Niv,Mod);
store256(p0,u0),store256(p1,u1),store256(p2,u2),store256(p3,u3);
}
rt=mul_upd_rt(rt,load256(info->rt3i+__builtin_ctzll(~k)),Mod);
}
store256(st_1+p,rt);
}
int tt=std::min(__builtin_ctzll(~(j>>_lg_itth))+_lg_itth,lgn);
for(idt L=_itth,l=L<<2;t<=tt;L=l,l<<=2,t+=2,++p){
if((j+_itth)==l){
if(shrk && l==n){
for(idt i=0;i<L;++i){
auto const p0=f+i,p1=p0+L,p2=p1+L,p3=p2+L;
const auto f2=load256(p2),f3=load256(p3),f0=load256(p0),f1=load256(p1);
const auto g3=mul_bsmfxd(Lsub32(f3,f2,Mod2),Img,ImgNiv,Mod),g2=add32(f2,f3,Mod2);
const auto g0=add32(f0,f1,Mod2),g1=sub32(f0,f1,Mod2);
const auto h0=add32(g0,g2,Mod2),h1=add32(g1,g3,Mod2);
const auto h2=sub32(g0,g2,Mod2),h3=sub32(g1,g3,Mod2);
const auto u0=shrk32(h0,Mod),u1=shrk32(h1,Mod);
const auto u2=shrk32(h2,Mod),u3=shrk32(h3,Mod);
store256(p0,u0),store256(p1,u1),store256(p2,u2),store256(p3,u3);
}
}
else{
for(idt i=0;i<L;++i){
auto const p0=f+i,p1=p0+L,p2=p1+L,p3=p2+L;
const auto f2=load256(p2),f3=load256(p3),f0=load256(p0),f1=load256(p1);
const auto g3=mul_bsmfxd(Lsub32(f3,f2,Mod2),Img,ImgNiv,Mod),g2=add32(f2,f3,Mod2);
const auto g0=add32(f0,f1,Mod2),g1=sub32(f0,f1,Mod2);
const auto h0=add32(g0,g2,Mod2),h1=add32(g1,g3,Mod2);
const auto h2=sub32(g0,g2,Mod2),h3=sub32(g1,g3,Mod2);
store256(p0,h0),store256(p1,h1),store256(p2,h2),store256(p3,h3);
}
}
}
else{
auto rt=load256(st_1+p);
const auto r1=_mm256_permutevar8x32_epi32(rt,id);
const auto r1Niv=_mm256_permutevar8x32_epi32(_mm256_mul_epu32(rt,Niv),id);
rt=mul_upd_rt(rt,load256(info->rt3i+__builtin_ctzll(~j>>t)),Mod);
const auto r2=_mm256_shuffle_epi32(r1,_MM_PERM_BBBB),r3=_mm256_shuffle_epi32(r1,_MM_PERM_DDDD);
const auto r2Niv=_mm256_shuffle_epi32(r1Niv,_MM_PERM_BBBB),r3Niv=_mm256_shuffle_epi32(r1Niv,_MM_PERM_DDDD);
store256(st_1+p,rt);
for(idt i=0;i<L;++i){
auto const p0=f+j+_itth-l+i,p1=p0+L,p2=p1+L,p3=p2+L;
const auto f0=load256(p0),f1=load256(p1),f2=load256(p2),f3=load256(p3);
const auto g0=add32(f0,f1,Mod2),g1=sub32(f0,f1,Mod2);
const auto g2=add32(f2,f3,Mod2),g3=mul_bsmfxd(Lsub32(f3,f2,Mod2),Img,ImgNiv,Mod);
const auto h0=Ladd32(g0,g2,Mod2),h1=Ladd32(g1,g3,Mod2);
const auto h2=Lsub32(g0,g2,Mod2),h3=Lsub32(g1,g3,Mod2);
const auto u0=shrk32(h0,Mod2),u1=mul_bsmfxd(h1,r1,r1Niv,Mod);
const auto u2=mul_bsmfxd(h2,r2,r2Niv,Mod),u3=mul_bsmfxd(h3,r3,r3Niv,Mod);
store256(p0,u0),store256(p1,u1),store256(p2,u2),store256(p3,u3);
}
}
}
}
if(shrk && nn==n && n<=_itth){
for(idt i=0;i<n;++i){
const auto f0=load256(f+i);
store256(f+i,shrk32(f0,Mod));
}
}
if(nn!=n){
for(idt i=0;i<nn;++i){
auto const p0=f+i,p1=f+nn+i;
const auto f0=load256(p0),f1=load256(p1);
const auto g0=add32(f0,f1,Mod2),g1=sub32(f0,f1,Mod2);
if constexpr(shrk){
const auto h0=shrk32(g0,Mod),h1=shrk32(g1,Mod);
store256(p0,h0),store256(p1,h1);
}
else{
store256(p0,g0),store256(p1,g1);
}
}
}
}
// Returns fx * f[0,8) * g[0,8) (mod x^8 - ww).
[[gnu::always_inline]] inline I256 convolve8(const I256*f,const I256*g,I256 ww,I256 fx,I256 Niv,I256 Mod,I256 Mod2){
const auto raa=load256(f),rbb=load256(g);
const auto taa=shrk32(raa,Mod2),bb=shrk32(mul_bsm(rbb,fx,Niv,Mod),Mod);
const auto aw=shrk32(mul_bsm(taa,ww,Niv,Mod),Mod);
const auto aa=shrk32(taa,Mod);
const auto awa=_mm256_permute2x128_si256(aa,aw,3);
const auto b0=_mm256_permute4x64_epi64(bb,0x00),b1=_mm256_shuffle_epi32(b0,_MM_PERM_CDAB);
const auto a0=aa,a1=_mm256_srli_epi64(a0,32);
const auto aw7=_mm256_alignr_epi8(aa,awa,12);
auto res00=_mm256_mul_epu32(a0,b0);
auto res01=_mm256_mul_epu32(a1,b0);
auto res10=_mm256_mul_epu32(aw7,b1);
auto res11=_mm256_mul_epu32(a0,b1);
const auto b2=_mm256_permute4x64_epi64(bb,0x55),b3=_mm256_shuffle_epi32(b2,_MM_PERM_CDAB);
const auto aw6=_mm256_alignr_epi8(aa,awa,8);
const auto aw5=_mm256_alignr_epi8(aa,awa,4);
res00=_mm256_add_epi64(res00,_mm256_mul_epu32(aw6,b2));
res01=_mm256_add_epi64(res01,_mm256_mul_epu32(aw7,b2));
res10=_mm256_add_epi64(res10,_mm256_mul_epu32(aw5,b3));
res11=_mm256_add_epi64(res11,_mm256_mul_epu32(aw6,b3));
const auto b4=_mm256_permute4x64_epi64(bb,0xaa),b5=_mm256_shuffle_epi32(b4,_MM_PERM_CDAB);
const auto aw3=_mm256_alignr_epi8(awa,aw,12);
res00=_mm256_add_epi64(res00,_mm256_mul_epu32(awa,b4));
res01=_mm256_add_epi64(res01,_mm256_mul_epu32(aw5,b4));
res10=_mm256_add_epi64(res10,_mm256_mul_epu32(aw3,b5));
res11=_mm256_add_epi64(res11,_mm256_mul_epu32(awa,b5));
const auto b6=_mm256_permute4x64_epi64(bb,0xff),b7=_mm256_shuffle_epi32(b6,_MM_PERM_CDAB);
const auto aw2=_mm256_alignr_epi8(awa,aw,8);
const auto aw1=_mm256_alignr_epi8(awa,aw,4);
res00=_mm256_add_epi64(res00,_mm256_mul_epu32(aw2,b6));
res01=_mm256_add_epi64(res01,_mm256_mul_epu32(aw3,b6));
res10=_mm256_add_epi64(res10,_mm256_mul_epu32(aw1,b7));
res11=_mm256_add_epi64(res11,_mm256_mul_epu32(aw2,b7));
res00=_mm256_add_epi64(res00,res10);
res01=_mm256_add_epi64(res01,res11);
return shrk32(reduce(res00,res01,Niv,Mod),Mod2);
}
inline void vector_convolution_direct(I256*f,const I256*g,idt lm,const FNTT32_info*const info){
u32 RR=info->one;
const auto mod=info->mod,niv=info->niv;
const auto Fx=_mm256_set1_epi32(mul_s((mod-((mod-1)>>(__builtin_ctzll(lm)))),info->r3,niv,mod));
const auto Niv=_mm256_set1_epi32(niv),Mod=_mm256_set1_epi32(mod),Mod2=_mm256_set1_epi32(info->mod2);
for(idt i=0;i<lm;++i){
store256(f+i,convolve8(f+i,g+i,_mm256_set1_epi32(RR),Fx,Niv,Mod,Mod2));
RR=mul(RR,info->RT1[__builtin_ctzll(~i)],niv,mod);
}
}
inline void vector_convolution_accumulate(I256*const result,const I256*const f,
const I256*const g,idt lm,
const FNTT32_info*const info){
u32 RR=info->one;
const auto mod=info->mod,niv=info->niv;
const auto Fx=_mm256_set1_epi32(mul_s((mod-((mod-1)>>(__builtin_ctzll(lm)))),info->r3,niv,mod));
const auto Niv=_mm256_set1_epi32(niv),Mod=_mm256_set1_epi32(mod),Mod2=_mm256_set1_epi32(info->mod2);
for(idt i=0;i<lm;++i){
const auto product=convolve8(f+i,g+i,_mm256_set1_epi32(RR),Fx,Niv,Mod,Mod2);
store256(result+i,add32(load256(result+i),product,Mod2));
RR=mul(RR,info->RT1[__builtin_ctzll(~i)],niv,mod);
}
}
} // namespace fast998_v2
} // namespace internal
} // namespace fps
} // namespace m1une
#endif // M1UNE_FPS_HAS_X86_SIMD
#line 24 "math/fps/convolution.hpp"
#ifdef M1UNE_FPS_HAS_X86_SIMD
#pragma GCC pop_options
#endif
#line 1 "math/modint.hpp"
#line 9 "math/modint.hpp"
namespace m1une {
namespace math {
template <uint32_t Modulus>
struct ModInt {
static_assert(0 < Modulus, "Modulus must be positive");
private:
uint32_t _v;
public:
static constexpr uint32_t mod() {
return Modulus;
}
static constexpr ModInt raw(uint32_t v) noexcept {
ModInt x;
x._v = v;
return x;
}
constexpr ModInt() noexcept : _v(0) {}
template <class Integer, std::enable_if_t<std::is_integral_v<Integer>, int> = 0>
constexpr ModInt(Integer v) noexcept {
if constexpr (std::is_signed_v<Integer>) {
int64_t x = static_cast<int64_t>(v) % static_cast<int64_t>(Modulus);
if (x < 0) x += Modulus;
_v = static_cast<uint32_t>(x);
} else {
_v = static_cast<uint32_t>(static_cast<uint64_t>(v) % Modulus);
}
}
constexpr uint32_t val() const noexcept {
return _v;
}
constexpr ModInt& operator++() noexcept {
_v++;
if (_v == Modulus) _v = 0;
return *this;
}
constexpr ModInt& operator--() noexcept {
if (_v == 0) _v = Modulus;
_v--;
return *this;
}
constexpr ModInt operator++(int) noexcept {
ModInt res = *this;
++*this;
return res;
}
constexpr ModInt operator--(int) noexcept {
ModInt res = *this;
--*this;
return res;
}
constexpr ModInt& operator+=(const ModInt& rhs) noexcept {
_v += rhs._v;
if (_v >= Modulus) _v -= Modulus;
return *this;
}
constexpr ModInt& operator-=(const ModInt& rhs) noexcept {
_v -= rhs._v;
if (_v >= Modulus) _v += Modulus;
return *this;
}
constexpr ModInt& operator*=(const ModInt& rhs) noexcept {
uint64_t z = _v;
z *= rhs._v;
_v = static_cast<uint32_t>(z % Modulus);
return *this;
}
constexpr ModInt& operator/=(const ModInt& rhs) noexcept {
return *this *= rhs.inv();
}
constexpr ModInt operator+(const ModInt& rhs) const noexcept {
return ModInt(*this) += rhs;
}
constexpr ModInt operator-(const ModInt& rhs) const noexcept {
return ModInt(*this) -= rhs;
}
constexpr ModInt operator*(const ModInt& rhs) const noexcept {
return ModInt(*this) *= rhs;
}
constexpr ModInt operator/(const ModInt& rhs) const noexcept {
return ModInt(*this) /= rhs;
}
constexpr bool operator==(const ModInt& rhs) const noexcept {
return _v == rhs._v;
}
constexpr bool operator!=(const ModInt& rhs) const noexcept {
return _v != rhs._v;
}
constexpr ModInt pow(long long n) const noexcept {
ModInt res = raw(1 % Modulus);
ModInt x = n < 0 ? inv() : *this;
uint64_t exponent = n < 0 ? uint64_t(-(n + 1)) + 1 : uint64_t(n);
while (exponent > 0) {
if (exponent & 1) res *= x;
x *= x;
exponent >>= 1;
}
return res;
}
constexpr ModInt inv() const noexcept {
int64_t a = _v, b = Modulus, u = 1, v = 0;
while (b) {
int64_t t = a / b;
a -= t * b;
std::swap(a, b);
u -= t * v;
std::swap(u, v);
}
assert(a == 1);
u %= Modulus;
if (u < 0) u += Modulus;
return raw(static_cast<uint32_t>(u));
}
friend std::ostream& operator<<(std::ostream& os, const ModInt& rhs) {
return os << rhs._v;
}
friend std::istream& operator>>(std::istream& is, ModInt& rhs) {
long long v;
is >> v;
rhs = ModInt(v);
return is;
}
};
using modint998244353 = ModInt<998244353>;
using modint1000000007 = ModInt<1000000007>;
template <int Id = 0>
struct DynamicModInt {
private:
uint32_t _v;
inline static uint32_t _mod = 1;
public:
static uint32_t mod() noexcept {
return _mod;
}
static void set_mod(uint32_t modulus) noexcept {
assert(modulus > 0);
assert(modulus <= uint32_t(1) << 31);
_mod = modulus;
}
static DynamicModInt raw(uint32_t v) noexcept {
assert(v < _mod);
DynamicModInt x;
x._v = v;
return x;
}
DynamicModInt() noexcept : _v(0) {}
template <class Integer, std::enable_if_t<std::is_integral_v<Integer>, int> = 0>
DynamicModInt(Integer v) noexcept {
if constexpr (std::is_signed_v<Integer>) {
int64_t x = static_cast<int64_t>(v) % static_cast<int64_t>(_mod);
if (x < 0) x += _mod;
_v = static_cast<uint32_t>(x);
} else {
_v = static_cast<uint32_t>(static_cast<uint64_t>(v) % _mod);
}
}
uint32_t val() const noexcept {
return _v;
}
DynamicModInt& operator++() noexcept {
_v++;
if (_v == _mod) _v = 0;
return *this;
}
DynamicModInt& operator--() noexcept {
if (_v == 0) _v = _mod;
_v--;
return *this;
}
DynamicModInt operator++(int) noexcept {
DynamicModInt result = *this;
++*this;
return result;
}
DynamicModInt operator--(int) noexcept {
DynamicModInt result = *this;
--*this;
return result;
}
DynamicModInt& operator+=(const DynamicModInt& rhs) noexcept {
_v += rhs._v;
if (_v >= _mod) _v -= _mod;
return *this;
}
DynamicModInt& operator-=(const DynamicModInt& rhs) noexcept {
_v -= rhs._v;
if (_v >= _mod) _v += _mod;
return *this;
}
DynamicModInt& operator*=(const DynamicModInt& rhs) noexcept {
_v = static_cast<uint32_t>(uint64_t(_v) * rhs._v % _mod);
return *this;
}
DynamicModInt& operator/=(const DynamicModInt& rhs) noexcept {
return *this *= rhs.inv();
}
DynamicModInt operator+(const DynamicModInt& rhs) const noexcept {
return DynamicModInt(*this) += rhs;
}
DynamicModInt operator-(const DynamicModInt& rhs) const noexcept {
return DynamicModInt(*this) -= rhs;
}
DynamicModInt operator*(const DynamicModInt& rhs) const noexcept {
return DynamicModInt(*this) *= rhs;
}
DynamicModInt operator/(const DynamicModInt& rhs) const noexcept {
return DynamicModInt(*this) /= rhs;
}
bool operator==(const DynamicModInt& rhs) const noexcept {
return _v == rhs._v;
}
bool operator!=(const DynamicModInt& rhs) const noexcept {
return _v != rhs._v;
}
DynamicModInt pow(long long exponent) const noexcept {
DynamicModInt result = raw(1 % _mod);
DynamicModInt base = exponent < 0 ? inv() : *this;
uint64_t magnitude =
exponent < 0 ? uint64_t(-(exponent + 1)) + 1 : uint64_t(exponent);
while (magnitude > 0) {
if (magnitude & 1) result *= base;
base *= base;
magnitude >>= 1;
}
return result;
}
DynamicModInt inv() const noexcept {
int64_t a = _v, b = _mod, u = 1, v = 0;
while (b) {
int64_t quotient = a / b;
a -= quotient * b;
std::swap(a, b);
u -= quotient * v;
std::swap(u, v);
}
assert(a == 1);
u %= _mod;
if (u < 0) u += _mod;
return raw(static_cast<uint32_t>(u));
}
friend std::ostream& operator<<(std::ostream& os, const DynamicModInt& rhs) {
return os << rhs._v;
}
friend std::istream& operator>>(std::istream& is, DynamicModInt& rhs) {
long long value;
is >> value;
rhs = DynamicModInt(value);
return is;
}
};
} // namespace math
} // namespace m1une
#line 29 "math/fps/convolution.hpp"
namespace m1une {
namespace fps {
namespace internal {
template <class Mint, class = void>
struct has_static_modulus : std::false_type {};
template <class Mint>
struct has_static_modulus<
Mint, std::void_t<decltype(std::integral_constant<uint32_t, Mint::mod()>{})>>
: std::true_type {};
constexpr uint32_t primitive_root_constexpr(uint32_t mod) {
if (mod == 2) return 1;
if (mod == 167772161) return 3;
if (mod == 469762049) return 3;
if (mod == 754974721) return 11;
if (mod == 998244353) return 3;
if (mod == 1224736769) return 3;
uint32_t divisors[32] = {};
int count = 0;
uint32_t x = mod - 1;
for (uint32_t p = 2; uint64_t(p) * p <= x; p++) {
if (x % p != 0) continue;
divisors[count++] = p;
while (x % p == 0) x /= p;
}
if (x > 1) divisors[count++] = x;
for (uint32_t g = 2;; g++) {
bool ok = true;
for (int i = 0; i < count; i++) {
uint64_t value = 1;
uint64_t base = g;
uint32_t exponent = (mod - 1) / divisors[i];
while (exponent > 0) {
if (exponent & 1) value = value * base % mod;
base = base * base % mod;
exponent >>= 1;
}
if (value == 1) {
ok = false;
break;
}
}
if (ok) return g;
}
}
constexpr int two_adic_order(uint32_t x) {
int result = 0;
while ((x & 1) == 0) {
x >>= 1;
result++;
}
return result;
}
template <class Mint>
struct NttRoots {
static constexpr int max_base = two_adic_order(Mint::mod() - 1);
std::array<Mint, max_base + 1> root;
std::array<Mint, max_base + 1> inverse_root;
std::array<Mint, max_base> rate;
std::array<Mint, max_base> inverse_rate;
std::array<Mint, max_base> rate_radix4;
std::array<Mint, max_base> inverse_rate_radix4;
NttRoots() {
constexpr uint32_t primitive_root = primitive_root_constexpr(Mint::mod());
for (int level = 1; level <= max_base; level++) {
root[level] = Mint(primitive_root).pow((Mint::mod() - 1) >> level);
inverse_root[level] = root[level].inv();
}
Mint product = 1;
Mint inverse_product = 1;
for (int i = 0; i + 1 < max_base; i++) {
rate[i] = root[i + 2] * product;
inverse_rate[i] = inverse_root[i + 2] * inverse_product;
product *= inverse_root[i + 2];
inverse_product *= root[i + 2];
}
product = 1;
inverse_product = 1;
for (int i = 0; i + 2 < max_base; i++) {
rate_radix4[i] = root[i + 3] * product;
inverse_rate_radix4[i] = inverse_root[i + 3] * inverse_product;
product *= inverse_root[i + 3];
inverse_product *= root[i + 3];
}
}
};
template <class Mint>
const NttRoots<Mint>& ntt_roots() {
static const NttRoots<Mint> roots;
return roots;
}
template <class Mint>
void ntt(std::vector<Mint>& a, bool inverse, bool normalize = true) {
const int n = int(a.size());
assert(n > 0 && (n & (n - 1)) == 0);
assert((Mint::mod() - 1) % uint32_t(n) == 0);
const auto& roots = ntt_roots<Mint>();
const int height = two_adic_order(uint32_t(n));
if (!inverse) {
int phase = 0;
while (phase < height) {
if (height - phase == 1) {
const int width = 1 << (height - phase - 1);
Mint twiddle = 1;
for (int block = 0; block < (1 << phase); block++) {
const int offset = block << (height - phase);
for (int i = 0; i < width; i++) {
const Mint left = a[offset + i];
const Mint right = a[offset + i + width] * twiddle;
a[offset + i] = left + right;
a[offset + i + width] = left - right;
}
if (block + 1 != (1 << phase))
twiddle *= roots.rate[__builtin_ctz(~uint32_t(block))];
}
phase++;
continue;
}
const int width = 1 << (height - phase - 2);
Mint twiddle = 1;
const Mint imaginary = roots.root[2];
for (int block = 0; block < (1 << phase); block++) {
const Mint twiddle2 = twiddle * twiddle;
const Mint twiddle3 = twiddle2 * twiddle;
const int offset = block << (height - phase);
for (int i = 0; i < width; i++) {
const uint64_t mod2 = uint64_t(Mint::mod()) * Mint::mod();
const uint64_t a0 = a[offset + i].val();
const uint64_t a1 = uint64_t(a[offset + i + width].val()) * twiddle.val();
const uint64_t a2 =
uint64_t(a[offset + i + 2 * width].val()) * twiddle2.val();
const uint64_t a3 =
uint64_t(a[offset + i + 3 * width].val()) * twiddle3.val();
const uint64_t a1na3i =
uint64_t(Mint(a1 + mod2 - a3).val()) * imaginary.val();
const uint64_t negative_a2 = mod2 - a2;
a[offset + i] = Mint(a0 + a2 + a1 + a3);
a[offset + i + width] = Mint(a0 + a2 + 2 * mod2 - a1 - a3);
a[offset + i + 2 * width] = Mint(a0 + negative_a2 + a1na3i);
a[offset + i + 3 * width] = Mint(a0 + negative_a2 + mod2 - a1na3i);
}
if (block + 1 != (1 << phase))
twiddle *= roots.rate_radix4[__builtin_ctz(~uint32_t(block))];
}
phase += 2;
}
} else {
int phase = height;
while (phase > 0) {
if (phase == 1) {
const int width = 1 << (height - phase);
Mint twiddle = 1;
for (int block = 0; block < (1 << (phase - 1)); block++) {
const int offset = block << (height - phase + 1);
for (int i = 0; i < width; i++) {
const Mint left = a[offset + i];
const Mint right = a[offset + i + width];
a[offset + i] = left + right;
a[offset + i + width] = (left - right) * twiddle;
}
if (block + 1 != (1 << (phase - 1)))
twiddle *= roots.inverse_rate[__builtin_ctz(~uint32_t(block))];
}
phase--;
continue;
}
const int width = 1 << (height - phase);
Mint twiddle = 1;
const Mint inverse_imaginary = roots.inverse_root[2];
for (int block = 0; block < (1 << (phase - 2)); block++) {
const Mint twiddle2 = twiddle * twiddle;
const Mint twiddle3 = twiddle2 * twiddle;
const int offset = block << (height - phase + 2);
for (int i = 0; i < width; i++) {
const uint64_t a0 = a[offset + i].val();
const uint64_t a1 = a[offset + i + width].val();
const uint64_t a2 = a[offset + i + 2 * width].val();
const uint64_t a3 = a[offset + i + 3 * width].val();
const uint64_t a2na3i =
uint64_t(Mint((Mint::mod() + a2 - a3) * inverse_imaginary.val()).val());
a[offset + i] = Mint(a0 + a1 + a2 + a3);
a[offset + i + width] =
Mint((a0 + Mint::mod() - a1 + a2na3i) * twiddle.val());
a[offset + i + 2 * width] = Mint(
(a0 + a1 + 2ULL * Mint::mod() - a2 - a3) * twiddle2.val());
a[offset + i + 3 * width] = Mint(
(a0 + Mint::mod() - a1 + Mint::mod() - a2na3i) * twiddle3.val());
}
if (block + 1 != (1 << (phase - 2)))
twiddle *= roots.inverse_rate_radix4[__builtin_ctz(~uint32_t(block))];
}
phase -= 2;
}
if (normalize) {
const Mint inverse_n = Mint(n).inv();
for (Mint& value : a) value *= inverse_n;
}
}
}
#ifdef M1UNE_FPS_HAS_X86_SIMD
#pragma GCC push_options
#pragma GCC target("avx2,bmi")
template <class Mint>
__attribute__((target("avx2,bmi"), hot))
std::vector<Mint> convolution_998244353_simd(const std::vector<Mint>& a,
const std::vector<Mint>& b) {
const int result_size = int(a.size() + b.size() - 1);
int n = 1;
while (n < result_size) n <<= 1;
const bool squaring = &a == &b;
auto* transformed_a = static_cast<uint32_t*>(
::operator new[](sizeof(uint32_t) * n, std::align_val_t(32)));
auto* transformed_b = squaring
? transformed_a
: static_cast<uint32_t*>(::operator new[](
sizeof(uint32_t) * n, std::align_val_t(32)));
if constexpr (std::is_same_v<Mint, math::ModInt<998244353>>) {
static_assert(sizeof(Mint) == sizeof(uint32_t) && std::is_trivially_copyable_v<Mint>);
std::memcpy(transformed_a, a.data(), sizeof(uint32_t) * a.size());
if (!squaring)
std::memcpy(transformed_b, b.data(), sizeof(uint32_t) * b.size());
} else {
for (int i = 0; i < int(a.size()); i++) transformed_a[i] = a[i].val();
if (!squaring)
for (int i = 0; i < int(b.size()); i++) transformed_b[i] = b[i].val();
}
std::memset(transformed_a + a.size(), 0, sizeof(uint32_t) * (n - a.size()));
if (!squaring)
std::memset(transformed_b + b.size(), 0, sizeof(uint32_t) * (n - b.size()));
static constexpr fast998_v2::FNTT32_info transform(998244353);
const std::size_t vector_size = std::size_t(n) >> 3;
fast998_v2::vector_dif(reinterpret_cast<__m256i*>(transformed_a), vector_size, &transform);
if (!squaring)
fast998_v2::vector_dif(reinterpret_cast<__m256i*>(transformed_b), vector_size,
&transform);
fast998_v2::vector_convolution_direct(
reinterpret_cast<__m256i*>(transformed_a),
reinterpret_cast<const __m256i*>(transformed_b), vector_size, &transform);
fast998_v2::vector_dit<true>(reinterpret_cast<__m256i*>(transformed_a), vector_size,
&transform);
std::vector<Mint> result(result_size);
for (int j = 0; j < result_size; j++) result[j] = Mint::raw(transformed_a[j]);
::operator delete[](transformed_a, std::align_val_t(32));
if (!squaring) ::operator delete[](transformed_b, std::align_val_t(32));
return result;
}
#pragma GCC pop_options
#endif
} // namespace internal
template <class Mint>
std::vector<Mint> convolution_naive(const std::vector<Mint>& a, const std::vector<Mint>& b) {
if (a.empty() || b.empty()) return {};
std::vector<Mint> result(a.size() + b.size() - 1);
if (a.size() < b.size()) {
for (int i = 0; i < int(a.size()); i++) {
for (int j = 0; j < int(b.size()); j++) result[i + j] += a[i] * b[j];
}
} else {
for (int j = 0; j < int(b.size()); j++) {
for (int i = 0; i < int(a.size()); i++) result[i + j] += a[i] * b[j];
}
}
return result;
}
template <class Mint>
std::vector<Mint> convolution_ntt(const std::vector<Mint>& a, const std::vector<Mint>& b) {
const int result_size = int(a.size() + b.size() - 1);
int n = 1;
while (n < result_size) n <<= 1;
assert((Mint::mod() - 1) % uint32_t(n) == 0);
#ifdef M1UNE_FPS_HAS_X86_SIMD
if constexpr (Mint::mod() == 998244353) {
if (n >= 64 && __builtin_cpu_supports("avx2"))
return internal::convolution_998244353_simd(a, b);
}
#endif
// Allocate the padded buffers directly. Constructing from the inputs and
// then resizing used to allocate and copy both large operands twice.
const bool squaring = &a == &b;
std::vector<Mint> fa(n);
std::copy(a.begin(), a.end(), fa.begin());
internal::ntt(fa, false);
const Mint inverse_n = Mint(n).inv();
if (squaring) {
for (int i = 0; i < n; i++) fa[i] *= fa[i] * inverse_n;
} else {
std::vector<Mint> fb(n);
std::copy(b.begin(), b.end(), fb.begin());
internal::ntt(fb, false);
for (int i = 0; i < n; i++) fa[i] *= fb[i] * inverse_n;
}
internal::ntt(fa, true, false);
fa.resize(result_size);
return fa;
}
namespace internal {
template <class Mint>
std::vector<Mint> convolution_998244353_blocked_scalar(const std::vector<Mint>& a,
const std::vector<Mint>& b,
int transform_size) {
assert(Mint::mod() == 998244353);
assert(transform_size >= 2 && (transform_size & (transform_size - 1)) == 0);
assert((Mint::mod() - 1) % uint32_t(transform_size) == 0);
const int block_size = transform_size / 2;
const int a_blocks = int((a.size() + block_size - 1) / block_size);
const int b_blocks = int((b.size() + block_size - 1) / block_size);
auto transform_blocks = [&](const std::vector<Mint>& values, int block_count) {
std::vector<std::vector<Mint>> blocks;
blocks.reserve(block_count);
for (int block = 0; block < block_count; block++) {
const int begin = block * block_size;
const int count = std::min(block_size, int(values.size()) - begin);
std::vector<Mint> transformed(transform_size);
std::copy_n(values.begin() + begin, count, transformed.begin());
ntt(transformed, false);
blocks.emplace_back(std::move(transformed));
}
return blocks;
};
std::vector<std::vector<Mint>> transformed_a = transform_blocks(a, a_blocks);
std::vector<std::vector<Mint>> transformed_b = transform_blocks(b, b_blocks);
const int result_size = int(a.size() + b.size() - 1);
std::vector<Mint> result(result_size);
std::vector<Mint> transformed_result(transform_size);
for (int diagonal = 0; diagonal < a_blocks + b_blocks - 1; diagonal++) {
std::fill(transformed_result.begin(), transformed_result.end(), Mint(0));
const int first_a = std::max(0, diagonal - (b_blocks - 1));
const int last_a = std::min(a_blocks - 1, diagonal);
for (int a_block = first_a; a_block <= last_a; a_block++) {
const int b_block = diagonal - a_block;
for (int i = 0; i < transform_size; i++)
transformed_result[i] +=
transformed_a[a_block][i] * transformed_b[b_block][i];
}
ntt(transformed_result, true);
const int output_offset = diagonal * block_size;
const int output_count = std::min(transform_size, result_size - output_offset);
for (int i = 0; i < output_count; i++)
result[output_offset + i] += transformed_result[i];
}
return result;
}
#ifdef M1UNE_FPS_HAS_X86_SIMD
class AlignedUint32Buffer {
private:
uint32_t* data_;
public:
explicit AlignedUint32Buffer(std::size_t size)
: data_(static_cast<uint32_t*>(
::operator new[](sizeof(uint32_t) * size, std::align_val_t(32)))) {}
AlignedUint32Buffer(const AlignedUint32Buffer&) = delete;
AlignedUint32Buffer& operator=(const AlignedUint32Buffer&) = delete;
AlignedUint32Buffer(AlignedUint32Buffer&& other) noexcept : data_(other.data_) {
other.data_ = nullptr;
}
AlignedUint32Buffer& operator=(AlignedUint32Buffer&& other) noexcept {
if (this == &other) return *this;
::operator delete[](data_, std::align_val_t(32));
data_ = other.data_;
other.data_ = nullptr;
return *this;
}
~AlignedUint32Buffer() {
::operator delete[](data_, std::align_val_t(32));
}
uint32_t* data() {
return data_;
}
const uint32_t* data() const {
return data_;
}
};
template <class Mint>
__attribute__((target("avx2,bmi"), hot))
std::vector<Mint> convolution_998244353_blocked_simd(const std::vector<Mint>& a,
const std::vector<Mint>& b,
int transform_size) {
assert(Mint::mod() == 998244353);
assert(transform_size >= 64 && (transform_size & (transform_size - 1)) == 0);
assert((Mint::mod() - 1) % uint32_t(transform_size) == 0);
const int block_size = transform_size / 2;
const int a_blocks = int((a.size() + block_size - 1) / block_size);
const int b_blocks = int((b.size() + block_size - 1) / block_size);
static constexpr fast998_v2::FNTT32_info transform(998244353);
const std::size_t vector_size = std::size_t(transform_size) / 8;
auto transform_blocks = [&](const std::vector<Mint>& values, int block_count) {
std::vector<AlignedUint32Buffer> blocks;
blocks.reserve(block_count);
for (int block = 0; block < block_count; block++) {
const int begin = block * block_size;
const int count = std::min(block_size, int(values.size()) - begin);
AlignedUint32Buffer transformed(transform_size);
if constexpr (std::is_same_v<Mint, math::ModInt<998244353>>) {
static_assert(sizeof(Mint) == sizeof(uint32_t) &&
std::is_trivially_copyable_v<Mint>);
std::memcpy(transformed.data(), values.data() + begin,
sizeof(uint32_t) * count);
} else {
for (int i = 0; i < count; i++)
transformed.data()[i] = values[begin + i].val();
}
std::memset(transformed.data() + count, 0,
sizeof(uint32_t) * (transform_size - count));
fast998_v2::vector_dif(reinterpret_cast<__m256i*>(transformed.data()),
vector_size, &transform);
blocks.emplace_back(std::move(transformed));
}
return blocks;
};
std::vector<AlignedUint32Buffer> transformed_a = transform_blocks(a, a_blocks);
std::vector<AlignedUint32Buffer> transformed_b = transform_blocks(b, b_blocks);
const int result_size = int(a.size() + b.size() - 1);
std::vector<Mint> result(result_size);
AlignedUint32Buffer transformed_result(transform_size);
for (int diagonal = 0; diagonal < a_blocks + b_blocks - 1; diagonal++) {
std::memset(transformed_result.data(), 0, sizeof(uint32_t) * transform_size);
const int first_a = std::max(0, diagonal - (b_blocks - 1));
const int last_a = std::min(a_blocks - 1, diagonal);
for (int a_block = first_a; a_block <= last_a; a_block++) {
const int b_block = diagonal - a_block;
fast998_v2::vector_convolution_accumulate(
reinterpret_cast<__m256i*>(transformed_result.data()),
reinterpret_cast<const __m256i*>(transformed_a[a_block].data()),
reinterpret_cast<const __m256i*>(transformed_b[b_block].data()),
vector_size, &transform);
}
fast998_v2::vector_dit<true>(
reinterpret_cast<__m256i*>(transformed_result.data()), vector_size,
&transform);
const int output_offset = diagonal * block_size;
const int output_count = std::min(transform_size, result_size - output_offset);
for (int i = 0; i < output_count; i++) {
uint32_t value = result[output_offset + i].val() + transformed_result.data()[i];
if (value >= Mint::mod()) value -= Mint::mod();
result[output_offset + i] = Mint::raw(value);
}
}
return result;
}
#endif
template <class Mint>
std::vector<Mint> convolution_998244353_blocked(const std::vector<Mint>& a,
const std::vector<Mint>& b,
int transform_size = 1 << 23) {
#ifdef M1UNE_FPS_HAS_X86_SIMD
if (transform_size >= 64 && __builtin_cpu_supports("avx2"))
return convolution_998244353_blocked_simd(a, b, transform_size);
#endif
return convolution_998244353_blocked_scalar(a, b, transform_size);
}
} // namespace internal
template <class Mint>
std::vector<Mint> convolution(const std::vector<Mint>& a, const std::vector<Mint>& b) {
if (a.empty() || b.empty()) return {};
if (std::min(a.size(), b.size()) <= 32) return convolution_naive(a, b);
const int result_size = int(a.size() + b.size() - 1);
int n = 1;
while (n < result_size) n <<= 1;
if constexpr (internal::has_static_modulus<Mint>::value) {
if constexpr (Mint::mod() == 998244353) {
if (n > (1 << 23))
return internal::convolution_998244353_blocked(a, b);
}
if ((Mint::mod() - 1) % uint32_t(n) == 0) return convolution_ntt(a, b);
}
using Mint1 = math::ModInt<167772161>;
using Mint2 = math::ModInt<469762049>;
using Mint3 = math::ModInt<754974721>;
assert(n <= (1 << 24));
[[maybe_unused]] const unsigned __int128 coefficient_bound =
static_cast<unsigned __int128>(std::min(a.size(), b.size())) * (Mint::mod() - 1) *
(Mint::mod() - 1);
[[maybe_unused]] const unsigned __int128 crt_modulus =
static_cast<unsigned __int128>(Mint1::mod()) * Mint2::mod() * Mint3::mod();
assert(coefficient_bound < crt_modulus);
auto converted_convolution = [&]<class OtherMint>() {
std::vector<OtherMint> converted_a(a.size());
std::vector<OtherMint> converted_b(b.size());
for (int i = 0; i < int(a.size()); i++) converted_a[i] = OtherMint(a[i].val());
for (int i = 0; i < int(b.size()); i++) converted_b[i] = OtherMint(b[i].val());
return convolution_ntt(converted_a, converted_b);
};
std::vector<Mint1> c1 = converted_convolution.template operator()<Mint1>();
std::vector<Mint2> c2 = converted_convolution.template operator()<Mint2>();
std::vector<Mint3> c3 = converted_convolution.template operator()<Mint3>();
static const uint64_t inverse_mod1_mod2 = Mint2(Mint1::mod()).inv().val();
static const uint64_t mod1_mod3 = Mint1::mod() % Mint3::mod();
static const uint64_t mod1_mod2_mod3 =
mod1_mod3 * (Mint2::mod() % Mint3::mod()) % Mint3::mod();
static const uint64_t inverse_mod1_mod2_mod3 = Mint3(uint32_t(mod1_mod2_mod3)).inv().val();
const uint64_t target_mod = Mint::mod();
const uint64_t mod1_target = Mint1::mod() % target_mod;
const uint64_t mod1_mod2_target = mod1_target * (Mint2::mod() % target_mod) % target_mod;
std::vector<Mint> result(result_size);
for (int i = 0; i < result_size; i++) {
const uint64_t r1 = c1[i].val();
const uint64_t r2 = c2[i].val();
const uint64_t r3 = c3[i].val();
const uint64_t first =
(r2 + Mint2::mod() - r1 % Mint2::mod()) % Mint2::mod() * inverse_mod1_mod2 %
Mint2::mod();
const uint64_t combined_mod3 =
(r1 % Mint3::mod() + mod1_mod3 * (first % Mint3::mod())) % Mint3::mod();
const uint64_t second =
(r3 + Mint3::mod() - combined_mod3) % Mint3::mod() * inverse_mod1_mod2_mod3 %
Mint3::mod();
uint64_t value = r1 % target_mod;
value = (value + mod1_target * (first % target_mod)) % target_mod;
value = (value + mod1_mod2_target * (second % target_mod)) % target_mod;
result[i] = Mint::raw(uint32_t(value));
}
return result;
}
} // namespace fps
} // namespace m1une
#ifdef M1UNE_FPS_HAS_X86_SIMD
#undef M1UNE_FPS_HAS_X86_SIMD
#endif
#line 27 "utilities/bigint.hpp"
namespace m1une {
namespace utilities {
struct BigInt {
static constexpr int BASE = 1000000000;
static constexpr int BASE_DIGITS = 9;
std::vector<int> a;
int sign;
BigInt() : sign(1) {}
BigInt(long long v) {
*this = v;
}
BigInt(const std::string& s) {
read(s);
}
BigInt& operator=(long long v) {
sign = 1;
unsigned long long magnitude = static_cast<unsigned long long>(v);
if (v < 0) {
sign = -1;
magnitude = 0 - magnitude;
}
a.clear();
for (; magnitude > 0; magnitude /= BASE) {
a.push_back(int(magnitude % BASE));
}
return *this;
}
BigInt& operator=(const std::string& s) {
read(s);
return *this;
}
void trim() {
while (!a.empty() && a.back() == 0) {
a.pop_back();
}
if (a.empty()) sign = 1;
}
void read(const std::string& s) {
sign = 1;
a.clear();
int pos = 0;
while (pos < (int)s.size() && (s[pos] == '-' || s[pos] == '+')) {
if (s[pos] == '-') sign = -1;
++pos;
}
a.reserve((int(s.size()) - pos + BASE_DIGITS - 1) / BASE_DIGITS);
for (int i = int(s.size()) - 1; i >= pos; i -= BASE_DIGITS) {
int x = 0;
for (int j = std::max(pos, i - BASE_DIGITS + 1); j <= i; ++j) {
x = x * 10 + (s[j] - '0');
}
a.push_back(x);
}
trim();
}
std::string to_string() const {
if (a.empty()) return "0";
static const auto digit_quads = [] {
std::array<char, 40000> digits{};
for (int value = 0; value < 10000; ++value) {
int current = value;
for (int index = 3; index >= 0; --index) {
digits[4 * value + index] = char('0' + current % 10);
current /= 10;
}
}
return digits;
}();
char leading[BASE_DIGITS];
const std::to_chars_result converted =
std::to_chars(leading, leading + BASE_DIGITS, a.back());
assert(converted.ec == std::errc());
const int leading_size = int(converted.ptr - leading);
std::string res((sign == -1) + leading_size +
(a.size() - 1) * BASE_DIGITS,
'0');
int offset = 0;
if (sign == -1) res[offset++] = '-';
std::copy(leading, converted.ptr, res.begin() + offset);
offset += leading_size;
for (int i = (int)a.size() - 2; i >= 0; --i) {
const unsigned value = unsigned(a[i]);
const unsigned leading_digit = value / 100000000;
const unsigned remaining = value - leading_digit * 100000000;
const unsigned middle = remaining / 10000;
const unsigned trailing = remaining - middle * 10000;
res[offset] = char('0' + leading_digit);
std::memcpy(res.data() + offset + 1, digit_quads.data() + 4 * middle, 4);
std::memcpy(res.data() + offset + 5, digit_quads.data() + 4 * trailing, 4);
offset += BASE_DIGITS;
}
return res;
}
bool is_zero() const {
return a.empty() || (a.size() == 1 && a[0] == 0);
}
BigInt operator-() const {
BigInt res = *this;
if (!is_zero()) res.sign = -sign;
return res;
}
BigInt abs() const {
BigInt res = *this;
res.sign = 1;
return res;
}
friend bool operator<(const BigInt& x, const BigInt& y) {
if (x.sign != y.sign) return x.sign < y.sign;
if (x.a.size() != y.a.size()) {
return (x.sign == 1) ? (x.a.size() < y.a.size()) : (x.a.size() > y.a.size());
}
for (int i = (int)x.a.size() - 1; i >= 0; --i) {
if (x.a[i] != y.a[i]) {
return (x.sign == 1) ? (x.a[i] < y.a[i]) : (x.a[i] > y.a[i]);
}
}
return false;
}
friend bool operator>(const BigInt& x, const BigInt& y) {
return y < x;
}
friend bool operator<=(const BigInt& x, const BigInt& y) {
return !(y < x);
}
friend bool operator>=(const BigInt& x, const BigInt& y) {
return !(x < y);
}
friend bool operator==(const BigInt& x, const BigInt& y) {
return x.sign == y.sign && x.a == y.a;
}
friend bool operator!=(const BigInt& x, const BigInt& y) {
return !(x == y);
}
BigInt& operator+=(const BigInt& other) {
if (other.is_zero()) return *this;
if (is_zero()) return *this = other;
if (sign != other.sign) {
const int comparison = magnitude_compare(a, other.a);
if (comparison == 0) {
a.clear();
sign = 1;
} else if (comparison > 0) {
subtract_magnitude_inplace(a, other.a);
} else {
std::vector<int> result = other.a;
subtract_magnitude_inplace(result, a);
a = std::move(result);
sign = other.sign;
}
return *this;
}
add_magnitude_inplace(a, other.a);
return *this;
}
BigInt& operator-=(const BigInt& other) {
if (other.is_zero()) return *this;
if (is_zero()) return *this = -other;
if (sign != other.sign) {
add_magnitude_inplace(a, other.a);
return *this;
}
const int comparison = magnitude_compare(a, other.a);
if (comparison == 0) {
a.clear();
sign = 1;
} else if (comparison > 0) {
subtract_magnitude_inplace(a, other.a);
} else {
std::vector<int> result = other.a;
subtract_magnitude_inplace(result, a);
a = std::move(result);
sign = -sign;
}
return *this;
}
BigInt& operator*=(int v) {
if (v == 0 || is_zero()) return *this = 0;
long long multiplier = v;
if (multiplier < 0) {
sign = -sign;
multiplier = -multiplier;
}
a.reserve(a.size() + 2);
long long carry = 0;
for (int i = 0; i < (int)a.size() || carry; ++i) {
if (i == (int)a.size()) a.push_back(0);
const long long cur = a[i] * multiplier + carry;
carry = cur / BASE;
a[i] = (int)(cur % BASE);
}
trim();
return *this;
}
private:
static constexpr int MULTIPLICATION_THRESHOLD = 128;
static constexpr int SQUARE_THRESHOLD = 224;
static constexpr int DIVISION_THRESHOLD = 64;
static constexpr int FFT_SPLIT = 1 << 15;
struct FftComplex {
double real;
double imaginary;
FftComplex operator+(const FftComplex& other) const {
return {real + other.real, imaginary + other.imaginary};
}
FftComplex operator-(const FftComplex& other) const {
return {real - other.real, imaginary - other.imaginary};
}
FftComplex operator*(const FftComplex& other) const {
return {real * other.real - imaginary * other.imaginary,
real * other.imaginary + imaginary * other.real};
}
FftComplex operator*(double scalar) const {
return {real * scalar, imaginary * scalar};
}
FftComplex conjugate() const {
return {real, -imaginary};
}
};
struct FftProducts {
FftComplex diagonal;
FftComplex cross;
};
static const std::vector<FftComplex>& fft_roots(int size) {
static std::vector<FftComplex> roots(2, FftComplex{1, 0});
if (int(roots.size()) < size) {
int length = int(roots.size());
roots.resize(size);
while (length < size) {
const long double angle = std::numbers::pi_v<long double> / length;
const long double step_real = std::cos(angle);
const long double step_imaginary = std::sin(angle);
for (int i = length; i < 2 * length; ++i) {
roots[i] = roots[i / 2];
if (i & 1) {
const long double real = roots[i].real;
const long double imaginary = roots[i].imaginary;
roots[i] = {
double(real * step_real - imaginary * step_imaginary),
double(real * step_imaginary + imaginary * step_real)};
}
}
length *= 2;
}
}
return roots;
}
#ifdef M1UNE_BIGINT_HAS_X86_SIMD
__attribute__((target("avx2,fma"), always_inline)) static inline __m256d
multiply_complex(__m256d value, __m256d root) {
const __m256d real = _mm256_movedup_pd(value);
const __m256d imaginary = _mm256_permute_pd(value, 0xf);
const __m256d swapped_root = _mm256_permute_pd(root, 0x5);
return _mm256_fmaddsub_pd(real, root,
_mm256_mul_pd(imaginary, swapped_root));
}
__attribute__((target("avx2,fma"), hot)) static void fft_simd(
FftComplex* values, int size) {
const std::vector<FftComplex>& roots = fft_roots(size);
for (int length = size / 2; length > 0; length /= 2) {
for (int offset = 0; offset < size; offset += 2 * length) {
int i = 0;
for (; i + 1 < length; i += 2) {
const __m256d even = _mm256_loadu_pd(reinterpret_cast<const double*>(
values + offset + i));
const __m256d odd = _mm256_loadu_pd(reinterpret_cast<const double*>(
values + offset + i + length));
const __m256d root = _mm256_loadu_pd(reinterpret_cast<const double*>(
roots.data() + length + i));
_mm256_storeu_pd(reinterpret_cast<double*>(values + offset + i),
_mm256_add_pd(even, odd));
_mm256_storeu_pd(
reinterpret_cast<double*>(values + offset + i + length),
multiply_complex(_mm256_sub_pd(even, odd), root));
}
for (; i < length; ++i) {
const FftComplex even = values[offset + i];
const FftComplex odd = values[offset + i + length];
values[offset + i] = even + odd;
values[offset + i + length] =
(even - odd) * roots[length + i];
}
}
}
}
__attribute__((target("avx2,fma"), hot)) static void inverse_fft_simd(
FftComplex* values, int size) {
const std::vector<FftComplex>& roots = fft_roots(size);
const __m256d conjugate_mask = _mm256_setr_pd(0.0, -0.0, 0.0, -0.0);
for (int length = 1; length < size; length *= 2) {
for (int offset = 0; offset < size; offset += 2 * length) {
int i = 0;
for (; i + 1 < length; i += 2) {
const __m256d even = _mm256_loadu_pd(reinterpret_cast<const double*>(
values + offset + i));
const __m256d value = _mm256_loadu_pd(reinterpret_cast<const double*>(
values + offset + i + length));
__m256d root = _mm256_loadu_pd(reinterpret_cast<const double*>(
roots.data() + length + i));
root = _mm256_xor_pd(root, conjugate_mask);
const __m256d odd = multiply_complex(value, root);
_mm256_storeu_pd(reinterpret_cast<double*>(values + offset + i),
_mm256_add_pd(even, odd));
_mm256_storeu_pd(
reinterpret_cast<double*>(values + offset + i + length),
_mm256_sub_pd(even, odd));
}
for (; i < length; ++i) {
const FftComplex even = values[offset + i];
const FftComplex value = values[offset + i + length];
const FftComplex root = roots[length + i];
const FftComplex odd = {
value.real * root.real + value.imaginary * root.imaginary,
value.imaginary * root.real - value.real * root.imaginary};
values[offset + i] = even + odd;
values[offset + i + length] = even - odd;
}
}
}
const __m256d inverse_size = _mm256_set1_pd(1.0 / double(size));
int i = 0;
for (; i + 1 < size; i += 2) {
const __m256d value = _mm256_loadu_pd(
reinterpret_cast<const double*>(values + i));
_mm256_storeu_pd(reinterpret_cast<double*>(values + i),
_mm256_mul_pd(value, inverse_size));
}
for (; i < size; ++i) {
values[i].real /= size;
values[i].imaginary /= size;
}
}
#endif
static void fft(FftComplex* values, int size) {
assert(size > 0 && (size & (size - 1)) == 0);
#ifdef M1UNE_BIGINT_HAS_X86_SIMD
fft_simd(values, size);
return;
#endif
const std::vector<FftComplex>& roots = fft_roots(size);
// Decimation in frequency leaves the spectrum bit-reversed. The inverse
// transform consumes that order directly, avoiding permutation passes.
for (int length = size / 2; length > 0; length /= 2) {
for (int offset = 0; offset < size; offset += 2 * length) {
for (int i = 0; i < length; ++i) {
const FftComplex even = values[offset + i];
const FftComplex odd = values[offset + i + length];
values[offset + i] = even + odd;
values[offset + i + length] =
(even - odd) * roots[length + i];
}
}
}
}
static void inverse_fft(FftComplex* values, int size) {
assert(size > 0 && (size & (size - 1)) == 0);
#ifdef M1UNE_BIGINT_HAS_X86_SIMD
inverse_fft_simd(values, size);
return;
#endif
const std::vector<FftComplex>& roots = fft_roots(size);
for (int length = 1; length < size; length *= 2) {
for (int offset = 0; offset < size; offset += 2 * length) {
for (int i = 0; i < length; ++i) {
const FftComplex even = values[offset + i];
const FftComplex value = values[offset + i + length];
const FftComplex root = roots[length + i];
const FftComplex odd = {
value.real * root.real + value.imaginary * root.imaginary,
value.imaginary * root.real - value.real * root.imaginary};
values[offset + i] = even + odd;
values[offset + i + length] = even - odd;
}
}
}
const double inverse_size = 1.0 / double(size);
for (int i = 0; i < size; ++i) {
values[i].real *= inverse_size;
values[i].imaginary *= inverse_size;
}
}
static FftComplex fft_low(const FftComplex& value,
const FftComplex& reflected) {
return (value + reflected) * 0.5;
}
static FftComplex fft_high(const FftComplex& value,
const FftComplex& reflected) {
const FftComplex difference = value - reflected;
return {difference.imaginary * 0.5, -difference.real * 0.5};
}
static void inverse_real_fft(FftComplex* values, int size) {
assert(size >= 2 && (size & (size - 1)) == 0);
const int half_size = size / 2;
const std::vector<FftComplex>& roots = fft_roots(size);
static std::vector<FftComplex> bit_reversed_roots;
static int cached_size = 0;
if (cached_size != size) {
bit_reversed_roots.resize(half_size);
std::vector<int> reversed(half_size);
const int shift = std::countr_zero(unsigned(half_size));
for (int i = 1; i < half_size; ++i) {
reversed[i] =
(reversed[i / 2] >> 1) | ((i & 1) << (shift - 1));
}
for (int i = 0; i < half_size; ++i) {
bit_reversed_roots[i] = roots[half_size + reversed[i]];
}
cached_size = size;
}
for (int i = 0; i < half_size; ++i) {
const FftComplex first = values[2 * i];
const FftComplex second = values[2 * i + 1];
const FftComplex even = (first + second) * 0.5;
const FftComplex odd =
((first - second) * 0.5) * bit_reversed_roots[i].conjugate();
values[i] = {even.real - odd.imaginary,
even.imaginary + odd.real};
}
inverse_fft(values, half_size);
}
static void trim_magnitude(std::vector<int>& value) {
while (!value.empty() && value.back() == 0) value.pop_back();
}
static bool magnitude_less(const std::vector<int>& lhs, const std::vector<int>& rhs) {
return magnitude_compare(lhs, rhs) < 0;
}
static int magnitude_compare(const std::vector<int>& lhs,
const std::vector<int>& rhs) {
if (lhs.size() != rhs.size()) return lhs.size() < rhs.size() ? -1 : 1;
for (int i = int(lhs.size()) - 1; i >= 0; --i) {
if (lhs[i] != rhs[i]) return lhs[i] < rhs[i] ? -1 : 1;
}
return 0;
}
static void add_magnitude_inplace(std::vector<int>& lhs,
const std::vector<int>& rhs) {
const int lhs_size = int(lhs.size());
const int rhs_size = int(rhs.size());
const int size = std::max(lhs_size, rhs_size);
if (lhs_size < rhs_size) lhs.resize(rhs_size);
int carry = 0;
int i = 0;
for (; i < rhs_size; ++i) {
const long long current = (long long)lhs[i] + rhs[i] + carry;
lhs[i] = int(current >= BASE ? current - BASE : current);
carry = current >= BASE;
}
while (i < size && carry) {
++lhs[i];
carry = lhs[i] == BASE;
if (carry) lhs[i] = 0;
++i;
}
if (carry) lhs.push_back(1);
}
static void subtract_magnitude_inplace(std::vector<int>& lhs,
const std::vector<int>& rhs) {
assert(!magnitude_less(lhs, rhs));
int borrow = 0;
for (int i = 0; i < int(rhs.size()) || borrow; ++i) {
int current = lhs[i] - borrow - (i < int(rhs.size()) ? rhs[i] : 0);
borrow = current < 0;
if (borrow) current += BASE;
lhs[i] = current;
}
assert(borrow == 0);
trim_magnitude(lhs);
}
static bool magnitude_less_equal(const std::vector<int>& lhs,
const std::vector<int>& rhs) {
return !magnitude_less(rhs, lhs);
}
static std::vector<int> add_magnitude(const std::vector<int>& lhs,
const std::vector<int>& rhs) {
std::vector<int> result(std::max(lhs.size(), rhs.size()) + 1);
for (int i = 0; i < int(result.size()) - 1; ++i) {
if (i < int(lhs.size())) result[i] += lhs[i];
if (i < int(rhs.size())) result[i] += rhs[i];
if (result[i] >= BASE) {
result[i] -= BASE;
result[i + 1]++;
}
}
trim_magnitude(result);
return result;
}
static std::vector<int> subtract_magnitude(const std::vector<int>& lhs,
const std::vector<int>& rhs) {
assert(!magnitude_less(lhs, rhs));
std::vector<int> result = lhs;
int borrow = 0;
for (int i = 0; i < int(result.size()); ++i) {
const long long current =
(long long)result[i] - borrow - (i < int(rhs.size()) ? rhs[i] : 0);
if (current < 0) {
result[i] = int(current + BASE);
borrow = 1;
} else {
result[i] = int(current);
borrow = 0;
}
}
assert(borrow == 0);
trim_magnitude(result);
return result;
}
static std::vector<int> multiply_naive(const std::vector<int>& lhs,
const std::vector<int>& rhs) {
if (lhs.empty() || rhs.empty()) return std::vector<int>();
std::vector<long long> product(lhs.size() + rhs.size());
constexpr long long REDUCTION = 4LL * BASE * BASE;
for (int i = 0; i < int(lhs.size()); ++i) {
for (int j = 0; j < int(rhs.size()); ++j) {
product[i + j] += (long long)lhs[i] * rhs[j];
if (product[i + j] >= REDUCTION) {
product[i + j] -= REDUCTION;
product[i + j + 1] += 4LL * BASE;
}
}
}
std::vector<int> result;
result.reserve(product.size() + 1);
long long carry = 0;
for (int i = 0; i < int(product.size()) || carry > 0; ++i) {
if (i < int(product.size())) carry += product[i];
result.push_back(int(carry % BASE));
carry /= BASE;
}
trim_magnitude(result);
return result;
}
static std::vector<int> square_naive(const std::vector<int>& value) {
if (value.empty()) return std::vector<int>();
std::vector<long long> product(2 * value.size());
constexpr long long REDUCTION = 4LL * BASE * BASE;
for (int i = 0; i < int(value.size()); ++i) {
product[2 * i] += (long long)value[i] * value[i];
if (product[2 * i] >= REDUCTION) {
product[2 * i] -= REDUCTION;
product[2 * i + 1] += 4LL * BASE;
}
for (int j = i + 1; j < int(value.size()); ++j) {
product[i + j] += 2LL * value[i] * value[j];
if (product[i + j] >= REDUCTION) {
product[i + j] -= REDUCTION;
product[i + j + 1] += 4LL * BASE;
}
}
}
std::vector<int> result;
result.reserve(product.size() + 1);
long long carry = 0;
for (int i = 0; i < int(product.size()) || carry > 0; ++i) {
if (i < int(product.size())) carry += product[i];
result.push_back(int(carry % BASE));
carry /= BASE;
}
trim_magnitude(result);
return result;
}
static std::vector<int> multiply_by_limb(const std::vector<int>& value,
int multiplier) {
assert(0 <= multiplier && multiplier < BASE);
if (value.empty() || multiplier == 0) return std::vector<int>();
std::vector<int> result;
result.reserve(value.size() + 1);
uint64_t carry = 0;
for (int limb : value) {
const uint64_t current = uint64_t(limb) * multiplier + carry;
result.push_back(int(current % BASE));
carry = current / BASE;
}
if (carry) result.push_back(int(carry));
return result;
}
static std::vector<int> multiply_ntt(const std::vector<int>& lhs,
const std::vector<int>& rhs) {
using Mint1 = math::ModInt<998244353>;
using Mint2 = math::ModInt<754974721>;
using Mint3 = math::ModInt<469762049>;
const int result_size = int(lhs.size() + rhs.size() - 1);
assert(result_size <= (1 << 24));
auto convolve = [&]<class Mint>() {
std::vector<Mint> x(lhs.begin(), lhs.end());
if (&lhs == &rhs) return fps::convolution(x, x);
std::vector<Mint> y(rhs.begin(), rhs.end());
return fps::convolution(x, y);
};
const std::vector<Mint1> residues1 = convolve.template operator()<Mint1>();
const std::vector<Mint2> residues2 = convolve.template operator()<Mint2>();
const std::vector<Mint3> residues3 = convolve.template operator()<Mint3>();
constexpr uint64_t MOD1 = Mint1::mod();
constexpr uint64_t MOD2 = Mint2::mod();
constexpr uint64_t MOD3 = Mint3::mod();
constexpr uint64_t MOD12 = MOD1 * MOD2;
const static uint64_t inverse_mod1_mod2 = Mint2(MOD1).inv().val();
const static uint64_t inverse_mod12_mod3 = Mint3(MOD12 % MOD3).inv().val();
[[maybe_unused]] const unsigned __int128 coefficient_bound =
static_cast<unsigned __int128>(std::min(lhs.size(), rhs.size())) *
(BASE - 1) * (BASE - 1);
[[maybe_unused]] constexpr unsigned __int128 CRT_MODULUS =
static_cast<unsigned __int128>(MOD12) * MOD3;
assert(coefficient_bound < CRT_MODULUS);
std::vector<int> result;
result.reserve(result_size + 2);
unsigned __int128 carry = 0;
for (int i = 0; i < result_size || carry > 0; ++i) {
if (i < result_size) {
const uint64_t first = residues1[i].val();
const uint64_t second = residues2[i].val();
const uint64_t third = residues3[i].val();
const uint64_t difference12 =
(second + MOD2 - first % MOD2) % MOD2;
const uint64_t quotient12 =
difference12 * inverse_mod1_mod2 % MOD2;
const uint64_t combined12 = first + MOD1 * quotient12;
const uint64_t difference3 =
(third + MOD3 - combined12 % MOD3) % MOD3;
const uint64_t quotient3 =
difference3 * inverse_mod12_mod3 % MOD3;
carry += combined12 +
static_cast<unsigned __int128>(MOD12) * quotient3;
}
result.push_back(int(carry % BASE));
carry /= BASE;
}
trim_magnitude(result);
return result;
}
struct ProductFingerprint {
uint32_t mod31;
uint32_t mod29;
};
template <int bits>
static uint32_t mersenne_reduce(uint64_t value) {
constexpr uint64_t MODULUS = (uint64_t(1) << bits) - 1;
value = (value & MODULUS) + (value >> bits);
value = (value & MODULUS) + (value >> bits);
if (value >= MODULUS) value -= MODULUS;
return uint32_t(value);
}
static ProductFingerprint magnitude_fingerprint(const std::vector<int>& value) {
ProductFingerprint result{0, 0};
for (int i = int(value.size()) - 1; i >= 0; --i) {
result.mod31 = mersenne_reduce<31>(uint64_t(result.mod31) * BASE + value[i]);
result.mod29 = mersenne_reduce<29>(uint64_t(result.mod29) * BASE + value[i]);
}
return result;
}
static bool product_matches(const std::vector<int>& lhs, const std::vector<int>& rhs,
const std::vector<int>& product) {
const ProductFingerprint lhs_value = magnitude_fingerprint(lhs);
const ProductFingerprint rhs_value = magnitude_fingerprint(rhs);
const ProductFingerprint actual = magnitude_fingerprint(product);
return actual.mod31 ==
mersenne_reduce<31>(uint64_t(lhs_value.mod31) * rhs_value.mod31) &&
actual.mod29 ==
mersenne_reduce<29>(uint64_t(lhs_value.mod29) * rhs_value.mod29);
}
static std::vector<int> multiply_fft(const std::vector<int>& lhs,
const std::vector<int>& rhs) {
const int result_size = int(lhs.size() + rhs.size() - 1);
const int transform_size = int(std::bit_ceil(unsigned(result_size)));
const unsigned __int128 coefficient_bound =
static_cast<unsigned __int128>(std::min(lhs.size(), rhs.size())) *
2 * (FFT_SPLIT - 1) * ((BASE - 1) / FFT_SPLIT);
if (transform_size > (1 << 20) || coefficient_bound >= (uint64_t(1) << 50)) {
return multiply_ntt(lhs, rhs);
}
std::unique_ptr<FftComplex[]> transformed_lhs(
new FftComplex[transform_size]);
for (int i = 0; i < int(lhs.size()); ++i) {
transformed_lhs[i] = {double(lhs[i] % FFT_SPLIT),
double(lhs[i] / FFT_SPLIT)};
}
std::fill(transformed_lhs.get() + lhs.size(),
transformed_lhs.get() + transform_size, FftComplex{0, 0});
fft(transformed_lhs.get(), transform_size);
const bool squaring = &lhs == &rhs;
std::unique_ptr<FftComplex[]> transformed_rhs(
new FftComplex[transform_size]);
if (!squaring) {
for (int i = 0; i < int(rhs.size()); ++i) {
transformed_rhs[i] = {double(rhs[i] % FFT_SPLIT),
double(rhs[i] / FFT_SPLIT)};
}
std::fill(transformed_rhs.get() + rhs.size(),
transformed_rhs.get() + transform_size, FftComplex{0, 0});
fft(transformed_rhs.get(), transform_size);
}
static std::vector<int> reflected_indices;
if (int(reflected_indices.size()) < transform_size) {
const int previous_size = int(reflected_indices.size());
reflected_indices.resize(transform_size);
for (int i = std::max(1, previous_size); i < transform_size; ++i) {
// Negating a frequency complements the bits below the highest
// set bit of its bit-reversed index.
reflected_indices[i] =
i ^ int(std::bit_floor(unsigned(i)) - 1);
}
}
auto calculate_products = [&](int index) {
const int opposite = reflected_indices[index];
const FftComplex lhs_reflected =
transformed_lhs[opposite].conjugate();
const FftComplex lhs_low =
fft_low(transformed_lhs[index], lhs_reflected);
const FftComplex lhs_high =
fft_high(transformed_lhs[index], lhs_reflected);
FftComplex rhs_low = lhs_low;
FftComplex rhs_high = lhs_high;
if (!squaring) {
const FftComplex rhs_reflected =
transformed_rhs[opposite].conjugate();
rhs_low = fft_low(transformed_rhs[index], rhs_reflected);
rhs_high = fft_high(transformed_rhs[index], rhs_reflected);
}
const FftComplex low_product = lhs_low * rhs_low;
const FftComplex high_product = lhs_high * rhs_high;
const FftComplex diagonal =
low_product + FftComplex{-high_product.imaginary,
high_product.real};
const FftComplex cross =
lhs_low * rhs_high + lhs_high * rhs_low;
return FftProducts{diagonal, cross};
};
for (int i = 0; i < transform_size; ++i) {
const int opposite = reflected_indices[i];
if (i > opposite) continue;
const FftProducts products = calculate_products(i);
FftProducts opposite_products = products;
if (i != opposite) opposite_products = calculate_products(opposite);
transformed_lhs[i] = products.diagonal;
transformed_lhs[opposite] = opposite_products.diagonal;
transformed_rhs[i] = products.cross;
transformed_rhs[opposite] = opposite_products.cross;
}
inverse_fft(transformed_lhs.get(), transform_size);
inverse_real_fft(transformed_rhs.get(), transform_size);
std::vector<int> result;
result.reserve(result_size + 2);
unsigned __int128 carry = 0;
for (int i = 0; i < result_size || carry > 0; ++i) {
if (i < result_size) {
const long long low = std::llround(transformed_lhs[i].real);
const long long high = std::llround(transformed_lhs[i].imaginary);
const FftComplex packed_cross = transformed_rhs[i / 2];
const long long cross = std::llround(
(i & 1) ? packed_cross.imaginary : packed_cross.real);
if (low < 0 || high < 0 || cross < 0) return multiply_ntt(lhs, rhs);
carry += low + static_cast<unsigned __int128>(cross) * FFT_SPLIT +
static_cast<unsigned __int128>(high) * FFT_SPLIT * FFT_SPLIT;
}
result.push_back(int(carry % BASE));
carry /= BASE;
}
trim_magnitude(result);
if (result.empty() || !product_matches(lhs, rhs, result)) {
return multiply_ntt(lhs, rhs);
}
return result;
}
static std::vector<int> multiply_magnitude(const std::vector<int>& lhs,
const std::vector<int>& rhs) {
if (lhs.empty() || rhs.empty()) return std::vector<int>();
if (lhs.size() == 1) return multiply_by_limb(rhs, lhs[0]);
if (rhs.size() == 1) return multiply_by_limb(lhs, rhs[0]);
if (&lhs == &rhs && lhs.size() <= SQUARE_THRESHOLD) {
return square_naive(lhs);
}
if (std::min(lhs.size(), rhs.size()) <= MULTIPLICATION_THRESHOLD) {
return multiply_naive(lhs, rhs);
}
return multiply_fft(lhs, rhs);
}
static std::pair<std::vector<int>, std::vector<int>> divide_by_limb(
const std::vector<int>& dividend, int divisor) {
assert(0 < divisor && divisor < BASE);
if (divisor == 1) {
return std::make_pair(dividend, std::vector<int>());
}
std::vector<int> quotient(dividend.size());
long long remainder = 0;
for (int i = int(dividend.size()) - 1; i >= 0; --i) {
const long long current = remainder * BASE + dividend[i];
quotient[i] = int(current / divisor);
remainder = current % divisor;
}
trim_magnitude(quotient);
std::vector<int> remainder_digits;
if (remainder != 0) remainder_digits.push_back(int(remainder));
return std::make_pair(std::move(quotient), std::move(remainder_digits));
}
static std::pair<std::vector<int>, std::vector<int>> divide_classical(
const std::vector<int>& dividend, const std::vector<int>& divisor) {
assert(!divisor.empty());
if (divisor.size() == 1) return divide_by_limb(dividend, divisor[0]);
if (magnitude_less(dividend, divisor)) {
return std::make_pair(std::vector<int>(), dividend);
}
const int normalization = BASE / (divisor.back() + 1);
std::vector<int> normalized_divisor(divisor.size());
uint64_t carry = 0;
for (int i = 0; i < int(divisor.size()); ++i) {
const uint64_t current = uint64_t(divisor[i]) * normalization + carry;
normalized_divisor[i] = int(current % BASE);
carry = current / BASE;
}
assert(carry == 0);
std::vector<int> normalized_dividend(dividend.size() + 1);
carry = 0;
for (int i = 0; i < int(dividend.size()); ++i) {
const uint64_t current = uint64_t(dividend[i]) * normalization + carry;
normalized_dividend[i] = int(current % BASE);
carry = current / BASE;
}
normalized_dividend[dividend.size()] = int(carry);
const int divisor_size = int(normalized_divisor.size());
const int quotient_size = int(dividend.size()) - divisor_size + 1;
const uint64_t leading_divisor = normalized_divisor.back();
const uint64_t second_divisor = normalized_divisor[divisor_size - 2];
std::vector<int> quotient(quotient_size);
for (int position = quotient_size - 1; position >= 0; --position) {
const uint64_t leading_dividend =
uint64_t(normalized_dividend[position + divisor_size]) * BASE +
normalized_dividend[position + divisor_size - 1];
uint64_t digit = leading_dividend / leading_divisor;
uint64_t remainder = leading_dividend % leading_divisor;
if (digit >= BASE) {
digit = BASE - 1;
remainder = leading_dividend - digit * leading_divisor;
}
while (remainder < BASE &&
digit * second_divisor >
remainder * BASE +
normalized_dividend[position + divisor_size - 2]) {
--digit;
remainder += leading_divisor;
}
uint64_t borrow = 0;
for (int i = 0; i < divisor_size; ++i) {
const uint64_t product =
digit * uint64_t(normalized_divisor[i]) + borrow;
const uint64_t low = product % BASE;
borrow = product / BASE;
if (uint64_t(normalized_dividend[position + i]) < low) {
normalized_dividend[position + i] =
int(uint64_t(normalized_dividend[position + i]) + BASE - low);
++borrow;
} else {
normalized_dividend[position + i] -= int(low);
}
}
long long top =
(long long)normalized_dividend[position + divisor_size] -
static_cast<long long>(borrow);
if (top < 0) {
--digit;
uint64_t add_carry = 0;
for (int i = 0; i < divisor_size; ++i) {
const uint64_t current =
uint64_t(normalized_dividend[position + i]) +
normalized_divisor[i] + add_carry;
normalized_dividend[position + i] = int(current % BASE);
add_carry = current / BASE;
}
top += add_carry;
}
assert(0 <= top && top < BASE);
normalized_dividend[position + divisor_size] = int(top);
quotient[position] = int(digit);
}
trim_magnitude(quotient);
std::vector<int> remainder(normalized_dividend.begin(),
normalized_dividend.begin() + divisor_size);
trim_magnitude(remainder);
std::pair<std::vector<int>, std::vector<int>> denormalized =
divide_by_limb(remainder, normalization);
assert(denormalized.second.empty());
return std::make_pair(std::move(quotient), std::move(denormalized.first));
}
static std::vector<int> reciprocal(const std::vector<int>& value, int degree) {
assert(!value.empty());
assert(BASE / 2 <= value.back() && value.back() < BASE);
assert(degree >= 0);
int precision = degree;
const int value_size = int(value.size());
while (precision > DIVISION_THRESHOLD) precision = (precision + 1) / 2;
std::vector<int> inverse(value_size + precision + 1);
inverse.back() = 1;
inverse = divide_classical(inverse, value).first;
while (precision < degree) {
std::vector<int> square = multiply_magnitude(inverse, inverse);
square.insert(square.begin(), 0);
const int copied = std::min(value_size, 2 * precision + 1);
const std::vector<int> leading(value.end() - copied, value.end());
// The original Newton formula right-aligns `leading` in 2p + 1
// limbs, then discards the same 2p + 1 low limbs of the product.
// Cancelling that shift avoids convolving a long zero prefix.
std::vector<int> correction = multiply_magnitude(square, leading);
assert(int(correction.size()) >= copied);
correction.erase(correction.begin(), correction.begin() + copied);
std::vector<int> shifted(precision + 1);
const std::vector<int> doubled = add_magnitude(inverse, inverse);
shifted.insert(shifted.end(), doubled.begin(), doubled.end());
inverse = subtract_magnitude(shifted, correction);
assert(!inverse.empty());
inverse.erase(inverse.begin());
precision *= 2;
}
assert(precision >= degree);
inverse.erase(inverse.begin(), inverse.begin() + precision - degree);
trim_magnitude(inverse);
return inverse;
}
static std::pair<std::vector<int>, std::vector<int>> divide_magnitude(
const std::vector<int>& dividend, const std::vector<int>& divisor) {
assert(!divisor.empty());
if (divisor.size() <= DIVISION_THRESHOLD ||
int(dividend.size()) - int(divisor.size()) <= DIVISION_THRESHOLD) {
return divide_classical(dividend, divisor);
}
if (dividend.size() > 2 * divisor.size()) {
const int size_ratio = int(dividend.size() / divisor.size());
const int block_multiplier = size_ratio >= 16 ? 3 : size_ratio >= 7 ? 2 : 1;
const int block_size = block_multiplier * int(divisor.size());
const int block_count =
(int(dividend.size()) + block_size - 1) / block_size;
const int normalization = BASE / (divisor.back() + 1);
const std::vector<int> normalized_divisor =
multiply_by_limb(divisor, normalization);
const int degree = block_size + 3;
const std::vector<int> inverse = reciprocal(normalized_divisor, degree);
const int discarded = int(divisor.size()) + degree;
auto divide_block = [&](const std::vector<int>& current) {
const std::vector<int> normalized_current =
multiply_by_limb(current, normalization);
std::vector<int> quotient_product =
multiply_magnitude(normalized_current, inverse);
std::vector<int> partial_quotient;
if (int(quotient_product.size()) > discarded) {
partial_quotient.assign(quotient_product.begin() + discarded,
quotient_product.end());
}
std::vector<int> product =
multiply_magnitude(normalized_divisor, partial_quotient);
while (magnitude_less(normalized_current, product)) {
partial_quotient =
subtract_magnitude(partial_quotient, std::vector<int>(1, 1));
product = subtract_magnitude(product, normalized_divisor);
}
std::vector<int> partial_remainder =
subtract_magnitude(normalized_current, product);
while (magnitude_less_equal(normalized_divisor, partial_remainder)) {
partial_quotient =
add_magnitude(partial_quotient, std::vector<int>(1, 1));
partial_remainder =
subtract_magnitude(partial_remainder, normalized_divisor);
}
trim_magnitude(partial_quotient);
trim_magnitude(partial_remainder);
std::pair<std::vector<int>, std::vector<int>> denormalized =
divide_by_limb(partial_remainder, normalization);
assert(denormalized.second.empty());
return std::make_pair(std::move(partial_quotient),
std::move(denormalized.first));
};
std::vector<int> quotient(dividend.size());
std::vector<int> remainder;
for (int block = block_count - 1; block >= 0; --block) {
const int begin = block * block_size;
const int end = std::min(begin + block_size, int(dividend.size()));
std::vector<int> current(dividend.begin() + begin,
dividend.begin() + end);
current.insert(current.end(), remainder.begin(), remainder.end());
trim_magnitude(current);
std::pair<std::vector<int>, std::vector<int>> partial =
divide_block(current);
assert(int(partial.first.size()) <= end - begin);
std::copy(partial.first.begin(), partial.first.end(),
quotient.begin() + begin);
remainder = std::move(partial.second);
}
trim_magnitude(quotient);
trim_magnitude(remainder);
return std::make_pair(std::move(quotient), std::move(remainder));
}
const int normalization = BASE / (divisor.back() + 1);
const std::vector<int> normalized_dividend =
multiply_magnitude(dividend, std::vector<int>(1, normalization));
const std::vector<int> normalized_divisor =
multiply_magnitude(divisor, std::vector<int>(1, normalization));
const int dividend_size = int(normalized_dividend.size());
const int divisor_size = int(normalized_divisor.size());
const int degree = dividend_size - divisor_size + 2;
const std::vector<int> inverse = reciprocal(normalized_divisor, degree);
std::vector<int> quotient = multiply_magnitude(normalized_dividend, inverse);
const int discarded = divisor_size + degree;
assert(discarded <= int(quotient.size()));
quotient.erase(quotient.begin(), quotient.begin() + discarded);
std::vector<int> product = multiply_magnitude(normalized_divisor, quotient);
while (magnitude_less(normalized_dividend, product)) {
quotient = subtract_magnitude(quotient, std::vector<int>(1, 1));
product = subtract_magnitude(product, normalized_divisor);
}
std::vector<int> remainder = subtract_magnitude(normalized_dividend, product);
while (magnitude_less_equal(normalized_divisor, remainder)) {
quotient = add_magnitude(quotient, std::vector<int>(1, 1));
remainder = subtract_magnitude(remainder, normalized_divisor);
}
trim_magnitude(quotient);
trim_magnitude(remainder);
std::pair<std::vector<int>, std::vector<int>> denormalized =
divide_by_limb(remainder, normalization);
assert(denormalized.second.empty());
return std::make_pair(std::move(quotient), std::move(denormalized.first));
}
public:
BigInt& operator*=(const BigInt& other) {
if (is_zero() || other.is_zero()) return *this = 0;
const int result_sign = sign * other.sign;
a = multiply_magnitude(a, other.a);
sign = result_sign;
trim();
return *this;
}
friend std::pair<BigInt, BigInt> divmod(const BigInt& a1, const BigInt& b1) {
if (b1.is_zero()) {
throw std::domain_error("BigInt division by zero");
}
std::pair<std::vector<int>, std::vector<int>> result =
divide_magnitude(a1.a, b1.a);
BigInt q, r;
q.a = std::move(result.first);
r.a = std::move(result.second);
q.sign = a1.sign * b1.sign;
r.sign = a1.sign;
q.trim();
r.trim();
return {q, r};
}
friend BigInt gcd(BigInt first, BigInt second) {
first.sign = 1;
second.sign = 1;
while (!second.is_zero()) {
first %= second;
std::swap(first, second);
}
return first;
}
BigInt& operator/=(const BigInt& other) {
return *this = divmod(*this, other).first;
}
BigInt& operator%=(const BigInt& other) {
return *this = divmod(*this, other).second;
}
friend BigInt operator+(BigInt x, const BigInt& y) {
return x += y;
}
friend BigInt operator-(BigInt x, const BigInt& y) {
return x -= y;
}
friend BigInt operator*(BigInt x, const BigInt& y) {
return x *= y;
}
friend BigInt operator/(BigInt x, const BigInt& y) {
return x /= y;
}
friend BigInt operator%(BigInt x, const BigInt& y) {
return x %= y;
}
friend std::ostream& operator<<(std::ostream& os, const BigInt& b) {
return os << b.to_string();
}
friend std::istream& operator>>(std::istream& is, BigInt& b) {
std::string s;
if (is >> s) b.read(s);
return is;
}
};
} // namespace utilities
} // namespace m1une
#ifdef M1UNE_BIGINT_HAS_X86_SIMD
#undef M1UNE_BIGINT_HAS_X86_SIMD
#endif
#line 5 "verify/math/rational.test.cpp"
#line 1 "utilities/fast_io.hpp"
#line 6 "utilities/fast_io.hpp"
#include <cerrno>
#line 8 "utilities/fast_io.hpp"
#include <cstddef>
#include <cstdio>
#include <cstdlib>
#line 13 "utilities/fast_io.hpp"
#include <iterator>
#line 15 "utilities/fast_io.hpp"
#include <sys/stat.h>
#line 18 "utilities/fast_io.hpp"
#include <unistd.h>
#line 20 "utilities/fast_io.hpp"
namespace m1une {
namespace utilities {
struct FastOutput;
namespace internal {
// Shared with the convenience helpers in template.hpp.
inline FastOutput* standard_output_instance = nullptr;
// Detect std::begin(x), std::end(x).
template <class T, class = void>
struct is_range : std::false_type {};
template <class T>
struct is_range<T, std::void_t<
decltype(std::begin(std::declval<T&>())),
decltype(std::end(std::declval<T&>()))
>> : std::true_type {};
template <class T>
inline constexpr bool is_range_v = is_range<T>::value;
template <class T>
using range_reference_t = decltype(*std::begin(std::declval<T&>()));
template <class T>
using range_value_t = std::remove_cv_t<std::remove_reference_t<range_reference_t<T>>>;
template <class T, class = void>
struct range_stored_value {
using type = range_value_t<T>;
};
template <class T>
struct range_stored_value<T, std::void_t<typename std::remove_cv_t<std::remove_reference_t<T>>::value_type>> {
using type = typename std::remove_cv_t<std::remove_reference_t<T>>::value_type;
};
template <class T>
using range_stored_value_t = typename range_stored_value<T>::type;
// Treat strings and C strings as scalar output objects, not as ranges.
template <class T>
struct is_char_array : std::false_type {};
template <class T, std::size_t N>
struct is_char_array<T[N]>
: std::bool_constant<std::is_same_v<std::remove_cv_t<T>, char>> {};
template <class T>
struct is_string_like
: std::bool_constant<
std::is_same_v<std::decay_t<T>, std::string>
|| std::is_same_v<std::decay_t<T>, const char*>
|| std::is_same_v<std::decay_t<T>, char*>
|| is_char_array<std::remove_reference_t<T>>::value
> {};
template <class T>
inline constexpr bool is_string_like_v = is_string_like<T>::value;
// ModInt-like type: x.val() is printable, and x can be assigned from long long.
template <class T, class = void>
struct has_val_method : std::false_type {};
template <class T>
struct has_val_method<T, std::void_t<decltype(std::declval<const T&>().val())>>
: std::true_type {};
template <class T>
inline constexpr bool has_val_method_v = has_val_method<T>::value;
template <class T, class = void>
struct has_static_mod_raw : std::false_type {};
template <class T>
struct has_static_mod_raw<
T, std::void_t<decltype(T::mod()), decltype(T::raw(std::declval<uint32_t>()))>>
: std::true_type {};
template <class T>
inline constexpr bool has_static_mod_raw_v = has_static_mod_raw<T>::value;
// libstdc++ before GCC 16 does not classify __int128 as an integral type in
// strict ISO modes such as -std=c++23. Keep the fast-I/O interface independent
// of that implementation detail.
template <class T>
inline constexpr bool is_integral_v =
std::is_integral_v<T>
|| std::is_same_v<std::remove_cv_t<T>, __int128_t>
|| std::is_same_v<std::remove_cv_t<T>, __uint128_t>;
template <class T>
inline constexpr bool is_signed_v =
std::is_signed_v<T>
|| std::is_same_v<std::remove_cv_t<T>, __int128_t>;
template <class T>
struct make_unsigned {
using type = std::make_unsigned_t<T>;
};
template <>
struct make_unsigned<__int128_t> {
using type = __uint128_t;
};
template <>
struct make_unsigned<__uint128_t> {
using type = __uint128_t;
};
template <class T>
using make_unsigned_t = typename make_unsigned<std::remove_cv_t<T>>::type;
} // namespace internal
struct FastInput {
static constexpr int buffer_size = 1 << 20;
private:
std::FILE* _stream;
char _buffer[buffer_size];
int _position;
int _length;
int _file_descriptor;
bool _streaming;
bool refill() {
_position = 0;
if (_streaming) {
ssize_t length;
do {
length = ::read(_file_descriptor, _buffer, buffer_size);
} while (length < 0 && errno == EINTR);
if (length <= 0) {
_length = 0;
return false;
}
_length = int(length);
} else {
_length = int(std::fread(_buffer, 1, buffer_size, _stream));
}
return _length != 0;
}
template <class T>
bool read_integer_from_stream(T& value) {
if (!skip_spaces()) return false;
int c = read_char_raw();
bool negative = false;
if (c == '-') {
negative = true;
c = read_char_raw();
}
if constexpr (internal::is_signed_v<T>) {
T result = 0;
while ('0' <= c && c <= '9') {
result = negative ? result * 10 - (c - '0')
: result * 10 + (c - '0');
c = read_char_raw();
}
value = result;
} else {
T result = 0;
while ('0' <= c && c <= '9') {
result = result * 10 + T(c - '0');
c = read_char_raw();
}
value = negative ? T(0) - result : result;
}
return true;
}
bool prepare_number() {
if (_length - _position >= 64) return true;
const int remaining = _length - _position;
if (remaining > 0) std::memmove(_buffer, _buffer + _position, remaining);
const int added = int(std::fread(_buffer + remaining, 1, buffer_size - remaining, _stream));
_position = 0;
_length = remaining + added;
if (_length < buffer_size) _buffer[_length] = '\0';
return _length != 0;
}
public:
explicit FastInput(std::FILE* stream = stdin)
: _stream(stream),
_position(0),
_length(0),
_file_descriptor(::fileno(stream)),
_streaming([&] {
struct stat status;
return _file_descriptor >= 0
&& ::fstat(_file_descriptor, &status) == 0
&& !S_ISREG(status.st_mode);
}()) {}
FastInput(const FastInput&) = delete;
FastInput& operator=(const FastInput&) = delete;
int read_char_raw() {
if (_position == _length && !refill()) return EOF;
return _buffer[_position++];
}
bool skip_spaces() {
int c = read_char_raw();
while (c != EOF && c <= ' ') c = read_char_raw();
if (c == EOF) return false;
--_position;
return true;
}
bool read(char& value) {
if (!skip_spaces()) return false;
value = char(read_char_raw());
return true;
}
bool read(std::string& value) {
if (!skip_spaces()) return false;
value.clear();
while (true) {
const int begin = _position;
while (_position < _length &&
static_cast<unsigned char>(_buffer[_position]) > ' ') {
++_position;
}
value.append(_buffer + begin, _position - begin);
if (_position < _length) {
++_position;
return true;
}
if (!refill()) return true;
}
}
bool read(bool& value) {
int x;
if (!read(x)) return false;
value = x != 0;
return true;
}
template <class T>
std::enable_if_t<
internal::is_integral_v<T>
&& !std::is_same_v<std::remove_cv_t<T>, bool>
&& !std::is_same_v<std::remove_cv_t<T>, char>,
bool
>
read(T& value) {
if (_streaming) return read_integer_from_stream(value);
if (!prepare_number()) return false;
int c = static_cast<unsigned char>(_buffer[_position++]);
while (c <= ' ') c = static_cast<unsigned char>(_buffer[_position++]);
bool negative = false;
if (c == '-') {
negative = true;
c = static_cast<unsigned char>(_buffer[_position++]);
}
if constexpr (internal::is_signed_v<T>) {
T result = 0;
while ('0' <= c && c <= '9') {
const int first = c - '0';
const int second = static_cast<unsigned char>(_buffer[_position]) - '0';
if (0 <= second && second <= 9) {
result = negative ? result * 100 - (first * 10 + second)
: result * 100 + (first * 10 + second);
++_position;
} else {
result = negative ? result * 10 - first : result * 10 + first;
}
c = static_cast<unsigned char>(_buffer[_position++]);
}
value = result;
} else {
T result = 0;
while ('0' <= c && c <= '9') {
const unsigned first = unsigned(c - '0');
const int second = static_cast<unsigned char>(_buffer[_position]) - '0';
if (0 <= second && second <= 9) {
result = result * 100 + T(first * 10 + unsigned(second));
++_position;
} else {
result = result * 10 + T(first);
}
c = static_cast<unsigned char>(_buffer[_position++]);
}
value = negative ? T(0) - result : result;
}
if (_position > _length) _position = _length;
return true;
}
template <class T>
std::enable_if_t<std::is_floating_point_v<T>, bool>
read(T& value) {
if (!skip_spaces()) return false;
int c = read_char_raw();
bool negative = false;
if (c == '-' || c == '+') {
negative = c == '-';
c = read_char_raw();
}
long double result = 0;
while ('0' <= c && c <= '9') {
result = result * 10 + (c - '0');
c = read_char_raw();
}
if (c == '.') {
long double place = 0.1L;
c = read_char_raw();
while ('0' <= c && c <= '9') {
result += (c - '0') * place;
place *= 0.1L;
c = read_char_raw();
}
}
if (c == 'e' || c == 'E') {
c = read_char_raw();
bool exponent_negative = false;
if (c == '-' || c == '+') {
exponent_negative = c == '-';
c = read_char_raw();
}
int exponent = 0;
while ('0' <= c && c <= '9') {
exponent = exponent * 10 + (c - '0');
c = read_char_raw();
}
long double scale = 1;
long double power = 10;
while (exponent > 0) {
if (exponent & 1) scale *= power;
power *= power;
exponent >>= 1;
}
result = exponent_negative ? result / scale : result * scale;
}
value = static_cast<T>(negative ? -result : result);
return true;
}
template <class T>
std::enable_if_t<
internal::has_val_method_v<T>
&& !internal::is_integral_v<T>
&& !internal::is_range_v<T>,
bool
>
read(T& value) {
long long x;
if (!read(x)) return false;
if constexpr (internal::has_static_mod_raw_v<T>) {
if (x >= 0 && uint64_t(x) < uint64_t(T::mod())) {
value = T::raw(uint32_t(x));
} else {
value = T(x);
}
} else {
value = T(x);
}
return true;
}
template <class First, class Second>
bool read(std::pair<First, Second>& value) {
if (!read(value.first)) return false;
return read(value.second);
}
template <class Range>
std::enable_if_t<
internal::is_range_v<Range>
&& !internal::is_string_like_v<Range>,
bool
>
read(Range& range) {
using StoredValue = internal::range_stored_value_t<Range>;
constexpr bool nested = internal::is_range_v<StoredValue>
&& !internal::is_string_like_v<StoredValue>;
for (auto&& value : range) {
if constexpr (std::is_same_v<StoredValue, bool> && !nested) {
bool x;
if (!read(x)) return false;
value = x;
} else {
if (!read(value)) return false;
}
}
return true;
}
template <class First, class Second, class... Rest>
bool read(First& first, Second& second, Rest&... rest) {
if (!read(first)) return false;
return read(second, rest...);
}
template <class T>
FastInput& operator>>(T& value) {
if (!read(value)) std::abort();
return *this;
}
};
struct FastOutput {
static constexpr int buffer_size = 1 << 20;
private:
inline static const auto digit_quads = [] {
std::array<char, 40000> result{};
for (int i = 0; i < 10000; i++) {
int value = i;
for (int j = 3; j >= 0; j--) {
result[4 * i + j] = char('0' + value % 10);
value /= 10;
}
}
return result;
}();
std::FILE* _stream;
char _buffer[buffer_size];
int _position;
int _precision;
std::chars_format _float_format;
char _range_separator;
std::string* _capture = nullptr;
template <class T>
std::string format_cell(const T& value) {
std::string result;
struct CaptureGuard {
std::string*& target;
std::string* previous;
~CaptureGuard() { target = previous; }
} guard{_capture, _capture};
_capture = &result;
write(value);
return result;
}
template <class Matrix>
void write_aligned_matrix(const Matrix& matrix) {
std::vector<std::vector<std::string>> rows;
std::vector<std::size_t> widths;
for (const auto& row : matrix) {
auto& cells = rows.emplace_back();
std::size_t column = 0;
for (const auto& value : row) {
cells.push_back(format_cell(value));
if (column == widths.size()) widths.push_back(0);
widths[column] = std::max(widths[column], cells.back().size());
++column;
}
}
bool first = true;
for (const auto& row : rows) {
if (!first) write_char('\n');
first = false;
for (std::size_t column = 0; column < row.size(); ++column) {
if (column != 0) write_char(_range_separator);
for (std::size_t padding = row[column].size();
padding < widths[column]; ++padding) {
write_char(' ');
}
write(row[column]);
}
}
}
public:
explicit FastOutput(std::FILE* stream = stdout)
: _stream(stream),
_position(0),
_precision(6),
_float_format(std::chars_format::general),
_range_separator(' ') {
if (_stream == stdout
&& internal::standard_output_instance == nullptr) {
internal::standard_output_instance = this;
}
}
FastOutput(const FastOutput&) = delete;
FastOutput& operator=(const FastOutput&) = delete;
~FastOutput() {
flush();
if (internal::standard_output_instance == this) {
internal::standard_output_instance = nullptr;
}
}
void flush() {
if (_position != 0) {
std::fwrite(_buffer, 1, _position, _stream);
_position = 0;
}
std::fflush(_stream);
}
void write_char(char c) {
if (_capture != nullptr) {
_capture->push_back(c);
return;
}
if (_position == buffer_size) flush();
_buffer[_position++] = c;
}
void write(const char* s) {
while (*s != '\0') write_char(*s++);
}
void write(const std::string& s) {
if (_capture != nullptr) {
_capture->append(s);
return;
}
std::size_t position = 0;
while (position < s.size()) {
if (_position == buffer_size) flush();
const std::size_t copied =
std::min<std::size_t>(buffer_size - _position, s.size() - position);
std::memcpy(_buffer + _position, s.data() + position, copied);
_position += int(copied);
position += copied;
}
}
void write(char c) {
write_char(c);
}
void write(bool value) {
write_char(value ? '1' : '0');
}
template <class T>
std::enable_if_t<std::is_floating_point_v<T>>
write(T value) {
char digits[128];
auto [end, error] = std::to_chars(
digits,
digits + sizeof(digits),
value,
_float_format,
_precision
);
if (error != std::errc()) std::abort();
for (const char* pointer = digits; pointer != end; pointer++) {
write_char(*pointer);
}
}
template <class T>
std::enable_if_t<
internal::is_integral_v<T>
&& !std::is_same_v<std::remove_cv_t<T>, bool>
&& !std::is_same_v<std::remove_cv_t<T>, char>
>
write(T value) {
using Raw = std::remove_cv_t<T>;
using Unsigned = internal::make_unsigned_t<Raw>;
Unsigned magnitude;
if constexpr (internal::is_signed_v<Raw>) {
if (value < 0) {
write_char('-');
magnitude = Unsigned(0) - Unsigned(value);
} else {
magnitude = Unsigned(value);
}
} else {
magnitude = value;
}
if (magnitude == 0) {
write_char('0');
return;
}
unsigned chunks[16];
int count = 0;
while (magnitude >= 10000) {
const Unsigned quotient = magnitude / 10000;
chunks[count++] = unsigned(magnitude - quotient * 10000);
magnitude = quotient;
}
if (_capture == nullptr && _position > buffer_size - 64) flush();
char captured[64];
char* const begin = _capture != nullptr ? captured : _buffer + _position;
char* destination = begin;
const unsigned leading = unsigned(magnitude);
const char* first = digit_quads.data() + 4 * leading;
int skip = leading < 10 ? 3 : leading < 100 ? 2 : leading < 1000 ? 1 : 0;
for (; skip < 4; skip++) *destination++ = first[skip];
while (count--) {
const char* digits = digit_quads.data() + 4 * chunks[count];
std::memcpy(destination, digits, 4);
destination += 4;
}
if (_capture != nullptr) {
_capture->append(begin, destination - begin);
} else {
_position += int(destination - begin);
}
}
template <class T>
std::enable_if_t<
internal::has_val_method_v<T>
&& !internal::is_integral_v<T>
&& !internal::is_range_v<T>
>
write(const T& value) {
write(value.val());
}
template <class First, class Second>
void write(const std::pair<First, Second>& value) {
write(value.first);
write_char(' ');
write(value.second);
}
template <class Range>
std::enable_if_t<
internal::is_range_v<Range>
&& !internal::is_string_like_v<Range>
>
write(const Range& range) {
using StoredValue = internal::range_stored_value_t<const Range>;
constexpr bool nested = internal::is_range_v<StoredValue>
&& !internal::is_string_like_v<StoredValue>;
bool first = true;
for (const auto& value : range) {
if (!first) write_char(nested ? '\n' : _range_separator);
first = false;
if constexpr (std::is_same_v<StoredValue, bool> && !nested) {
write(static_cast<bool>(value));
} else {
write(value);
}
}
}
template <class First, class... Rest>
void print(const First& first, const Rest&... rest) {
write(first);
((write_char(' '), write(rest)), ...);
}
void println() {
write_char('\n');
}
void set_precision(int precision) {
_precision = precision;
}
void set_fixed(int precision = 6) {
_float_format = std::chars_format::fixed;
_precision = precision;
}
void set_general(int precision = 6) {
_float_format = std::chars_format::general;
_precision = precision;
}
void set_range_separator(char separator) {
_range_separator = separator;
}
template <class Matrix>
void write_aligned(const Matrix& matrix) {
using Row = internal::range_stored_value_t<const Matrix>;
using Cell = internal::range_stored_value_t<const Row>;
static_assert(internal::is_range_v<Row> && !internal::is_string_like_v<Row>,
"write_aligned requires a two-dimensional range");
static_assert(!internal::is_range_v<Cell> || internal::is_string_like_v<Cell>,
"write_aligned requires scalar cells");
write_aligned_matrix(matrix);
}
template <class Matrix>
void println_aligned(const Matrix& matrix) {
write_aligned(matrix);
write_char('\n');
}
template <class... Args>
void println(const Args&... args) {
print(args...);
write_char('\n');
}
template <class T>
FastOutput& operator<<(const T& value) {
write(value);
return *this;
}
};
} // namespace utilities
} // namespace m1une
#line 12 "verify/math/rational.test.cpp"
namespace {
using Fraction = m1une::math::Rational<long long>;
using BigInt = m1une::utilities::BigInt;
using BigFraction = m1une::math::Rational<BigInt>;
void test_fixed() {
constexpr Fraction zero;
constexpr Fraction half(2, 4);
constexpr Fraction negative(3, -6);
static_assert(zero.numerator() == 0);
static_assert(zero.denominator() == 1);
static_assert(half == Fraction(1, 2));
static_assert(negative == Fraction(-1, 2));
static_assert(half + negative == 0);
static_assert(Fraction(2, 3) + Fraction(5, 6) == Fraction(3, 2));
static_assert(Fraction(2, 3) * Fraction(9, 4) == Fraction(3, 2));
static_assert(Fraction(2, 3) / Fraction(4, 9) == Fraction(3, 2));
static_assert(Fraction(-7, 3).floor() == -3);
static_assert(Fraction(-7, 3).ceil() == -2);
static_assert(Fraction(-7, 3).trunc() == -2);
static_assert(Fraction(1, 3) < Fraction(1, 2));
static_assert(m1une::math::abs(Fraction(-3, 4)) == Fraction(3, 4));
[[maybe_unused]] Fraction large(
std::numeric_limits<long long>::max(),
std::numeric_limits<long long>::max() - 1
);
assert(large > 1);
std::stringstream stream;
stream << Fraction(-6, 8) << ' ' << Fraction(5);
assert(stream.str() == "-3/4 5");
Fraction first;
Fraction second;
stream.seekg(0);
stream >> first >> second;
assert(first == Fraction(-3, 4));
assert(second == 5);
}
void test_randomized() {
std::uint64_t state = 1301;
auto random = [&state]() {
state ^= state << 7;
state ^= state >> 9;
return state;
};
for (int trial = 0; trial < 100000; ++trial) {
long long a = static_cast<long long>(random() % 2001) - 1000;
long long b = 1 + static_cast<long long>(random() % 1000);
long long c = static_cast<long long>(random() % 2001) - 1000;
long long d = 1 + static_cast<long long>(random() % 1000);
Fraction first(a, b);
Fraction second(c, d);
[[maybe_unused]] __int128_t left = __int128_t(a) * d;
[[maybe_unused]] __int128_t right = __int128_t(c) * b;
assert((first <=> second) == (left <=> right));
[[maybe_unused]] Fraction sum = first + second;
assert(
__int128_t(sum.numerator()) * b * d
== (__int128_t(a) * d + __int128_t(c) * b)
* sum.denominator()
);
[[maybe_unused]] Fraction product = first * second;
assert(
__int128_t(product.numerator()) * b * d
== __int128_t(a) * c * product.denominator()
);
if (c != 0) {
[[maybe_unused]] Fraction quotient = first / second;
assert(
__int128_t(quotient.numerator()) * b * c
== __int128_t(a) * d * quotient.denominator()
);
}
}
}
void test_bigint() {
BigInt power = 1;
for (int i = 0; i < 100; ++i) power *= 10;
BigFraction reduced(power * 6, power * 8);
assert(reduced == BigFraction(3, 4));
assert(reduced.numerator() == 3);
assert(reduced.denominator() == 4);
BigFraction large(power, 3);
BigFraction inverse(3, power);
assert(large * inverse == 1);
assert(large > inverse);
assert(large + BigFraction(1, 3) == BigFraction(power + 1, 3));
assert(large / large == 1);
BigFraction negative(-(power + 1), power);
assert(negative.floor() == -2);
assert(negative.ceil() == -1);
assert(negative.trunc() == -1);
assert(negative.sign() == -1);
assert(abs(negative) == -negative);
BigFraction integer = 5;
assert(integer + 1 == 6);
assert(BigFraction(1, 2).to_long_double() == 0.5L);
long double near_two_thirds =
BigFraction(power * 2 + 1, power * 3 + 1).to_long_double();
assert(0.66L < near_two_thirds && near_two_thirds < 0.67L);
std::stringstream stream;
stream << large;
BigFraction parsed;
stream >> parsed;
assert(parsed == large);
std::uint64_t state = 2309;
auto random = [&state]() {
state ^= state << 7;
state ^= state >> 9;
return state;
};
auto assert_same = [](const BigFraction& big, const Fraction& small) {
assert(big.numerator().to_string() == std::to_string(small.numerator()));
assert(big.denominator().to_string() == std::to_string(small.denominator()));
};
for (int trial = 0; trial < 10000; ++trial) {
long long a = static_cast<long long>(random() % 2001) - 1000;
long long b = 1 + static_cast<long long>(random() % 1000);
long long c = static_cast<long long>(random() % 2001) - 1000;
long long d = 1 + static_cast<long long>(random() % 1000);
BigFraction big_first = BigFraction(BigInt(a), BigInt(b));
BigFraction big_second = BigFraction(BigInt(c), BigInt(d));
Fraction small_first(a, b);
Fraction small_second(c, d);
assert_same(big_first + big_second, small_first + small_second);
assert_same(big_first - big_second, small_first - small_second);
assert_same(big_first * big_second, small_first * small_second);
assert((big_first <=> big_second) == (small_first <=> small_second));
if (c != 0) {
assert_same(big_first / big_second, small_first / small_second);
}
}
}
} // namespace
int main() {
m1une::utilities::FastInput fast_input;
m1une::utilities::FastOutput fast_output;
test_fixed();
test_randomized();
test_bigint();
long long a, b;
fast_input >> a >> b;
fast_output << a + b << '\n';
}
verify/math/rational.test.cpp