#ifdef DEBUG #define _GLIBCXX_DEBUG #define DEXPR( LL , BOUND , VALUE , DEBUG_VALUE ) CEXPR( LL , BOUND , DEBUG_VALUE ) #define CERR( ANSWER ) cerr << ANSWER << endl; #define LIBRARY_SEARCH if( LibrarySearch() != 0 ){ QUIT; }; #else #pragma GCC optimize ( "O3" ) #pragma GCC optimize( "unroll-loops" ) #pragma GCC target ( "sse4.2,fma,avx2,popcnt,lzcnt,bmi2" ) #define DEXPR( LL , BOUND , VALUE , DEBUG_VALUE ) CEXPR( LL , BOUND , VALUE ) #define CERR( ANSWER ) #define LIBRARY_SEARCH #endif #include using namespace std; using uint = unsigned int; using ll = long long; using ull = unsigned long long; #define ATT __attribute__( ( target( "sse4.2,fma,avx2,popcnt,lzcnt,bmi2" ) ) ) #define TYPE_OF( VAR ) decay_t #define UNTIE ios_base::sync_with_stdio( false ); cin.tie( nullptr ) #define CEXPR( LL , BOUND , VALUE ) constexpr LL BOUND = VALUE #define CIN( LL , A ) LL A; cin >> A #define ASSERT( A , MIN , MAX ) assert( ( MIN ) <= A && A <= ( MAX ) ) #define CIN_ASSERT( A , MIN , MAX ) CIN( TYPE_OF( MAX ) , A ); ASSERT( A , MIN , MAX ) #define SET_ASSERT( A , MIN , MAX ) cin >> A; ASSERT( A , MIN , MAX ) #define GETLINE( A ) string A; getline( cin , A ) #define GETLINE_SEPARATE( A , SEPARATOR ) string A; getline( cin , A , SEPARATOR ) #define FOR( VAR , INITIAL , FINAL_PLUS_ONE ) for( TYPE_OF( FINAL_PLUS_ONE ) VAR = INITIAL ; VAR < FINAL_PLUS_ONE ; VAR ++ ) #define FOREQ( VAR , INITIAL , FINAL ) for( TYPE_OF( FINAL ) VAR = INITIAL ; VAR <= FINAL ; VAR ++ ) #define FOREQINV( VAR , INITIAL , FINAL ) for( TYPE_OF( INITIAL ) VAR = INITIAL ; VAR >= FINAL ; VAR -- ) #define FOR_ITR( ARRAY , ITR , END ) for( auto ITR = ARRAY .begin() , END = ARRAY .end() ; ITR != END ; ITR ++ ) #define REPEAT( HOW_MANY_TIMES ) FOR( VARIABLE_FOR_REPEAT , 0 , HOW_MANY_TIMES ) #define QUIT return 0 #define COUT( ANSWER ) cout << ( ANSWER ) << "\n" #define RETURN( ANSWER ) COUT( ANSWER ); QUIT #define SET_PRECISION( PRECISION ) cout << fixed << setprecision( PRECISION ) #define DOUBLE( PRECISION , ANSWER ) SET_PRECISION << ( ANSWER ) << "\n"; QUIT // 圧縮用 #define TE template #define TY typename #define US using #define ST static #define IN inline #define CL class #define PU public #define OP operator #define CE constexpr #define CO const #define NE noexcept #define RE return #define WH while #define VO void #define VE vector #define LI list #define BE begin #define EN end #define SZ size #define MO move #define TH this #define CRI CO int& #define CRUI CO uint& #define CRL CO ll& #define DEFINITION_OF_GET_FOR_ZETA_TRANSFORM( MU ) \ list sub = Sub( t ); \ sub.push_front( t ); \ U answer = Zero(); \ \ while( ! sub.empty() ){ \ \ const T& t_sub = sub.front(); \ answer = Sum( answer , Prod( m_val[e_inv( t_sub )] , MU( t_sub , t ) ) ); \ sub.pop_front(); \ \ } \ \ return answer; \ #define DEFINITION_OF_INVERSE_IMAGE_SUM_FOR_ZETA_TRANSFORM( MU ) \ const T t = f_inv_max( s ); \ list sub = r( s ); \ U answer = Zero(); \ \ while( ! sub.empty() ){ \ \ const T& t_sub = f_inv_max( sub.front() ); \ answer = Sum( answer , Prod( m_val[e_inv( t_sub )] , MU( t_sub , t ) ) ); \ sub.pop_front(); \ \ } \ \ return answer; \ template class ZetaTransformBody { private: int m_length; map m_memory; vector m_memory_inv; protected: int m_size; U m_val[size_max]; public: // 仮想継承用にダミーのデフォルトコンストラクタを設定。 inline ZetaTransformBody( const int& size = 0 ); inline void Add( const T& t , const U& u ); inline void AddAll( const U& u ); inline ZetaTransformBody& operator+=( const ZetaTransformBody& a ); inline void MultiplyAll( const U& u ); inline ZetaTransformBody& operator*=( const ZetaTransformBody& a ); // 子クラスの半順序のメビウス関数のデフォルトの再帰式を使うため、 // 再帰深度が浅い場合にしか使えない。 inline U Get( const T& t ); // muは子クラスの半順序のメビウス関数。 template inline U Get( const T& t ); inline const U& InitialSegmentSum( const T& t ); // f_inv_maxは定義域にr(s)を含む部分写像T->Sであり、要件 // (1) Sは半順序集合である。 // (2) r(s)は重複を持たない。(よって集合とみなす) // (3) sはr(s)の最大元である。 // (4) 像f_inv_max(s)の要素を上界に持つTの要素全体Rへのfの制限f_Rは順序保存写像R->r(s)である。 // (5) r(s)の任意の要素tに対しf_inv_max(t)が逆像f_R^{-1}(t)の最大元である。 // を満たすfが存在する場合にのみ以下の2つをサポート。 // f( t ) = sを満たすRの要素t全体を渡る総和取得。 template r(const S&)> inline U InverseImageSum( const S& s ); template r(const S&) , int mu(const T&,const T&)> inline U InverseImageSum( const S& s ); // f( t ) <= sを満たすRの要素t全体を渡る総和取得。(結果的にrは使わないが要件上はrの存在が必要) template inline const U& IntervalInverseImageSum( const S& s ); virtual T e( const int& i ); virtual int e_inv( const T& t ); private: virtual const U& Zero() const = 0; virtual U Sum( const U& u0 , const U& u1 ) const = 0; virtual U Prod( const U& u0 , const U& u1 ) const = 0; virtual list Sub( const T& t ) const = 0; virtual list Sup( const T& t ) const = 0; virtual int Moevius( const T& t0 , const T& t1 ) = 0; }; template class SemiRingForZetaTransform : virtual public ZetaTransformBody { public: inline SemiRingForZetaTransform( const int& dummy ); private: inline const U& Zero() const; inline U Sum( const U& u0 , const U& u1 ) const; inline U Prod( const U& u0 , const U& u1 ) const; }; template E(const T&) , list E_inv(const T&) , typename U , int size_max> class PartiallyOrderedSetForZetaTransform : virtual public ZetaTransformBody { private: map > m_moevius; inline list Sub( const T& t ) const; inline list Sup( const T& t ) const; // デフォルトの再帰式であるため、再帰深度が浅い場合にしか使えない。 inline int Moevius( const T& t0 , const T& t1 ); }; template class EnumerationForZetaTransform : virtual public ZetaTransformBody { private: inline int e( const int& i ); inline int e_inv( const int& t ); }; // 入力の範囲内で要件 // (1) Eがint上の半順序<のグラフでE_invが(int,>)のグラフである。 // を満たす場合にのみ以下をサポート。 // 0による初期化O(size_max) // ゼータ変換前の配列による初期化O(始切片のサイズの総和) // 一点加算O(終切片[t,∞)のサイズ) // 全体加算O(size) // 各点加法O(size) // 全体乗算O(size) // (int,<)がjoin半束である場合の畳み込み乗法O(size) // 一点取得O(始切片(-∞,t]のサイズ×メビウス関数の計算量) // 区間和取得O(1) // 逆像和取得O(始切片(-∞,f_inv_max(r_max)]のサイズ×メビウス関数の計算量) // 区間逆像和取得O(1) template E(const int&) , list E_inv(const int&) , int size_max> class ZetaTransform : public PartiallyOrderedSetForZetaTransform , public EnumerationForZetaTransform { public: inline ZetaTransform( const int& size ); inline ZetaTransform( const int& size , const ll ( &a )[size_max] ); private: inline const ll& Zero() const; inline ll Sum( const ll& u0 , const ll& u1 ) const; inline ll Prod( const ll& u0 , const ll& u1 ) const; }; // 入力の範囲内で要件 // (1) 2^digit <= size_max // (2) (U,a_U:U^2->U,z_U:1->U,m_U:U^2->U)が単位的環である。 // を満たす場合にのみ以下をサポート。 // z_U()による初期化O(size_max) // ゼータ変換前の配列による初期化O(digit 2^digit)(可換加法モノイド性を使う) // 一点加算O(2^digit)(可換加法モノイド性を使う) // 全体加算O(2^digit)(可換加法モノイド性だけでも実装できるが単位的半環性を使う) // 各点加法O(2^digit)(加法モノイド性を使う) // 全体乗算O(2^digit)(半環性を使う) // 畳み込み乗法O(2^digit)(半環性を使う) // 一点取得O(2^digit)(単位的環性を使う。愚直と同じオーダー) // 多点取得O(digit 2^digit)(単位的環性を使う) // 区間和取得O(1)(可換加法モノイド性を使う) // 逆像和取得O(始切片(-∞,f_inv_max(r_max)]のサイズ)(単位的環性を使う) // 区間逆像和取得O(1)(可換加法モノイド性を使う) template class FastZetaTransform : public SemiRingForZetaTransform , public EnumerationForZetaTransform { private: int m_digit; public: inline FastZetaTransform( const int& digit ); inline FastZetaTransform( const int& digit , const U ( &a )[size_max] ); inline FastZetaTransform& operator+=( const FastZetaTransform& a ); inline FastZetaTransform& operator*=( const FastZetaTransform& a ); inline void FastMoeviusTransform( U ( &a )[size_max] ); private: inline list Sub( const int& t ) const; inline list Sup( const int& t ) const; inline int Moevius( const int& t0 , const int& t1 ); }; // 入力の範囲内で要件 // (1) EがT上の半順序<のグラフでE_invが(T,>)のグラフである。 // (2) (U,a_U:U^2->U,z_U:1->U,m_U:U^2->U)が単位的環である。 // を満たす場合にのみ以下をサポート。 // z_U()による初期化O(size_max) // 一点加算O(終切片[t,∞)のサイズ×log_2(size))(可換加法モノイド性を使う) // 全体加算O(size)(可換加法モノイド性だけでも実装できるが単位的半環性を使う) // 各点加法O(size)(加法モノイド性を使う) // 全体乗算O(size)(半環性を使う) // (T,<)がjoin半束である場合の畳み込み乗法O(size)(半環性を使う) // 一点取得O(始切片(-∞,t]のサイズ×メビウス関数の計算量×log_2(size))(単位的環性を使う) // 区間和取得O(log_2(size))(可換加法モノイド性を使う) // 逆像和取得O(始切片(-∞,f_inv_max(r_max)]のサイズ×メビウス関数の計算量×log_2(size))(単位的環性を使う) // 区間逆像和取得O(log_2(size))(可換加法モノイド性を使う) template E(const T&) , list E_inv(const T&) , typename U , U a_U(const U&,const U&) , const U& z_U() , U m_U(const U&,const U&) , int size_max> class MemorisationZetaTransform : public SemiRingForZetaTransform , public PartiallyOrderedSetForZetaTransform { public: inline MemorisationZetaTransform( const int& size ); }; // 入力の範囲内で要件 // (1) EがT上の半順序<のグラフでE_invが(T,>)のグラフである。 // (2) (U,a_U:U^2->U,z_U:1->U,m_U:U^2->U)が単位的環である。 // (3) enum_T:int->Tとenum_T_inv:int->Tが互いに逆写像である。 // を満たす場合にのみ以下をサポート。 // z_U()による初期化O(size_max) // 一点加算O(終切片[t,∞)のサイズ)(可換加法モノイド性を使う) // 全体加算O(size)(可換加法モノイド性だけでも実装できるが単位的半環性を使う) // 各点加法O(size)(加法モノイド性を使う) // 全体乗算O(size)(半環性を使う) // (T,<)がjoin半束である場合の畳み込み乗法O(size)(半環性を使う) // 一点取得O(始切片(-∞,t]のサイズ×メビウス関数の計算量)(単位的環性を使う) // 区間和取得O(1)(可換加法モノイド性を使う) // 逆像和取得O(始切片(-∞,f_inv_max(r_max)]のサイズ×メビウス関数の計算量)(単位的環性を使う) // 区間逆像和取得O(1)(可換加法モノイド性を使う) template E(const T&) , list E_inv(const T&) , typename U , U a_U(const U&,const U&) , const U& z_U() , U m_U(const U&,const U&) , int size_max , T enum_T(const int&) , int enum_T_inv(const T&)> class EnumerationZetaTransform : public SemiRingForZetaTransform , public PartiallyOrderedSetForZetaTransform { public: inline EnumerationZetaTransform( const int& size ); private: inline T e( const int& i ); inline int e_inv( const T& t ); }; template inline ZetaTransformBody::ZetaTransformBody( const int& size ) : m_length() , m_memory() , m_memory_inv() , m_size( size ) , m_val() {} template inline SemiRingForZetaTransform::SemiRingForZetaTransform( const int& dummy ) : ZetaTransformBody( dummy ) { using base = ZetaTransformBody; const U& zero = z_U(); if( base::m_val[0] != zero ){ for( int i = 0 ; i < base::m_size ; i++ ){ base::m_val[i] = zero; } } } template E(const int&) , list E_inv(const int&) , int size_max> inline ZetaTransform::ZetaTransform( const int& size ) : ZetaTransformBody( size ) , PartiallyOrderedSetForZetaTransform() {} template E(const int&) , list E_inv(const int&) , int size_max> inline ZetaTransform::ZetaTransform( const int& size , const ll ( &a )[size_max] ) : ZetaTransformBody( size ) , PartiallyOrderedSetForZetaTransform() , EnumerationForZetaTransform() { using base = ZetaTransformBody; for( int i = 0 ; i < base::m_size ; i++ ){ ll& m_val_i = base::m_val[i]; list sub_i = E( i ); while( ! sub_i.empty() ){ m_val_i += a[sub_i.front()]; sub_i.pop_front(); } } } template inline FastZetaTransform::FastZetaTransform( const int& digit ) : ZetaTransformBody( size ) , SemiRingForZetaTransform( size ) , EnumerationForZetaTransform() {} template inline FastZetaTransform::FastZetaTransform( const int& digit , const U ( &a )[size_max] ) : ZetaTransformBody( 1 << digit ) , SemiRingForZetaTransform( size ) , m_digit( digit ) { using base = ZetaTransformBody; for( int i = 0 ; i < base::m_size ; i++ ){ base::m_val[i] = a[i]; } int power = 1; for( int d = 0 ; d < m_digit ; d++ , power <<= 1 ){ for( int i = 0 ; i < base::m_size ; i++ ){ if( ( i & power ) != 0 ){ U& m_val_i = base::m_val[i]; m_val_i = a_U( m_val_i , base::m_val[i ^ power] ); } } } } template E(const T&) , list E_inv(const T&) , typename U , U a_U(const U&,const U&) , const U& z_U() , U m_U(const U&,const U&) , int size_max> inline MemorisationZetaTransform::MemorisationZetaTransform( const int& size ) : ZetaTransformBody( size ) , SemiRingForZetaTransform( size ) , PartiallyOrderedSetForZetaTransform() {} template E(const T&) , list E_inv(const T&) , typename U , U a_U(const U&,const U&) , const U& z_U() , U m_U(const U&,const U&) , int size_max , T enum_T(const int&) , int enum_T_inv(const T&)> inline EnumerationZetaTransform::EnumerationZetaTransform( const int& size ) : ZetaTransformBody( size ) , SemiRingForZetaTransform( size ) , PartiallyOrderedSetForZetaTransform() {} template inline void ZetaTransformBody::Add( const T& t , const U& u ) { list sup = Sup( t ); sup.push_front( t ); while( ! sup.empty() ){ U& m_val_i = m_val[e_inv( sup.front() )]; m_val_i = Sum( m_val_i , u ); sup.pop_front(); } return; } template inline void ZetaTransformBody::AddAll( const U& u ) { for( int i = 0 ; i < m_size ; i++ ){ U& m_val_i = m_val[i]; m_val_i = Sum( m_val_i , Prod( Sub( e( i ) ).size() , u ) ); } return; } template inline ZetaTransformBody& ZetaTransformBody::operator+=( const ZetaTransformBody& a ) { for( int i = 0 ; i < m_size ; i++ ){ U& m_val_i = m_val[i]; m_val_i = Sum( m_val_i , a.m_val[i] ); } return *this; } template inline void ZetaTransformBody::MultiplyAll( const U& u ) { for( int i = 0 ; i < m_size ; i++ ){ U& m_val_i = m_val[i]; m_val_i = Prod( m_val_i , u ); } return; } template inline ZetaTransformBody& ZetaTransformBody::operator*=( const ZetaTransformBody& a ) { for( int i = 0 ; i < m_size ; i++ ){ U& m_val_i = m_val[i]; m_val_i = Prod( m_val_i , a.m_val[i] ); } return *this; } template inline U ZetaTransformBody::Get( const T& t ) { DEFINITION_OF_GET_FOR_ZETA_TRANSFORM( Moevius ); } template template inline U ZetaTransformBody::Get( const T& t ) { DEFINITION_OF_GET_FOR_ZETA_TRANSFORM( mu ); } template inline const U& ZetaTransformBody::InitialSegmentSum( const T& t ) { return m_val[e_inv( t )]; } template template r(const S&)> inline U ZetaTransformBody::InverseImageSum( const S& s ) { DEFINITION_OF_INVERSE_IMAGE_SUM_FOR_ZETA_TRANSFORM( Moevius ); } template template r(const S&) , int mu(const T&,const T&)> inline U ZetaTransformBody::InverseImageSum( const S& s ) { DEFINITION_OF_INVERSE_IMAGE_SUM_FOR_ZETA_TRANSFORM( mu ); } template template inline const U& ZetaTransformBody::IntervalInverseImageSum( const S& s ) { return m_val[e_inv( f_inv_max( s ) )]; } template inline void FastZetaTransform::FastMoeviusTransform( U ( &a )[size_max] ) { using base = ZetaTransformBody; for( int i = 0 ; i < base::m_size ; i++ ){ a[i] = base::m_val[i]; } int power = 1; for( int d = 0 ; d < m_digit ; d++ , power <<= 1 ){ for( int i = 0 ; i < base::m_size ; i++ ){ if( ( i & power ) != 0 ){ U& a_i = a[i]; a_i = a_U( a_i , m_U( -1 , a[i ^ power] ) ); } } } } template T ZetaTransformBody::e( const int& i ) { assert( i < m_length ); return m_memory_inv[i]; } template inline int EnumerationForZetaTransform::e( const int& i ) { return i; } template E(const T&) , list E_inv(const T&) , typename U , U a_U(const U&,const U&) , const U& z_U() , U m_U(const U&,const U&) , int size_max , T enum_T(const int&) , int enum_T_inv(const T&)> inline T EnumerationZetaTransform::e( const int& i ) { return enum_T( i ); } template int ZetaTransformBody::e_inv( const T& t ) { if( m_memory.count( t ) == 0 ){ assert( m_length < m_size ); m_memory_inv.push_back( t ); return m_memory[t] = m_length++; } return m_memory[t]; } template inline int EnumerationForZetaTransform::e_inv( const int& t ) { return t; } template E(const T&) , list E_inv(const T&) , typename U , U a_U(const U&,const U&) , const U& z_U() , U m_U(const U&,const U&) , int size_max , T enum_T(const int&) , int enum_T_inv(const T&)> inline int EnumerationZetaTransform::e_inv( const T& t ) { return enum_T_inv( t ); } template inline const U& SemiRingForZetaTransform::Zero() const { return z_U(); } template E(const int&) , list E_inv(const int&) , int size_max> inline const ll& ZetaTransform::Zero() const { static const ll zero = 0; return zero; } template inline U SemiRingForZetaTransform::Sum( const U& u0 , const U& u1 ) const { return a_U( u0 , u1 ); } template E(const int&) , list E_inv(const int&) , int size_max> inline ll ZetaTransform::Sum( const ll& u0 , const ll& u1 ) const { return u0 + u1; } template inline U SemiRingForZetaTransform::Prod( const U& u0 , const U& u1 ) const { return m_U( u0 , u1 ); } template E(const int&) , list E_inv(const int&) , int size_max> inline ll ZetaTransform::Prod( const ll& u0 , const ll& u1 ) const { return u0 * u1; } template E(const T&) , list E_inv(const T&) , typename U , int size_max> inline list PartiallyOrderedSetForZetaTransform::Sub( const T& t ) const { return E( t ); } template inline list FastZetaTransform::Sub( const int& t ) const { list sub{}; sub.push_back( t ); int sub_size = 1; int power = 1; for( int d = 0 ; d < m_digit ; d++ , power <<= 1 ){ if( ( t & power ) != 0 ){ auto itr = sub.begin(); for( int i = 0 ; i < sub_size ; i++ , itr++ ){ sub.push_back( *itr ^ power ); } sub_size <<= 1; } } return sub; } template E(const T&) , list E_inv(const T&) , typename U , int size_max> inline list PartiallyOrderedSetForZetaTransform::Sup( const T& t ) const { return E_inv( t ); } template inline list FastZetaTransform::Sup( const int& t ) const { list sup{}; sup.push_back( t ); int sup_size = 1; int power = 1; for( int d = 0 ; d < m_digit ; d++ , power <<= 1 ){ if( ( t & power ) == 0 ){ auto itr = sup.begin(); for( int i = 0 ; i < sup_size ; i++ , itr++ ){ sup.push_back( *itr | power ); } sup_size <<= 1; } } return sup; } template E(const T&) , list E_inv(const T&) , typename U , int size_max> inline int PartiallyOrderedSetForZetaTransform::Moevius( const T& t0 , const T& t1 ) { using base = ZetaTransformBody; const int i = base::e_inv( t0 ); const int j = base::e_inv( t1 ); map& moevius_t0 = m_moevius[i]; bool found = moevius_t0.count( j ) == 1; int& answer = moevius_t0[j]; if( ! found ){ if( i == j ){ answer = 1; } else { list sub = E( t1 ); while( ! sub.empty() ){ answer -= Moevius( t0 , sub.front() ); sub.pop_front(); } } } return answer; } template int FastZetaTransform::Moevius( const int& t0 , const int& t1 ) { int t = t1 ^ t0; int count = 0; while( t != 0 ){ t -= t & -t; count++; } return ( count & 1 ) == 0 ? 1 : -1; } inline DEXPR( int , bound_Ai , 1000000 , 100 ); inline list E( const int& i ) { cerr << "ここが呼ばれるのは意図通りではありません。" << endl; abort(); list answer{}; int n = 0; while( n <= bound_Ai ){ answer.push_back( n ); n += i; } return answer; } template class PrimeEnumeration { public: INT m_val[length_max]; int m_length; inline constexpr PrimeEnumeration(); }; template inline constexpr PrimeEnumeration::PrimeEnumeration() : m_val() , m_length( 0 ) { bool is_comp[val_limit] = {}; for( INT i = 2 ; i < val_limit ; i++ ){ if( is_comp[i] == false ){ INT j = i; while( ( j += i ) < val_limit ){ is_comp[j] = true; } m_val[m_length++] = i; if( m_length >= length_max ){ break; } } } } template void SetPrimeFactorisation( const PrimeEnumeration& prime , const INT& n , vector& P , vector& exponent ) { INT n_copy = n; int i = 0; while( i < prime.m_length ){ const INT& p = prime.m_val[i]; if( p * p > n_copy ){ break; } if( n_copy % p == 0 ){ P.push_back( p ); exponent.push_back( 1 ); INT& exponent_back = exponent.back(); n_copy /= p; while( n_copy % p == 0 ){ exponent_back++; n_copy /= p; } } i++; } if( n_copy != 1 ){ P.push_back( n_copy ); exponent.push_back( 1 ); } return; } template void SetPrimeFactorisation( const INT& n , vector& P , vector& exponent , vector& P_power ) { INT n_copy = n; INT p = 2; if( n_copy % p == 0 ){ P.push_back( p ); exponent.push_back( 1 ); P_power.push_back( p ); INT& exponent_back = exponent.back(); INT& P_power_back = P_power.back(); n_copy /= p; while( n_copy % p == 0 ){ exponent_back++; P_power_back *= p; n_copy /= p; } } p++; while( p * p <= n_copy ){ if( n_copy % p == 0 ){ P.push_back( p ); exponent.push_back( 1 ); P_power.push_back( p ); INT& exponent_back = exponent.back(); INT& P_power_back = P_power.back(); n_copy /= p; while( n_copy % p == 0 ){ exponent_back++; P_power_back *= p; n_copy /= p; } } p += 2; } if( n_copy != 1 ){ P.push_back( n_copy ); exponent.push_back( 1 ); } return; } // sqrt( n )以下の最大の素数 < val_limit の時のみサポート。 template list EnumerateDivisor( const PrimeEnumeration& prime , INT n ) noexcept { vector P{}; vector exponent{}; SetPrimeFactorisation( prime , n , P , exponent ); const int length = P.size(); list divisor{}; divisor.push_back( 1 ); auto begin = divisor.begin() , end = divisor.end(); for( int i = 0 ; i < length ; i++ ){ const INT& P_i = P[i]; const int& exponent_i = exponent[i]; list temp{}; INT power = 1; for( int e = 1 ; e <= exponent_i ; e++ ){ power *= P_i; for( auto itr = begin ; itr != end ; itr++ ){ temp.push_back( *itr * power ); } } while( ! temp.empty() ){ divisor.push_back( temp.front() ); temp.pop_front(); } } return divisor; } inline CEXPR( int , sqrt_bound_Ai , 1000 ); inline constexpr PrimeEnumeration prime{}; inline list E_inv( const int& i ) { list answer = EnumerateDivisor( prime , i ); answer.pop_back(); return answer; } int Ai; inline int id( const int& i ) { return i; } inline list r( const int& n ) { vector P{}; vector exponent{}; SetPrimeFactorisation( prime , Ai , P , exponent ); const int length = P.size(); list sqfr_divisor{}; sqfr_divisor.push_back( 1 ); auto begin = sqfr_divisor.begin() , end = sqfr_divisor.end(); for( int i = 0 ; i < length ; i++ ){ const int& P_i = P[i]; list temp{}; for( auto itr = begin ; itr != end ; itr++ ){ temp.push_back( *itr * P_i ); } while( ! temp.empty() ){ sqfr_divisor.push_back( temp.front() ); temp.pop_front(); } } return sqfr_divisor; } template int MoeviusFunction( const PrimeEnumeration& prime , INT n ) noexcept { static int memory[bound_Ai + 1] = {}; static bool uninitialised = true; if( uninitialised ){ for( int i = 0 ; i < bound_Ai + 1 ; i++ ){ memory[i] = -2; } uninitialised = false; } int& answer = memory[n]; if( answer == -2 ){ answer = 1; int i = 0; while( i < prime.m_length && n > 1 ){ const int& p = prime.m_val[i]; if( n % p == 0 ){ n /= p; if( n % p == 0 ){ return answer = 0; } answer *= -1; } i++; } n == 1 ? answer : answer *= -1; } return answer; } inline int mu( const int& t0 , const int& t1 ) { return MoeviusFunction( prime , t0 / t1 ); } US ull = unsigned long long;IN CEXPR(uint,P,998244353);TE IN CE INT& RS(INT& n)NE{RE n < 0?((((++n)*= -1)%= M)*= -1)+= M - 1:n %= M;}TE IN CE uint& RS(uint& n)NE{RE n %= M;}TE IN CE ull& RS(ull& n)NE{RE n %= M;}TE IN CE INT& RSP(INT& n)NE{CE CO uint trunc = (1 << 23)- 1;INT n_u = n >> 23;n &= trunc;INT n_uq = (n_u / 7)/ 17;n_u -= n_uq * 119;n += n_u << 23;RE n < n_uq?n += P - n_uq:n -= n_uq;}TE <> IN CE ull& RS(ull& n)NE{CE CO ull Pull = P;CE CO ull Pull2 = (Pull - 1)* (Pull - 1);RE RSP(n > Pull2?n -= Pull2:n);}TE IN CE INT RS(INT&& n)NE{RE MO(RS(n));}TE IN CE INT RS(CO INT& n)NE{RE RS(INT(n));} #define SFINAE_FOR_MOD(DEFAULT)TY T,enable_if_t >::value>* DEFAULT #define DC_OF_CM_FOR_MOD(FUNC)IN bool OP FUNC(CO Mod& n)CO NE #define DC_OF_AR_FOR_MOD(FUNC)IN Mod OP FUNC(CO Mod& n)CO NE;TE IN Mod OP FUNC(T&& n)CO NE; #define DF_OF_CM_FOR_MOD(FUNC)TE IN bool Mod::OP FUNC(CO Mod& n)CO NE{RE m_n FUNC n.m_n;} #define DF_OF_AR_FOR_MOD(FUNC,FORMULA)TE IN Mod Mod::OP FUNC(CO Mod& n)CO NE{RE MO(Mod(*TH)FUNC ## = n);}TE TE IN Mod Mod::OP FUNC(T&& n)CO NE{RE FORMULA;}TE IN Mod OP FUNC(T&& n0,CO Mod& n1)NE{RE MO(Mod(forward(n0))FUNC ## = n1);} TE CL Mod{PU:uint m_n;IN CE Mod()NE;IN CE Mod(CO Mod& n)NE;IN CE Mod(Mod& n)NE;IN CE Mod(Mod&& n)NE;TE IN CE Mod(CO T& n)NE;TE IN CE Mod(T& n)NE;TE IN CE Mod(T&& n)NE;IN CE Mod& OP=(CO Mod& n)NE;IN CE Mod& OP=(Mod&& n)NE;IN CE Mod& OP+=(CO Mod& n)NE;IN CE Mod& OP-=(CO Mod& n)NE;IN CE Mod& OP*=(CO Mod& n)NE;IN Mod& OP/=(CO Mod& n);IN CE Mod& OP<<=(int n)NE;IN CE Mod& OP>>=(int n)NE;IN CE Mod& OP++()NE;IN CE Mod OP++(int)NE;IN CE Mod& OP--()NE;IN CE Mod OP--(int)NE;DC_OF_CM_FOR_MOD(==);DC_OF_CM_FOR_MOD(!=);DC_OF_CM_FOR_MOD(<);DC_OF_CM_FOR_MOD(<=);DC_OF_CM_FOR_MOD(>);DC_OF_CM_FOR_MOD(>=);DC_OF_AR_FOR_MOD(+);DC_OF_AR_FOR_MOD(-);DC_OF_AR_FOR_MOD(*);DC_OF_AR_FOR_MOD(/);IN CE Mod OP<<(int n)CO NE;IN CE Mod OP>>(int n)CO NE;IN CE Mod OP-()CO NE;IN CE Mod& SignInvert()NE;IN CE Mod& Double()NE;IN CE Mod& Halve()NE;IN Mod& Invert();TE IN CE Mod& PositivePW(T&& EX)NE;TE IN CE Mod& NonNegativePW(T&& EX)NE;TE IN CE Mod& PW(T&& EX);IN CE VO swap(Mod& n)NE;IN CE CO uint& RP()CO NE;ST IN CE Mod DeRP(CO uint& n)NE;ST IN CE uint& Normalise(uint& n)NE;ST IN CO Mod& Inverse(CO uint& n)NE;ST IN CO Mod& Factorial(CO uint& n)NE;ST IN CO Mod& FactorialInverse(CO uint& n)NE;ST IN CO Mod& zero()NE;ST IN CO Mod& one()NE;TE IN CE Mod& Ref(T&& n)NE;}; #define SFINAE_FOR_MN(DEFAULT)TY T,enable_if_t,decay_t >::value>* DEFAULT #define DC_OF_AR_FOR_MN(FUNC)IN MN OP FUNC(CO MN& n)CO NE;TE IN MN OP FUNC(T&& n)CO NE; #define DF_OF_CM_FOR_MN(FUNC)TE IN bool MN::OP FUNC(CO MN& n)CO NE{RE m_n FUNC n.m_n;} #define DF_OF_AR_FOR_MN(FUNC,FORMULA)TE IN MN MN::OP FUNC(CO MN& n)CO NE{RE MO(MN(*TH)FUNC ## = n);}TE TE IN MN MN::OP FUNC(T&& n)CO NE{RE FORMULA;}TE IN MN OP FUNC(T&& n0,CO MN& n1)NE{RE MO(MN(forward(n0))FUNC ## = n1);} TE CL MN :PU Mod{PU:IN CE MN()NE;IN CE MN(CO MN& n)NE;IN CE MN(MN& n)NE;IN CE MN(MN&& n)NE;TE IN CE MN(CO T& n)NE;TE IN CE MN(T&& n)NE;IN CE MN& OP=(CO MN& n)NE;IN CE MN& OP=(MN&& n)NE;IN CE MN& OP+=(CO MN& n)NE;IN CE MN& OP-=(CO MN& n)NE;IN CE MN& OP*=(CO MN& n)NE;IN MN& OP/=(CO MN& n);IN CE MN& OP<<=(int n)NE;IN CE MN& OP>>=(int n)NE;IN CE MN& OP++()NE;IN CE MN OP++(int)NE;IN CE MN& OP--()NE;IN CE MN OP--(int)NE;DC_OF_AR_FOR_MN(+);DC_OF_AR_FOR_MN(-);DC_OF_AR_FOR_MN(*);DC_OF_AR_FOR_MN(/);IN CE MN OP<<(int n)CO NE;IN CE MN OP>>(int n)CO NE;IN CE MN OP-()CO NE;IN CE MN& SignInvert()NE;IN CE MN& Double()NE;IN CE MN& Halve()NE;IN CE MN& Invert();TE IN CE MN& PositivePW(T&& EX)NE;TE IN CE MN& NonNegativePW(T&& EX)NE;TE IN CE MN& PW(T&& EX);IN CE uint RP()CO NE;IN CE Mod Reduce()CO NE;ST IN CE MN DeRP(CO uint& n)NE;ST IN CO MN& Formise(CO uint& n)NE;ST IN CO MN& Inverse(CO uint& n)NE;ST IN CO MN& Factorial(CO uint& n)NE;ST IN CO MN& FactorialInverse(CO uint& n)NE;ST IN CO MN& zero()NE;ST IN CO MN& one()NE;ST IN CE uint Form(CO uint& n)NE;ST IN CE ull& Reduction(ull& n)NE;ST IN CE ull& ReducedMU(ull& n,CO uint& m)NE;ST IN CE uint MU(CO uint& n0,CO uint& n1)NE;ST IN CE uint BaseSquareTruncation(uint& n)NE;TE IN CE MN& Ref(T&& n)NE;};TE IN CE MN Twice(CO MN& n)NE;TE IN CE MN Half(CO MN& n)NE;TE IN CE MN Inverse(CO MN& n);TE IN CE MN PW(CO MN& n,CO T& EX);TE IN CE MN<2> PW(CO MN<2>& n,CO T& p);TE IN CE T Square(CO T& t);TE <> IN CE MN<2> Square >(CO MN<2>& t);TE IN CE VO swap(MN& n0,MN& n1)NE;TE IN string to_string(CO MN& n)NE;TE IN basic_ostream& OP<<(basic_ostream& os,CO MN& n); TE CL COantsForMod{PU:COantsForMod()= delete;ST CE CO bool g_even = ((M & 1)== 0);ST CE CO uint g_memory_bound = 1000000;ST CE CO uint g_memory_LE = M < g_memory_bound?M:g_memory_bound;ST IN CE ull MNBasePW(ull&& EX)NE;ST CE uint g_M_minus = M - 1;ST CE uint g_M_minus_2 = M - 2;ST CE uint g_M_minus_2_neg = 2 - M;ST CE CO int g_MN_digit = 32;ST CE CO ull g_MN_base = ull(1)<< g_MN_digit;ST CE CO uint g_MN_base_minus = uint(g_MN_base - 1);ST CE CO uint g_MN_digit_half = (g_MN_digit + 1)>> 1;ST CE CO uint g_MN_base_sqrt_minus = (1 << g_MN_digit_half)- 1;ST CE CO uint g_MN_M_neg_inverse = uint((g_MN_base - MNBasePW((ull(1)<< (g_MN_digit - 1))- 1))& g_MN_base_minus);ST CE CO uint g_MN_base_mod = uint(g_MN_base % M);ST CE CO uint g_MN_base_square_mod = uint(((g_MN_base % M)* (g_MN_base % M))% M);};TE IN CE ull COantsForMod::MNBasePW(ull&& EX)NE{ull prod = 1;ull PW = M;WH(EX != 0){(EX & 1)== 1?(prod *= PW)&= g_MN_base_minus:prod;EX >>= 1;(PW *= PW)&= g_MN_base_minus;}RE prod;} #include #define SET_VE_32_128_FOR_SIMD(UINT,VE_NAME,SCALAR0,SCALAR1,SCALAR2,SCALAR3)CE CO UINT VE_NAME ## _copy[4] ={SCALAR0,SCALAR1,SCALAR2,SCALAR3};ST CO __m128i v_ ## VE_NAME = _mm_load_si128((__m128i*)VE_NAME ##_copy) #define SET_VE_64_128_FOR_SIMD(UINT,VE_NAME,SCALAR0,SCALAR1)CE CO UINT VE_NAME ## _copy[2] ={SCALAR0,SCALAR1};ST CO __m128i v_ ## VE_NAME = _mm_load_si128((__m128i*)VE_NAME ##_copy) #define SET_VE_64_256_FOR_SIMD(ULL,VE_NAME,SCALAR0,SCALAR1,SCALAR2,SCALAR3)CE CO ULL VE_NAME ## _copy[4] ={SCALAR0,SCALAR1,SCALAR2,SCALAR3};ST CO __m256i v_ ## VE_NAME = _mm256_load_si256((__m256i*)VE_NAME ##_copy) #define SET_CO_VE_32_128_FOR_SIMD(UINT,VE_NAME,SCALAR)SET_VE_32_128_FOR_SIMD(UINT,VE_NAME,SCALAR,SCALAR,SCALAR,SCALAR) #define SET_CO_VE_64_128_FOR_SIMD(ULL,VE_NAME,SCALAR)SET_VE_64_128_FOR_SIMD(ULL,VE_NAME,SCALAR,SCALAR) #define SET_CO_VE_64_256_FOR_SIMD(ULL,VE_NAME,SCALAR)SET_VE_64_256_FOR_SIMD(ULL,VE_NAME,SCALAR,SCALAR,SCALAR,SCALAR) TE CL COantsForSIMDForMod{PU:COantsForSIMDForMod()= delete;ST IN CO __m128i& v_M()NE;ST IN CO __m128i& v_Mull()NE;ST IN CO __m128i& v_M_minus()NE;ST IN CO __m128i& v_M_neg_inverse()NE;ST IN CO __m128i& v_digitull()NE;};TE IN CO __m128i& COantsForSIMDForMod::v_M()NE{SET_CO_VE_32_128_FOR_SIMD(uint,M,M);RE v_M;}TE IN CO __m128i& COantsForSIMDForMod::v_Mull()NE{SET_CO_VE_64_128_FOR_SIMD(ull,Mull,M);RE v_Mull;}TE IN CO __m128i& COantsForSIMDForMod::v_M_minus()NE{SET_CO_VE_32_128_FOR_SIMD(uint,M_minus,M - 1);RE v_M_minus;}TE IN CO __m128i& COantsForSIMDForMod::v_M_neg_inverse()NE{SET_CO_VE_32_128_FOR_SIMD(uint,M_neg_inverse,COantsForMod::g_MN_M_neg_inverse);RE v_M_neg_inverse;}TE IN CO __m128i& COantsForSIMDForMod::v_digitull()NE{SET_CO_VE_64_128_FOR_SIMD(ull,digitull,COantsForMod::g_MN_digit);RE v_digitull;}TE IN __m128i& SIMD_RS_32_128(__m128i& v)NE{CO __m128i& v_M = COantsForSIMDForMod::v_M();RE v -= v_M * _mm_cmpgt_epi32(v,v_M);}TE IN __m128i& SIMD_RS_64_128(__m128i& v)NE{ull v_copy[2];_mm_store_si128((__m128i*)v_copy,v);for(uint i = 0;i < 2;i++){ull& v_copy_i = v_copy[i];v_copy_i = (v_copy_i < M?0:M);}RE v -= _mm_load_si128((__m128i*)v_copy);}TE IN __m256i& SIMD_RS_64_256(__m256i& v)NE{ull v_copy[4];_mm256_store_si256((__m256i*)v_copy,v);for(uint i = 0;i < 4;i++){ull& v_copy_i = v_copy[i];v_copy_i = (v_copy_i < M?0:M);}RE v -= _mm256_load_si256((__m256i*)v_copy);}IN CE int SIMD_Shuffle(CRI a0,CRI a1,CRI a2,CRI a3)NE{RE (a0 << (0 << 1))+ (a1 << (1 << 1))+ (a2 << (2 << 1))+ (a3 << (3 << 1));}TE IN VO SIMD_Addition_32_64(CO Mod& a0,CO Mod& a1,CO Mod& b0,CO Mod& b1,Mod& c0,Mod& c1)NE{uint a_copy[4] ={a0.m_n,a1.m_n,0,0};uint b_copy[4] ={b0.m_n,b1.m_n,0,0};__m128i v_a = _mm_load_si128((__m128i*)a_copy);v_a += _mm_load_si128((__m128i*)b_copy);ST CO __m128i& v_M_minus = COantsForSIMDForMod::v_M_minus();ST CO __m128i& v_M = COantsForSIMDForMod::v_M();v_a += _mm_cmpgt_epi32(v_a,v_M_minus)& v_M;_mm_store_si128((__m128i*)a_copy,v_a);c0.m_n = MO(a_copy[0]);c1.m_n = MO(a_copy[1]);RE;}TE IN VO SIMD_Addition_32_128(CO Mod& a0,CO Mod& a1,CO Mod& a2,CO Mod& a3,CO Mod& b0,CO Mod& b1,CO Mod& b2,CO Mod& b3,Mod& c0,Mod& c1,Mod& c2,Mod& c3)NE{uint a_copy[4] ={a0.m_n,a1.m_n,a2.m_n,a3.m_n};uint b_copy[4] ={b0.m_n,b1.m_n,b2.m_n,b3.m_n};__m128i v_a = _mm_load_si128((__m128i*)a_copy)+ _mm_load_si128((__m128i*)b_copy);_mm_store_si128((__m128i*)a_copy,v_a);for(uint i = 0;i < 4;i++){b_copy[i] = a_copy[i] < M?0:M;}v_a -= _mm_load_si128((__m128i*)b_copy);_mm_store_si128((__m128i*)a_copy,v_a);c0.m_n = MO(a_copy[0]);c1.m_n = MO(a_copy[1]);c2.m_n = MO(a_copy[2]);c3.m_n = MO(a_copy[3]);RE;}TE IN VO SIMD_Substracition_32_64(CO Mod& a0,CO Mod& a1,CO Mod& b0,CO Mod& b1,Mod& c0,Mod& c1)NE{uint a_copy[4] ={a0.m_n,a1.m_n,0,0};uint b_copy[4] ={b0.m_n,b1.m_n,0,0};__m128i v_a = _mm_load_si128((__m128i*)a_copy);__m128i v_b = _mm_load_si128((__m128i*)b_copy);_mm_store_si128((__m128i*)a_copy,v_a);for(uint i = 0;i < 2;i++){b_copy[i] = a_copy[i] < b_copy[i]?M:0;}(v_a += _mm_load_si128((__m128i*)b_copy))-= v_b;_mm_store_si128((__m128i*)a_copy,v_a);c0.m_n = MO(a_copy[0]);c1.m_n = MO(a_copy[1]);RE;}TE IN VO SIMD_Subtraction_32_128(CO Mod& a0,CO Mod& a1,CO Mod& a2,CO Mod& a3,CO Mod& b0,CO Mod& b1,CO Mod& b2,CO Mod& b3,Mod& c0,Mod& c1,Mod& c2,Mod& c3)NE{uint a_copy[4] ={a0.m_n,a1.m_n,a2.m_n,a3.m_n};uint b_copy[4] ={b0.m_n,b1.m_n,b2.m_n,b3.m_n};__m128i v_a = _mm_load_si128((__m128i*)a_copy);__m128i v_b = _mm_load_si128((__m128i*)b_copy);_mm_store_si128((__m128i*)a_copy,v_a);for(uint i = 0;i < 4;i++){b_copy[i] = a_copy[i] < b_copy[i]?M:0;}(v_a += _mm_load_si128((__m128i*)b_copy))-= v_b;_mm_store_si128((__m128i*)a_copy,v_a);c0.m_n = MO(a_copy[0]);c1.m_n = MO(a_copy[1]);c2.m_n = MO(a_copy[2]);c3.m_n = MO(a_copy[3]);RE;} US MP = Mod

