No.1035 Color Box

yukicoder.me

$ M $ 種類のペンキをすべて使わなくても良い場合は、 $ N $ 個の箱がそれぞれが $ M $ 個の選択肢を持っているため $M ^{N}$ 通りとなります。

ここで問題となってしまうのが、$ M $ 種類のペンキのうち1回も使われないものが存在することです。

そこで、

$$ (0種類以上のペンキを使わない) - (1種類以上のペンキを使わない) + (2種類以上のペンキを使わない) $$

といった具合で包除原理を適用します。

$ i $ 種類以上のペンキを使わない場合の数は、

  • $ M $ 種類のペンキのうち使用する $ M - i $ 種類を選ぶ $ \,_{M} C_{M - i} $
  • $ N $ 個の箱の色を $ M - i $ 種類のペンキのうちから選ぶ $ (M - i)^{N} $

の積をとることで、

$$ _{M}C_{M - i} \times (M - i) ^ {N} $$

となります。したがって、全体の場合の数は

$$ \sum _{i = 0}^{M - 1} (-1)^{i} \times _{M}C_{M - i} \times (M - i) ^ {N} $$

と求めることが出来ます。 $ M - i $種類以上のペンキを使用しないとすると、

$$ \sum _{i = 1}^{M} (-1)^{M - i} \times _{M}C_{i} \times i ^ {N} $$

とすることもできます。

#include <atcoder/all>
using namespace std;
using mint = atcoder::modint1000000007;

int main(){
    int n, m;
    cin >> n >> m;
    
    vector<mint> fact(m + 1), inv(m + 1);
    fact[0] = 1;
    for(int i = 1; i <= m; i++) fact[i] = i * fact[i - 1];
    inv[m] = 1 / fact[m];
    for(int i = m; i >= 1; i--) inv[i - 1] = i * inv[i];
    
    mint ans;
    for(int i = 1; i <= m; i++){
        ans += (m - i & 1 ? -1 : 1) * inv[m - i] * inv[i] * mint(i).pow(n);
    }
    
    cout << (ans * fact[m]).val() << '\n';
}

yukicoder No.2164 Equal Balls

yukicoder.me

解説というよりかは、ACするまでに考えたこととそれからの高速化について考えたことのメモです。
具体的には、私の初AC時は${\rm O}(NM + M ^ {3} \log M)$ だったのですが、そこから ${\rm O} ( NM + M ^ {2} \log^{2} M ) $への高速化を行いました。
また、計算量を考えるうえで $ M = {\rm max} A = {\rm max} B $ を暗黙のうちに用いています。

ACするまでに考えたこと

まず、Aの和、Bの和をそれぞれ持ったDPを考えたくなりますが、$N = 10 ^ 5 $ のため計算量的に無理そうだと悟ります。そこで、(Aの和) $-$ (Bの和)といった変換を行うことで、和について持つkeyは1つで良くなりそうだと気づきますが、これでもまだ厳しそうです(実際は直前 $ M $ 個の情報全部持たなければいけないためおそらくkeyは1つでは無理)。以降、 $s = $ (Aの和) $-$ (Bの和)とします。
DPの遷移について考えてみると、$M - 1$ 項目までは$s$ はなんでもよいですが、 $ M $ 項目以降は直前の $ M $ 個の和が完全に一致している($ s = 0 $ である)必要があることに気づきます。そして、1項目が $ s = 1 $ だった場合、$ M + 1 $に遷移する際は直前 $ M $ 個の $ s $ が1減ることになります。すなわち、 $ M $ 個の和について常に $ s = 0 $ が成り立っていなければならないので、この均衡を保つためには、 $ M + 1 $ でも $ s = 1 $ とならなければいけないことが分かります。

したがって、要素番号 $ i $ の $ i \bmod M $ が等しいものについてはすべてにおいて $ s $ が等しいものを選択し、すべての余剰 $ 0 \leq x \leq M - 1 $についての $ s $ の和が0となっているものを数え上げればいいことになります。これならば、マイナスを考慮しても状態が$ 300 \times 600 $程度に抑えられそうなため、遷移に $ {\rm O} (M) $ 程度かけても間に合いそうだと考えました。

方針がたったため、各要素 $ i $ について$\bmod M $ でまとめます。すなわち、余剰が同じものについて $ s $ が同じものについての積をとります。しかし、ここで少し困りました。 $ A[i] $の中から $ a $ 個選ぶ、 $ B[i] $ の中から $ b $ 個選ぶ場合、$ s = a - b $ となり、それぞれ $ A[i] $ と $ B[i] $ についての二項係数の積で求めることができますが、同じ$s $ の中でも($s = 1$を例に挙げると、$(a, b) = (1, 0), (2, 1), (3, 2)$などのように) 複数のものが同じ$s$に集められる場合があることに気がつきます。愚直にすべての積をとってしまうと、$ {\rm O} (M ^ 2)$かかってしまい、状態数を考えると$ {\rm O} (NM ^ 2) $まで膨れ上がるため、高速化を考える必要がありました。まず、思いついたものとしてそれぞれを$ {\rm O} (M \log M) $で畳み込むことですが、それでも $ M = 600 $程度の $ {\rm O} (N M \log M) $ が間に合うとは思えません。そこで、DEGwerさんの数え上げpdfに二項係数の畳み込みみたいなの載ってないかなと思い、見に行くと30ページの公式集に見つかります。

