rep マクロでハマった

注意: この記事は C++ 初心者が書いたものです.

前置き

最近競プロで C++ を使い始めたのですが,まずは人の真似からということで便利マクロの代表格であろう rep マクロを導入していました.rep マクロにはいろいろなバリエーションがあると思いますが,私はこれを次のように define していました.

#define rep(i, n) for (int i = 0; i < (n); ++i)

 0\leq i \lt n を満たす整数  i を昇順に走査するマクロです.そして,この派生として  l\leq i\lt r を満たす整数  i を昇順に走査する replr なるマクロも定義しており,その定義は次のようなものでした.

#define replr(i, l, r) for (int i = (l); i < (r); ++i)

既に予想できた人もいるかもしれないんですが,i の型を int で固定していたことが原因でバグを埋め込んでしまいました.というのも,例えば lrlong long であり,int に収まらない値であった場合に意図通りの動作をしてくれません.

rep(i, n)nint に収まらないことなんてそうそうないのでまだいいんですが,replr の方は int に収まらないことがたまにあるので直します.

解決策

常に long long にしておけば良いのはそうなんですが,不必要に大きなサイズの型を使うのは個人的に嫌なので,定義を次のように変更してみました.

#define replr(i, l, r) for (decltype(r) i = 0; i < (r); ++i)

decltyper の型を取得しています.しかし,これではうまくいかないことがあります.例えば,次のような場合です.

#include <iostream>

constexpr int N = 10;

int main() {
    for (decltype(N) i = 0; i < N; ++i) {
        std::cout << i << std::endl;
    }
}

Nconst なので decltype(N) i = 0; の宣言で i の型は const int となり,++i の部分で値を変更しようとしているためコンパイルエラーとなります.

そこで,もう一工夫して次のように定義を変更してみます.

#include <type_traits>

#define replr(i, l, r) for (std::remove_const_t<decltype(r)> i = l; i < r; ++i) {

これでめでたく動きました.

おわりに

私は C++ のことを何も分かっていないので,「こうした方がいい」,「これはマズイ」みたいなことがあれば教えてください.

第一種スターリング数の末尾項を計算する

こどふぉブログ https://codeforces.com/blog/entry/47764 の最後で紹介されているテクニックです.第一種スターリング数  \begin{bmatrix}N\\ N-i\end{bmatrix}\bmod P を全ての  i=0,\ldots,K に対して  O(K ^ 2\log N +\log P) で計算できます.

解説

以下, a _ {N, i}:=\begin{bmatrix}N\\ N-i\end{bmatrix} とおきます.

 \displaystyle
\mathcal{S}_i(l,r)=\{S\subset\{l,\ldots,r-1\}\mid |S|=i\}

とすると,任意の非負整数  n に対して

 \displaystyle
\begin{align}
a_{2n,i}
&=\sum_{S\in\mathcal{S}_i(0,2n)}\prod_{v\in S}v\\
&=\sum_{j=0}^i\left(\sum_{S\in\mathcal{S}_{i-j}(0,n)}\prod_{v\in S}v\right)\left(\sum_{S\in\mathcal{S}_j(n,2n)}\prod_{v\in S}v\right)\\
&=\sum_{j=0}^ia_{n,i-j}\left(\sum_{S\in\mathcal{S}_j(0,n)}\prod_{v\in S}(n+v)\right)\\
&=\sum_{j=0}^ia_{n,i-j}\sum_{k=0}^j\begin{pmatrix}n-k\\ j-k\end{pmatrix}n^{j-k}\left(\sum_{S\in\mathcal{S}_k(0,n)}\prod_{v\in S}v\right)\\
&=\sum_{j=0}^ia_{n,i-j}\sum_{k=0}^j\begin{pmatrix}n-k\\ j-k\end{pmatrix}n^{j-k}a_{n,k}\\
a_{2n+1,i}&=\begin{cases}
1&(i=0)\\
a_{2n,i}+2n\cdot a_{2n,i-1}&(i\gt 0)\
\end{cases}
\end{align}

が成り立つのでダブリング的に  a _ n を求めていくことが出来ます.

このままでは  a _ n から  a _ {2n} を求めるのに  \Theta(K ^ 3) 掛かってしまいそうですが, \displaystyle b _ j := \sum _ {k=0} ^ j\begin{pmatrix}n-k \\ j - k \end{pmatrix} n ^ {j - k} a _ {n, k} i に依存しないことに注目します.すると,これは  O(K ^ 2) で前計算できるので,全体としても  O(K ^ 2) a _ {2n} を求めることができます.

ここで, n が大きい場合は階乗の前計算が出来ないため二項係数の計算が少し難しくなる点に注意です. b _ j に集める形 ( j\rightarrow k 順のループ) で書こうとすると  n-k の値が動いて厄介なので, b _ j に配る形 ( k\rightarrow j 順のループ) で書くと楽だと思います.

以上を実装したのが次になります. i=1,\ldots,K \mathrm{mod}\ P における乗法逆元の計算をサボっているので計算量は  \Theta(K ^ 2\log N + K\log P) となっていますが,実用上は誤差レベルだと思います.

折り畳み

実装例 (C++)

注意:  k\lt P を仮定しています.これが満たされない場合,逆元が存在しないことがあるので二項係数の計算がうまくいきません.

/**
 * return:
 *   vector<modint> v s.t. v[i] = S1[n,n-i] for i=0,...,k, where S1 is the stirling number of the first kind (unsigned).
 * constraints:
 * - 0 <= n <= 10^18
 * - 0 <= k <= 5000
 * - k < mod
 */
template <typename modint>
std::vector<modint> last_terms_of_stirling_numbers_of_the_first_kind(unsigned long long n, size_t k) {
    std::vector<modint> invs(k + 1);
    for (size_t i = 1; i <= k; ++i) invs[i] = modint(1) / i;
    std::vector<modint> a(k + 1, 0);
    a[0] = 1;
    size_t bit_length = 0;
    while (n >> bit_length != 0) ++bit_length;
    unsigned long long m = 0;
    for (size_t bit = bit_length; bit --> 0;) {
        std::vector<modint> b(k + 1, 0);
        for (size_t j = 0; j <= k; ++j) {
            modint tmp = 1;
            for (size_t i = j; i <= k; ++i) {
                b[i] += a[j] * tmp;
                tmp *= (m - i) * invs[i - j + 1] * m;
            }
        }
        for (size_t i = k + 1; i --> 0;) {
            modint sum = 0;
            for (size_t j = 0; j <= i; ++j) sum += a[j] * b[i - j];
            a[i] = sum;
        }
        m <<= 1;
        if ((n >> bit) & 1) {
            for (size_t i = k; i > 0; --i) a[i] += m * a[i - 1];
            ++m;
        }
    }
    return a;
}

verify: https://codeforces.com/contest/1516/submission/114464936


以上が元の記事に書かれていたことを私なりに再構築して説明した部分になります.そして,ここからは蛇足です.

 a _ {2n, i} を求める式を眺めていると,FFT で高速化出来そうな感じがします.実際,

 \displaystyle
\begin{align}
a_{2n,i}
&=\sum_{j=0}^ia_{n,i-j}\sum_{k=0}^j\begin{pmatrix}n-k\\ j-k\end{pmatrix}n^{j-k}a_{n,k}\\
&=\sum_{j=0}^ia_{n,i-j}\cdot\dfrac{1}{(n-j)!}\sum_{k=0}^j\dfrac{n^{j-k}}{(j-k)!}\cdot a_{n,k}\cdot(n-k)!\\
&=\sum_{j=0}^ia_{n,i-j}\cdot\dfrac{[x^j]f_n g_n}{(n-j)!}\quad\left(\begin{aligned}
f_n:=&\sum_{k=0}^K \dfrac{n^k}{k!}x^k\\
g_n:=&\sum_{k=0}^K a_{n,k}\cdot(n-k)!\cdot x^k
\end{aligned}\right)\\
&=[x^i]a_nh_n\quad\left(h_n:=\sum_{j=0}^K\dfrac{[x^j]f_ng_n}{(n-j)!}x^j\right)
\end{align}

より  a _ {2n} = a _ n h _ n が成り立ち,2 回の畳み込みで  a _ {2n} が求まりそうです.しかし,式をよく見ると, g _ n h _ n の係数を求める際に大きな値 (最大で  N/2 程度) の階乗が必要になっています.

これは簡単には求まりませんが, O(\sqrt{P}(\log P) ^ 2) あるいは  O(\sqrt{P}\log P) で求めるアルゴリズムが存在します 1.従って,これを用いれば全体  O((K\log K + \sqrt{P}\log P)\log N) でこの問題を解くことが出来ます.

ただし, n P 以上となることがある場合は階乗の逆元が存在しない可能性があるためこの方法は使えません

おわりに

私は E8 さんの競プロ典型 90 問の 1 問 005 - Restricted Digits(★7) で DP テーブルをダブリング的に求める高速化手法を知りました.そして,丁度良いタイミングで同じ考え方を使った高速化の例を見つけたので書いてみました.正直思いつける気はしませんが...


  1. 解説記事があったんですが,消えてしまった上に私は他の解説記事を知らないのでリンクを張れません…

下に凸な数列同士の (min, +) convolution

より制約を弱くして一方のみが下に凸だとした場合も SMAWK Algorithm を使えば畳み込む数列の長さの和に対して線形時間で求まるんですが,両方ともが下に凸という制約を課した場合はかなり簡単に,かつ定数倍もかなり軽いアルゴリズムで求まることが分かったのでこの記事を書きました.

問題