;US MNP = MN

;TE IN CE uint MN::Form(CO uint& n)NE{ull n_copy = n;RE uint(MO(Reduction(n_copy *= COantsForMod::g_MN_base_square_mod)));}TE IN CE ull& MN::Reduction(ull& n)NE{ull n_sub = n & COantsForMod::g_MN_base_minus;RE ((n += ((n_sub *= COantsForMod::g_MN_M_neg_inverse)&= COantsForMod::g_MN_base_minus)*= M)>>= COantsForMod::g_MN_digit)< M?n:n -= M;}TE IN CE ull& MN::ReducedMU(ull& n,CO uint& m)NE{RE Reduction(n *= m);}TE IN CE uint MN::MU(CO uint& n0,CO uint& n1)NE{ull n0_copy = n0;RE uint(MO(ReducedMU(ReducedMU(n0_copy,n1),COantsForMod::g_MN_base_square_mod)));}TE IN CE uint MN::BaseSquareTruncation(uint& n)NE{CO uint n_u = n >> COantsForMod::g_MN_digit_half;n &= COantsForMod::g_MN_base_sqrt_minus;RE n_u;}TE IN CE MN::MN()NE:Mod(){static_assert(! COantsForMod::g_even);}TE IN CE MN::MN(CO MN& n)NE:Mod(n){}TE IN CE MN::MN(MN& n)NE:Mod(n){}TE IN CE MN::MN(MN&& n)NE:Mod(MO(n)){}TE TE IN CE MN::MN(CO T& n)NE:Mod(n){static_assert(! COantsForMod::g_even);Mod::m_n = Form(Mod::m_n);}TE TE IN CE MN::MN(T&& n)NE:Mod(forward(n)){static_assert(! COantsForMod::g_even);Mod::m_n = Form(Mod::m_n);}TE IN CE MN& MN::OP=(CO MN& n)NE{RE Ref(Mod::OP=(n));}TE IN CE MN& MN::OP=(MN&& n)NE{RE Ref(Mod::OP=(MO(n)));}TE IN CE MN& MN::OP+=(CO MN& n)NE{RE Ref(Mod::OP+=(n));}TE IN CE MN& MN::OP-=(CO MN& n)NE{RE Ref(Mod::OP-=(n));}TE IN CE MN& MN::OP*=(CO MN& n)NE{ull m_n_copy = Mod::m_n;RE Ref(Mod::m_n = MO(ReducedMU(m_n_copy,n.m_n)));}TE IN MN& MN::OP/=(CO MN& n){RE OP*=(MN(n).Invert());}TE IN CE MN& MN::OP<<=(int n)NE{RE Ref(Mod::OP<<=(n));}TE IN CE MN& MN::OP>>=(int n)NE{RE Ref(Mod::OP>>=(n));}TE IN CE MN& MN::OP++()NE{RE Ref(Mod::Normalise(Mod::m_n += COantsForMod::g_MN_base_mod));}TE IN CE MN MN::OP++(int)NE{MN n{*TH};OP++();RE n;}TE IN CE MN& MN::OP--()NE{RE Ref(Mod::m_n < COantsForMod::g_MN_base_mod?((Mod::m_n += M)-= COantsForMod::g_MN_base_mod):Mod::m_n -= COantsForMod::g_MN_base_mod);}TE IN CE MN MN::OP--(int)NE{MN n{*TH};OP--();RE n;}DF_OF_AR_FOR_MN(+,MN(forward(n))+= *TH);DF_OF_AR_FOR_MN(-,MN(forward(n)).SignInvert()+= *TH);DF_OF_AR_FOR_MN(*,MN(forward(n))*= *TH);DF_OF_AR_FOR_MN(/,MN(forward(n)).Invert()*= *TH);TE IN CE MN MN::OP<<(int n)CO NE{RE MO(MN(*TH)<<= n);}TE IN CE MN MN::OP>>(int n)CO NE{RE MO(MN(*TH)>>= n);}TE IN CE MN MN::OP-()CO NE{RE MO(MN(*TH).SignInvert());}TE IN CE MN& MN::SignInvert()NE{RE Ref(Mod::m_n > 0?Mod::m_n = M - Mod::m_n:Mod::m_n);}TE IN CE MN& MN::Double()NE{RE Ref(Mod::Double());}TE IN CE MN& MN::Halve()NE{RE Ref(Mod::Halve());}TE IN CE MN& MN::Invert(){assert(Mod::m_n > 0);RE PositivePW(uint(COantsForMod::g_M_minus_2));}TE TE IN CE MN& MN::PositivePW(T&& EX)NE{MN PW{*TH};(--EX)%= COantsForMod::g_M_minus_2;WH(EX != 0){(EX & 1)== 1?OP*=(PW):*TH;EX >>= 1;PW *= PW;}RE *TH;}TE TE IN CE MN& MN::NonNegativePW(T&& EX)NE{RE EX == 0?Ref(Mod::m_n = 1):PositivePW(forward(EX));}TE TE IN CE MN& MN::PW(T&& EX){bool neg = EX < 0;assert(!(neg && Mod::m_n == 0));RE neg?PositivePW(forward(EX *= COantsForMod::g_M_minus_2_neg)):NonNegativePW(forward(EX));}TE IN CE uint MN::RP()CO NE{ull m_n_copy = Mod::m_n;RE MO(Reduction(m_n_copy));}TE IN CE Mod MN::Reduce()CO NE{ull m_n_copy = Mod::m_n;RE Mod::DeRP(MO(Reduction(m_n_copy)));}TE IN CE MN MN::DeRP(CO uint& n)NE{RE MN(Mod::DeRP(n));}TE IN CO MN& MN::Formise(CO uint& n)NE{ST MN memory[COantsForMod::g_memory_LE] ={zero(),one()};ST uint LE_curr = 2;WH(LE_curr <= n){memory[LE_curr] = DeRP(LE_curr);LE_curr++;}RE memory[n];}TE IN CO MN& MN::Inverse(CO uint& n)NE{ST MN memory[COantsForMod::g_memory_LE] ={zero(),one()};ST uint LE_curr = 2;WH(LE_curr <= n){memory[LE_curr] = MN(Mod::Inverse(LE_curr));LE_curr++;}RE memory[n];}TE IN CO MN& MN::Factorial(CO uint& n)NE{ST MN memory[COantsForMod::g_memory_LE] ={one(),one()};ST uint LE_curr = 2;ST MN val_curr{one()};MN val_last{one()};WH(LE_curr <= n){memory[LE_curr++] = val_curr *= ++val_last;}RE memory[n];}TE IN CO MN& MN::FactorialInverse(CO uint& n)NE{ST MN memory[COantsForMod::g_memory_LE] ={one(),one()};ST uint LE_curr = 2;ST MN val_curr{one()};MN val_last{one()};WH(LE_curr <= n){memory[LE_curr] = val_curr *= Inverse(LE_curr);LE_curr++;}RE memory[n];}TE IN CO MN& MN::zero()NE{ST CE CO MN z{};RE z;}TE IN CO MN& MN::one()NE{ST CE CO MN o{DeRP(1)};RE o;}TE TE IN CE MN& MN::Ref(T&& n)NE{RE *TH;}TE IN CE MN Twice(CO MN& n)NE{RE MO(MN(n).Double());}TE IN CE MN Half(CO MN& n)NE{RE MO(MN(n).Halve());}TE IN CE MN Inverse(CO MN& n){RE MO(MN(n).Invert());}TE IN CE MN PW(CO MN& n,CO T& EX){RE MO(MN(n).PW(T(EX)));}TE IN CE VO swap(MN& n0,MN& n1)NE{n0.swap(n1);}TE IN string to_string(CO MN& n)NE{RE to_string(n.RP())+ " + MZ";}TE IN basic_ostream& OP<<(basic_ostream& os,CO MN& n){RE os << n.RP();} TE IN CE Mod::Mod()NE:m_n(){}TE IN CE Mod::Mod(CO Mod& n)NE:m_n(n.m_n){}TE IN CE Mod::Mod(Mod& n)NE:m_n(n.m_n){}TE IN CE Mod::Mod(Mod&& n)NE:m_n(MO(n.m_n)){}TE TE IN CE Mod::Mod(CO T& n)NE:m_n(RS(n)){}TE TE IN CE Mod::Mod(T& n)NE:m_n(RS(decay_t(n))){}TE TE IN CE Mod::Mod(T&& n)NE:m_n(RS(forward(n))){}TE IN CE Mod& Mod::OP=(CO Mod& n)NE{RE Ref(m_n = n.m_n);}TE IN CE Mod& Mod::OP=(Mod&& n)NE{RE Ref(m_n = MO(n.m_n));}TE IN CE Mod& Mod::OP+=(CO Mod& n)NE{RE Ref(Normalise(m_n += n.m_n));}TE IN CE Mod& Mod::OP-=(CO Mod& n)NE{RE Ref(m_n < n.m_n?(m_n += M)-= n.m_n:m_n -= n.m_n);}TE IN CE Mod& Mod::OP*=(CO Mod& n)NE{RE Ref(m_n = COantsForMod::g_even?RS(ull(m_n)* n.m_n):MN::MU(m_n,n.m_n));}TE <> IN CE MP& MP::OP*=(CO MP& n)NE{ull m_n_copy = m_n;RE Ref(m_n = MO((m_n_copy *= n.m_n)< P?m_n_copy:RSP(m_n_copy)));}TE IN Mod& Mod::OP/=(CO Mod& n){RE OP*=(Mod(n).Invert());}TE IN CE Mod& Mod::OP<<=(int n)NE{WH(n-- > 0){Normalise(m_n <<= 1);}RE *TH;}TE IN CE Mod& Mod::OP>>=(int n)NE{WH(n-- > 0){((m_n & 1)== 0?m_n:m_n += M)>>= 1;}RE *TH;}TE IN CE Mod& Mod::OP++()NE{RE Ref(m_n < COantsForMod::g_M_minus?++m_n:m_n = 0);}TE IN CE Mod Mod::OP++(int)NE{Mod n{*TH};OP++();RE n;}TE IN CE Mod& Mod::OP--()NE{RE Ref(m_n == 0?m_n = COantsForMod::g_M_minus:--m_n);}TE IN CE Mod Mod::OP--(int)NE{Mod n{*TH};OP--();RE n;}DF_OF_CM_FOR_MOD(==);DF_OF_CM_FOR_MOD(!=);DF_OF_CM_FOR_MOD(>);DF_OF_CM_FOR_MOD(>=);DF_OF_CM_FOR_MOD(<);DF_OF_CM_FOR_MOD(<=);DF_OF_AR_FOR_MOD(+,Mod(forward(n))+= *TH);DF_OF_AR_FOR_MOD(-,Mod(forward(n)).SignInvert()+= *TH);DF_OF_AR_FOR_MOD(*,Mod(forward(n))*= *TH);DF_OF_AR_FOR_MOD(/,Mod(forward(n)).Invert()*= *TH);TE IN CE Mod Mod::OP<<(int n)CO NE{RE MO(Mod(*TH)<<= n);}TE IN CE Mod Mod::OP>>(int n)CO NE{RE MO(Mod(*TH)>>= n);}TE IN CE Mod Mod::OP-()CO NE{RE MO(Mod(*TH).SignInvert());}TE IN CE Mod& Mod::SignInvert()NE{RE Ref(m_n > 0?m_n = M - m_n:m_n);}TE IN CE Mod& Mod::Double()NE{RE Ref(Normalise(m_n <<= 1));}TE IN CE Mod& Mod::Halve()NE{RE Ref(((m_n & 1)== 0?m_n:m_n += M)>>= 1);}TE IN Mod& Mod::Invert(){assert(m_n > 0);uint m_n_neg;RE m_n < COantsForMod::g_memory_LE?Ref(m_n = Inverse(m_n).m_n):(m_n_neg = M - m_n < COantsForMod::g_memory_LE)?Ref(m_n = M - Inverse(m_n_neg).m_n):PositivePW(uint(COantsForMod::g_M_minus_2));}TE <> IN Mod<2>& Mod<2>::Invert(){assert(m_n > 0);RE *TH;}TE TE IN CE Mod& Mod::PositivePW(T&& EX)NE{Mod PW{*TH};EX--;WH(EX != 0){(EX & 1)== 1?OP*=(PW):*TH;EX >>= 1;PW *= PW;}RE *TH;}TE <> TE IN CE Mod<2>& Mod<2>::PositivePW(T&& EX)NE{RE *TH;}TE TE IN CE Mod& Mod::NonNegativePW(T&& EX)NE{RE EX == 0?Ref(m_n = 1):Ref(PositivePW(forward(EX)));}TE TE IN CE Mod& Mod::PW(T&& EX){bool neg = EX < 0;assert(!(neg && m_n == 0));neg?EX *= COantsForMod::g_M_minus_2_neg:EX;RE m_n == 0?*TH:(EX %= COantsForMod::g_M_minus)== 0?Ref(m_n = 1):PositivePW(forward(EX));}TE IN CO Mod& Mod::Inverse(CO uint& n)NE{ST Mod memory[COantsForMod::g_memory_LE] ={zero(),one()};ST uint LE_curr = 2;WH(LE_curr <= n){memory[LE_curr].m_n = M - MN::MU(memory[M % LE_curr].m_n,M / LE_curr);LE_curr++;}RE memory[n];}TE IN CO Mod& Mod::Factorial(CO uint& n)NE{ST Mod memory[COantsForMod::g_memory_LE] ={one(),one()};ST uint LE_curr = 2;WH(LE_curr <= n){memory[LE_curr] = MN::Factorial(LE_curr).Reduce();LE_curr++;}RE memory[n];}TE IN CO Mod& Mod::FactorialInverse(CO uint& n)NE{ST Mod memory[COantsForMod::g_memory_LE] ={one(),one()};ST uint LE_curr = 2;WH(LE_curr <= n){memory[LE_curr] = MN::FactorialInverse(LE_curr).Reduce();LE_curr++;}RE memory[n];}TE IN CE VO Mod::swap(Mod& n)NE{std::swap(m_n,n.m_n);}TE IN CE CO uint& Mod::RP()CO NE{RE m_n;}TE IN CE Mod Mod::DeRP(CO uint& n)NE{Mod n_copy{};n_copy.m_n = n;RE n_copy;}TE IN CE uint& Mod::Normalise(uint& n)NE{RE n < M?n:n -= M;}TE IN CO Mod& Mod::zero()NE{ST CE CO Mod z{};RE z;}TE IN CO Mod& Mod::one()NE{ST CE CO Mod o{DeRP(1)};RE o;}TE TE IN CE Mod& Mod::Ref(T&& n)NE{RE *TH;}TE IN CE Mod Twice(CO Mod& n)NE{RE MO(Mod(n).Double());}TE IN CE Mod Half(CO Mod& n)NE{RE MO(Mod(n).Halve());}TE IN Mod Inverse(CO Mod& n){RE MO(Mod(n).Invert());}TE IN CE Mod Inverse_COrexpr(CO uint& n)NE{RE MO(Mod::DeRP(RS(n)).NonNegativePW(M - 2));}TE IN CE Mod PW(CO Mod& n,CO T& EX){RE MO(Mod(n).PW(T(EX)));}TE IN CE VO swap(Mod& n0,Mod& n1)NE{n0.swap(n1);}TE IN string to_string(CO Mod& n)NE{RE to_string(n.RP())+ " + MZ";}TE IN basic_ostream& OP<<(basic_ostream& os,CO Mod& n){RE os << n.RP();} // using MOD = MP; using MOD = MNP; inline MOD Sum( const MOD& u0 , const MOD& u1 ) { return u0 + u1; }; inline const MOD& Zero() { return MOD::zero(); }; inline MOD Prod( const MOD& u0 , const MOD& u1 ) { return u0 * u1; }; int main() { UNTIE; // LIBRARY_SEARCH; // CEXPR( int , bound_T , 200000 ); // CIN_ASSERT( T , 1 , bound_T ); CEXPR( int , bound_N , 200000 ); // 0が5個 // CEXPR( ll , bound_N , 1000000000 ); // 0が9個 // CEXPR( ll , bound_N , 1000000000000000000 ); // 0が18個 // CEXPR( int , bound_M , 100000 ); // 0が5個 // CEXPR( ll , bound_M , 1000000000 ); // 0が9個 // CEXPR( ll , bound_M , 1000000000000000000 ); // 0が18個 CIN_ASSERT( N , 1 , bound_N ); EnumerationZetaTransform zt{ bound_Ai + 1 }; REPEAT( N ){ SET_ASSERT( Ai , 1 , bound_Ai ); zt.Add( Ai , 1 + zt.InitialSegmentSum( 1 ) - zt.InverseImageSum( 1 ) ); } RETURN( zt.InitialSegmentSum( 1 ) ); }