drive.google.com

$$ \sum_{i = 0}^{k} \binom{n}{i} \binom{m}{k -i } = \binom{n + m}{k} $$

らしいです。$ (1 + x) ^ {N}$の係数が二項係数になるので、

$$ (1 + x) ^ {A[i]} (1 + x) ^ {B[i]} = (1 + x) ^ {A[i] + B[i]} $$

から自明だった...。ただ、$a + b$ではなく,$a - b$で使用しているためこの公式ができるか確信が持てなかったのですが,パスカルの三角形は左右対称的で$B[i]$から$b $ 個選ぶ場合の数は、$B[i]$から $b$ 個選ばないものを選ぶ場合の数に等しいことからこの公式が使えそうだと思いいたります。 このようにして、各余剰について$s$に対する二項係数の積を求めるパートは$ {\rm O} (NM )$でできました。

続いて、最後の和の数え上げパートですが、ここではMの和を情報に持つ必要があるため、状態数が $ {\rm O} (M^{2}) $ 、遷移が $ {\rm O} (M)$ 、遷移回数が $ {\rm O} ( M ) $ 、全体が$ {\rm O} (M ^ {4})$の計算量になり、ここでも高速化が必要になりました。ここの高速化は畳み込みしか思いつかなかったのもありますが、今回は長さ $ M ^ 2 $ 程度のものを $ M $ 回畳み込めばよいため、間に合いそうだと思いました。よって、計算量 $ {\rm O} ( M^{3} \log M ) $ 程度で間に合いました。

高速化について

無事ACできたのですが、初提出の実行時間は4041 msで本当に想定解なのか疑わしいです...。解説を見ようにもこの日はACとっても見れない仕様になっているようでした。
最終的には、以下の点を改良することで1000 msを切る程度まで実行時間を改良することができました。

  • 二項係数の計算を階乗前計算から和の漸化式による計算を切り替える(パスカルの三角形からの計算)
  • $ \bmod M $ が等しいものをまとめるパートで毎回余剰を求めないようにする
  • $ \bmod M $ が異なるものについて全体の畳み込みを行う際にサイズが近いもので畳み込みされるようにする

前2つはなるべく $ \bmod $ の計算回数を減らそうというものなのでいいとして、最後のものについて軽く触れます。

これは、端的に言えば、先頭から畳み込み演算を行っていたものを、並列的な畳み込み演算に切り替えたことになります。配列のサイズを $ N $ とすると、どちらも $N - 1 $ 回の計算であることには変わりないですが、畳み込みのサイズが後者の方が抑えられそうな気がしました。

この発想に至った経緯について触れます。畳み込みの計算量は配列 $ A $ のサイズを $ | A | $ 、配列 $ B $ のサイズを $ | B | $ とすると、 $ {\rm O} ( ( |A| + |B| ) \log ( |A| + |B| ) ) $になります。$ \log $の計算量評価は難しいため、ひとまず $\log$ なしで考えてみると畳み込みの累積演算は以下の問題に帰着されるかと思います。

$ N $ 個の要素からなる配列 $ A_i $ が与えられます。要素$i$、$j$について合体させるには、コストが $ A_i + A_j $ かかり、合体後のサイズは$ A_i + A_j $ となります。最終的に配列の要素が1つになるまで合体させるのにかかる最小コストはいくつですか?

これは、AtCoderのEducational DP Contest / DP まとめコンテスト 「N - Slimes 」の隣接選択の条件が外れただけのものになります。

atcoder.jp

隣接選択の条件が外れたより近い問題としてはAtCoder Beginner Contest 252の「F - Bread」やプログラミングコンテストチャレンジブック第2版49ページの問題などがります。いずれもなるべく小さいサイズのもの2つ選んでマージするという貪欲法によって最適解が得られました。以上の経験から、先頭から累積的に計算するよりも並列的に計算した方が計算量が良くなるのではないかと予想しました。

atcoder.jp

