#ifdef DEBUG #define _GLIBCXX_DEBUG #define UNTIE ios_base::sync_with_stdio( false ); cin.tie( nullptr ); signal( SIGABRT , &AlertAbort ) #define DEXPR( LL , BOUND , VALUE , DEBUG_VALUE ) CEXPR( LL , BOUND , DEBUG_VALUE ) #define CERR( ANSWER ) cerr << ANSWER << endl; #define COUT( ANSWER ) cout << ANSWER << endl #define ASSERT( A , MIN , MAX ) CERR( "ASSERTチェック: " << ( MIN ) << ( ( MIN ) <= A ? "<=" : ">" ) << A << ( A <= ( MAX ) ? "<=" : ">" ) << ( MAX ) ); assert( ( MIN ) <= A && A <= ( MAX ) ) #define LIBRARY_SEARCH if( LibrarySearch() != 0 ){ QUIT; }; #define START_WATCH( PROCESS_NAME ) StartWatch( PROCESS_NAME ) #define STOP_WATCH( HOW_MANY_TIMES ) StopWatch( HOW_MANY_TIMES ) #else #pragma GCC optimize ( "O3" ) #pragma GCC optimize( "unroll-loops" ) #pragma GCC target ( "sse4.2,fma,avx2,popcnt,lzcnt,bmi2" ) #define UNTIE ios_base::sync_with_stdio( false ); cin.tie( nullptr ) #define DEXPR( LL , BOUND , VALUE , DEBUG_VALUE ) CEXPR( LL , BOUND , VALUE ) #define CERR( ANSWER ) #define COUT( ANSWER ) cout << ANSWER << "\n" #define ASSERT( A , MIN , MAX ) assert( ( MIN ) <= A && A <= ( MAX ) ) #define LIBRARY_SEARCH #define START_WATCH( PROCESS_NAME ) #define STOP_WATCH( HOW_MANY_TIMES ) #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 CEXPR( LL , BOUND , VALUE ) constexpr LL BOUND = VALUE #define CIN( LL , A ) LL A; cin >> A #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 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 QUIT return 0 #define SET_PRECISION( DECIMAL_DIGITS ) cout << fixed << setprecision( DECIMAL_DIGITS_ ) #define RETURN( ANSWER ) COUT( ( ANSWER ) ); QUIT #ifdef DEBUG inline void AlertAbort( int n ) { CERR( "abort関数が呼ばれました。assertマクロのメッセージが出力されていない場合はオーバーフローの有無を確認をしてください。" ); } void StartWatch( const string& process_name = "nothing" ); void StopWatch( const int& how_many_times = 1 ); #endif template inline T Absolute( const T& a ){ return a > 0 ? a : -a; } template inline T Residue( const T& a , const T& p ){ return a >= 0 ? a % p : ( a % p ) + p; } #define POWER( ANSWER , ARGUMENT , EXPONENT ) \ static_assert( ! is_same::value && ! is_same::value ); \ TYPE_OF( ARGUMENT ) ANSWER{ 1 }; \ { \ TYPE_OF( ARGUMENT ) ARGUMENT_FOR_SQUARE_FOR_POWER = ( ARGUMENT ); \ TYPE_OF( EXPONENT ) EXPONENT_FOR_SQUARE_FOR_POWER = ( EXPONENT ); \ while( EXPONENT_FOR_SQUARE_FOR_POWER != 0 ){ \ if( EXPONENT_FOR_SQUARE_FOR_POWER % 2 == 1 ){ \ ANSWER *= ARGUMENT_FOR_SQUARE_FOR_POWER; \ } \ ARGUMENT_FOR_SQUARE_FOR_POWER *= ARGUMENT_FOR_SQUARE_FOR_POWER; \ EXPONENT_FOR_SQUARE_FOR_POWER /= 2; \ } \ } \ #define POWER_MOD( ANSWER , ARGUMENT , EXPONENT , MODULO ) \ ll ANSWER{ 1 }; \ { \ ll ARGUMENT_FOR_SQUARE_FOR_POWER = ( MODULO + ( ( ARGUMENT ) % MODULO ) ) % MODULO; \ TYPE_OF( EXPONENT ) EXPONENT_FOR_SQUARE_FOR_POWER = ( EXPONENT ); \ while( EXPONENT_FOR_SQUARE_FOR_POWER != 0 ){ \ if( EXPONENT_FOR_SQUARE_FOR_POWER % 2 == 1 ){ \ ANSWER = ( ANSWER * ARGUMENT_FOR_SQUARE_FOR_POWER ) % MODULO; \ } \ ARGUMENT_FOR_SQUARE_FOR_POWER = ( ARGUMENT_FOR_SQUARE_FOR_POWER * ARGUMENT_FOR_SQUARE_FOR_POWER ) % MODULO; \ EXPONENT_FOR_SQUARE_FOR_POWER /= 2; \ } \ } \ #define FACTORIAL_MOD( ANSWER , ANSWER_INV , INVERSE , MAX_INDEX , CONSTEXPR_LENGTH , MODULO ) \ static ll ANSWER[CONSTEXPR_LENGTH]; \ static ll ANSWER_INV[CONSTEXPR_LENGTH]; \ static ll INVERSE[CONSTEXPR_LENGTH]; \ { \ ll VARIABLE_FOR_PRODUCT_FOR_FACTORIAL = 1; \ ANSWER[0] = VARIABLE_FOR_PRODUCT_FOR_FACTORIAL; \ FOREQ( i , 1 , MAX_INDEX ){ \ ANSWER[i] = ( VARIABLE_FOR_PRODUCT_FOR_FACTORIAL *= i ) %= MODULO; \ } \ ANSWER_INV[0] = ANSWER_INV[1] = INVERSE[1] = VARIABLE_FOR_PRODUCT_FOR_FACTORIAL = 1; \ FOREQ( i , 2 , MAX_INDEX ){ \ ANSWER_INV[i] = ( VARIABLE_FOR_PRODUCT_FOR_FACTORIAL *= INVERSE[i] = MODULO - ( ( ( MODULO / i ) * INVERSE[MODULO % i] ) % MODULO ) ) %= MODULO; \ } \ } \ // 二分探索テンプレート // EXPRESSIONがANSWERの広義単調関数の時、EXPRESSION >= TARGETの整数解を格納。 #define BS( ANSWER , MINIMUM , MAXIMUM , EXPRESSION , DESIRED_INEQUALITY , TARGET , INEQUALITY_FOR_CHECK , UPDATE_U , UPDATE_L , UPDATE_ANSWER ) \ static_assert( ! is_same::value && ! is_same::value ); \ ll ANSWER = MINIMUM; \ if( MINIMUM <= MAXIMUM ){ \ ll VARIABLE_FOR_BINARY_SEARCH_L = MINIMUM; \ ll VARIABLE_FOR_BINARY_SEARCH_U = MAXIMUM; \ ANSWER = ( VARIABLE_FOR_BINARY_SEARCH_L + VARIABLE_FOR_BINARY_SEARCH_U ) / 2; \ ll VARIABLE_FOR_DIFFERENCE_FOR_BINARY_SEARCH; \ while( VARIABLE_FOR_BINARY_SEARCH_L != VARIABLE_FOR_BINARY_SEARCH_U ){ \ VARIABLE_FOR_DIFFERENCE_FOR_BINARY_SEARCH = ( EXPRESSION ) - ( TARGET ); \ CERR( "二分探索中: " << VARIABLE_FOR_BINARY_SEARCH_L << "<=" << ANSWER << "<=" << VARIABLE_FOR_BINARY_SEARCH_U << ":" << EXPRESSION << "-" << TARGET << "=" << VARIABLE_FOR_DIFFERENCE_FOR_BINARY_SEARCH ); \ if( VARIABLE_FOR_DIFFERENCE_FOR_BINARY_SEARCH INEQUALITY_FOR_CHECK 0 ){ \ VARIABLE_FOR_BINARY_SEARCH_U = UPDATE_U; \ } else { \ VARIABLE_FOR_BINARY_SEARCH_L = UPDATE_L; \ } \ ANSWER = UPDATE_ANSWER; \ } \ CERR( "二分探索終了: " << VARIABLE_FOR_BINARY_SEARCH_L << "<=" << ANSWER << "<=" << VARIABLE_FOR_BINARY_SEARCH_U << ":" << EXPRESSION << ( EXPRESSION > TARGET ? ">" : EXPRESSION < TARGET ? "<" : "=" ) << TARGET ); \ CERR( ( EXPRESSION DESIRED_INEQUALITY TARGET ? "二分探索成功" : "二分探索失敗" ) ); \ assert( EXPRESSION DESIRED_INEQUALITY TARGET ); \ } else { \ CERR( "二分探索失敗: " << MINIMUM << ">" << MAXIMUM ); \ assert( MINIMUM <= MAXIMUM ); \ } \ // 単調増加の時にEXPRESSION >= TARGETの最小解を格納。 #define BS1( ANSWER , MINIMUM , MAXIMUM , EXPRESSION , TARGET ) \ BS( ANSWER , MINIMUM , MAXIMUM , EXPRESSION , >= , TARGET , >= , ANSWER , ANSWER + 1 , ( VARIABLE_FOR_BINARY_SEARCH_L + VARIABLE_FOR_BINARY_SEARCH_U ) / 2 ) \ // 単調増加の時にEXPRESSION <= TARGETの最大解を格納。 #define BS2( ANSWER , MINIMUM , MAXIMUM , EXPRESSION , TARGET ) \ BS( ANSWER , MINIMUM , MAXIMUM , EXPRESSION , <= , TARGET , > , ANSWER - 1 , ANSWER , ( VARIABLE_FOR_BINARY_SEARCH_L + 1 + VARIABLE_FOR_BINARY_SEARCH_U ) / 2 ) \ // 単調減少の時にEXPRESSION >= TARGETの最大解を格納。 #define BS3( ANSWER , MINIMUM , MAXIMUM , EXPRESSION , TARGET ) \ BS( ANSWER , MINIMUM , MAXIMUM , EXPRESSION , >= , TARGET , < , ANSWER - 1 , ANSWER , ( VARIABLE_FOR_BINARY_SEARCH_L + 1 + VARIABLE_FOR_BINARY_SEARCH_U ) / 2 ) \ // 単調減少の時にEXPRESSION <= TARGETの最小解を格納。 #define BS4( ANSWER , MINIMUM , MAXIMUM , EXPRESSION , TARGET ) \ BS( ANSWER , MINIMUM , MAXIMUM , EXPRESSION , <= , TARGET , <= , ANSWER , ANSWER + 1 , ( VARIABLE_FOR_BINARY_SEARCH_L + VARIABLE_FOR_BINARY_SEARCH_U ) / 2 ) \ // 圧縮用 #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 ASK_DETAILS( ... ) \ CERR( "問題の区分は以下の中で何番に該当しますか?" ); \ problems = { __VA_ARGS__ }; \ problems_size = problems.size(); \ FOR( i , 0 , problems_size ){ \ CERR( i << ": " << problems[i] ); \ } \ cin >> num; \ CERR( "" ); \ num_temp = 0; \ if( num < 0 || num >= problems_size ){ \ CERR( "返答は" << problems_size - 1 << "以下の非負整数にしてください。" ); \ CERR( "終了します。" ); \ CERR( "" ); \ return -1; \ } \ int LibrarySearch( int num = -1 ) { vector problems{}; int problems_size = 13; int num_temp = 0; string reply{}; if( num == -1 ){ CERR( "ライブラリーを探索しますか?[y/n]" ); cin >> reply; if( reply == "n" ){ CERR( "ライブラリーを探索せずに続行します。" ); CERR( "" ); return 0; } else if( reply != "y" ){ CERR( "y/nのいずれかで答えてください。" ); CERR( "終了します。" ); CERR( "" ); return -1; } CERR( "" ); CERR( "ライブラリーを探索します。" ); ASK_DETAILS( "数や関数に関する問題。" , "配列に関する問題。" , "文字列に関する問題。" , "順列に関する問題。" , "矩形領域に関する問題。" , "グラフに関する問題。" , "部分和問題。" , "確率/期待値に関する問題。" , "ゲームに関する問題。" , "論理に関する問題。" , "半順序集合に関する問題。" , "関数適用に関する問題。" , "構築問題。" ); } else { CERR( "" ); } CEXPR( int , num_graph , 5 ); CEXPR( int , num_subsequence_sum , 6 ); CEXPR( int , num_game , 8 ); if( num == num_temp++ ){ CERR( "入力は1つの数か、1つの数と法を表す数ですか?[y/n]" ); cin >> reply; CERR( "" ); if( reply == "y" ){ CERR( "まずは小さい入力の場合を愚直に計算し、OEISで検索しましょう。" ); CERR( "https://oeis.org/?language=japanese" ); CERR( "" ); CERR( "次に出力の定義と等価な式を考察しましょう。" ); CERR( "- 単調ならば、冪乗や階乗" ); CERR( "- 定義にp進法が使われていれば、各種探索アルゴリズム" ); CERR( "- 入力が素数に近い場合に規則性があれば、p進付値、p進法、" ); CERR( " オイラー関数、約数の個数など" ); CERR( "を検討しましょう。" ); } else if( reply == "n" ){ ASK_DETAILS( "求解問題。" , "最大化問題。" , "数え上げ問題。" , "その他。" ); if( num == num_temp++ ){ CERR( "まず" ); CERR( "- 単調関数は二分探索" ); CERR( "- 可微分関数はニュートン法" ); CERR( "- 一次関数は掃き出し法" ); CERR( "を検討しましょう。" ); CERR( "" ); CERR( "それらが間に合わない時は" ); CERR( "- 探索範囲を絞れそうな場合は全探策" ); CERR( "- 多変数の合成関数で表せる場合は半分全列挙" ); CERR( "を検討しましょう。" ); } else if( num == num_temp++ ){ CERR( "- 凸関数は三分探索" ); CERR( "- 可微分関数は微分のニュートン法" ); CERR( "を検討しましょう。" ); } else if( num == num_temp++ ){ CERR( "- 変数の対称性があれば大小関係を制限した全探策" ); CERR( "- 何らかの約数となるなど動く範囲が狭い変数があればそれらを決め打った全探策" ); CERR( "- 多変数の合成関数で表せる場合は半分全列挙" ); CERR( "を検討しましょう。" ); } else if( num == num_temp++ ){ CERR( "このケースのライブラリー探索は不完全です。" ); } } else { CERR( "y/nのいずれかで答えてください。" ); CERR( "終了します。" ); CERR( "" ); return -1; } CERR( "" ); CERR( "マルチテストケースの場合は以下の前計算を検討しましょう;" ); CERR( "素数列挙、約数列挙、オイラー関数の値列挙、サブゴールとなる関係式を満たす解列挙。" ); } else if( num == num_temp++ ){ ASK_DETAILS( "区間処理問題。" , "最大化問題。" , "最長部分列問題。" , "数え上げ問題。" , "部分和問題。" ); if( num == num_temp++ ){ ASK_DETAILS( "可換群構造+を使う問題。" , "可換羃等モノイド構造∨を使う問題。" , "モノイド構造*を使う問題。" , "非結合的マグマ構造*を使う問題。" , "集合へのマグマ作用(*,\\cdot)を使う問題。" , "モノイドへのマグマ作用(+,\\cdot)を使う問題。" , "定数とのmaxを取った値の区間和取得を使う問題。" ); if( num == num_temp++ ){ CERR( "- 区間加算/区間取得が必要ならば可換群BIT" ); CERR( " \\Mathematics\\SetTheory\\DirectProduct\\AffineSpace\\BIT\\Template" ); CERR( "- 一点代入/一点加算/区間取得が必要ならば可換群平方分割" ); CERR( " \\Mathematics\\SetTheory\\DirectProduct\\AffineSpace\\SqrtDecomposition\\Template" ); CERR( "- 区間以外の領域で加算/全更新後の一点取得が必要ならば階差数列" ); CERR( " \\Mathematics\\SetTheory\\DirectProduct\\Tree\\DifferenceSeqeuence" ); CERR( "を検討しましょう。" ); } else if( num == num_temp++ ){ CERR( "- 一点代入/区間加算/一点取得/区間取得が必要ならば可換羃等モノイドBIT" ); CERR( " \\Mathematics\\SetTheory\\DirectProduct\\AffineSpace\\BIT\\IntervalMax\\Template" ); CERR( "を検討しましょう。" ); } else if( num == num_temp++ ){ CERR( "- 一点代入/区間取得が必要ならばモノイドBIT" ); CERR( " \\Mathematics\\SetTheory\\DirectProduct\\AffineSpace\\BIT\\Template\\Monoid" ); CERR( "- 一点加算/区間加算/一点取得/区間取得が必要ならばモノイド平方分割" ); CERR( " \\Mathematics\\SetTheory\\DirectProduct\\AffineSpace\\SqrtDecomposition\\Template\\Monoid" ); CERR( "- 一点代入/一点取得/区間取得が必要ならばモノイドセグメント木" ); CERR( " \\Mathematics\\SetTheory\\DirectProduct\\AffineSpace\\SegmentTree" ); CERR( "を検討しましょう。" ); } else if( num == num_temp++ ){ CERR( "- 写像のコード化" ); CERR( " \\Mathematics\\Function\\Encoder" ); CERR( "によりモノイドに帰着させることを検討しましょう。" ); } else if( num == num_temp++ ){ CERR( "- 一点作用/区間作用/一点取得が必要ならば双対平方分割" ); CERR( " \\Mathematics\\SetTheory\\DirectProduct\\AffineSpace\\SqrtDecomposition\\Template\\Dual" ); CERR( "を検討しましょう。" ); } else if( num == num_temp++ ){ CERR( "- 区間代入/区間作用/区間加算/一点取得/区間取得が必要な場合は遅延評価平方分割" ); CERR( " \\Mathematics\\SetTheory\\DirectProduct\\AffineSpace\\SqrtDecomposition\\Template\\LazyEvaluation" ); CERR( "を検討しましょう。" ); } else if( num == num_temp++ ){ CERR( "maxで全体更新でなく区間更新をする場合の汎用的な解法は分かりません。" ); CERR( "例えば区間が包含について単調でmaxを取る値も単調であれば全体更新と" ); CERR( "同様の処理ができます。状況に応じた考察をしましょう。" ); CERR( "" ); CERR( "maxで全体更新をする場合、maxを取る値は単調である場合に帰着できます。" ); CERR( "maxで全体更新をしない場合、つまりただmaxの区間和を処理するだけの場合、" ); CERR( "クエリの順番を入れ替えることができるので、単調な全体更新に帰着できます。" ); CERR( "従って以下では単調な全体更新の問題を考えます。" ); CERR( "" ); CERR( "maxを取る定数を変数化し、元の値との大小を表す{0,1}値の係数を考えます。" ); CERR( "すると区間作用前後の値は統一的にその係数と変数を使って表せます。" ); CERR( "配列の各成分の係数の値が変化するイベントとクエリをソートして管理し、" ); CERR( "クエリがイベントを跨ぐたびに係数を更新することを検討しましょう。" ); CERR( "" ); CERR( "例えばクエリB_qに対するmax(A_i,B_q)の区間和は、" ); CERR( "- 優先度つきキューA'={(A_i,i)|i}(構築O(N log N))" ); CERR( "- (B_q,q)_qをソートしたB'(構築O(Q log Q))" ); CERR( "- 長さNの数列C=(0,...,0)(構築O(N))" ); CERR( "を用意し、B'を前から探索して順に各クエリ(B_q,q)を処理します。" ); CERR( "具体的にはA'を前から探索して順にA_iCL PrimeEnumeration{PU:INT m_val[LE_max];int m_LE;CE PrimeEnumeration();};TE CE PrimeEnumeration::PrimeEnumeration():m_val(),m_LE(0){bool is_comp[val_limit] ={};for(INT i = 2;i < val_limit;i++){if(is_comp[i] == false){INT j = i;WH((j += i) < val_limit){is_comp[j] = true;}m_val[m_LE++] = i;if(m_LE >= LE_max){break;}}}}TE VO SetPrimeFactorisation(CO PrimeEnumeration& prime,CO INT& n,VE& P,VE& EX){INT n_copy = n;int i = 0;WH(i < prime.m_LE){CO INT& p = prime.m_val[i];if(p * p > n_copy){break;}if(n_copy % p == 0){P.push_back(p);EX.push_back(1);INT& 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;} #define EULER_FUNCTION EX.clear();INT AN = n;CO INT SZ = P.SZ();for(INT i = 0;i < SZ;i++){AN -= AN / P[i];}RE AN; TE INT EulerFunction(CO PrimeEnumeration& prime,CO INT& n,VE& P,VE& EX){SetPrimeFactorisation(prime,n,P,EX);EULER_FUNCTION;}TE IN INT EulerFunction(CO PrimeEnumeration& prime,CO INT& n){VE P{};VE EX{};RE EulerFunction(prime,n,P,EX);}TE VO MemoriseEulerFunction(CO PrimeEnumeration& prime,INT2(&memory)[SZ]){ST INT quotient[SZ];for(int n = 1;n < SZ;n++){memory[n] = quotient[n] = n;}for(int i = 0;i < prime.m_LE;i++){CO INT& p_i = prime.m_val[i];int n = 0;WH((n += p_i)< SZ){INT2& memory_n = memory[n];INT& quotient_n = quotient[n];memory_n -= memory_n / p_i;WH((quotient_n /= p_i)% p_i == 0){}}}for(int n = val_limit;n < SZ;n++){CO INT& quotient_n = quotient[n];if(quotient_n != 1){INT2& memory_n = memory[n];memory_n -= memory_n / quotient_n;}}RE;} 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.m_val[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(); \ \ 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]; vector m_val; 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; }; 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] ); inline ZetaTransform( const int& size , const vector&a ); 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 ZetaTransformBody::ZetaTransformBody( const int& size ) : m_length() , m_memory() , m_memory_inv() , m_size( size ) , m_val(size_max) {} 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] ) : template E(const int&) , list E_inv(const int&) , int size_max> inline ZetaTransform::ZetaTransform( const int& size , const vector& a ) : ZetaTransformBody( size ) , PartiallyOrderedSetForZetaTransform() , EnumerationForZetaTransform() { using base = ZetaTransformBody; 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] ) : ZetaTransformBody( 1 << digit ) , SemiRingForZetaTransform( size ) { 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::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; } template int TwoAryMoeviusFunction( const PrimeEnumeration& prime , const int& t0 , const int& t1 ) { // constexpr int size_max = val_limit * val_limit; constexpr int size_max = 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.m_LE ; i++ ){ const INT& p_i = prime.m_val[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 ]; } inline PrimeEnumeration prime{}; inline CEXPR( int , bound_N , 300000 ); // 0が5個 list divisor[bound_N+1] = {}; inline list E(const int& i ) { return divisor[i]; } inline list E_inv(const int&) { return list(); } // 使わない inline int mu( const int& n , const int& m ) { return TwoAryMoeviusFunction( prime , n , m ); } ll f( const int& N , const ll& N2 , const ll& m ) { // ll A[bound_N + 1]; // vector A( bound_N + 1 ); // FOREQ( i , 1 , N ){ // A[i] = m * i / N2; // } // ZetaTransform zt{ N + 1 , A }; START_WATCH( "construct zt" ); vector zt( bound_N + 1 ); FOREQ( i , 1 , N ){ zt[i] = m * i / N2; } FOR( i , 0 , prime.m_LE ){ const int& p_i = prime.m_val[i]; FOREQINV( n , N / p_i , 0 ){ zt[n * p_i] -= zt[n]; } } STOP_WATCH(); START_WATCH( "compute answer" ); ll answer = 0; FOREQ( i , 1 , N ){ // answer += f.Get( i ); answer += zt[i]; } STOP_WATCH(); return answer; } int main() { UNTIE; LIBRARY_SEARCH; CEXPR( int , bound_T , 15000 ); CIN_ASSERT( T , 1 , bound_T ); static ll euler_sum[bound_N + 1]; MemoriseEulerFunction( prime , euler_sum ); FOREQ( i , 2 , bound_N ){ euler_sum[i] += euler_sum[i-1]; } MemoriseEnumerateDivisor( divisor ); FOREQ( i , 1 , bound_N ){ divisor[i].pop_back(); } REPEAT( T ){ CIN_ASSERT( N , 1 , bound_N ); ll N2 = ll( N ) * N; CIN_ASSERT( K , 1 , N2 ); ll& M = euler_sum[N]; if( K == M ){ COUT( "1/1" ); continue; } if( K >= M * 2 ){ COUT( -1 ); continue; } bool over = K > M; if( over ){ K = M * 2 - K; } BS1( m , 0 , N2 - 1 , f( N , N2 , m ) , K ); FOREQ( q , 1 , N ){ int p = q * m / N2; if( p > q * ( m - 1 ) / N2 ){ if( over ){ swap( p , q ); } COUT( p << "/" << q ); break; } } } // CEXPR( int , bound_N , 1000000000 ); // 0が9個 // CEXPR( ll , bound_N , 1000000000000000000 ); // 0が18個 // CEXPR( int , bound_M , 100000 ); // 0が5個 // // CEXPR( int , bound_M , 1000000000 ); // 0が9個 // // CEXPR( ll , bound_M , 1000000000000000000 ); // 0が18個 // CIN_ASSERT( M , 0 , bound_M ); // CEXPR( ll , P , 998244353 ); // CEXPR( ll , P , 1000000007 ); // CEXPR( int , bound_Q , 100000 ); // CIN_ASSERT( Q , 1 , bound_Q ); // REPEAT( Q ){ // COUT( N ); // } // RETURN( N ); QUIT; }