長さ  N の数列  A=(A _ 0,\ldots,A _ {N-1}) と長さ  M の数列  B=(B _ 0,\ldots, B _ {M-1}) が与えられます.このとき,すべての  0\leq k\leq N + M - 2 に対して

 \displaystyle
C _ k = \min _ {i + j = k} A _ i + B _ j

を求めてください.ただし, A,B はともに下に凸であるとします.つまり,次の 2 つの条件を満たします.

条件:

  • 全ての  1\leq i\leq N-2 に対して  A _ i - A _ {i - 1} \leq A _ {i + 1} - A _ i
  • 全ての  1\leq j\leq M-2 に対して  B _ j - B _ {j - 1} \leq B _ {j + 1} - B _ j

解法

 k に対して, i + j = k の下で  A _ i + B _ j を最小化する  i, j をそれぞれ  i _ k, j _ k と表す.ただし,そのような  i, j の組が複数ある場合は,適当に一つだけ選ぶことにする.

まず, k=0 の場合を考える.このときは,明らかに  i _ k = j _ k = 0 である.

続いて, k\gt 0 の場合を考える.このとき,実は  (i _ k, j _ k) = (i _ {k - 1}, j _ {k - 1} + 1) または  (i _ k, j _ k) = (i _ {k - 1} + 1, j _ {k - 1}) の 2 通りを試すだけでよい.これは,以下のような場合分けから従う.

  •  i \geq i _ {k - 1} + 1 のとき
 \displaystyle
\begin{aligned}
A _ i + B _ {k - i} 
&= (A _ {i - 1} + B _ {k - i}) + (A _ i - A _ {i - 1}) \\
&\geq (A _ {i _ {k - 1}} + B _ {j _ {k - 1}}) + (A _ {i _ {k - 1} + 1} - A _ {i _ {k - 1}}) \\
&=A_{i_{k-1}+1}+B_{j_{k-1}}
\end{aligned}
  •  j \geq j _ {k - 1} + 1 のとき
 \displaystyle
\begin{aligned}
A _ {k - j} + B _ j 
&= (A _ {k - j} + B _ {j - 1}) + (B _ j - B _ {j - 1}) \\
&\geq (A _ {i _ {k - 1}} + B _ {j _ {k - 1}}) + (B _ {j _ {k - 1} + 1} - B _ {j _ {k - 1}}) \\
&=A_{i_{k-1}}+B_{j_{k-1} + 1}
\end{aligned}

この事実を用いることで, k=0,\ldots,N+M-2 に対して  C _ k \Theta (N + M) 時間で求めることが出来る.

実装

std::vector<int> convolve(const std::vector<int> &a, const std::vector<int> &b) {
    constexpr int INF = std::numeric_limits<int>::max();
    int n = a.size();
    int m = b.size();
    if (n == 0 or m == 0) return {};
    std::vector<int> c(n + m - 1);
    c[0] = a[0] + b[0];
    int i = 0, j = 0;
    for (int k = 1; k <= n + m - 2; ++k) {
        int x = i + 1 == n ? INF : a[i + 1] + b[j];
        int y = j + 1 == m ? INF : a[i] + b[j + 1];
        if (x <= y) {
            c[k] = x;
            ++i;
        } else {
            c[k] = y;
            ++j;
        }
    }
    return c;
}

おまけ (片方のみが下に凸な場合)

以下, A は任意の数列, B は下に凸な数列としても一般性を失いません (逆なら swap すればよいので).

 (M+N-1) \times N 行列  P を次の図のように構築します.ただし, P を陽に持った時点で計算量が  O( (N+M) M) となってしまうので,陰に持ちます (つまり,必要になった時に計算する).

f:id:suisen_kyopro:20210527201819p:plain

このとき,結論から言えば, P は totally monotone であることが示せます.従って,SMAWK Algorithm を用いることですべての  k に対する  \underset{i}{\arg \min}\, P _ {k, i} \Theta(N + M) で求めることができ,同時に  C _ k も求まります.SMAWK Algorithm に関しては,参考文献 [1] が分かりやすいので,詳細は省かせて頂きます (不等号の誤植があった気はしますが).

SMAWK Algorithm の説明を省いた代わりに,ここでは  P が totally monotone であることを確認しておきましょう.これは,任意の 2\times 2 部分行列が monotone であることが言えればよいです.

 \displaystyle
\begin{pmatrix}
a & b\\
c & d
\end{pmatrix}

ここで,上のような  2\times 2 行列に対して,次のいずれかを満たす場合に monotone であるといいます.

  1.  c\gt d
  2.  c=d かつ  a\leq b
  3.  c\lt d かつ  a\lt b

つまり, 0\leq k _ 1\lt k _ 2\lt N + M - 1,\ 0\leq i _ 1\lt i _ 2\lt N を満たす  k _ 1, k _ 2, i _ 1, i _ 2 であって,部分行列  \displaystyle
\begin{pmatrix}
P _ {i _ 1, k _ 1}\;\, P _ {i _ 2, k _ 1} \\
P _ {i _ 1, k _ 2}\;\, P _ {i _ 2, k _ 2}
\end{pmatrix}
に関して

  1.  P _ {i _ 1, k _ 2}=P _ {i _ 2, k _ 2} かつ  P _ {i _ 1, k _ 1}\gt P _ {i _ 2, k _ 1}
  2.  P _ {i _ 1, k _ 2}\lt P _ {i _ 2, k _ 2} かつ  P _ {i _ 1, k _ 1}\geq P _ {i _ 2, k _ 1}

となるようなものが存在しないことを言えればよいです.まず, \infty に対して適切な値を入れておくことで, \infty を含むような  2\times 2 部分行列を monotone にできます 1.従って,以下では部分行列が  \infty を含まない場合 ( B が配列外参照を起こさない場合) を考えます.

 B が配列外参照を起こさないという仮定の下では,部分行列は次のように表すことが出来ます.

 \displaystyle
\begin{pmatrix}
P _ {i _ 1, k _ 1}& P _ {i _ 2, k _ 1} \\
P _ {i _ 1, k _ 2}& P _ {i _ 2, k _ 2}
\end{pmatrix}=
\begin{pmatrix}
A _ {i _ 1} + B _ {k _ 1 - i _ 1} & A _ {i _ 2} + B _ {k _ 1 - i _ 2} \\
A _ {i _ 1} + B _ {k _ 2 - i _ 1} & A _ {i _ 2} + B _ {k _ 2 - i _ 2}
\end{pmatrix}

では,以下のいずれかが成立するという仮定を置いて,矛盾を導きます.

  1.  P _ {i _ 1, k _ 2}=P _ {i _ 2, k _ 2} かつ  P _ {i _ 1, k _ 1}\gt P _ {i _ 2, k _ 1}
  2.  P _ {i _ 1, k _ 2}\lt P _ {i _ 2, k _ 2} かつ  P _ {i _ 1, k _ 1}\geq P _ {i _ 2, k _ 1}

方針としては, i _ 1\lt i _ 2,\ k _ 1\lt k _ 2 に注意しながら, B が下に凸であるという性質を活かせるように式変形を頑張るだけです.

 \displaystyle
\begin{aligned}
& \begin{cases}
P _ {i _ 1, k _ 1}\gt P _ {i _ 2, k _ 1}\\
P _ {i _ 1, k _ 2}=P _ {i _ 2, k _ 2}
\end{cases}\quad \text{or}\quad
\begin{cases}
P _ {i _ 1, k _ 1}\geq P _ {i _ 2, k _ 1}\\
P _ {i _ 1, k _ 2}\lt P _ {i _ 2, k _ 2}
\end{cases} \\
\Rightarrow &
\begin{cases}
A _ {i _ 1} + B _ {k _ 2 - i _ 1}\gt A _ {i _ 2} + B _ {k _ 2 - i _ 2}\\
A _ {i _ 1} + B _ {k _ 1 - i _ 1}\leq A _ {i _ 2} + B _ {k _ 1 - i _ 2}
\end{cases}\quad \text{or}\quad
\begin{cases}
A _ {i _ 1} + B _ {k _ 2 - i _ 1}\geq A _ {i _ 2} + B _ {k _ 2 - i _ 2}\\
A _ {i _ 1} + B _ {k _ 1 - i _ 1}\lt A _ {i _ 2} + B _ {k _ 1 - i _ 2}
\end{cases} \\
\Rightarrow &
A _ {i _ 1} + B _ {k _ 2 - i _ 1} - B _ {k _ 2 - i _ 2} \lt A _ {i _ 1} + B _ {k _ 1 - i _ 1} - B _ {k _ 1 - i _ 2} \\
\Rightarrow &
B _ {k _ 2 - i _ 1} - B _ {k _ 2 - i _ 2}\lt B _ {k _ 1 - i _ 1} - B _ {k _ 1 - i _ 2} \\
\Rightarrow &
\begin{cases}
B _ {k _ 2 - i _ 1} - B _ {k _ 2 - i _ 2}\lt B _ {k _ 1 - i _ 1} - B _ {k _ 1 - i _ 2} \\
B _ {k _ 2 - i _ 1} - B _ {k _ 1 - i _ 1}\lt B _ {k _ 2 - i _ 2} - B _ {k _ 1 - i _ 2} \\
\end{cases} \\
\Rightarrow &
\begin{cases}
\displaystyle
\sum _ {j = k _ 2 - i _ 2} ^ {k _ 2 - i _ 1 - 1} B _ {j + 1} - B _ {j}
\lt
\sum _ {j = k _ 1 - i _ 2} ^ {k _ 1 - i _ 1 - 1} B _ {j + 1} - B _ {j} \\
\displaystyle
\sum _ {j = k _ 1 - i _ 1} ^ {k _ 2 - i _ 1 - 1} B _ {j + 1} - B _ {j}
\lt
\sum _ {j = k _ 1 - i _ 2} ^ {k _ 2 - i _ 2 - 1} B _ {j + 1} - B _ {j} \\
\end{cases} \\
\Rightarrow &
\begin{cases}
(i _ 2 - i _ 1)(B _ {k _ 2 - i _ 2 + 1} - B _ {k _ 2 - i _ 2})
\lt
(i _ 2 - i _ 1)(B _ {k _ 1 - i _ 1 + 1} - B _ {k _ 1 - i _ 1})\\
\displaystyle
(k _ 2 - k _ 1)(B _ {k _ 1 - i _ 1 + 1} - B _ {k _ 1 - i _ 1})
\lt
(k _ 2 - k _ 1)(B _ {k _ 2 - i _ 2 + 1} - B _ {k _ 2 - i _ 2})
\end{cases} \\
\Rightarrow &
\begin{cases}
B _ {k _ 2 - i _ 2 + 1} - B _ {k _ 2 - i _ 2}
\lt
B _ {k _ 1 - i _ 1 + 1} - B _ {k _ 1 - i _ 1}\\
\displaystyle
B _ {k _ 1 - i _ 1 + 1} - B _ {k _ 1 - i _ 1}
\lt
B _ {k _ 2 - i _ 2 + 1} - B _ {k _ 2 - i _ 2}
\end{cases} \\
\end{aligned} \\