実際に評価してみると、8個の要素からなる配列を先頭からマージしていく場合、

  • $ a_1 a_2 a_3 a_4 a_5 a_6 a_7 a_8 $  
    合計コスト $0$
  • $ ( a_1 + a_2 ) a_3 a_4 a_5 a_6 a_7 a_8 $  
    合計コスト $a_1+a_2$
  • $ ( ( a_1+a_2)+a_3) a_4 a_5 a_6 a_7 a_8 $  
    合計コスト $2(a_1+a_2) + a_3$
  • $ ( ( ( a_1+a_2)+a_3)+a_4) a_5 a_6 a_7 a_8 $  
    合計コスト $3(a_1+a_2) + 2a_3 + a_4$
  • $ ( ( ( ( a_1+a_2)+a_3)+a_4)+a_5) a_6 a_7 a_8 $  
    合計コスト $4(a_1+a_2) + 3a_3 + 2a_4 + a_5$
  • $ ( ( ( ( ( a_1+a_2)+a_3)+a_4)+a_5)+a_6) a_7 a_8 $  
    合計コスト $5(a_1+a_2) + 4a_3 + 3a_4 + 2a_5 + a_6$
  • $ ( ( ( ( ( ( a_1+a_2)+a_3)+a_4)+a_5)+a_6)+a_7) a_8 $  
    合計コスト $6(a_1+a_2) + 5a_3 + 4a_4 + 3a_5 + 2a_6 + a_7$
  • $ ( ( ( ( ( ( ( a_1+a_2)+a_3)+a_4)+a_5)+a_6)+a_7)+a_8) $  
    合計コスト $7(a_1+a_2) + 6a_3 + 5a_4 + 4a_5 + 3a_6 + 2 a_7 + a_8$

のようになります。配列 $ A $ のサイズを $ M $ 、 $ A $ の要素をすべて $ 1 $ とすると、合計コストは $ M ^ 2 ( M + 1 ) / 2 $ 程度になります。($ A $ の要素が$ 1 $の場合は、$ M (M + 1) / 2$ 程度)

8個の要素からなる配列を並列的にマージしていく場合、

  • $ a_1 a_2 a_3 a_4 a_5 a_6 a_7 a_8 $  
    合計コスト $0$
  • $ (a_1+a_2) (a_3+a_4) (a_5+a_6) (a_7+a_8) $  
    合計コスト $a_1 + a_2 +a_3 + a_4 + a_5 + a_6 + a_7 + a_8$
  • $ ( ( a_1+a_2)+(a_3+a_4)) ( ( a_5+a_6)+(a_7+a_8)) $  
    合計コスト $2(a_1 + a_2 +a_3 + a_4 + a_5 + a_6 + a_7 + a_8)$
  • $ ( ( ( a_1+a_2)+(a_3+a_4))+ ( ( a_5+a_6)+(a_7+a_8))) $  
    合計コスト $3(a_1 + a_2 +a_3 + a_4 + a_5 + a_6 + a_7 + a_8)$

となります。配列 $ A $ のサイズを $ M $ 、 $ A $ の要素をすべて $ M $ とすると、合計コストは $ M^ 2 \log M $ 程度になります。($ A $ の要素が$ 1 $の場合は、$ M \log M $ 程度)

若干疑わしかったため、念のためEDPC - TのACコードに

8
1 1 1 1 1 1 1 1

を投げてみます。 すると無事 $ 8 \times \log_2 8 = 8 \times 3 = 24 $ が返ってきました。

畳み込みの計算量を $ {\rm O} (M_1 \log M_2 ) $ とすると、ここでは $ M_1 $ の部分を考えたことになります。 $ M_2 $ の部分は $ M ^ {2} $ になっても定数倍を考慮すると、$ \log M $ に落ちます。したがって、先頭から累積的にマージさせると$ {\rm O} (NM + M ^ {3} \log M) $ であったのに対し、並列的にマージすることによって計算量を $ {\rm O}( NM + M^{2} \log ^{2} M ) $ まで落とすことができました。 $ NM $ は二項係数の積を求めるパートに依存する計算量です。
C++の標準ライブラリでは、std::accumulateが先頭から累積的にマージするもので、C++17以降で使えるstd::reduceが並列的に計算を行うものになります。C++に馴染みのない方へ簡単に説明すると、accumulate的なマージがデータ構造でいうと累積和にあたり、reduce的なマージがデータ構造でいうとセグメント木にあたるかと思います。

実際に、accumulateとreduceを使って実装してみるとaccumulateからreduceに切り替えることによって、2948 msから996 msへとなり2秒近く早くなりました。こんなことによってオーダーレベルで改善されてるとは思わなかった...。writer様やtester様は添え字をstd::queueに持ってマージしているようです。速い解法の人はみんなstd::queue、std::deque、std::priority_queueなどを使ったいたので、強い人には常識なのかな...?

std::accmulateを用いた実装 (2948 ms)

#include <atcoder/all>
using namespace std;
using namespace atcoder;
using mint = modint998244353;

