#ifdef INCLUDE_MAIN DEXPR( int , bound_N , 1000000 , 100 ); // 0が6個 DEXPR( int , sqrt_bound_N , 1000 , 10 ); // 0が3個 DECLARATION_OF_FUNCTIONS_FOR_DIVISOR_ZETA_TRANSFORM( sqrt_bound_N ); ll N; DEFINITION_OF_FUNCTIONS_FOR_DIVISOR_ZETA_TRANSFORM( N ); inline void Solve() { // // 数 // // DEXPR( ll , bound_N , 100000 , 100 ); // 0が5個 // // // DEXPR( ll , bound_N , 1000000000 , 100 ); // 0が9個 // // // DEXPR( ll , bound_N , 1000000000000000000 , 100 ); // 0が18個 // // CEXPR( TYPE_OF( bound_N ) , bound_M , bound_N ); // // // DEXPR( ll , bound_M , 100000 , 100 ); // 0が5個 // // // DEXPR( ll , bound_M , 1000000000 , 100 ); // 0が9個 // // // DEXPR( ll , bound_M , 1000000000000000000 , 100 ); // 0が18個 SET_LL( N ); ll a[bound_N+1]; a[0] = 0; FOREQ( i , 1 , N ){ a[i] = ( N / i ) * ( N / i - 1 ) / 2; } ZetaTransform zt{ N + 1 , a , true }; // // CIN( ll , M ); // // CIN( ll , N , M , K ); // // CIN_ASSERT( N , 1 , bound_N ); // ランダムテスト用。上限のデフォルト値は10^5。 // // CIN_ASSERT( M , 1 , bound_M ); // ランダムテスト用。上限のデフォルト値は10^5。 ll answer = N * ( N - 1 ) / 2 - zt.Get( 1 ); FOREQ( i , 2 , N ){ answer -= N / i - 1; } // // MP answer = 0; // // auto answer = Answer( N , M , K ); // // COUT( answer ); RETURN( answer ); // // 文字列 // CIN( string , S ); // // CIN( string , T ); // ll answer = 0; // // MP answer = 0; // // auto answer = Answer( N , M , K ); // // COUT( answer ); // RETURN( answer ); // // 配列 // CIN_A( ll , A , N ); // // CIN_A( ll , B , N ); // // vector A( N ); // // vector B( N ); // // ll A[bound_N]; // 関数(コンストラクタ)の引数に使う。長さのデフォルト値は10^5。 // // ll B[bound_N]; // 関数(コンストラクタ)の引数に使う。長さのデフォルト値は10^5。 // // FOR( i , 0 , N ){ // // cin >> A[i] >> B[i]; // // } // ll answer = 0; // // MP answer = 0; // // auto answer = Answer( N , M , K ); // // COUT( answer ); // // COUT_A( A , N ); // RETURN( answer ); // // 順列 // vector A( N ); // vector A_inv( N ); // FOR( i , 0 , N ){ // cin >> A[i]; // A_inv[--A[i]] = i; // } // ll answer = 0; // // MP answer = 0; // // auto answer = Answer( N , M , K ); // // COUT( answer ); // // COUT_A( A , N ); // RETURN( answer ); // // グラフ // e.resize( N ); // // e.resize( N ); // FOR( j , 0 , M ){ // CIN_ASSERT( uj , 1 , N ); // CIN_ASSERT( vj , 1 , N ); // uj--; // vj--; // e[uj].push_back( vj ); // e[vj].push_back( uj ); // // CIN( ll , wj ); // // e[uj].push_back( { vj , wj } ); // // e[vj].push_back( { uj , wj } ); // } // ll answer = 0; // // MP answer = 0; // // auto answer = Answer( N , M , K ); // // COUT( answer ); // // COUT_A( A , N ); // RETURN( answer ); // // 座標圧縮や単一クエリタイプなどのための入力格納 // vector > data( M ); // FOR( j , 0 , M ){ // CIN( ll , x , y , z ); // data[j] = { x , y , z }; // } // ll answer = 0; // // MP answer = 0; // // auto answer = Answer( N , M , K ); // // COUT( answer ); // // COUT_A( A , N ); // RETURN( answer ); // // 一般のクエリ // CIN( int , Q ); // // DEXPR( int , bound_Q , 100000 , 100 ); // 基本不要。 // // CIN_ASSERT( Q , 1 , bound_Q ); // 基本不要。 // // vector > query( Q ); // // vector > query( Q ); // FOR( q , 0 , Q ){ // CIN( int , type ); // if( type == 1 ){ // CIN( int , x , y ); // // query[q] = { type , x , y }; // } else if( type == 2 ){ // CIN( int , x , y ); // // query[q] = { type , x , y }; // } else { // CIN( int , x , y ); // // query[q] = { type , x , y }; // } // // CIN( int , x , y ); // // // query[q] = { x , y }; // } // // sort( query , query + Q ); // // FOR( q , 0 , Q ){ // // auto& [x,y] = query[q]; // // // auto& [type,x,y] = query[q]; // // } // auto answer = Answer( N , M , K ); // // COUT( answer ); // // COUT_A( A , N ); // RETURN( answer ); // // グリッド // // DEXPR( int , bound_H , 2000 , 30 ); // // // DEXPR( int , bound_H , 100000 , 10 ); // 0が5個 // // // CEXPR( int , bound_H , 1000000000 ); // 0が9個 // // // CEXPR( int , bound_W , bound_H ); // // static_assert( ll( bound_H ) * bound_W < ll( 1 ) << 31 ); // // CEXPR( int , bound_HW , bound_H * bound_W ); // // // CEXPR( int , bound_HW , 100000 ); // 0が5個 // // // CEXPR( int , bound_HW , 1000000 ); // 0が6個 // cin >> H >> W; // // SET_ASSERT( H , 1 , bound_H ); // ランダムテスト用。上限のデフォルト値は2*10^3。 // // SET_ASSERT( W , 1 , bound_W ); // ランダムテスト用。上限のデフォルト値は2*10^3。 // H_minus = H - 1; // W_minus = W - 1; // HW = H * W; // // assert( HW <= bound_HW ); // 基本不要。上限のデフォルト値は4*10^6。 // vector S( H ); // FOR( i , 0 , H ){ // cin >> S[i]; // // SetEdgeOnGrid( S[i] , i , e ); // // SetWallOnGrid( S[i] , i , non_wall ); // } // // {h,w}へデコード: EnumHW( v ) // // {h,w}をコード: EnumHW_inv( h , w ); // // (i,j)->(k,h)の方向番号を取得: DirectionNumberOnGrid( i , j , k , h ); // // v->wの方向番号を取得: DirectionNumberOnGrid( v , w ); // // 方向番号の反転U<->D、R<->L: ReverseDirectionNumberOnGrid( n ); // ll answer = 0; // // MP answer = 0; // // auto answer = Answer( N , M , K ); // // COUT( answer ); // // COUT_A( A , N ); // RETURN( answer ); } REPEAT_MAIN(1); #else // INCLUDE_MAIN #ifdef INCLUDE_SUB ll Naive( int N , int M , int K ) { ll answer = N + M + K; return answer; } ll Answer( ll N , ll M , ll K ) { // START_WATCH; ll answer = N + M + K; // // TLに準じる乱択や全探索。デフォルトの猶予は100.0[ms]。 // CEXPR( double , TL , 2000.0 ); // while( CHECK_WATCH( TL ) ){ // } return answer; } inline void Experiment() { // CEXPR( int , bound , 10 ); // FOREQ( N , 0 , bound ){ // FOREQ( M , 0 , bound ){ // FOREQ( K , 0 , bound ){ // COUT( N , M , K , ":" , Naive( N , M , K ) ); // } // } // // cout << Naive( N ) << ",\n"[N==bound]; // } } inline void SmallTest() { // CEXPR( int , bound , 10 ); // FOREQ( N , 0 , bound ){ // FOREQ( M , 0 , bound ){ // FOREQ( K , 0 , bound ){ // COMPARE( N , M , K ); // } // } // // COMPARE( N ); // } } #define INCLUDE_MAIN #include __FILE__ #else // INCLUDE_SUB #ifdef INCLUDE_LIBRARY /* C-x 3 C-x o C-x C-fによるファイル操作用 BFS: c:/Users/user/Documents/Programming/Mathematics/Geometry/Graph/BreadthFirstSearch/compress.txt CoordinateCompress: c:/Users/user/Documents/Programming/Mathematics/SetTheory/DirectProduct/CoordinateCompress/compress.txt DFSOnTree c:/Users/user/Documents/Programming/Mathematics/Geometry/Graph/DepthFirstSearch/Tree/a.hpp Divisor: c:/Users/user/Documents/Programming/Mathematics/Arithmetic/Prime/Divisor/compress.txt Polynomial c:/Users/user/Documents/Programming/Mathematics/Polynomial/compress.txt UnionFind c:/Users/user/Documents/Programming/Utility/VLTree/UnionFindForest/compress.txt */ // VVV 常設でないライブラリは以下に挿入する。 TE CL PrimeEnumeration{PU:bool m_is_composite[val_limit];INT m_val[LE_max];int m_LE;CE PrimeEnumeration();CE CO INT& OP[](CRI n) CO;CE CO INT& Get(CRI n) CO;CE CO bool& IsComposite(CRI i) CO;CE CRI LE() CO NE;}; TE CE PrimeEnumeration::PrimeEnumeration():m_is_composite(),m_val(),m_LE(0){for(INT i = 2;i < val_limit;i++){if(! m_is_composite[i]){INT j = i;WH((j += i)< val_limit){m_is_composite[j] = true;}m_val[m_LE++] = i;if(m_LE >= LE_max){break;}}}}TE CE CO INT& PrimeEnumeration::OP[](CRI n)CO{assert(n < m_LE);RE m_val[n];}TE CE CO INT& PrimeEnumeration::Get(CRI n)CO{RE OP[](n);}TE CE CO bool& PrimeEnumeration::IsComposite(CRI i)CO{assert(i < val_limit);RE m_is_composite[i];}TE CE CRI PrimeEnumeration::LE()CO NE{RE m_LE;} TE VO SetPrimeFactorisation(CO PrimeEnumeration& prime,CO INT1& n,VE& P,VE& EX){INT1 n_copy = n;int i = 0;WH(i < prime.m_LE){CO INT2& p = prime[i];if(p * p > n_copy){break;}if(n_copy % p == 0){P.push_back(p);EX.push_back(1);INT3& EX_back = EX.back();n_copy /= p;WH(n_copy % p == 0){EX_back++;n_copy /= p;}}i++;}if(n_copy != 1){P.push_back(n_copy);EX.push_back(1);}RE;} TE INT CountDivisor(CO PrimeEnumeration& prime,INT n) NE{VE P{};VE EX{};SetPrimeFactorisation(prime,n,P,EX);P.clear();CO int LE = EX.SZ();INT AN = 1;for(int i = 0;i < LE;i++){AN *= EX[i] + 1;}RE AN;} TE LI EnumerateDivisor(CO PrimeEnumeration& prime,INT n) NE{VE P{};VE EX{};SetPrimeFactorisation(prime,n,P,EX);CO int LE = P.SZ();LI divisor{};divisor.push_back(1);auto BE = divisor.BE(),EN = divisor.EN();for(int i = 0;i < LE;i++){CO INT& P_i = P[i];CRI EX_i = EX[i];LI temp{};INT PW = 1;for(int e = 1;e <= EX_i;e++){PW *= P_i;for(auto IT = BE;IT != EN;IT++){temp.push_back(*IT * PW);}}WH(! temp.empty()){divisor.push_back(temp.front());temp.pop_front();}}RE divisor;} TE VO MemoriseEnumerateDivisor(LI(&memory)[SZ_max])NE{for(int d = 1;d < SZ_max;d++){int n = 0;WH((n += d)< SZ_max){memory[n].push_back(d);}}RE;} TE int MoeviusFunction(CO PrimeEnumeration& prime,INT n)NE{int AN = 1;int i = 0;WH(i < prime.m_LE && n > 1){CRI p = prime[i++];if(n % p == 0){if((n /= p)% p == 0){RE 0;}AN *= -1;}}RE n == 1?AN:AN *= -1;} #define DEFINITION_OF_GET_FOR_ZETA_TRANSFORM( MU ) \ list sub = Sub( t ); \ sub.push_front( t ); \ U answer = Zero(); \ CERR( "Get for" , t ); \ \ while( ! sub.empty() ){ \ \ const T& t_sub = sub.front(); \ answer = Sum( answer , Prod( m_val[e_inv( t_sub )] , MU( t_sub , t ) ) ); \ CERR( "t_sub ==" , t_sub , ", answer +=" , Prod( m_val[e_inv( t_sub )] , MU( t_sub , t ) ) ); \ sub.pop_front(); \ \ } \ \ CERR( "answer ==" , answer ); \ 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; \ // Tの各要素tに対し、t以下の要素を渡る総和を管理。 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& InitialSegmentInverseImageSum( 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; }; // E_invはAdd(privateにはSup)にのみ使用。 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,>)のグラフである。 // を満たす場合にのみ以下をサポート。ただしE_invはAdd(privateにはSup)にのみ使用。 // 0による初期化O(size_max) // ゼータ変換前の配列による初期化O(size_max+始切片のサイズの総和) // ゼータ変換後の配列による初期化O(size_max) // 一点加算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] , const bool& transformed ); 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(size_max + digit 2^digit)(可換加法モノイド性を使う) // ゼータ変換後の配列による初期化O(size_max) // 一点加算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] , const bool& transformed ); 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)が単位的環である。 // を満たす場合にのみ以下をサポート。ただしE_invはAdd(privateにはSup)にのみ使用。 // 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が互いに逆写像である。 // を満たす場合にのみ以下をサポート。ただしE_invはAdd(privateにはSup)にのみ使用。 // 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] , const bool& transformed ) : ZetaTransformBody( size ) , PartiallyOrderedSetForZetaTransform() , EnumerationForZetaTransform() { using base = ZetaTransformBody; if( transformed ){ for( int i = 0 ; i < base::m_size ; i++ ){ base::m_val[i] = a[i]; } } else { for( int i = 0 ; i < base::m_size ; i++ ){ ll& m_val_i = base::m_val[i] = a[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() , m_digit( digit ) {} template inline FastZetaTransform::FastZetaTransform( const int& digit , const U ( &a )[size_max] , const bool& transformed ) : ZetaTransformBody( 1 << digit ) , SemiRingForZetaTransform( size ) { using base = ZetaTransformBody; for( int i = 0 ; i < base::m_size ; i++ ){ base::m_val[i] = a[i]; } if( ! transformed ){ 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::InitialSegmentInverseImageSum( 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() ){ cerr << "メビウス関数のデフォルトの再帰式が呼ばれています。" << endl; cerr << "このエラー出力回数が多い場合、再帰深度が深過ぎる可能性があります。" << endl; 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; } #define DECLARATION_OF_FUNCTIONS_FOR_DIVISOR_ZETA_TRANSFORM( CONSTEXPR_SQRT_BOUND_NUM ) \ constexpr PrimeEnumeration PRIME_FOR_DIVISOR_ZETA_FUNCTION{}; \ inline int DivisorMoeviusFunction( const int& t0 , const int& t1 ); \ inline int MultipleMoeviusFunction( const int& t0 , const int& t1 ); \ inline list DivisorEdge( const int& t ); \ inline list MultipleEdge( const int& t ); \ #define DEFINITION_OF_FUNCTIONS_FOR_DIVISOR_ZETA_TRANSFORM( NUM ) \ inline int DivisorMoeviusFunction( const int& t0 , const int& t1 ) { return TwoAryMoeviusFunction( PRIME_FOR_DIVISOR_ZETA_FUNCTION , t0 , t1 ); } \ inline int MultipleMoeviusFunction( const int& t0 , const int& t1 ) { return TwoAryMoeviusFunction( PRIME_FOR_DIVISOR_ZETA_FUNCTION , t1 , t0 ); } \ inline list DivisorEdge( const int& t ) { return DivisorEdgeForZetaTransform( PRIME_FOR_DIVISOR_ZETA_FUNCTION , NUM , t ); } \ inline list MultipleEdge( const int& t ) { return MultipleEdgeForZetaTransform( NUM , t ); } \ // 古典的な1変数メビウス関数(../../Arithmetic/Prime/Divisor/a_Body.hpp)の2変数化。 // 最初のみO(val_lim log val_lim)で後はO(1)。 // これを用いてmu(const int&,const int&)を // - t0 <= t1がt1 % t0 == 0で与えられる場合はこのままprimeを渡す。 // - t0 <= t1がt0 % t1 == 0で与えられる場合は変数を反対にしてprimeを渡す。 // と定義してZetaTansformBodyのGetに渡す。 template int TwoAryMoeviusFunction( const PrimeEnumeration& prime , const int& t0 , const int& t1 ); template inline list DivisorEdgeForZetaTransform( const PrimeEnumeration& prime , const int& t_max , const int& t ); inline list MultipleEdgeForZetaTransform( const int& t_max , const int& t ); template int TwoAryMoeviusFunction( const PrimeEnumeration& prime , const int& t0 , const int& t1 ) { constexpr int size_max = val_limit * val_limit; static int memory[size_max]; static bool uninitialised = true; if( uninitialised ){ vector rest( size_max ); for( int n = 0 ; n < size_max ; n++ ){ memory[n] = 1; rest[n] = n; } for( int i = 0 ; i < prime.LE() ; i++ ){ const INT& p_i = prime[i]; int n = 0; while( ( n += p_i ) < size_max ){ int& rest_n = rest[n] /= p_i; memory[n] *= ( rest_n % p_i == 0 ? 0 : -1 ); } } for( int n = val_limit ; n < size_max ; n++ ){ if( rest[n] != 1 ){ memory[n] *= -1; } } uninitialised = false; } assert( t1 % t0 == 0 ); return memory[ t1 / t0 ]; } template inline list DivisorEdgeForZetaTransform( const PrimeEnumeration& prime , const int& t_max , const int& t ) { assert( 0 <= t && t <= t_max ); list answer{}; if( t != 0 ){ answer = EnumerateDivisor( prime , t ); answer.pop_back(); } return answer; } inline list MultipleEdgeForZetaTransform( const int& t_max , const int& t ) { assert( 0 <= t && t <= t_max ); list answer{}; if( t != 0 ){ int temp = t; while( ( temp += t ) <= t_max ){ answer.push_back( temp ); } } return answer; } // AAA 常設でないライブラリは以上に挿入する。 #define INCLUDE_SUB #include __FILE__ #else // INCLUDE_LIBRARY // #define REACTIVE // #define USE_GETLINE #ifdef DEBUG #define _GLIBCXX_DEBUG #define REPEAT_MAIN( BOUND ) START_MAIN; signal( SIGABRT , &AlertAbort ); AutoCheck( exec_mode , use_getline ); if( exec_mode == sample_debug_mode || exec_mode == submission_debug_mode || exec_mode == library_search_mode ){ return 0; } else if( exec_mode == experiment_mode ){ Experiment(); return 0; } else if( exec_mode == small_test_mode ){ SmallTest(); return 0; }; DEXPR( int , bound_test_case_num , BOUND , min( BOUND , 100 ) ); int test_case_num = 1; if( exec_mode == solve_mode ){ if constexpr( bound_test_case_num > 1 ){ SET_ASSERT( test_case_num , 1 , bound_test_case_num ); } } else if( exec_mode == random_test_mode ){ CERR( "ランダムテストを行う回数を指定してください。" ); SET_LL( test_case_num ); } FINISH_MAIN #define DEXPR( LL , BOUND , VALUE , DEBUG_VALUE ) CEXPR( LL , BOUND , DEBUG_VALUE ) #define ASSERT( A , MIN , MAX ) CERR( "ASSERTチェック: " , ( MIN ) , ( ( MIN ) <= A ? "<=" : ">" ) , A , ( A <= ( MAX ) ? "<=" : ">" ) , ( MAX ) ); assert( ( MIN ) <= A && A <= ( MAX ) ) #define SET_ASSERT( A , MIN , MAX ) if( exec_mode == solve_mode ){ SET_LL( A ); ASSERT( A , MIN , MAX ); } else if( exec_mode == random_test_mode ){ CERR( #A , " = " , ( A = GetRand( MIN , MAX ) ) ); } else { assert( false ); } #define SOLVE_ONLY static_assert( __FUNCTION__[0] == 'S' ) #define CERR( ... ) VariadicCout( cerr , __VA_ARGS__ ) << endl #define COUT( ... ) VariadicCout( cout << "出力: " , __VA_ARGS__ ) << endl #define CERR_A( A , N ) OUTPUT_ARRAY( cerr , A , N ) << endl #define COUT_A( A , N ) cout << "出力: "; OUTPUT_ARRAY( cout , A , N ) << endl #define CERR_ITR( A ) OUTPUT_ITR( cerr , A ) << endl #define COUT_ITR( A ) cout << "出力: "; OUTPUT_ITR( cout , A ) << endl #else #pragma GCC optimize ( "O3" ) #pragma GCC optimize ( "unroll-loops" ) #pragma GCC target ( "sse4.2,fma,avx2,popcnt,lzcnt,bmi2" ) #define REPEAT_MAIN( BOUND ) START_MAIN; CEXPR( int , bound_test_case_num , BOUND ); int test_case_num = 1; if constexpr( bound_test_case_num > 1 ){ SET_ASSERT( test_case_num , 1 , bound_test_case_num ); } FINISH_MAIN #define DEXPR( LL , BOUND , VALUE , DEBUG_VALUE ) CEXPR( LL , BOUND , VALUE ) #define ASSERT( A , MIN , MAX ) assert( ( MIN ) <= A && A <= ( MAX ) ) #define SET_ASSERT( A , MIN , MAX ) SET_LL( A ); ASSERT( A , MIN , MAX ) #define SOLVE_ONLY #define CERR( ... ) #define COUT( ... ) VariadicCout( cout , __VA_ARGS__ ) << ENDL #define CERR_A( A , N ) #define COUT_A( A , N ) OUTPUT_ARRAY( cout , A , N ) << ENDL #define CERR_ITR( A ) #define COUT_ITR( A ) OUTPUT_ITR( cout , A ) << ENDL #endif #ifdef REACTIVE #define ENDL endl #else #define ENDL "\n" #endif #ifdef USE_GETLINE #define SET_LL( A ) { GETLINE( A ## _str ); A = stoll( A_str ); } #define GETLINE_SEPARATE( SEPARATOR , ... ) SOLVE_ONLY; string __VA_ARGS__; VariadicGetline( cin , SEPARATOR , __VA_ARGS__ ) #define GETLINE( ... ) SOLVE_ONLY; GETLINE_SEPARATE( '\n' , __VA_ARGS__ ) #else #define SET_LL( A ) cin >> A #define CIN( LL , ... ) SOLVE_ONLY; LL __VA_ARGS__; VariadicCin( cin , __VA_ARGS__ ) #define SET_A( A , N ) SOLVE_ONLY; FOR( VARIABLE_FOR_CIN_A , 0 , N ){ cin >> A[VARIABLE_FOR_CIN_A]; } #define CIN_A( LL , A , N ) vector A( N ); SET_A( A , N ); #endif #include using namespace std; using uint = unsigned int; using ll = long long; using ull = unsigned long long; using ld = long double; using lld = __float128; template using T2 = pair; template using T3 = tuple; template using T4 = tuple; using path = pair; #define ATT __attribute__( ( target( "sse4.2,fma,avx2,popcnt,lzcnt,bmi2" ) ) ) #define START_MAIN int main(){ ios_base::sync_with_stdio( false ); cin.tie( nullptr ) #define FINISH_MAIN REPEAT( test_case_num ){ if constexpr( bound_test_case_num > 1 ){ CERR( "testcase " , VARIABLE_FOR_REPEAT_test_case_num , ":" ); } Solve(); CERR( "" ); } } #define START_WATCH chrono::system_clock::time_point watch = chrono::system_clock::now() #define CURRENT_TIME static_cast( chrono::duration_cast( chrono::system_clock::now() - watch ).count() / 1000.0 ) #define CHECK_WATCH( TL_MS ) ( CURRENT_TIME < TL_MS - 100.0 ) #define TYPE_OF( VAR ) decay_t #define CEXPR( LL , BOUND , VALUE ) constexpr LL BOUND = VALUE #define CIN_ASSERT( A , MIN , MAX ) TYPE_OF( MAX ) A; SET_ASSERT( A , MIN , MAX ) #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 AUTO_ITR( ARRAY ) auto itr_ ## ARRAY = ARRAY .begin() , end_ ## ARRAY = ARRAY .end() #define FOR_ITR( ARRAY ) for( AUTO_ITR( ARRAY ) , itr = itr_ ## ARRAY ; itr_ ## ARRAY != end_ ## ARRAY ; itr_ ## ARRAY ++ , itr++ ) #define REPEAT( HOW_MANY_TIMES ) FOR( VARIABLE_FOR_REPEAT_ ## HOW_MANY_TIMES , 0 , HOW_MANY_TIMES ) #define SET_PRECISION( DECIMAL_DIGITS ) cout << fixed << setprecision( DECIMAL_DIGITS ) #define OUTPUT_ARRAY( OS , A , N ) FOR( VARIABLE_FOR_OUTPUT_ARRAY , 0 , N ){ OS << A[VARIABLE_FOR_OUTPUT_ARRAY] << (VARIABLE_FOR_OUTPUT_ARRAY==N-1?"":" "); } OS #define OUTPUT_ITR( OS , A ) { auto ITERATOR_FOR_OUTPUT_ITR = A.begin() , END_FOR_OUTPUT_ITR = A.end(); bool VARIABLE_FOR_OUTPUT_ITR = ITERATOR_FOR_COUT_ITR != END_FOR_COUT_ITR; while( VARIABLE_FOR_OUTPUT_ITR ){ OS << *ITERATOR_FOR_COUT_ITR; ( VARIABLE_FOR_OUTPUT_ITR = ++ITERATOR_FOR_COUT_ITR != END_FOR_COUT_ITR ) ? OS : OS << " "; } } OS #define RETURN( ... ) SOLVE_ONLY; COUT( __VA_ARGS__ ); return #define COMPARE( ... ) auto naive = Naive( __VA_ARGS__ ); auto answer = Answer( __VA_ARGS__ ); bool match = naive == answer; COUT( #__VA_ARGS__ , ":" , naive , match ? "==" : "!=" , answer ); if( !match ){ return; } // 入出力用 template inline basic_istream& VariadicCin( basic_istream& is ) { return is; } template inline basic_istream& VariadicCin( basic_istream& is , Arg& arg , ARGS&... args ) { return VariadicCin( is >> arg , args... ); } template inline basic_istream& VariadicGetline( basic_istream& is , const char& separator ) { return is; } template inline basic_istream& VariadicGetline( basic_istream& is , const char& separator , Arg& arg , ARGS&... args ) { return VariadicGetline( getline( is , arg , separator ) , separator , args... ); } template inline basic_ostream& VariadicCout( basic_ostream& os , const Arg& arg ) { return os << arg; } template inline basic_ostream& VariadicCout( basic_ostream& os , const Arg1& arg1 , const Arg2& arg2 , const ARGS&... args ) { return VariadicCout( os << arg1 << " " , arg2 , args... ); } // デバッグ用 #ifdef DEBUG inline void AlertAbort( int n ) { CERR( "abort関数が呼ばれました。assertマクロのメッセージが出力されていない場合はオーバーフローの有無を確認をしてください。" ); } void AutoCheck( int& exec_mode , const bool& use_getline ); inline void Solve(); inline void Experiment(); inline void SmallTest(); inline void RandomTest(); ll GetRand( const ll& Rand_min , const ll& Rand_max ); int exec_mode; CEXPR( int , solve_mode , 0 ); CEXPR( int , sample_debug_mode , 1 ); CEXPR( int , submission_debug_mode , 2 ); CEXPR( int , library_search_mode , 3 ); CEXPR( int , experiment_mode , 4 ); CEXPR( int , small_test_mode , 5 ); CEXPR( int , random_test_mode , 6 ); #ifdef USE_GETLINE CEXPR( bool , use_getline , true ); #else CEXPR( bool , use_getline , false ); #endif #endif // 圧縮用 #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& // VVV 常設ライブラリは以下に挿入する。 // AAA 常設ライブラリは以上に挿入する。 #define INCLUDE_LIBRARY #include __FILE__ #endif // INCLUDE_LIBRARY #endif // INCLUDE_SUB #endif // INCLUDE_MAIN