以上より,矛盾が示せたので  P は totally monotone であることが分かりました.

参考文献

[1] http://web.cs.unlv.edu/larmore/Courses/CSC477/monge.pdf


  1. 右上の区画では  i に対して単調増加な値を,左下の区画では  i に対して単調減少な値を入れるとよいです

京セラプログラミングコンテスト2021(AtCoder Beginner Contest 200)

なんかコンテストが始まった瞬間 Chrome 君がフリーズしてめちゃくちゃ焦った.

A - Century

答えは  \left\lfloor\dfrac{N-1}{100}\right\rfloor+1

実装 (Python)

print((int(input()) - 1) // 100 + 1)

B - 200th ABC-200

毎回愚直に操作してよい.多くの言語は % 演算子で剰余を取ることが出来る.

実装

N, K = map(int, input().split())
for _  in range(K):
    if N % 200:
        N = N * 1000 + 200
    else:
        N //= 200
print(N)

C - Ringo's Favorite Numbers 2

 A _ i - A _ j 200 の倍数という条件は  A _ i \equiv A _ j\pmod{200} と言い換えられるので, \mathrm{cnt} [ i ] = \# \{j\mid A _ j \bmod 200 = i \} を計算しておく.すると,答えは

 \displaystyle
\sum_{j=0}^{199}\frac{\mathrm{cnt} [ i ](\mathrm{cnt} [ i ] - 1)}{2}

と表される.

実装 (Python)

N = int(input())
A = list(map(int, input().split()))
cnt = [0] * 200
for v in A:
    cnt[v % 200] += 1
ans = 0
for c in cnt:
    ans += c * (c - 1) // 2
print(ans)

D - Happy Birthday! 2

私の解法は after_contest で落ちました...

E - Patisserie ABC 2

 k 個の非負整数  A _ 1,\ldots, A _ k であって, 0\leq A _ i\lt N\;(i=1,\ldots,k) かつ和が  M であるようなものの個数は

 \displaystyle
[x^M](1+\cdots x^{N-1})^k=[x^M]\left(\frac{1-x^N}{1-x}\right)^k

と表せる.今回は  k=1,2,3 に対してこの多項式を ( \mathrm{mod}\; x ^ {3N} で) 予め計算しておくことで, i+j+k\to i\to j\to k の順に全体  O(N) で確定していくことができる.

 1 - x ^ N を掛けるのは各係数  c _ i に対して  c _ i\leftarrow c _ i - c _ {i - N} とすればよく, 1 - x で割るのは  c _ i\leftarrow c _ 0 + \cdots + c _ i とすればよい.in-place に計算を行う場合は,前者は次数の降順に,後者は次数の昇順に計算を行うようにすると良い.

なお,数学をすることで各多項式の係数の prefix sum あるいは suffix sum は  O(1) で求められるので,二分探索により  O(\log N) で解くことも可能である.

実装 (Python)

N, K = map(int, input().split())
f1 = [0] * (3 * N + 1)
f1[0] = 1
f1[N] = -1
for i in range(3 * N):
    f1[i + 1] += f1[i]
f2 = f1[:]
for i in range(3 * N):
    f2[i + 1] += f2[i]
for i in range(3 * N, N - 1, -1):
    f2[i] -= f2[i - N]
f3 = f2[:]
for i in range(3 * N):
    f3[i + 1] += f3[i]
for i in range(3 * N, N - 1, -1):
    f3[i] -= f3[i - N]
for x in range(3 * N - 2):
    if K > f3[x]:
        K -= f3[x]
        continue
    for i in range(min(N, x + 1)):
        if K > f2[x - i]:
            K -= f2[x - i]
            continue
        for j in range(min(N, x - i + 1)):
            if K > f1[x - i - j]:
                K -= f1[x - i - j]
                continue
            k = x - i - j
            print(i + 1, j + 1, k + 1)
            exit(0)

F - Minflip Summation

まず,?0 または 1 に決めた後の問題の答えを考える.文字列を円環として見る (即ち,先頭の文字と末尾の文字も隣接しているとみなす) ことにすると,「全て同じ文字である」という条件は,「隣接する文字の xor を取ってできる文字列が全て  0 である」と言い換えられる.つまり,元の文字列  S に対して隣接 xor を取った文字列  T を以下のように定めたとき, T の全ての項が  0 となるまでに必要な最小の操作回数を考えればよい.

 \displaystyle
T_i=S_{i-1}\oplus S_{i}\quad(i=0,\ldots,|S|-1)

ここで, \oplus は xor 演算, S _ {-1} = S _ {|S| - 1} とする.

 S に対する区間 flip 操作は, T において,任意の  2 項を選択して flip する操作となる.従って,最小の操作回数は  T _ i = 1 となる  i の個数を  2 で割ったものに等しい. T _ i = 1 となる  i の個数が奇数だと困ったことになるが,実は必ず偶数になることが証明できる.厳密な証明は置いておくことにしてお気持ちを説明すると,円環  S S _ 0 を始点に時計回りに順番に見ていった場合に,仮に変化点が奇数個であれば,最後に  S _ 0 に戻ってきたときに  S _ 0 とは異なる文字を見ているはずなので矛盾する,という寸法である.

? を決め打ちした場合の答えが分かったので,元の問題を考える. S の隣接項の関係のみが重要であるということが分かったので, S _ {i - 1}, S _ i の答えへの寄与を各々考えると,次のようになる.ただし, q S _ i = ? となる  i の個数である.

 S _ {i - 1}  S _ i 寄与
 0  0  0
 0  1  k\cdot 2^{kq}
 0 ?  k\cdot 2^{kq-1}
 1  0  k\cdot 2^{kq}
 1  1  0
 1 ?  k\cdot 2^{kq-1}
?  0  k\cdot 2^{kq-1}
?  1  k\cdot 2^{kq-1}
? ?  k\cdot 2^{kq-1}

全てを説明していると冗長になるので最後のパターンだけ説明する. S _ {i - 1} = S _ i = ? のとき,答えに正の寄与をもたらすのは  S _ {i - 1} = 0, S _ i = 1 および  S _ {i - 1} = 1, S _ i = 0 のケースであり,他の  kq - 2 個の ? の決め方はそれぞれ  2 ^ {kq-2} 通りずつある.従って,この場合の寄与は  2\cdot 2 ^ {kq-2} = 2 ^ {kq - 1} と計算できる.表で  k が掛かっているのは,文字列が  k 回繰り返されるためである.以上より,表に従って寄与の和を計算し,最後に  2 で割ったものが答えである...

と言いたいところだが,実はコーナーケースが存在する. S = ?,  K = 1 の場合,上で説明した  S _ {i - 1} = S _ i = ? での議論が破綻する.というのは,最終的にできる円環において,先頭と末尾が同一項を指してしまうためである.このとき,先頭と末尾が異なるような ? への割り当ては存在しないので,答えへの寄与は  0 となるはずである.しかし,表中の式に従って寄与を計算すると  1 となり間違った答えを出してしまう.同様に  S= 0 S= 1 のケースでも先頭と末尾の項が同一項を指すことになるが,この場合はもともと寄与が  0 であったので不具合は起こらない.

 K\gt 1 であったり,あるいは  |S|\gt 1 であれば最終的にできる円環の長さが  2 以上となるのでこのようなことは起こらない.

実装 (Python)

以下は  P := 10 ^ 9 + 7 として全体  O(|S| + \log K + \log P) で動作する.

P = 10 ** 9 + 7

def solve(S: str, k: int):
    if S == '?' and k == 1:
        return 0

    ans = 0
    for prv, cur in zip(S[-1] + S, S):
        if prv == '?' or cur == '?':
            ans += 1
        elif prv != cur:
            ans += 2
    return ans * k * pow(2, k * S.count('?') - 2, P) % P

if __name__ == '__main__':
    S = input()
    k = int(input())
    print(solve(S, k))

短くしてみた (135 byte).

S,K=input(),int(input())
print((K>(S=='?'))*sum([1,2*(b!=a)][b<'?'>a]for a,b in zip(S,S[-1]+S))*K*pow(2,S.count('?')*K-2,P:=10**9+7)%P)

ZONeエナジー プログラミングコンテスト “HELLO SPACE”

だめだめのだめだめでした,残念.

A - UFO襲来

各添え字  i に対して愚直にチェックする.100 点だけど実質的に for 文が必要でびっくりした.

実装 (Python)

Python だと部分文字列の切り出しが簡単.

S = input()
T = "ZONe"
ans = 0
for i in range(len(S) - 4 + 1):
    ans += S[i:i + 4] == T
print(ans)

B - 友好の印

二分探索をやった (200 点とは).高さ  X を決め打つと,これが ok である必要十分条件は次の図のように考えられる.

f:id:suisen_kyopro:20210501223303p:plain

実装 (Python)

N, D, H = map(int, input().split())
p = []
for _ in range(N):
    d, h = map(int, input().split())
    p.append((d, h))
p.sort()
l = 0
r = H
for _ in range(50):
    X = (l + r) / 2
    ok = True
    for d, h in p:
        ok &= h <= (H - X) * d / D + X
    if ok:
        r = t
    else:
        l = t
print(f'{r:.10f}')

C - MAD TEAM

またまた二分探索.以下はスキルの種類数を  M として  X を達成できるかを考える. A _ i \geq X となる  i の集合を  S _ 1 とする.他も同様にして  S _ 2,\ldots, S _ M まで定める.すると,この問題は次のように言い換えられる.

問題:「 T=\{i,j,k\}\;(1\leq i\lt j\lt k\leq N) なる  T であって, S _ m\cap T\neq \emptyset\;(m=1.\ldots,M) を満たすものは存在するか?」

以降はこの判定問題を解くことを考える.まず,各  i=1,\ldots,N に対して  U _ i = \{m\mid i\in S _ i \} を計算する.ここで, 1\leq i\lt j\leq N を固定すると,問題は更に次のように言い換えられる.

問題:「 i,j が与えられる. U _ i \cup U _ j\cup U _ k = \{1,\ldots, M\} を満たす  k は存在するか?」

従って,すべての  S\subset\{1,\ldots,M\} に対して,

 \displaystyle
\mathrm{exists}[S]=\begin{cases}
\mathrm{True}&\text{if $\ \exists k\ $ s.t. $\ S\subset U _ k$}\\
\mathrm{False}&\text{otherwise}
\end{cases}

を予め計算しておけば  i,j を固定したときの問題に高速に答えることが出来る.

 U _ i の計算は  O(N M) \mathrm{exists} の計算は  O(2 ^ M M) で可能なので,結局全体  O( (N ^ 2 + NM + 2 ^ M M)\min(\max A, \max B, \max C, \max D, \max E) ) でこの問題を解くことが出来た.

実装 (C++)

// マクロなどは省略
int main() {
    int n;
    read(n);
    vec<vec<int>> params(n, vec<int>(5));
    rep(i, n) read(params[i]);
    int l = 0, r = 1000000001;
    while (r - l > 1) {
        int x = (l + r) >> 1;
        vec<int> ss(n);
        vec<bool> ex(1 << 5, false);
        rep(i, n) {
            int s = 0;
            rep(j, 5) {
                int t = (params[i][j] >= x);
                s |= t << j;
            }
            ss[i] = s;
            ex[s] = true;
        }
        rrep(s, 1 << 5) {
            rep(j, 5) {
                if ((s >> j) & 1) continue;
                int t = s | (1 << j);
                ex[s] = ex[s] or ex[t];
            }
        }
        bool ok = false;
        rep(i, n) {
            repl(j, i + 1, n) {
                int s = ((1 << 5) - 1) ^ (ss[i] | ss[j]);
                if (ex[s]) {
                    ok = true;
                    break;
                }
            }
            if (ok) break;
        }
        if (ok) {
            l = x;
        } else {
            r = x;
        }
    }
    print(l);
    return 0;
}

D - 宇宙人からのメッセージ

反転クエリは今反転しているかどうかのフラグを持つことで対処する.そして,反転されている場合は文字列の先頭に,反転されていない場合は文字列の末尾に文字を追加すればよい.このような操作を高速に行えるデータ構造として Deque が存在する.ここで,Deque は先頭/末尾への要素の追加/削除がならし定数時間で行えるデータ構造である.

なお,本問では連続する文字を削除する必要がある.これは,追加する文字が,その時点の文字列における追加しようとしている側の端の文字と同じ文字であった場合に追加の代わりに削除を行うようにすればよい (日本語が下手ですみません) .

実装 (C++)

// マクロなどは省略
int main() {
    string s;
    read(s);
    deque<char> t;
    bool rev = false;
    for (char c : s) {
        if (c == 'R') {
            rev = not rev;
        } else {
            if (rev) {
                if (t.size() and t.front() == c) {
                    t.pop_front();
                } else {
                    t.push_front(c);
                }
            } else {
                if (t.size() and t.back() == c) {
                    t.pop_back();
                } else {
                    t.push_back(c);
                }
            }
        }
    }
    while (t.size()) {
        if (rev) {
            cout << t.back();
            t.pop_back();
        } else {
            cout << t.front();
            t.pop_front();
        }
    }
    cout << '\n';
    return 0;
}

E - 潜入

ただの最短経路問題に見えるが, (r,\ c) から  (r-i,\ c) への移動についてすべての辺を張ってしまうと辺の数が最大で  R ^ 2 C 本程度になってしまい,TL に間に合わない可能性が高い (間に合ったという人もいるみたいですが).

 r を減らす方向の移動は,初めに  1 のコストを払った後, r 1 減るたびにコストを  1 だけ払うと考えられる.これは,例えば「入場料」を支払って電車に乗るようなものである.この考え方に従って,今入場しているかどうかで頂点の状態を場合分けをすると,次の図のように辺を張ることが出来る.

f:id:suisen_kyopro:20210501235035p:plain

実装 (Java)

// ライブラリなどは省略
public static void solve(ExtendedScanner sc, FastPrintStream pw) {
    int n = sc.nextInt();
    int m = sc.nextInt();
    Digraph<SimpleEdge> g = new Digraph<>(2 * n * m);
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m - 1; ++j) {
            int a = sc.nextInt();
            g.addEdge(new SimpleEdge(i * m + j, i * m + j + 1, a));
            g.addEdge(new SimpleEdge(i * m + j + 1, i * m + j, a));
        }
    }
    for (int i = 0; i < n - 1; ++i) {
        for (int j = 0; j < m; ++j) {
            int b = sc.nextInt();
            g.addEdge(new SimpleEdge(i * m + j, (i + 1) * m + j, b));
            g.addEdge(new SimpleEdge(n * m + (i + 1) * m + j, n * m + i * m + j, 1));
        }
    }
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j) {
            g.addEdge(new SimpleEdge(i * m + j, n * m + i * m + j, 1));
            g.addEdge(new SimpleEdge(n * m + i * m + j, i * m + j, 0));
        }
    }
    Dijkstra<SimpleEdge> dij = new Dijkstra<>(g, 0);
    pw.println(dij.distance(n * m - 1));
}