int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);
    int N, M, r = 600;
    cin >> N >> M;
    vector<int> A(N), B(N);
    for(auto &&v:A)cin >> v;
    for(auto &&v:B)cin >> v;

    vector<vector<mint>> comb(r + 1), tb(M, vector<mint>(r + 1, 1));
    for(int i = 1; i <= r; i++){
        comb[i].resize(i + 1);
        comb[i][0] = comb[i][i] = 1;
        for(int j = 1; j < i; j++) comb[i][j] = comb[i - 1][j - 1] + comb[i - 1][j];
    }

    for(int i = 0, rem = 0; i < N; i++){
        for(int j = -B[i]; j <= A[i]; j++) tb[rem][j + 300] *= comb[A[i] + B[i]][B[i] + j];
        for(int j = A[i] + 1; j <= 300; j++) tb[rem][j + 300] = 0;
        for(int j = B[i] + 1; j <= 300; j++) tb[rem][-j + 300] = 0;
        if(++rem >= M)rem -= M;
    }

    cout << accumulate(tb.begin(), tb.end(), vector<mint>(1, 1), [](vector<mint> lhs, vector<mint> rhs){
        return convolution(lhs, rhs);
    })[300 * M].val() << '\n';
}

std::reduceを用いた実装 (996 ms)

#include <atcoder/all>
using namespace std;
using namespace atcoder;
using mint = modint998244353;

int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);
    int N, M, r = 600;
    cin >> N >> M;
    vector<int> A(N), B(N);
    for(auto &&v:A)cin >> v;
    for(auto &&v:B)cin >> v;

    vector<vector<mint>> comb(r + 1), tb(M, vector<mint>(r + 1, 1));
    for(int i = 1; i <= r; i++){
        comb[i].resize(i + 1);
        comb[i][0] = comb[i][i] = 1;
        for(int j = 1; j < i; j++) comb[i][j] = comb[i - 1][j - 1] + comb[i - 1][j];
    }

    for(int i = 0, rem = 0; i < N; i++){
        for(int j = -B[i]; j <= A[i]; j++) tb[rem][j + 300] *= comb[A[i] + B[i]][B[i] + j];
        for(int j = A[i] + 1; j <= 300; j++) tb[rem][j + 300] = 0;
        for(int j = B[i] + 1; j <= 300; j++) tb[rem][-j + 300] = 0;
        if(++rem >= M)rem -= M;
    }

    cout << reduce(tb.begin(), tb.end(), vector<mint>(1, 1), [](vector<mint> lhs, vector<mint> rhs){
        return convolution(lhs, rhs);
    })[300 * M].val() << '\n';
}

yukicoder No.2149 Vanitas Vanitatum

yukicoder.me

はじめに

フック長の公式を初めて触る人が本問題の解説を自分なりに解釈しただけのものになります。表現や厳密性が怪しいところもあるかと思いますが、ご了承下さい。

問題概要

広義単調数列$A_1$、$A_2$、$\cdots$、$A_N$に対して、広義単調が保たれたまま以下の操作を行う。

    $A_i \geq 2$となる添え字$i$を選び2を引く。
    $A_i = A_{i + 1} \geq 1$となる添え字$i$を選び、それぞれ1を引く。

すべての要素を$A_i = 0$にするのに、何通りの操作があるか。

解法

はじめに、数列の「増加」、「変化なし」を01列に対応させます(増加を1、変化なしを0とする)。入力例として

6
2 2 3 5 6 6

を挙げると、

のようになります。それぞれの操作について見ていくと、1つの要素を選んで$-2$を行う操作は、01文字列のうち「110」を「011」に変化させることに相当します。

一方、隣接する2つの要素を選んで$-1$を行う操作は「100」を「001」に変化させることに相当します。

これらの操作を繰り返して、ソート完了済みの状態「000000111111」にするのに何通りの操作があるかが本問題の答えとなります。

「110」→「011」、「100」→「001」の操作はそれぞれ1つ飛ばしの「10」列を選んで「01」に変化させる操作になります。したがって、01列の偶数要素と奇数要素で分類させると、隣接する要素についてのみ考えることができ、見通しがよくなります。

ここで、偶奇分類後の01列について、それぞれの「0」の前にある「1」の個数を箱に対応させた図を書くと、上図のようになります。隣接する「10」を選んで「01」にして、ソート完了済みの状態まで操作を行うことから、箱の合計は転倒数に一致します(0の前の1の個数の和をとっているのでそれはそう...)。
奇数列と偶数列の01文字列の転倒数の合計が全体で行う操作回数の合計になりますので、奇数列の転倒数を$M_1$、偶数列の転倒数を$M_2$とすると、$M_1 + M_2$の操作回数のうち奇数列と偶数列のどちらで操作を行うかを選ぶことになりますので、本問題の答えは、 $$ \binom{M_1 + M_2}{M_1} \times (奇数列をソート完了済みにする操作方法の数)\times (偶数列をソート完了済みにする操作方法の数) $$ となります。

したがって、本問題は

01列が与えられます。隣接する要素を入れ換える操作を最小回数行って、ソート完了済みにするのに何通りの操作がありますか?

といった問題へと言い換えられました。最小回数である必要があるのは、「10」を「01」に変換する転倒数を1減らす操作しか許されておらず、「00」、「01」、「11」を入れ換えることは禁じられた操作を行っていることになるためです。