F - 出会いと別れ

 S:=\{0,\ldots,N-1\}\backslash \{A_1,\ldots,A_{M}\}=\{B_1,\ldots,B_K\}\;(K=N-M) と定義する.

頂点  0 から頂点  v へとたどり着けるための必要十分条件は,ある整数  0\leq p\leq K と長さ  p の整数列  1\leq i _ 1\lt i _ 2\lt\cdots\lt i _ p\leq K であって, v=B _ {i _ 1}\oplus\cdots\oplus B _ {i _ p} を満たすものが存在することである.ここで,\oplus は bit 毎の排他的論理和を取る演算である.つまり, B _ 1,\ldots,B _ K の線形結合によってすべての  0\leq v\lt N を表示できるか?という判定問題を考えればよいことになる.あとは, \log_2N 個の基底を取り出し (取り出せなかったら -1 を出力して終了する), N 頂点  N\log _ 2 N 辺のグラフで適当な全域木を取ってやればよい.

実装 (Java)

public static void solve(ExtendedScanner sc, FastPrintStream pw) {
    int n = sc.nextInt();
    int m = sc.nextInt();
    boolean[] ok = new boolean[n];
    Arrays.fill(ok, true);
    for (int i = 0; i < m; ++i) {
        int a = sc.nextInt();
        ok[a] = false;
    }
    IntArrayList bases = new IntArrayList();
    IntArrayList xs = new IntArrayList();
    for (int i = 0; i < n; ++i) {
        if (!ok[i]) continue;
        int e = i;
        PrimitiveIterator.OfInt iter = bases.iterator();
        while (iter.hasNext()) e = Math.min(e, e ^ iter.nextInt());
        if (e != 0) {
            bases.add(e);
            xs.add(i);
        }
    }
    if (bases.size() != Integer.numberOfTrailingZeros(n)) {
        pw.println(-1);
        return;
    }
    boolean[] vis = new boolean[n];
    IntDeque dq = new IntDeque(n);
    vis[0] = true;
    dq.addLast(0);
    while (dq.size() > 0) {
        int u = dq.removeFirst();
        for (var iter = xs.iterator(); iter.hasNext();) {
            int x = iter.nextInt();
            int v = u ^ x;
            if (!vis[v]) {
                vis[v] = true;
                dq.addLast(v);
                pw.print(u).print(' ').println(v);
            }
        }
    }
}

第六回 アルゴリズム実技検定 (バーチャル参加)

177:09 で完走.自己最速なので嬉しい.

f:id:suisen_kyopro:20210425084931p:plain

A - POSTal Code

正規表現が楽?

実装 (Python)

import re
if re.compile(r'[0-9]{3}-[0-9]{4}').match(input()):
    print('Yes')
else:
    print('No')

B - PASTal Code

o の位置を探すのが楽だと思う.

実装 (Python)

X = input().find('o') // 4 + 1
print(X if X else 'none')