左から$i$番目に積みあがった箱の数に着目すると、これは「左から$i$番目の$0$の前にある$1$の個数」に相当します。すなわち、それぞれの$0$について何回操作を行えるかを表したものになります。
これらの箱に

    同じ行にあるものでは左にあるものほど小さい
    同じ列にあるものでは上にあるものほど小さい

といった制約のもので、何回目にその$0$についての操作を行うかの番号を書いていき、全体を埋める場合の数が操作方法の場合の数に一致します。(仮にもし、同じ行について左にあるものよりも先に小さい数を書いてしまったら、「00」のswapによって、0の追い越しが発生してしまい操作回数は最小とはならないことになります。)

話を戻すと、この左上ほど小さい数が来るように数を書き込んだ最終的な状態の数を効率よく求めることができるのがフック長の公式です。(原義のヤング図形では単調減少数列を上から並べたものに対して左上ほど小さい数を書き込むことになりますが、今回は後ろから操作を決め打っていると考えると同じものを数えていることになります。)
フック長の公式は、次元(最終的な状態数)を$d_{\lambda}$、全体の箱の数を$n$、$i$行目、$j$列目のフック長を${\rm hook}(i, j)$とすると、以下の式で表されます。

$$ d_{\lambda} = \frac{n!}{\prod {\rm hook}(i, j)} $$

箱の数$n$の階乗をフック長の総積で割ったものになります。
フック長とは、この図で言うと左にある箱の数、上にある箱の数に自分自身を加えたものになります。「101100」で表される箱について、左から2番目、下から1番目にある箱では左に1つ、上に2つ、自分自身1つでフック長は$4$となります。

これをすべての箱に対して行い、フック長の公式に代入すると、

$$ d_{\lambda} = \frac{7!}{1\times 1\times 2\times 4\times 2\times 3\times 5} = 21 $$

となり、操作番号を書き込む、最終的な状態の数は21と求められました。すなわち、奇数列の「101100」に対して最小回数の隣接swapを行いソート完了済みの状態にするには21通りの方法があると言えます。

同様に、偶数列の「100110」に対して操作方法の数をフック長の公式から求めると、

$$ d_{\lambda} = \frac{5!}{1\times 2 \times 1 \times 2 \times 5} = 6 $$ で6通りの方法があります。

したがって、奇数列の箱の数は7個、偶数列の箱の数は5個でしたので、

6
2 2 3 5 6 6

に対する操作の数は、

$$ \binom{7 + 5}{7}\times 21 \times 6 = 792 \times 21 \times 6 = 99792 $$ 通りとなります。

余談になりますが、AtCoderの「Judge System Update Test Contest 202004」のC問題でヤング図形に対する数の書き込み方の個数を求める問題が出題されていたようです。(私が初見で解いたときは、順列全探索で間に合う制約だったので特に意識しませんでした。)

atcoder.jp

同じWriterさんが2ヶ月ほど前にyukicoderに出題なさっていた「No.2092 Conjugation」なんかもフック長の計算に役立ちそうです。

yukicoder.me

(私は本問題のフック長の計算の際に、上に積まれている箱の数ではなく下に積まれている箱の数を数えてしまいサンプルは合ったのに通らないといった事態が発生しました。もし同じ事態に直面している方がいらっしゃったらこれらの類題で確認してみると良いかもしれないです。)


続いて、0通りとなる場合、すなわち、操作を行ってもすべての要素を$0$とできない場合について考えます。これは、結論から言うと、奇数列、偶数列に対してそれぞれ操作を1回行う度に数列の合計が2減ることになりますので、 $ 2\times(M_1 + M_2) \neq \sum A_i $ のときに、0通りとなります。

私はこのことについて、すべての数列の要素の和が奇数の場合、0通りになるかと考えました。もとの01列を考えると「110」→「011」、「100」→「001」となり、全体の転倒数が2減ることになります。したがって、もとの01列における0の前にある1の個数は$A_i$に一致しますので、$\sum A_i \bmod 2 \neq 1$ が条件ではないかと考えました。
しかし、「1 2 3 4」→「10101010」のように、$A_i$ の要素の合計が偶数であり、奇数列と偶数列それぞれソート済みであるのにも関わらず、0と1の数に偏りがあるため「00001111」の最終状態までもっていくことができないものがあります。したがって、全体で操作できる回数から考える必要がありました。

実装例

実装には、AtCoder Libraryのmodint998244353を用いています。
また、二項係数の分母とフック長の公式の分子で両者の階乗が打ち消されるため、二項係数のための階乗の前計算等は行っていません。

#include <atcoder/all>
using namespace std;
using mint = atcoder::modint998244353;