Pythonfind 関数は見つからなかった場合に  -1 が返るので,この場合  X=0 となる.そのほかの場合は必ず  X\gt 0 になっている.

C - 携帯電話の購入

 A _ i Bset で持って,intersection の要素数 Q 以上のものを数えるとよい.

実装 (Python)

N, M = map(int, input().split())
A = []
for i in range(N):
   K, *X = map(int, input().split())
   A.append(set(X))
P, Q = map(int, input().split())
B = set(map(int, input().split()))
print(sum(Q <= len(B & X) for X in A))

Python では & 演算子set 同士の intersection を取ることが出来て楽に書ける.

D - K項足し算

ここから計算量を意識しないといけなくなる.愚直にやると  \Theta(N ^ 2) となって厳しいので,累積和で高速化する.

つまり,

 \displaystyle
S_i=\sum_{j=0}^{i-1}A_j

とおくと,元の数列の区間 A _ l + \cdots + A _ {r - 1}

 \displaystyle
A_l+\cdots+A_{r-1}=\left(\sum_{j=0}^{r-1}A_j\right)-\left(\sum_{j=0}^{l-1}A_j\right)=S_r-S_l

より  O(1) で計算できる.従って,全体  \Theta(N) でこの問題を解くことが出来る.

実装 (Python)

N, K, *A = map(int, open(0).read().split())
S = [0] + A
for i in range(N):
    S[i + 1] += S[i]
for i in range(N - K + 1):
    print(S[i + K] - S[i])

E - 前から3番目

実装が辛いことに...

deque を知っていますか?という問題.deque は先頭および末尾要素の追加と削除がならし  O(1) でできるデータ構造.

実装 (Python)

from collections import deque
N = int(input())
A = deque()
S = input()
for i, c in enumerate(S):
    if 'L' == c:
        A.appendleft(i + 1)
    if 'R' == c:
        A.append(i + 1)
    if 'A' == c:
        if len(A) <= 0:
            print("ERROR")
        else:
            print(A.popleft())
    if 'B' == c:
        if len(A) <= 1:
            print("ERROR")
        else:
            v = A.popleft()
            print(A.popleft())
            A.appendleft(v)
    if 'C' == c:
        if len(A) <= 2:
            print("ERROR")
        else:
            u = A.popleft()
            v = A.popleft()
            print(A.popleft())
            A.appendleft(v)
            A.appendleft(u)
    if 'D' == c:
        if len(A) <= 0:
            print("ERROR")
        else:
            print(A.pop())
    if 'E' == c:
        if len(A) <= 1:
            print("ERROR")
        else:
            v = A.pop()
            print(A.pop())
            A.append(v)
    if 'F' == c:
        if len(A) <= 2:
            print("ERROR")
        else:
            u = A.pop()
            v = A.pop()
            print(A.pop())
            A.append(v)
            A.append(u)

F - 安全装置

場合分けを正しくするのが少し大変.現時点で何秒間連続で  L 以上の負荷がかかっているかを表す変数  C を持っておく.はじめ, C = 0 である.

(可能な場合は) 各タスクの処理を開始する時点で必ず  C\lt T を満たすという仮定を置いて考えると良い.場合分けで破滅しそうな時はこういう仮定をちゃんと意識するのが結構大事だと思う (1WA したので偉そうに言えない).

  1.  B _ i\lt L の場合:  C = 0 にリセットされ,答えに  A _ i を加算するだけでよい.
  2.  B _ i \geq L かつ  C + A _ i \lt T の場合:  C\leftarrow C + A _ i とし,答えに  A _ i を加算するだけでよい.
  3.  B _ i \geq L かつ  C + A _ i = T の場合: ちょうどタスク終了時に安全装置が起動するので,タスクのやり直しは発生しない.従って,答えに  A _ i + X を加算して  C = 0 にリセットすればよい.
  4.  B _ i \geq L かつ  C + A _ i \gt T の場合: タスクを  T - C\lt A _ i 秒処理した時点で安全装置が起動するので,タスクのやり直しが発生する.従って,答えに  T-C+A _ i +X を加算し, C\leftarrow A _ i とすればよい.ただし, A _ i \gt T であればタスクのやり直しが永遠に失敗するためここで forever を出力して終了する必要がある.また, A _ i = T であればタスク終了時に再び安全装置が起動するので答えに更に  X を加算し, C = 0 にリセットする必要がある.

ここまで場合分けが出来れば,あとはそのまま実装すれば良い.

実装 (Python)

N, L, T, X = map(int, input().split())
C = 0
ans = 0
for i in range(N):
    A, B = map(int, input().split())
    if B < L:
        C = 0
        ans += A
    elif C + A < T:
        ans += A
        C += A
    elif C + A == T:
        ans += A + X
        C = 0
    else:
        ans += T - C + X + A
        C = A
        if C > T:
            print('forever')
            exit(0)
        elif C == T:
            ans += X
            C = 0
print(ans)

G - 一日一歩

突然難易度が上がった気がする.

動かないことが許されるので,一度到達可能になったらあとで到達不可能になることはない.そこで,頂点  1 から徐々に陣地を広げていくイメージで少しずつ辺を追加していくと考える.

「まだ使っていない辺のうち,到達可能な頂点の集合の要素を少なくとも一方の端点として持つもの」の集合を優先度付きキューに入れて管理すると,追加する  X _ i 以下の辺を効率的に探すことができる.優先度付きキューには合計で高々  M 回 push され,合計で高々  M しか pop されないので,優先度付きキューの操作にかかる計算量は全体  O(M\log M) と評価される.

あとは,頂点  1 と連結な頂点の個数が高速に求まればよく,これは例えば Union Find 木によりクエリ当たり  O(\alpha(N)) ( \alphaアッカーマン関数逆関数) で計算できる.