int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);
    int n;
    cin >> n;
    vector<int> A(n + 1);

    //01文字列の生成
    string str;
    for(int i = 1; i <= n; i++){
        cin >> A[i];
        str += string(A[i] - A[i - 1], '1') + "0";
    }

    //奇数要素と偶数要素についてそれぞれヤング図形を生成
    //s0, s1にはヤング図形の箱の合計(転倒数)が保持されている
    int cnt = 0, s0 = 0, s1 = 0;
    vector<int> odd, even;
    for(int i = 0; i < str.size(); i += 2){
        if(str[i] == '0'){
            odd.push_back(cnt);
            s0 += cnt;
        }else cnt++;
    }
    cnt = 0;
    for(int i = 1; i < str.size(); i += 2){
        if(str[i] == '0'){
            even.push_back(cnt);
            s1 += cnt;
        }else cnt++;
    }

    //条件の確認
    if(2 * (s0 + s1) != accumulate(A.begin(), A.end(), 0)){
        cout << 0 << '\n';
        return 0;
    }

    //フック長の積を計算(odd_v, even_v)に格納
    mint odd_v = 1, even_v = 1;
    vector<int> odd_cnt(1001), even_cnt(1001);
    for(int i = 0; i < odd.size(); i++){
        for(int j = 0; j < odd[i]; j++){
            odd_v *= odd[i] - j + odd_cnt[j];
            odd_cnt[j]++;
        }
    }
    for(int i = 0; i < even.size(); i++){
        for(int j = 0; j < even[i]; j++){
            even_v *= even[i] - j + even_cnt[j];
            even_cnt[j]++;
        }
    }

    mint ans = 1 / (odd_v * even_v);
    for(int i = 1; i <= s0 + s1; i++)ans *= i;

    cout << ans.val() << '\n';
}

幾何ライブラリの整備

式変形が大変だったりするのでメモ用。

2つの直線の交点

位置ベクトル$(s_1, t_1)$、$(s_2, t_2)$、方向ベクトル$(a_1, b_1)$、$(a_2, b_2)$としたときの交点座標。 方程式に直すと、 $$ y - t_1 = \frac{b_1}{a_1} (x - s_1) $$ $$ y - t_2 = \frac{b_2}{a_2} (x - s_2) $$ となる。これらの交点を求める。 $y$ について第1式と第2式をつなげる $$ \frac{b_1}{a_1} (x - s_1) + t_1 = \frac{b_2}{a_2} (x - s_2) + t_2 $$ $$ a_2b_1(x - s_1) + a_1a_2t_1 = a_1b_2(x - s_2) + a_1a_2t_2 $$ $$ -(a_1b_2 - a_2 b_1)x = (a_2b_1s_1-a_1b_2s_2) + a_1a_2(t_2 - t_1) $$ $$ x = \frac{a_1b_2s_2 - a_2 b_1 s_1 + a_1a_2(t_1 - t_2)}{a_1 b_2 - a_2 b_1} $$ これを第1式に代入してみる。 $$ y = \frac{b_1}{a_1} \left( \frac{a_1b_2s_2 - a_2 b_1 s_1 + a_1a_2(t_1 - t_2)}{a_1 b_2 - a_2 b_1} - s_1 \right) + t_1 $$ $$ y = \frac{b_1}{a_1} \left( \frac{a_1b_2s_2 - a_2 b_1 s_1 + a_1a_2(t_1 - t_2) - s_1({a_1 b_2 - a_2 b_1})}{a_1 b_2 - a_2 b_1} \right) + t_1 $$ $$ y = \frac{b_1}{a_1} \left( \frac{a_1b_2s_2 - a_2 b_1 s_1 + a_1a_2(t_1 - t_2) - a_1 b_2s_1 + a_2 b_1s_1}{a_1 b_2 - a_2 b_1} \right) + t_1 $$ $$ y = \frac{b_1}{a_1} \left( \frac{a_1b_2s_2 - a_1 b_2s_1 + a_1a_2(t_1 - t_2)}{a_1 b_2 - a_2 b_1} \right) + t_1 $$ $$ y =\frac{b_1b_2s_2 - b_1 b_2s_1 + b_1a_2(t_1 - t_2)}{a_1 b_2 - a_2 b_1} + t_1 $$ $$ y =\frac{b_1b_2(s_2 - s_1) + b_1a_2(t_1 - t_2) + t_1 (a_1 b_2 - a_2 b_1)}{a_1 b_2 - a_2 b_1} $$ $$ y =\frac{b_1b_2(s_2 - s_1) + a_1 b_2 t_1 - a_2 b_1t_2}{a_1 b_2 - a_2 b_1} $$

yukicoder No.119 旅行のツアーの問題

yukicoder.me

解法

燃やす埋める問題と呼ばれるものです。最大流(最小カット)を使います。
・$x$番目の国に行かない
・$x$番目の国に行くけど、ツアーに行かない
・$x$番目の国に行き、ツアーにも行く
という3択があります。まず、 M 個の条件がない場合は、以下のように頂点を倍加して、辺を張ることでこの3択を実現できます。
ここでは、$s$を始点、$t$を終点とします。

今回は、
①ツアーに行く満足度$A$
②ツアーに行く満足度+ツアーに行かないけど国に行く満足度$A+B$
③ツアーに行かないけど国に行く$B$
という順番で上流から辺を張りました。順番は任意で構いませんが、M個の条件を追加する際に真ん中が$A+B$の容量の辺だと都合が良いと思います。答えは$A+B$を持っていて、そこから流量の分だけ損失があると考えます。
それぞれがカットされた場合を考えます。

①$A$がカットされた(すなわち流量$A$が流れた)
$A$がカットされた(すなわち流量$A$が流れた)場合は、答えは$A+B-A=B$になります。すなわち、国に行くけどツアーに行かないのが最適ということになります。

②$A+B$がカットされた(すなわち流量$A+B$が流れた)
$A+B$がカットされた(すなわち流量$A+B$が流れた)場合は、答えは$A+B-(A+B)=0$になります。すなわち、国に行かないことが最適となります。(M個の条件がない場合ではこの選択肢は取り得ないことが分かるかと思います。)

③$B$がカットされた(すなわち流量$B$が流れた)
$B$がカットされた(すなわち流量$B$が流れた)場合は、答えは$A+B-B=A$になります。すなわち、国に行きツアーにも参加することが最適となります。

言い換えると、それぞれの辺をカットすることは以下の選択肢を選んだことに相当することが分かるかと思います。 ①国に行くけどツアーに行かない
②国に行かない
③国に行きツアーにも行く

そして、M個の条件を考えます。$u$<$v$の条件が成り立っており、$u$、$v$両方の国を訪れる場合、$u$の国でツアーを利用したら$v$の国でもツアーを利用しなければならないことにします。
結論から言ってしまうと、以下のように$u'$から$v$に容量が十分に大きいINFの辺を張ることで達成できます。 仮に、国$u$に行きツアーにも行く選択をした後、国$v$に行くけどツアーに行かない選択をした場合、容量が$B_u$と$A_v$の辺がカットされますが、$s \to u\to u' \to v \to v' \to t$の経路でまだ流量を流すことができますので、最適でないことが分かります。 これはすなわち、国$u$で③の選択肢を取ったならば、国$v$では②か③の選択肢を選ばなければいけないといったような状況に相当します。 今回は、国$u$で③の選択肢を取ったならば条件先の国$v$では②か③の2択になるといった具合でしたが、$u'\to v'$などに容量INFの辺を張ると③を選んだら条件先では③の一択になるなどといった応用もできます。容量INFの辺の始点よりも下流にある辺をカットする場合は、容量INFの辺の終点よりも下流の辺をカットしなければならないということが一般的に言えます。
逆に、国$v$で①の選択肢を取ったらならば、国$u$では①か②の2択になるといった見方もできます。容量INFの辺の終点よりも上流にある辺をカットする場合は、容量INFの辺の始点よりも上流の辺をカットしなければならないということが一般的に言えます。 辺の向きと始点と終点を逆にしたものも流量が同じになることが分かるかと思いますが、これも国$u$でINFの辺の終点より上流の選択肢(③)を選んだならば、国$v$ではINFの辺の始点よりも上流の選択肢(②か③)を選ばなければいけないと解釈することができます。(最初に、①、②、③の辺の順序は②が真ん中ならば任意で良いと書かせて頂いたのもこのためです。)

まとめると、
・INFの辺の始点よりも下流の辺をカットしたならば、INFの辺の終点よりも下流の辺をカットしなければならない
・INFの辺の終点よりも上流の辺をカットしたならば、INFの辺の始点よりも上流の辺をカットしなければならない
といったことを考えることで、M個の条件を追加することができました。
これらが同じことを指しているのは以下の2つが対偶関係となっているためです。
・INFの辺の始点よりも下流の辺をカットするならば、INFの辺の終点よりも下流の辺をカットする
・INFの辺の終点よりも下流の辺をカットしないならば、INFの辺の始点よりも下流の辺をカットしない

今回の問題で言うと、以下の2つが対偶関係にあることが分かるかと思います。
・国$u$でツアーに行ったならば、国$v$に行くときはツアーに行く
・国$v$に行くときにツアーに行かなかったならば、国$u$でツアーに行ってはいけない

実装例

実装にはAtCoder Libraryのmax_flowを用いています。

#include<atcoder/all>
using namespace std;

int main(){
    int n, m, u, v;
    cin >> n;
    atcoder::mf_graph<int> g(2 * n + 2);
    int s = 2 * n, t = s + 1, ans = 0;
    for(int i = 0; i < n; i++){
        cin >> u >> v;
        ans += u + v;
        g.add_edge(s, i, u);
        g.add_edge(i, i + n, u + v);
        g.add_edge(i + n, t, v);
    }
    cin >> m;
    for(int i = 0; i < m; i++){
        cin >> u >> v;
        g.add_edge(u + n, v, 1 << 30);
    }
    ans -= g.flow(s, t);
    cout << ans << endl;
}

yukicoder No.62 リベリオン(Extra)

yukicoder.me

解法の概要

中国余剰定理を用いた多倍長整数を必要としない解法です。

解法