以上より,全体  O(N+Q+M(\log M + \alpha(N)) でこの問題を解くことが出来た.

実装 (Java)

一日に 2 歩歩くことをうまく禁止する必要があるので注意.

自作のクラスが色々出てきますが,詳細は省略します (私の Java コードは大体 20000 Byte とかで,全部貼ってしまうとえらいことになるのでご容赦を...).

int n = sc.nextInt();
int m = sc.nextInt();
int q = sc.nextInt();
// 無向グラフ
Graph<SimpleEdge> g = new Graph<>(n);
for (int i = 0; i < m; ++i) {
    int u = sc.nextInt() - 1;
    int v = sc.nextInt() - 1;
    int w = sc.nextInt();
    g.addEdge(new SimpleEdge(u, v, w));
}
// Union Find 木
DisjointSetUnion dsu = new DisjointSetUnion(n);
PriorityQueue<SimpleEdge> pq = new PriorityQueue<>();
for (var e : g.getEdges(0)) pq.add(e);
// i 日目に追加する辺を一時的に保存するための List
ArrayList<SimpleEdge> l = new ArrayList<>();
for (int i = 0; i < q; ++i) {
    int x = sc.nextInt();
    while (pq.size() > 0) {
        var e = pq.peek();
        if (e.cost > x) break;
        pq.poll();
        int u = e.from;
        int v = e.to;
        if (dsu.same(0, v)) {
            int tmp = u; u = v; v = tmp;
        }
        if (dsu.same(0, v)) continue;
        dsu.merge(0, v);
        l.addAll(g.getEdges(v)); // ここで pq に追加すると 1 日に 2 歩歩くことを許してしまうので注意
    }
    pw.println(dsu.size(0));
    pq.addAll(l);
    l.clear();
}

H - カンカンマート

缶切りの条件がややこしいので,買う缶切りの個数を全探索できないかを考える.

以下, S _ 0 = \{i\mid T _ i = 0\},\ S _ 1 = \{i\mid T _ i = 1\} とし,列  A,\ B A = (P _ i) _ {i \in S _ 0},\ B = (P _ i) _ {i \in S _ 1} で定義する.

缶切りの個数  X を固定したとき, B からは最大  KX 個の要素を選ぶことができる.一方, A からはいくつでも要素を選んでよい.選択可能な要素の (多重) 集合が決まっている場合,明らかに小さい順に貪欲に選んでいくのが最適である.従って,予め  B を昇順にソートしておくと, B _ 1, \ldots, B _ { \min (KX, |B| ) } が選択可能であるとしてよい.

以上より,多重集合および下位  M 個の和を効率的に管理できるデータ構造があればこの問題を解くことが出来ることが分かった.これは, 2 本の優先度付きキューを持ち,片方に下位  M 要素を,もう片方にその他の要素を入れることにすれば実現できる.

実装 (Java)

LongPrioritySum が多重集合の下位  M 要素の和を管理するクラスである.

int n = sc.nextInt();
int m = sc.nextInt();
int k = sc.nextInt();
long q = sc.nextInt();
LongPrioritySum ps = new LongPrioritySum(m, false);
IntArrayList t = new IntArrayList();
for (int i = 0; i < n; ++i) {
    int p = sc.nextInt();
    if (sc.nextInt() == 0) {
        ps.add(p);
    } else {
        t.add(p);
    }
}
t.sort();
int cnt = ps.size();
long ans = cnt >= m ? ps.sum() : Long.MAX_VALUE;
int max = (n + k - 1) / k;

for (int x = 1; x <= max; ++x) {
    int l = k * (x - 1);
    int r = Math.min(l + k, t.size());
    for (int j = l; j < r; ++j) {
        ps.add(t.get(j));
    }
    if (ps.size() >= m) {
        ans = Math.min(ans, q * x + ps.sum());
    }
}
pw.println(ans);

今回は LongPrioritySum が本質的なので,実装を示しておく (LongPrioeityQueue 標準ライブラリには存在しないので注意).

public final class LongPrioritySum {
    private final boolean isDescending;
    private final LongPriorityQueue topk;
    private final LongPriorityQueue other;
    private int k;
    private long sumK = 0;
    public LongPrioritySum(int k, boolean descending) {
        this.isDescending = descending;
        this.topk = new LongPriorityQueue(!descending);
        this.other = new LongPriorityQueue(descending);
        this.k = k;
    }
    public LongPrioritySum(int k) {
        this(k, false);
    }
    public int  getK()      {return k;}
    public void setK(int k) {this.k = k;}
    public void incrK()     {this.k++;}
    public void decrK()     {this.k--;}
    public void addK(int v) {this.k += v;}
    public void subK(int v) {this.k -= v;}
    public void add(long v) {
        if (topk.size() == 0 || topk.getFirst() < v ^ isDescending) {
            other.add(v);
        } else {
            topk.add(v);
            sumK += v;
        }
    }
    public long sum() {
        int size = topk.size();
        if (size > k) {
            int d = size - k;
            while (d --> 0) {
                long v = topk.removeFirst();
                sumK -= v;
                other.add(v);
            }
        } else if (size < k){
            int d = k - size;
            while (d --> 0 && other.size() > 0) {
                long v = other.removeFirst();
                sumK += v;
                topk.add(v);
            }
        }
        return sumK;
    }
    public int size() {
        return topk.size() + other.size();
    }
}

I - 魚釣り

典型的な DP の問題.

 \displaystyle
\mathrm{dp0}[c][i][j]=\text{区画 $(i,\ j)$ に至るまでに $c$ 回釣っており,かつ区画 $(i,\ j)$ で魚を釣っていない場合の最大値}
 \displaystyle
\mathrm{dp1}[c][i][j]=\text{区画 $(i,\ j)$ に至るまでに $c$ 回釣っており,かつ区画 $(i,\ j)$ で魚を釣っている場合の最大値}

と定義する. 2 つのテーブルに分けたのは,同じ区画で高々  1 回しか釣れないという制約を考慮するためである.遷移は実装を見るのが一番早いと思う.

実装 (Java)

int h = sc.nextInt();
int w = sc.nextInt();
long[][] a = sc.longs(h, w);
long[][][] dp0 = new long[h + w][h][w];
long[][][] dp1 = new long[h + w][h][w];
for (int c = 1; c < h + w; ++c) {
    for (int i = 0; i < h; ++i) for (int j = 0; j < w; ++j) {
        // 区画 (i, j) で釣る場合の遷移
        dp1[c][i][j] = dp0[c - 1][i][j] + a[i][j];
        if (i > 0) { // 区画 (i - 1, j) から移動してくる場合の遷移
            dp0[c][i][j] = Math.max(dp0[c][i][j], dp0[c][i - 1][j]);
            dp0[c][i][j] = Math.max(dp0[c][i][j], dp1[c][i - 1][j]);
        }
        if (j > 0) { // 区画 (i, j - 1) から移動してくる場合の遷移
            dp0[c][i][j] = Math.max(dp0[c][i][j], dp0[c][i][j - 1]);
            dp0[c][i][j] = Math.max(dp0[c][i][j], dp1[c][i][j - 1]);
        }
    }
    // k = c の時の答えは max(dp0[c][h - 1][w - 1], dp1[c][h - 1][w - 1])
    pw.println(Math.max(dp0[c][h - 1][w - 1], dp1[c][h - 1][w - 1]));
}

J - ポイントとコストと

まず, p によらず  (Y _ i - C) ^ 2 の和は一定であるので別で計算しておく.すると,あとは

 \displaystyle
\min_p\sum_{i=1}^N (X_i-p)^2

を求める問題となる.これはかなり有名な問題で, \displaystyle\underset{p}{\mathrm{argmin}}\sum _ {i=1} ^ N (X _ i - p) ^ 2 = \dfrac{\sum _ {i = 1} ^ N X _ i}{N} という結果が知られている.従って,あとは決めた  p を用いて計算をすればよい.

ただ,この結果を知らなくても,展開すれば  p 2 次式になるので平方完成すれば最適な  p の値を求めることが出来る.

 \displaystyle
\begin{align}
\sum_{i=1}^N (X_i-p)^2
&=Np^2-\left(\sum_{i=1}^NX_i\right)p+\sum_{i=1}^NX_i^2\\
&=N\left(p-\dfrac{\sum_{i=1}^NX_i}{N}\right)^2+C\quad (C\equiv \mathrm{const})
\end{align}

あるいは最小値を持つことを仮定し, p微分して微分係数 0 となる位置を計算してもよい.

実装 (Python)

N, C = map(int, input().split())
ans = sx = 0
xs = []
for _ in range(N):
    x, y = map(int, input().split())
    sx += x
    ans += (y - C) ** 2
    xs.append(x)
mx = sx / N
for x in xs:
    ans += (x - mx) ** 2
print(ans)

K - 共通クーポン

 V _ i = U _ i - P _ i とする.次を求めたい.

 \displaystyle
\max _ {S\subset\{1,\ldots,N\}} \sum_{i\in S} V_i+20\times\left\lfloor\dfrac{\sum_{i\in S}P_i}{100}\right\rfloor

厄介なのは切り捨て処理なので, \left(\sum_{i\in S}P_i\right)\bmod 100 を状態に持って次のように DP テーブルを定義する.

 \displaystyle
\mathrm{dp}[i][m]=\max_{S\subset\{1,\ldots,i\}}\left\{V_i+20\times\left\lfloor\dfrac{\sum_{i\in S}P_i}{100}\right\rfloor\;\middle|\; \sum_{i\in S}P_i\equiv m\;\left(\mathrm{mod}\;100\right)\right\}

これも遷移は実装を見た方が早いと思う.

実装 (Python)

N = int(input())
items = []
for _ in range(N):
    P, U = map(int, input().split())
    V = U - P
    items.append((V, P))
inf = 1000000000
dp = [[-inf] * 100 for _ in range(N + 1)]
dp[0][0] = 0
for i, (v, p) in enumerate(items):
    for j in range(100):
        q, r = divmod(j + p, 100)
        # 商品を買う場合の遷移
        dp[i + 1][r] = max(dp[i + 1][r], dp[i][j] + v + q * 20)
        # 商品を買わない場合の遷移
        dp[i + 1][j] = max(dp[i + 1][j], dp[i][j])
print(max(dp[N]))

L - 都市計画

使う環状交差点を  2 ^ M 通り全探索し, N 個のタワーと使うと決めた環状交差点を頂点とし,辺を適切に張ったグラフの最小全域木を求めればよい.辺の重みは  O(1) で計算できるので,全体  O(2 ^ M (N + M) ^ 2)

と言うのは簡単だけど,位置関係によって距離の計算が変わるので面倒くさい...

manabitimes.jp

ここを見てください (丸投げ)

実装 (Java)

static long d2(long dx, long dy) {
    return dx * dx + dy * dy;
}

static double d(long dx, long dy) {
    return Math.sqrt(d2(dx, dy));
}

static class Edge implements Comparable<Edge> {
    int u, v;
    double c;
    Edge(int u, int v, double c) {
        this.u = u; this.v = v; this.c = c;
    }
    @Override
    public int compareTo(Edge o) {
        return Double.compare(c, o.c);
    }
}

public static void solve(ExtendedScanner sc, FastPrintStream pw) {
    int n = sc.nextInt();
    int m = sc.nextInt();
    int[] tx = new int[n];
    int[] ty = new int[n];
    for (int i = 0; i < n; ++i) {
        tx[i] = sc.nextInt();
        ty[i] = sc.nextInt();
    }
    int[] cx = new int[m];
    int[] cy = new int[m];
    int[] cr = new int[m];
    for (int j = 0; j < m; ++j) {
        cx[j] = sc.nextInt();
        cy[j] = sc.nextInt();
        cr[j] = sc.nextInt();
    }
    int[][] bits = BitUtil.bits(m);
    double ans = Double.MAX_VALUE;
    for (int s = 0; s < 1 << m; ++s) {
        PriorityQueue<Edge> pq = new PriorityQueue<>();
        int l = bits[s].length;
        for (int i = 0; i < n; ++i) {
            for (int j = i + 1; j < n; ++j) {
                pq.add(new Edge(i, j, d(tx[i] - tx[j], ty[i] - ty[j])));
            }
            for (int j = 0; j < l; ++j) {
                int p = bits[s][j];
                pq.add(new Edge(i, n + j, Math.abs(d(tx[i] - cx[p], ty[i] - cy[p]) - cr[p])));
            }
        }
        for (int i = 0; i < l; ++i) {
            for (int j = i + 1; j < l; ++j) {
                int p = bits[s][i], q = bits[s][j];
                int x1 = cx[p], y1 = cy[p], r1 = cr[p];
                int x2 = cx[q], y2 = cy[q], r2 = cr[q];
                if (r1 < r2) {
                    int tmp;
                    tmp = x1; x1 = x2; x2 = tmp;
                    tmp = y1; y1 = y2; y2 = tmp;
                    tmp = r1; r1 = r2; r2 = tmp;
                }
                double dc = d(x1 - x2, y1 - y2);
                double c;
                if (dc >= r1 + r2) {
                    c = dc - (r1 + r2);
                } else if (dc >= r1 - r2) {
                    c = 0;
                } else {
                    c = r1 - r2 - dc;
                }
                pq.add(new Edge(n + i, n + j, c));
            }
        }
        DisjointSetUnion dsu = new DisjointSetUnion(n + l);
        double sum = 0;
        while (dsu.size(0) != n + l) {
            var e = pq.poll();
            int u = e.u, v = e.v;
            if (dsu.same(u, v)) continue;
            dsu.merge(u, v);
            sum += e.c;
        }
        ans = Math.min(ans, sum);
    }
    pw.println(ans);
}

M - 等しい数

今回の問題ではこれが一番難しかった.

まず, A に現れる値 (クエリも含める) を座圧して  0,\ldots,M-1 と番号を振りなおす.各クエリの時点における  A _ i の値を必要に応じて知れるようにするため,区間代入を載せた双対セグ木 (/遅延セグ木) を用意しておく.そして,各  v に対して集合  S _ v = \{i \mid A _ i = v\} を管理することを考える. S _ v はそのまま set などで持つのではなく,区間の集合として持つ (いわゆる区間を set で管理するテク) ことにすると,愚直にやっても実はいい感じに計算量が償却される.

具体的には,次のような操作を行う.図中の left が指す値を取得するために双対セグ木が必要となっている.

f:id:suisen_kyopro:20210425082542p:plain
図 1. クエリ処理前
f:id:suisen_kyopro:20210425082611p:plain
図 2-1. クエリ  (L,\ R,\ x)=(1,\ 11,\ 2) の処理
f:id:suisen_kyopro:20210425082716p:plain
図 2-2. クエリ  (L,\ R,\ x)=(1,\ 11,\ 2) の処理
f:id:suisen_kyopro:20210425082730p:plain
図 2-3. クエリ  (L,\ R,\ x)=(1,\ 11,\ 2) の処理
f:id:suisen_kyopro:20210425082745p:plain
図 2-4. クエリ  (L,\ R,\ x)=(1,\ 11,\ 2) の処理
f:id:suisen_kyopro:20210425082756p:plain
図 2-5. クエリ  (L,\ R,\ x)=(1,\ 11,\ 2) の処理
f:id:suisen_kyopro:20210425082517p:plain
図 3. クエリ処理後

区間が set に挿入される回数を考える.まず,初期状態に一致させるために  N 回挿入される.その後に挿入される回数は,クエリの個数  Q に一致する.また,区間が部分的に削除される回数 (つまり,図 2-1, 2-5 のような状態になる回数) は,各クエリに対して高々  2 回しか起こらない.

区間が set から削除される回数は挿入される回数と部分的に削除される回数の和で抑えられるので,結局 set の操作回数は全体で  O(N+Q) 回である.従って,set パートの計算量は全体で  O( (N+Q)\log(N+Q) ) である.

また,双対セグ木に対するクエリ数は区間代入の回数と区間削除の回数の和であり,全体  O((N+Q)\log N) である.

以上より,全体  O( (N+Q)\log(N+Q) ) でこの問題が解けた.

実装 (Java)

int n = sc.nextInt();
int[] a = sc.ints(n);
int q = sc.nextInt();
int[] vals = Arrays.copyOf(a, n + q);
IntPair3[] qs = new IntPair3[q];
for (int i = 0; i < q; ++i) {
    int l = sc.nextInt() - 1;
    int r = sc.nextInt();
    int v = sc.nextInt();
    qs[i] = new IntPair3(l, r, v);
    vals[n + i] = v;
}
IntCompress cmp = new IntCompress(vals, false);
for (int i = 0; i < n; ++i) {
    a[i] = cmp.compress(a[i]);
}
for (int i = 0; i < q; ++i) {
    qs[i].trd = cmp.compress(qs[i].trd);
}
int m = cmp.compressedSize();
long[] num = new long[m];
IntDualSegmentTree seg = new IntDualSegmentTree(
    a, (f, x) -> f, (f, g) -> f, -1
);
IntRangeSet[] rs = new IntRangeSet[m];
for (int i = 0; i < m; ++i) rs[i] = new IntRangeSet();
for (int i = 0; i < n; ++i) {
    rs[a[i]].add(i);
    ++num[a[i]];
}
long ans = 0;
for (long c : num) {
    ans += c * (c - 1) / 2;
}
for (var query : qs) {
    int l = query.fst, r = query.snd;
    int v = query.trd;
    int cur = l;
    while (cur < r) {
        int u = seg.get(cur);
        int nxt = rs[u].getInclusingRange(cur).getRight();
        long rem;
        if (nxt < r) {
            rem = rs[u].removeAll(cur, nxt);
            cur = nxt + 1;
        } else {
            rem = rs[u].removeAll(cur, r - 1);
            cur = r;
        }
        long x = num[u], y = num[u] - rem;
        ans -= x * (x - 1) / 2 - y * (y - 1) / 2;
        num[u] -= rem;
    }
    long add = rs[v].addAll(l, r - 1);
    long x = num[v] + add;
    long y = num[v];
    ans += x * (x - 1) / 2 - y * (y - 1) / 2;
    num[v] = x;
    seg.apply(l, r, v);
    pw.println(ans);
}

N - 活動

活動の順序さえ決まれば簡単な DP で求まるので,順序を一意に定められないかを考える.体力が  X のときに,活動  i,\ j を行うとする.このとき,

  •  i,\ j の順で活動した場合: 得られる得点は  X\cdot a _ i + (X - b _ i)\cdot a _ j
  •  j,\ i の順で活動した場合: 得られる得点は  X\cdot a _ j + (X - b _ j)\cdot a _ i

である. i,\ j の順で活動する方がよいのは  X\cdot a _ i + (X - b _ i)\cdot a _ j\gt X\cdot a _ j + (X - b _ j)\cdot a _ i が成り立つ場合で,これを整理すると次のようになる.

 \displaystyle
\begin{align}
& X\cdot a _ i + (X - b _ i)\cdot a _ j\gt X\cdot a _ j + (X - b _ j)\cdot a _ i\\
\Leftrightarrow
& -b_ia_j\gt -b_ja_i\\
\Leftrightarrow
& \dfrac{a_i}{b_i}\gt\dfrac{a_j}{b_j}\quad(\because b_i,\ b_j\gt 0)\\
\end{align}

従って, a _ i/b _ i の大きい順に使うかどうかを決めていけばよいことが分かった.あとは, a _ i/b _ i の降順に活動をソートし,

 \displaystyle
\mathrm{dp}[i][j]=\text{活動$i$までで残りの体力が$j$であるような選び方における得点の最大値}

として DP をすればよい.

実装 (Java)

体力が負になることを許容する点に注意しながら実装する必要がある.負の index に対応するために添え字に下駄を履かせることもあるが,今回は直接終了状態に遷移するのが最適なのでその必要はない.

int n = sc.nextInt();
int h = sc.nextInt();
IntPair[] p = new IntPair[n];
for (int i = 0; i < n; ++i) {
    int a = sc.nextInt();
    int b = sc.nextInt();
    p[i] = new IntPair(a, b);
}
Arrays.sort(p, (p1, p2) -> {
    long a1 = p1.fst, b1 = p1.snd;
    long a2 = p2.fst, b2 = p2.snd;
    return Long.signum(a2 * b1 - a1 * b2);
});
long ans = 0;
long[][] dp = new long[n + 1][h + 1];
dp[0][h] = 0;
for (int i = 0; i < n; ++i) {
    long gain = p[i].fst;
    int damage = p[i].snd;
    for (int j = 0; j <= h; ++j) {
        dp[i + 1][j] = Math.max(dp[i + 1][j], dp[i][j]);
        if (j >= damage) {
            dp[i + 1][j - damage] = Math.max(dp[i + 1][j - damage], dp[i][j] + j * gain);
        } else { // 一旦体力が負になるとあとは活動しないのが最適なので,直接答えに遷移する.
            ans = Math.max(ans, dp[i][j] + j * gain);
        }
    }
}
pw.println(Math.max(ans, LongArrays.max(dp[n])));

O - 最短距離クエリ

制約が明らかに変なのでこれを活かす.入力を全域木とその他の  K:=M-N+1\leq 11 本の辺に分けておく.

木の場合,ダブリングなどにより LCA O(\log N) で求めることができるので,各頂点の深さを持つことで任意の  u,\ v に対し距離  \mathrm{dist}(u,\ v) O(\log N) で計算することが出来る.

各クエリ  u,\ v に対し, u,\ v および  K 本の余剰辺の端点以外の頂点を潰して考える.このとき高々  2K+2\leq 24 頂点のグラフになるので,LCA による距離の計算によって適切に辺の重みを設定すれば,Dijkstra 法などによりクエリを  O(K ^ 2 \log N) で処理できる.ダブリングにより LCA を求める場合,全体の計算量は  O(N\log N + Q K ^ 2 \log N) となる.

実装 (Java)

int n = sc.nextInt();
int m = sc.nextInt();
DisjointSetUnion dsu = new DisjointSetUnion(n);
int k = m - (n - 1);
IntPair[] edges = new IntPair[k]; 
TreeBuilder tb = new TreeBuilder(n);
for (int i = 0, j = 0; i < m; ++i) {
    int u = sc.nextInt() - 1;
    int v = sc.nextInt() - 1;
    if (dsu.same(u, v)) {
        edges[j++] = new IntPair(u, v);
    } else {
        dsu.merge(u, v);
        tb.addEdge(u, v);
    }
}
TreeSet<Integer> need = new TreeSet<>();
int[] ids = new int[n];
Tree t = tb.build();
EulerTourLCA lca = new EulerTourLCA(t);
int q = sc.nextInt();
while (q --> 0) {
    int u = sc.nextInt() - 1;
    int v = sc.nextInt() - 1;
    need.add(u);
    need.add(v);
    for (var e : edges) {
        need.add(e.fst);
        need.add(e.snd);
    }
    int c = need.size();
    int[] vs = new int[c];
    for (int i = 0; i < c; ++i) {
        vs[i] = need.pollFirst();
        ids[vs[i]] = i; 
    }
    long[][] g = new long[c][c];
    for (long[] r : g) Arrays.fill(r, Dijkstra2.UNREACHABLE);
    for (int i = 0; i < c; ++i) {
        g[i][i] = 0;
        for (int j = i + 1; j < c; ++j) {
            int x = vs[i], y = vs[j];
            g[i][j] = g[j][i] = lca.dist(x, y);
        }
    }
    for (var e : edges) {
        int x = e.fst, y = e.snd;
        int i = ids[x], j = ids[y];
        g[i][j] = g[j][i] = 1;
    }
    long ans = new Dijkstra2(g, ids[u]).getDistance(ids[v]);
    pw.println(ans);
}

AtCoder Beginner Contest 199(Sponsored by Panasonic)

(D 以外全部を解くのにかかった時間) < (D を解くのにかかった時間) になってしまった...

A - Square Inequality

YNeos を書こうとして頭と手が止まった.普通に書いた方がはやいね...

A, B, C = map(int, input().split())
print("YNeos"[A ** 2 + B ** 2 >= C ** 2::2])

B - Intersection

つまりは \displaystyle \max _ {1\leq i\leq N} A _ i \leq x \leq \min _ {1\leq j\leq N} B _ j を満たす  x を数える問題. \max A \gt \min B のケースに注意.

N = int(input())
L = max(map(int, input().split()))
R = min(map(int, input().split()))
print(max(0, R - L + 1))

C - IPFL

type 2 のクエリで実際に swap する代わりに,swap されているかどうかのフラグを持つ.もし swap されていれば,添え字  i (0-indexed) は次のように読み替えられる.

  •  i\lt N のとき:  i+N
  •  i\geq N のとき:  i-N

この変換法則を用いて type 1 のクエリを処理する.フラグが true の場合,最後の出力時に swap し忘れると間違えるので注意.

下に示すコードでは rev が swap のフラグとなっている.

N = int(input())
S = [c for c in input()]
Q = int(input())
rev = False
for _ in range(Q):
    t, a, b = map(int, input().split())
    a -= 1
    b -= 1
    if t == 1:
        if rev:
            if a < N:
                a += N
            else:
                a -= N
            if b < N:
                b += N
            else:
                b -= N
        S[a], S[b] = S[b], S[a]
    else:
        rev = not rev
if rev:
    S[:N], S[N:] = S[N:], S[:N]
print(''.join(S))

D - RGB Coloring 2

やばいやつ.

そのまま全探索をすると  O(M 3 ^ N) で厳しいので,少し工夫して探索する必要がある.ある頂点  u の色が決まっている場合, u に隣接する頂点  v の色は  u の色以外の  2 通りしか存在しない.

従って,各連結成分に対して始点  u を取って色を決め, u からの距離が近い順であったり,あるいは  u を根とする (連結成分に対する) 全域木をとって DFS/BFS の行きがけ順であったり,適切な順番に色を決めることにするとよい.そうすると,サイズ  C の連結成分に対する答えは  O((C + M) 2 ^ C) 程度で求めることが出来る.

また,各連結成分は独立に数えられるので,各連結成分に対する答えを求めて最後にかけ合わせればよい.

以下は Java による AC コード.DisjointSetUnion クラスや Graph などの実装は省略している.

public static void solve(ExtendedScanner sc, FastPrintStream pw) {
    int n = sc.nextInt();
    int m = sc.nextInt();
    int[] us = new int[m];
    int[] vs = new int[m];
    DisjointSetUnion dsu = new DisjointSetUnion(n);
    for (int i = 0; i < m; ++i) {
        us[i] = sc.nextInt() - 1;
        vs[i] = sc.nextInt() - 1;
        dsu.merge(us[i], vs[i]);
    }
    int[][] groups = dsu.groups(); // 連結成分を取得
    long ans = 1;
    for (int[] group : groups) {
        boolean[] con = new boolean[n]; // 連結成分に属するかどうか
        int k = group.length;
        Graph<SimpleEdge> g = new Graph<>(k);
        for (int id = 0; id < k; ++id) {
            con[group[id]] = true;
        }
        for (int i = 0; i < m; ++i) {
            int u = us[i];
            int v = vs[i];
            if (con[u] && con[v]) {
                int ui = ids[u];
                int vi = ids[v];
                g.addEdge(new SimpleEdge(ui, vi));
            }
        }
        ans *= solve(g);
    }
    pw.println(ans);
}

static long solve(Graph<SimpleEdge> g) {
    int n = g.getV();
    IntArrayList pre = new IntArrayList(); // BFS で訪問した順番
    int[] par = new int[n]; // 全域木における親
    Arrays.fill(par, -1);
    boolean[] vis = new boolean[n];
    vis[0] = true;
    IntDeque dq = new IntDeque(n);
    dq.addLast(0);
    while (dq.size() > 0) {
        int u = dq.removeLast();
        pre.add(u);
        for (var e : g.getEdges(u)) {
            int v = e.to;
            if (!vis[v]) {
                vis[v] = true;
                par[v] = u;
                dq.addLast(v);
            }
        }
    }
    int[] color = new int[n];
    long ans = 0;
    Out:for (int i = 0; i < 1 << n - 1; ++i) {
        Arrays.fill(color, 1, n, -1);
        for (int j = 1; j < n; ++j) { // 色を BFS の訪問順に決定していく.
            int v = pre.get(j);
            color[v] = (color[par[v]] + 1 + ((i >> (v - 1)) & 1)) % 3;
        }
        for (int u = 0; u < n; ++u) {
            for (var e : g.getEdges(u)) {
                int v = e.to;
                if (color[u] == color[v]) {
                    continue Out;
                }
            }
        }
        ++ans;
    }
    return ans * 3;
}

E - Permutation

典型的な bit DP で,確実に D よりは易しいと思う.

 \displaystyle
\mathrm{dp}[S]=\text{制約に違反せず $\{a_1,\ldots,a_{|S|}\}=S$ となるような $a_1,\ldots,a_{|S|}$ の決め方}

と定義する.遷移に関しては,最後に追加する  x\in S を全探索して,実際に制約に違反していないかを確かめればよい.

public static void solve(ExtendedScanner sc, FastPrintStream pw) {
    int n = sc.nextInt();
    int m = sc.nextInt();
    int[] xs = new int[m];
    int[] ys = new int[m];
    int[] zs = new int[m];
    for (int i = 0; i < m; ++i) {
        xs[i] = sc.nextInt();
        ys[i] = sc.nextInt();
        zs[i] = sc.nextInt();
    }
    int[][] bits = BitUtil .bits(n); // bits[i] := i の 2 進数表記において立っている bit の位置を記録した配列
    long[] dp = new long[1 << n];
    dp[0] = 1;
    for (int s = 1; s < 1 << n; ++s) {
        int[] cnt = new int[n + 1]; // cnt[i] = i 未満の要素の個数
        for (int x : bits[s]) {
            ++cnt[x + 1];
        }
        for (int i = 0; i < n; ++i) {
            cnt[i + 1] += cnt[i];
        }
        int k = bits[s].length; // k=|S|
        for (int x : bits[s]) { // 追加する要素を全探索
            int t = s ^ (1 << x);
            boolean ok = true;
            for (int j = 0; j < m; ++j) {
                if (xs[j] >= k) { // Xj<k なる j に関する制約は満たされていることが保証されている
                    ok &= cnt[ys[j]] <= zs[j];
                }
            }
            if (ok) dp[s] += dp[t];
        }
    }
    pw.println(dp[(1 << n) - 1]);
}

F - Graph Smoothing

一回の選択で辺  (u,\ v) が選択されたときの頂点  i への寄与を考える.

  •  i \in \{u,v\} のとき:  \dfrac{A _ u + A _ v}{2M}
  •  i \not\in \{u,v\} のとき:  \dfrac{A _ i}{M}

従って, N\times N 行列  B の要素を全て  0 で初期化し,各辺  (u,\ v) に対して

  •  B _ {u, u} \leftarrow B _ {u, u} + \dfrac{1}{2M}
  •  B _ {u, v} \leftarrow B _ {u, v} + \dfrac{1}{2M}
  •  B _ {v, u} \leftarrow B _ {v, u} + \dfrac{1}{2M}
  •  B _ {v, v} \leftarrow B _ {v, v} + \dfrac{1}{2M}
  •  B _ {i, i} \leftarrow B _ {i, i} + \dfrac{1}{M}\quad (i\not\in\{u,v\})

で更新し,行列累乗で  B ^ K を求めてベクトル  A に掛けると答えのベクトルが求まる.

// mod 関連の演算をやるクラス
static final ModArithmetic ma = ModArithmetic1000000007.INSTANCE;
// mod M 上の行列演算をやるクラス
static final ModMatrixFactory mmf = new ModMatrixFactory(ma);

public static void solve(ExtendedScanner sc, FastPrintStream pw) {
    int n = sc.nextInt();
    int m = sc.nextInt();
    int k = sc.nextInt();
    long im = ma.inv(m);
    long i2 = ma.inv(2);
    long[] a = sc.longs(n);
    long[][] mat = new long[n][n];
    for (int i = 0; i < m; ++i) {
        int u = sc.nextInt() - 1;
        int v = sc.nextInt() - 1;
        for (int j = 0; j < n; ++j) {
            if (j == u) {
                mat[j][j] += ma.mul(i2, im);
                mat[j][v] += ma.mul(i2, im);
            } else if (j == v) {
                mat[j][j] += ma.mul(i2, im);
                mat[j][u] += ma.mul(i2, im);
            } else {
                mat[j][j] += im;
            }
        }
    }
    long[] x = mmf.create(mat).pow(k).mul(a);
    for (long e : x) {
        pw.println(e);
    }
}