密室の幅を$W$、奥行きを$H$、時間を$D$秒間止め、まみの座標を$ (M_x, M_y) $、ほむらの座標を$(H_x, H_y)$ 、弾丸の一秒間に進む速度を$ (V_x,V_y) $とします。
まず、$(V_x,V_y)$について最大公約数をとります。例えば、1秒後に$(2, 4)$という速度であった場合、$(1, 2)$は0.5秒後に到達することになります。整数で扱うには都合が悪いため、速度を$(V_x/{\rm gcd}(V_x,V_y), V_y/{\rm gcd}(V_x,V_y))$とすることで整数座標には整数時間にしか訪れないという変換ができます。当然速度が${\rm gcd}(V_x,V_y)$分の1となるため、時間$D$を$D*{\rm gcd}(V_x,V_y)$としてあげます。
以下、$V_x$,$V_y$が登場した際は $V_x/{\rm gcd}(V_x, V_y)$、$V_y/{\rm gcd}(V_x, V_y)$を指すものとします。

この密室を無限に広い空間に変換します。これは、$(2W, 2H)$の空間が周期的に現れていると解釈することで可能です。($(W, H)$を横断する度に鏡のように反転しているというイメージです。)
$(2W, 2H)$の中に$(W, H)$の領域が4つあって、それらが鏡のような対称性をもっていると考えると、まみは${\rm mod}(2W, 2H)$として、$(M_x, M_y)$、$(-M_x, M_y)$、$(M_x, -M_y)$、$(-M_x, -M_y)$の4点にいることが分かるかと思います。

まみが$(M_x, M_y)$にいる場合を考えます。時間を$t$とすると、弾丸の座標とまみの座標が一致する条件は $$ M_x = V_x t + H_x \,\,({\rm mod}\,2W) $$ $$ M_y = V_yt + H_y \,\,({\rm mod}\,2H) $$ となります。時刻$t$について変形すると、 $$ t = \frac{M_x - H_x}{V_x}\,\,({\rm mod}\,2W) = \frac{M_y - H_y}{V_y}\,\,({\rm mod}\,2H) $$ となります。この際、$V_x$と$2W$、$V_y$と$2H$は互いに素にしておく必要があります(互いに素でないと逆元を求めることができない)ので、それぞれ、${\rm gcd}(V_x, 2W)$と${\rm gcd}(V_y, 2H)$をとっておき、それぞれ
$$ V_x \leftarrow \frac{V_x}{{\rm gcd}(V_x, 2W)} m_x = \frac{M_x - H_x}{{\rm gcd}(V_x, 2W)} D_x = \frac{2W }{ {\rm gcd}(V_x, 2W)} $$ $$ V_y \leftarrow \frac{V_y}{{\rm gcd}(V_y, 2H)} m_y = \frac{M_y - H_y}{{\rm gcd}(V_y, 2H)} D_y = \frac{2W}{ {\rm gcd}(V_y, 2H)} $$ としておきます。$M_x - H_x$が${\rm gcd}(V_x, 2W)$で割り切れない場合や$M_y - H_y$が${\rm gcd}(V_y, 2H)$で割り切れない場合は逆元が存在しないため解が存在しないことになります。
これらをもとの式に当てはめると、 $$ t = \frac{m_x}{V_x}\,\,({\rm mod}\,D_x) = \frac{m_y}{V_y}\,\,({\rm mod}\,D_y) $$ となりますので、この$t$ を中国余剰定理で求めることにより、$t \leq D \times {\rm gcd}(V_x, V_y) $となっているかを判定すればこの問題を解くことができます。他のまみの座標$(-M_x, M_y)$、$(M_x, -M_y)$、$(-M_x, -M_y)$についても同様にして解くことができます。

実装例

実装にはAtCoder Libraryのinv_modとcrtを用いています。

#include<atcoder/all>
using namespace std;
using namespace atcoder;

int main(){
    int t;
    cin >> t;
    while(t--){
        long long W, H, D, Mx, My, Hx, Hy, Vx, Vy, g, gx, gy, Dx, Dy, mx, my;
        cin >> W >> H >> D >> Mx >> My >> Hx >> Hy >> Vx >> Vy;
        g = gcd(Vx, Vy);
        Vx /= g, Vy /= g, D *= g;
        gx = gcd(Vx, 2 * W), gy = gcd(Vy, 2 * H);
        Vx /= gx, Dx = 2 * W / gx;
        Vy /= gy, Dy = 2 * H / gy;
        bool ans = false;
        for(int i = 0; i < 4; i++){
            mx = (i & 1? 1: -1) * Mx - Hx;
            my = (i & 2? 1: -1) * My - Hy;
            if(mx % gx != 0 || my % gy != 0)continue;
            mx /= gx, my /= gy;
            mx *= inv_mod(Vx, Dx);
            my *= inv_mod(Vy, Dy);
            auto T = crt({mx, my}, {Dx, Dy});
            if(T.second == 0)continue;
            if(T.first <= D)ans = true;
        }
        cout << (ans ? "Hit" : "Miss") << '\n';
    }
}