Educational Codeforces Round 104 F - Ones

問題

リンク

正整数 $n$ が与えられる。$1$ や $11$、$111$ などの $1$ のみからなる数の加減算によって $n$ を作る (例: $24=11+11+1+1, 102=111-11+1+1$)。式中に現れる $1$ の個数の最小値を求めよ。

解法

$\displaystyle \underbrace {11 \cdots 1} _ { \text { $k$ 個 } } = \dfrac {10 ^ k - 1 } { 9 }$ を使う個数を $a_k$ とおく。ただし、正の数として使う場合は $+1$ 個、負の数として使う場合は $-1$ 個としてカウントする。このとき、求めるべき値は次のように表される。

$$\min\left\{\sum _ { k = 1 } ^ { \infty } k \cdot | a _ k | \;\middle|\; \sum _ { k = 1 } ^ { \infty } (10 ^ k - 1) a _ k = 9 n\right\}$$

$\displaystyle S\coloneqq \sum _ { k = 1 } ^ \infty a_k$ を固定したときの答えを $A _ S$ とおくと、

$$A_S = \min\left\{\sum _ { k = 1 } ^ { \infty } k \cdot | a _ k | \;\middle|\; \sum _ { k = 1 } ^ { \infty } 10 ^ k \cdot a _ k = 9 n + S\land \sum _ { k = 1 } ^ \infty a_k = S \right\}$$

であり、求める値は $\min_S A_S$ である。

$a_k\neq 0$ となる最大の $k$ は $\lceil \log _ {10} |9n+S| \rceil + 1$ 以下としてよく、また、$| a _ k | \lt 10$ の場合のみを考えても十分なので $S$ は $|S| \leq 9 \cdot (\lceil\log _ {10} |9n+S| \rceil + 1)$ を満たす範囲を調べれば十分である。

$A_S$ は桁 DP で (上から何桁決めたか ($=i$)、上から $i$ 桁だけ見たときの目標との差分、上から $i$ 桁の $a_k$ の和) を持てば求めることができる。差分は高々 1 としてよいので、DP テーブルの大きさは $O( ( \log n) ^ 2 )$、遷移は $O(1)$ となる。

桁 DP を $O(\log n)$ 回行うので、全体の計算量は $O((\log n) ^ 3)$ となる。ただし、様々な部分に定数倍として $10$ や $20$ といった大きな係数が付くので注意。

感想

制約が本当は $10 ^ {50}$ なのを $2 ^ {50}$ と読み違えて long long を使ってて永遠に合わなかったのが悲しい...

$S$ を固定して解いた方が頭が壊れないと思って分離したけど、$S$ を固定せずに DP をして適切に枝刈りをすれば $O((\log n) ^ 2 \log \log n)$ とかになると思います (map でサボってますが、こんなイメージです: 提出)。

ICPC 2019-2020 North-Western Russia Regional Contest C - Cross-Stitch

前置き

公式解説がロシア語なので、自分の解法を書き留めておきます。

問題

ここ

解法

裏面の縫い目の個数は決まっている (表面の縫い目の個数よりちょうど $1$ だけ少ない) ので、裏面のすべての縫い目の長さを $1$ に出来るならば、明らかにそれが最適。実は、これは常に達成可能である。

$(h+1)\times (w+1)$ 頂点のグラフを考え、各 X に対して下図のように有向辺を張る。ここで、X が偶数行目にあるか、あるいは奇数行目にあるかで辺の張り方が異なることに注意する。

f:id:suisen_kyopro:20211116191859p:plain

上の図で、すべての頂点において入次数と出次数が等しい。また、制約から、不要な頂点を取り除いたグラフは連結なのでオイラー閉路が存在する。従って、任意のオイラー閉路を一つ選んで裏面の縫い目を一つ取り除くことで所望の解を得る。時間計算量は $O(hw)$。

実装 (C++)

提出リンク

折り畳み

オイラー路を求めるやつを実際に書くのは初めてで、何かを参照したわけでもないので使用実績はこの問題だけです。何が言いたいかというと、他の問題で使うとバグってて大変なことになるかもしれません。

#include <iostream>

#include <algorithm>
#include <cassert>
#include <optional>
#include <vector>

namespace suisen {
    struct DirectedEulerianGraph {
        DirectedEulerianGraph() {}
        DirectedEulerianGraph(int n) : n(n), g(n), in_deg(n, 0) {}

        void add_edge(int u, int v) {
            g[u].push_back(v);
            ++in_deg[v];
        }

        std::optional<std::vector<int>> euler_circuit(int start = 0) {
            std::size_t edge_num = 0;
            std::vector<std::vector<bool>> used(n);
            for (int i = 0; i < n; ++i) {
                const int sz = g[i].size();
                edge_num += sz;
                used[i].resize(sz, false);
                if (in_deg[i] != sz) return std::nullopt;
            }
            std::vector<int> res;
            auto dfs = [&](auto dfs, int u) -> void {
                for (std::size_t i = 0; i < g[u].size(); ++i) {
                    if (used[u][i]) continue;
                    const int v = g[u][i];
                    used[u][i] = true;
                    dfs(dfs, v);
                }
                res.push_back(u);
            };
            dfs(dfs, start);
            std::reverse(res.begin(), res.end());
            if (res.size() != edge_num + 1) return std::nullopt;
            return res;
        }
        std::optional<std::vector<int>> eulerian_trail() {
            int s = -1, t = -1, invalid = -1;
            for (int i = 0; i < n; ++i) {
                int out_deg = g[i].size();
                if (out_deg == in_deg[i] + 1) {
                    (s < 0 ? s : invalid) = i;
                } else if (out_deg == in_deg[i] - 1) {
                    (t < 0 ? t : invalid) = i;
                } else if (out_deg != in_deg[i]) {
                    invalid = i;
                }
            }
            if (s < 0 or t < 0 or invalid >= 0) return std::nullopt;
            add_edge(t, s);
            auto res = *euler_circuit(s);
            res.pop_back();
            // remove edge (t, s)
            g[t].pop_back();
            if (res.back() == t) return res;
            auto is_ts_edge = [&](int u, int v) {
                return u == t and v == s;
            };
            std::rotate(res.begin(), std::adjacent_find(res.begin(), res.end(), is_ts_edge) + 1, res.end());
            return res;
        }
        
        const std::vector<int>& operator[](int u) const {
            return g[u];
        }
    private:
        int n;
        std::vector<std::vector<int>> g;
        std::vector<int> in_deg;
    };
}

using suisen::DirectedEulerianGraph;

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int w, h;
    std::cin >> w >> h;
    std::vector<std::string> s(h);
    for (auto &row : s) std::cin >> row;

    DirectedEulerianGraph g((w + 1) * (h + 1));
    auto id = [&](int i, int j) {
        return i * (w + 1) + j;
    };
    int si = -1, sj = -1;
    for (int i = 0; i < h; ++i) {
        for (int j = 0; j < w; ++j) {
            if (s[i][j] != 'X') continue;
            if (i & 1) {
                si = i + 1, sj = j;
                g.add_edge(id(i + 1, j + 1), id(i, j));
                g.add_edge(id(i + 1, j), id(i, j + 1));
                g.add_edge(id(i, j), id(i + 1, j));
                g.add_edge(id(i, j + 1), id(i + 1, j + 1));
            } else {
                si = i, sj = j;
                g.add_edge(id(i, j), id(i + 1, j + 1));
                g.add_edge(id(i, j + 1), id(i + 1, j));
                g.add_edge(id(i + 1, j), id(i, j));
                g.add_edge(id(i + 1, j + 1), id(i, j + 1));
            }
        }
    }
    auto ans = *g.euler_circuit(id(si, sj));
    ans.pop_back();
    std::cout << ans.size() - 1 << '\n';
    for (int ij : ans) {
        int i = ij / (w + 1), j = ij % (w + 1);
        std::cout << j << ' ' << i << '\n';
    }
    return 0;
}

感想

明らかに L の方が難しいと思うんだけど、何故かあんまり解かれてなくてバチャ中は見逃してしまった...

AOJ No.2446 Enumeration

問題概要

リンク先 を見てください

解法

制約的に、bit DP をしたい。

$$ \begin{aligned} \mathrm{dp}(S) :=\; & 1 \text{ 以上 } M \text{ 以下の整数 } x \text{ であって、}\\ &\text{ある } i\in S \text{ が存在して } A_i \text{ が } x \text{ を割り切るようなものの個数} \end{aligned} $$

がすべての $S\subseteq\{1,\ldots,N\}$ に対して求まればこの問題が解ける。

集合 $S$ に対して、任意の $i\in \{1,\ldots,N\} - S$ を一つ選んで添加することを考える。即ち、$\mathrm{dp}(S)$ が既知として、$\mathrm{dp}(S\cup\{i\})$ の値を求める。

$\mathrm{dp}(S\cup\{i\})$ と $\mathrm{dp}(S)$ の差分は、以下を全て満たすような整数 $x$ の個数である。

$$ \begin{aligned} & 1\leq x\leq M, \\ & A_i \mid x, \\ & A_j \nmid x\quad \forall j\in S. \end{aligned} $$

これは、2 つめの条件を考慮して $x=x'A_i$ とおけば、

$$ \begin{aligned} & 1\leq x'\leq \left\lfloor M/A_i \right\rfloor, \\ & A_j/\gcd(A_i, A_j) \nmid x'\quad \forall j\in S, \end{aligned} $$

を満たす整数 $x'$ の個数と言い換えられる。従って、包除原理により一つの $\mathrm{dp}(S\cup\{i\})$ を $\Theta(N\cdot 2 ^ {|S|})$ で求められる。具体的には、

$$ \begin{aligned} \mathrm{dp}(S\cup\{i\}) &= \mathrm{dp}(S)+\sum _ {T\subseteq S} (-1)^{|T|} \left\lfloor \dfrac{\left\lfloor \dfrac{M}{A_i}\right\rfloor}{\underset{j \in T}{\mathrm{lcm}} \dfrac{A_j}{\gcd(A_i, A_j)}} \right\rfloor \\ &=\mathrm{dp}(S)+\sum _ {T\subseteq S} (-1)^{|T|} \left\lfloor \dfrac{M}{\underset{j \in T\cup\{i\}}{\mathrm{lcm}}A_j} \right\rfloor \end{aligned} $$

とすればよい。

しかし、このままでは間に合わないので高速化する必要がある。計算をもう少し進めると、

$$ \begin{aligned} \mathrm{dp}(S) &=\mathrm{dp}(\emptyset)+\sum _ {\emptyset\neq T\subseteq S} (-1)^{|T|+1} \left\lfloor \dfrac{M}{\underset{j \in T}{\mathrm{lcm}}\ A_j} \right\rfloor\\ &=M-\sum _ {T\subseteq S} (-1)^{|T|}\left\lfloor \dfrac{M}{\underset{j \in T}{\mathrm{lcm}}\ A_j} \right\rfloor \end{aligned} $$

を得る。高速ゼータ変換を用いることで、すべての $S\subseteq \{1,\ldots,N\}$ に対する $\mathrm{dp}(S)$ を $\Theta(N\cdot 2 ^ N\cdot \log M)$ で求めることができる。なお、$\mathrm{lcm}$ の計算も同時に bit DP で求めることで $\Theta((N + \log M) \cdot 2 ^ N)$ へと計算量を落とすことができる。

Tenka1 Programmer Contest F - ModularPowerEquation!!

問題リンク

atcoder.jp

問題概要

正整数 $ A, M $ が与えられるので、

$$ A ^ x \equiv x \pmod M $$

を満たす $2\times 10 ^ {18}$ 以下の正整数 $x$ を一つ求めよ。

制約

  • $1\leq A\leq 10 ^ 9$
  • $1\leq M\leq 10 ^ 9$

解法

Step 1.

$\gcd(A,M)=1$ ならば $ A ^ { \phi (M) } \equiv 1 \pmod M $ なので、まずはこの場合に帰着したい。

まず、$1\leq x\leq \lfloor \log _ 2 M\rfloor$ の場合を愚直に試し、答えが見つかれば直ちに出力して終了する。

以下、$x\gt\lfloor\log _ 2 M\rfloor$ とする。

このとき、$g\coloneqq\gcd(A ^ {\text{(十分大きな整数)}}, M)$ に対して答えとなる $x$ は $\gcd(A ^ {x},M)=g$ を満たす。従って、$g\mid A ^ x$ および $ g \mid M $ が成り立つから、

$$ \begin{aligned} & A ^ x \equiv x \pmod M \\ \Rightarrow & \exists k \in \mathbb{Z} \; \text{ s.t. } \; A ^ x - x = k M \\ \Rightarrow & \exists k \in \mathbb{Z} \; \text{ s.t. } \; \frac{x}{g} = \frac{A ^ x}{g} - k \cdot \frac{M}{g} \in \mathbb{Z} \\ \Rightarrow & g \mid x \end{aligned} $$

が必要である。$x=gx'$ を満たす正整数 $x'$ および $ M = g M' $ をみたす正整数 $M'$ を導入し、満たすべき条件式を

$$ \begin{cases} g x ' \gt \lfloor \log _ 2 M \rfloor \\ A ^ {g x '} \equiv g x ' \pmod {M '} \end{cases} $$

と書き換える。ここで、一つ目の条件が必要なことに注意する。ただし、求めた答えが小さすぎた場合は条件を満たすようになるまで周期を足すだけでよいので、以降はこの条件を考えない。

さて、この時点で $\gcd(A,M')=1$ となり、扱いやすくなった。加えて、$\gcd(g,M')=1$ も成立していることに注意する。

Step 2.

$ B \coloneqq A ^ g $ とすれば、条件式は

$$ B ^ {x '} \equiv g x ' \pmod{M '} $$

と書ける。$\gcd(A,M')=1$ なので $\gcd(B,M')=1$ であり、つまり左辺の周期は $\phi(M')$ の約数となる。そこで、$x'=q\cdot \phi(M')+r$ を満たす整数 $q,r$ を導入すると

$$ B ^ r \equiv g \cdot ( q \cdot \phi ( M ' ) + r ) \pmod { M ' } $$

となる。$r=0$ としてみると、

$$1\equiv g\cdot q\cdot \phi(M')\pmod{M'}$$

となる。いま $\gcd(g,M')=1$ なので、$\bmod M'$ において $g$ の乗法逆元は必ず存在する。一方、$\phi(M')$ に関しては $\gcd(M',\phi(M'))=1$ とは限らないが、仮にこれが成り立つとすれば、

$$q\equiv\bigl(g\cdot\phi(M')\bigr)^{-1}\pmod{M'}$$

を満たす $q$ を選ぶことで条件を満たす $x'$ を求めることが出来る (先述の通り、この時点で答えが小さすぎた場合は適切に数周期分だけ足す必要がある)。

残ったのは、$\gcd(M',\phi(M'))\neq 1$ の場合である。

Step 3.

条件式を $q,r$ で分離して

$$ g ^ { - 1 } \cdot B ^ r - r \equiv q \cdot \phi ( M ' ) $$

とする。見やすさのため、左辺を一時的に $a(r)$ と表記する。つまり、

$$a(r)\equiv q\cdot \phi(M')\pmod {M'}$$

である。$g'\coloneqq \gcd(M',\phi(M'))\gt 1$ とおけば、上式より

$$ \begin{aligned} & a(r) \equiv q \cdot \phi ( M ' ) \pmod { M ' } \\ \Rightarrow & \exists k \in \mathbb{Z} \; \text{ s.t. } \; a(r) - q \cdot \phi ( M ' ) = k M ' \\ \Rightarrow & \exists k \in \mathbb{Z} \; \text{ s.t. } \; \frac{ a ( r ) }{ g ' } = q \cdot \frac { \phi ( M ' ) }{ g ' }+k \cdot \frac{ M ' } { g ' } \in \mathbb{Z} \\ \Rightarrow& g ' \mid a(r) \end{aligned} $$

が必要である。$\gcd\left(\dfrac{\phi(M')}{g'},\dfrac{M'}{g'}\right)=1$ であるから、$g'\mid a(r)$ を満たす $r$ が一つ求まれば、$q$ は

$$q\equiv\frac{a(r)}{g'}\cdot\left(\frac{\phi(M')}{g'}\right)^{-1}\pmod{\frac{M'}{g'}}$$

を満たすものを選べばよく、この問題が解けたことになる。ここで、$g'\mid a(r)$ を満たす $r$ はどのようなものであるかを考えると

$$ \begin{aligned} & g'\mid a(r)\\ \Leftrightarrow & a(r)\equiv 0 \pmod{g'}\\ \Leftrightarrow & g ^ {-1}\cdot B ^ r\equiv r\pmod{g'} \end{aligned} $$

となり、元の問題と非常によく似た構造が現れる。少し異なるのは左辺に係数が付いている点であるが、この係数を加えても今までと全く同様の議論が行えるため、初めから係数付きの問題を解いており、元の問題はその係数が $1$ であるような特殊ケースだと思えば自然な再帰構造となる。

再帰的に問題を解けばよいことが分かったので、終端条件を確認する。$g'\lt\phi(M')\lt M'$ および $ g' \mid M' \mid M $ から、法 $ M $ が必ず $ M $ 以外の約数へと移っていくことが分かる。従って、最終的には $1$ となるか、あるいは途中で Step 1 や Step 2 で確認したような部分的に解けるケースに帰着されることで終端する。従って $M=1$ のケースだけ処理すれば良く、このケースは容易に解くことが出来る。

以上より、TL に間に合うかと値の上限制約を考えなければ、この問題を解くことが出来た。

上限制約

合同式の周期は $\mathrm{lcm}(\phi(M),M)\leq 10^{18}$ の約数なので、計算途中に適切に mod を取ればこの制約を破らないような解が構成できる。

計算量解析

まず、Step 1 で $1\leq x\leq\lfloor \log_2 M\rfloor$ を試すパートは愚直にやっても $O(\log M\log \log M)$ である。Step 2 や Step 3 では定数回の $\gcd$ およびモジュラ逆数の計算を行うが、これらは $O(\log M)$ である。Step 2, 3 で $\phi(M')$ が必要となるが、これは $N$ の素因数分解表示

$$N=\prod_{p}p^{q_p}$$

に対して

$$\phi(N)=\prod_{p} p^{q_p}\left(1-\dfrac{1}{p}\right)$$

として求められることを利用すれば $O(\sqrt{M})$ で求めることが出来る。以上より、再帰の 1 段分の計算量は $O(\sqrt{M})$ である。

再帰の段数は、$ M $ が毎回半分以下になるので $O(\log M)$ で押さえることができる。

$$ \sum _ {i=0} ^ {\infty}\frac{1}{2 ^ i}+\frac{1}{2 ^ i\sqrt{2}}=2+\sqrt{2}\lt\infty$$

より、全体でも $O(\sqrt{M})$ であり、十分高速である。

所感

銀 diff なので当然難しいし、解説はかなり天才的でびっくりした。ただどこかで同じような議論を見たことがある気がしていて、不動点典型 (?) だったりするのかな。

自分の解法と公式解説で共通してるのはある程度 $ x $ が大きいことを課せば再帰できるってことで、これは離散対数問題を任意 mod で解く時にも使うから他の問題でも使うことがあるかもしれないと思ったり。

第6回 ドワンゴからの挑戦状 予選 E - Span Covering の Θ(X^2((log X)^2+log N)+N) 解法

(追記 2022/04/29) : 計算量を  \Theta (X ^ 2 ( ( \log X ) ^ 2 + \log N ) + N ) に改善できたので、追記しました。

問題

atcoder.jp

解法

包除原理の適用を考える。具体的には、

 \displaystyle
A_i:= \text{区間 $[i-1,i)$ が被覆されていないような置き方の集合}

とおき、非負整数  k に対して集合  S_k

 \displaystyle
S_k:=\{(a_0,\ldots,a_k)\mid a_0+\cdots+a_k=X+1,\;\text{$a_i$ は正整数} \}

で定め、

 \displaystyle
\begin{aligned}
\mathrm{ans}
&=\sum_{k=0}^X (-1)^k \sum_{(a_0,\ldots,a_k)\in S_k}
\left\vert \bigcap_{j=0}^k A_{a_0+\cdots+a_j}\right\vert
\\
&=\sum_{k=0}^X (-1)^k \sum_{(a_0,\ldots,a_k)\in S_k}\prod_{i=1}^N\sum_{j=0}^k\max\{0,a_j-L_i\}\\
\end{aligned}

により答えを求める。

 X=15,k=4,(a _ 0,a _ 1,a _ 2,a _ 3,a _ 4)=(3,4,1,5,3)\in S_k の例

 \displaystyle
\begin{aligned}
\sum_{j=0}^k\max\{0,a_j-L_i\}&=\sum_{j\in\{j'\mid a_{j'}\geq L_i\}}(a_j-L_i)\\
&=\left(\sum_{j\in\{j'\mid a_{j'}\geq L_i\}}a_j\right)-L_i\cdot \vert \{j'\mid a_{j'}\geq L_i\}\vert
\end{aligned}

なので、 x に対して  c := \vert \{ j' \mid a _ {j' } \geq x\} \vert および  \displaystyle s := \sum _ { j \in \{ j' \mid a _ {j'} \geq x\}} a_j の値が分かれば

 \displaystyle
\begin{aligned}
\sum_{j=0}^k\max\{0,a_j-x\}&=s-cx
\end{aligned}

として計算することが出来る。そして、これは  a_j\geq x となる  a_j を全て決めてしまえば 以降変化することが無い

以上を踏まえて、次のような DP を考える。

 \displaystyle
\mathrm{dp} [x] [c] [s] := \sum _ { \overset { \scriptstyle a _ 0 + \cdots + a _ {c-1} = s , }{a _ 0, \ldots , a _ {c-1} \geq x} }\prod _ {i \in \{ i' \mid L_{i'} \geq x\} } \sum _ {j=0} ^ {c-1} \max \{0, a _ j - L _ i \}

このテーブルを用れば、

 \displaystyle
\sum_{k=0}^{X}(-1)^k\cdot \mathrm{dp}[1][k+1][X+1]

として答えを求めることが出来る。

 \mathrm{dp} テーブルの計算を考える。初期化は容易であり、

 \displaystyle
\mathrm{dp}[X+2][c][s]=\begin{cases}
1&\text{if $c=s=0$ }\\
0&\text{otherwise}
\end{cases}

とすればよい。更新は、 x の降順に求めていくことにすれば、列  a x を挿入する個数  p で場合分けをし、以下のように計算できる。

 \displaystyle
\mathrm{dp}[x][c][s]=(s-cx)^{\mathrm{cnt}[x]}\sum_{p=0}^{\min\{c,\lfloor s/x\rfloor\}} \binom{c}{p}\cdot \mathrm{dp}[x+1][c-p][s-px]

ここで、 \mathrm{cnt}[x]:= \vert\{i\mid L_{i}=x\}\vert とおいた。二項係数および  \mathrm{cnt} (s-cx) ^ {\mathrm{cnt}[x]} の取得が  O(1) でできる仮定すれば、調和級数の部分和の解析を用いることで、愚直に更新しても全体  \Theta(X ^ 3\log X) でこのテーブルのすべての値を求めることが出来ることが分かる。

このままでも十分高速であるが、少しの工夫でさらに計算量を落とすことが出来る。 c について  c\leq \left\lfloor\dfrac{X+1}{x}\right\rfloor の部分しか見なくてよいこと、 s について  s\geq cx の部分しか見なくてよいことから、べき乗等が  O(1) で求まるという仮定の下で、テーブル更新の計算量は

 \displaystyle
\begin{aligned}
\sum_{x=1}^{X+1}\sum_{c=0}^{\lfloor(X+1)/x\rfloor}\sum_{s=cx}^{X+1} (c+1)&=
\sum_{x=1}^{X+1}\sum_{c=0}^{\lfloor(X+1)/x\rfloor}(X+2-cx)(c+1)\\
&\leq (X+2)\sum_{x=1}^{X+1}\sum_{c=0}^{\lfloor(X+1)/x\rfloor}(c+1)\\
&\leq (X+2)\sum_{x=1}^{X+1}\sum_{c=0}^{\lfloor(X+1)/x\rfloor}\left(\left\lfloor\dfrac{X+1}{x}\right\rfloor+1\right)\\
&\leq (X+2)\sum_{x=1}^{X+1}\left(\dfrac{X+1}{x}+1\right)^2\\
&= (X+2)\sum_{x=1}^{X+1}\left( \dfrac{X^2}{x^2} + O(X) \right) \\
&= O(X^3) \quad\left(\because\sum_{i=1}^\infty \frac{1}{i^2}=\frac{\pi^2}{6}\lt \infty \right)\\
\end{aligned}

となり  \log を落とすことが出来る。また、以下に示すようにこの和は  \Omega(X ^ 3) であるから、即ち  \Theta(X ^ 3) である。

証明

 \displaystyle
\begin{aligned}
\sum_{x=1}^{X+1}\sum_{c=0}^{\lfloor(X+1)/x\rfloor}\sum_{s=cx}^{X+1} (c+1)&\geq 
\sum_{c=0}^{X+1}(X+2-c)(c+1) \quad(x=1\; \text{の部分だけを取り出した})\\
&=\dfrac{1}{6}(X+2)(X+3)(X+4)\\
&\geq \dfrac{1}{6} X ^ 3
\end{aligned}

より、示された。(証明終)

前計算に関して、二項係数に  \Theta(X ^ 2) \mathrm{cnt} \Theta(N)、べき乗に  \Theta(NX) だけかけるとすると、全体  \Theta(NX+X ^ 3) でこの問題を解くことが出来た。

なお、 \mathrm{cnt} が入力として与えられるような場合で  \mathrm{cnt}[x] が非常に大きくなりうるときは、各  x に対して  a ^ {\mathrm{cnt}[x]}\;(0\leq a\leq X+1) \Theta(X ^ 2 \log C) (ただし、 C:= \max_x \mathrm{cnt}[x] とした) だけかけて前計算するなどしてべき乗を高速に求められるようにすれば、全体  \Theta(X ^ 3 + X ^ 2 \log C) となる。なお、都度べき乗を計算しても  \Theta(X ^ 3 + X ^ 2 \log X \log C) である。

実装

提出 (C++, 504ms)

更なる高速化について (2022/04/29 追記)

以下の dp において、各  x=1,\ldots,X+1 に対して  \displaystyle 0 \leq c \leq \left\lfloor \dfrac{X + 1}{x} \right\rfloor の部分だけを見ればよいということであった。

 \displaystyle
\begin{aligned}
\mathrm{dp}[X+2][c][s]&=\begin{cases}
1&\text{if $c=s=0$ }\\
0&\text{otherwise}
\end{cases}, \\

\mathrm{dp}[x][c][s] &= (s-cx)^{\mathrm{cnt}[x]}\sum_{p=0}^{\min\{c,\lfloor s/x\rfloor\}} \binom{c}{p}\cdot \mathrm{dp}[x+1][c-p][s-px].
\end{aligned}

状態数は  \Theta(X ^ 2 \log X) であるが、遷移に  \Theta(c) だけ掛かっていたためテーブル更新の計算量が全体で  \Theta(X ^ 3) となっていた。ここでは、遷移の計算量を削減することで計算量を  \Theta(X ^ 2 (\log X) ^ 2) に落とす方法について述べる。即ち、元の問題は  \Theta (X ^ 2 ( ( \log X ) ^ 2 + \log N ) + N ) 時間で、入力で  \mathrm{cnt} が与えられる場合の問題は、 \mathrm{cnt} の最大値を  C として  \Theta (X ^ 2 ( ( \log X ) ^ 2 + \log C ) ) 時間で解くことが出来る。

dp の更新式を変形をすると、次のようになる。

 \displaystyle
\dfrac{\mathrm{dp}[x][c][s]}{c!} = (s-cx)^{\mathrm{cnt}[x]}\sum_{p=0}^{\min\{c,\lfloor s/x\rfloor\}} \dfrac{1}{p!} \cdot \dfrac{\mathrm{dp}[x+1][c-p][s-px]}{(c-p)!}.

ここで  \displaystyle \mathrm{dp'}[x][c][s]:=\dfrac{\mathrm{dp}[x][c][s]}{c!} と定義すれば、 \mathrm{dp'} は次のように計算できる。

 \displaystyle
\begin{aligned}
\mathrm{dp'}[X+2][c][s]&=\begin{cases}
1&\text{if $c=s=0$ }\\
0&\text{otherwise}
\end{cases}, \\
\mathrm{dp'}[x][c][s] &= (s-cx)^{\mathrm{cnt}[x]}\sum_{p=0}^{\min\{c,\lfloor s/x\rfloor\}} \dfrac{1}{p!} \cdot \mathrm{dp'}[x+1][c-p][s-px].
\end{aligned}

 x を固定し、 V[c][s] := \mathrm{dp'}[x+1][c][s] とおく。各  c, s に対する

 \displaystyle
\sum _ {p=0} ^ {\min\{c,\lfloor s/x\rfloor\}} \dfrac{1}{p!} \cdot V[c-p][s-px]

が高速に求まればよいが、 s = ax + b\; (0\leq b \lt x) と表すと、 求める和は以下のように書ける。

 \displaystyle
S _ b [c][a] := \sum_{p=0}^{\min\{c,a\}} \dfrac{1}{p!} \cdot V[c-p][(a-p)x+b]

従って、多項式  f および  g _ b\; (0\leq b\lt x) を以下のように定義すれば、

 \displaystyle
\begin{aligned}
f(t, u) &:= \sum _ {n = 0} ^ \infty \dfrac{1}{n!} t ^ n u ^ n \\
g _ b (t, u) &:= \sum _ {n = 0} ^ { \lfloor (X + 1) / x \rfloor } \sum _ {m = 0} ^ { \lfloor (X + 1 - b) / x \rfloor } V[n][mx + b] t ^ n u ^ m
\end{aligned}

 S _ b [c][a] = [t ^ c u ^ a] f \ast g _ b として計算される。 b を固定すると  f \ast g _ b高速フーリエ変換などにより  \Theta\left(\dfrac{X ^ 2}{ x ^ 2 } \log X\right) 時間で計算できるので、全ての  b に対して  f \ast g _ b を計算するのは、 \Theta\left(\dfrac{X ^ 2}{ x } \log X\right) 時間で可能である。従って、 \mathrm{dp}' 全体を求めるのにかかる計算量は  \Theta(X ^ 2 ( (\log X) ^ 2 + \log C) ) となる。

上記の解法のボトルネックは畳み込みの計算であり、定数倍は軽くない。 f の形が特殊なことを利用すると、定数倍を多少改善できる。

 \displaystyle
\begin{aligned}
S' _ {b, d} [c] &:= \sum_{p=0}^{c} \dfrac{1}{p!} \cdot V[c-p][(c-p+d)x+b], \\
f' (t) &:= \sum _ {n = 0} ^ \infty \dfrac{1}{n!} t ^ n \\
g' _ {b, d} (t) &:= \sum _ {n = 0} ^ { \lfloor (X + 1 - b) / x \rfloor } V[n][(n+d)x + b] t ^ n
\end{aligned}

として (記号が紛らわしいが、 '微分ではない)、 0 \leq b\lt x,\; 0\leq d\leq \lfloor (X+1-b)/x \rfloor を固定すれば、 S' _ {b, d} [c] = [t ^ c] f' \ast g' _ {b, d} として計算される。従って、 \lfloor (X + 1) / x \rfloor 次程度の多項式の畳み込みを  X+1 回以下だけ行うことで全ての  S' _ {b, d} [c] を計算することが出来る。これは、 \lfloor (X + 1) ^ 2 / x ^ 2 \rfloor 次程度の多項式の畳み込みを  x 回行うよりも定数倍が良い。

なお、 0 \leq d としてよい理由は、 s \lt cx ならば  \mathrm{dp}'[x][c][s] = 0 が成り立つためである。

実装

提出 (C++, 440ms)

 \mathrm{mod}\;10 ^ 9 + 7 での畳み込みということもあって遅いですね。

AtCoder Beginner Contest 213 H - Stroll

前置き

公式解説とは異なる  O(N ^ 3 T\log T) 解法を説明します。

問題リンク

atcoder.jp

解法

まず、形式的冪級数  f _ i を次のように定めます。

 \displaystyle
f_i=\sum_{d=0}^\infty \text{($d$ キロメートル歩いて頂点 $i$ にいる場合の数)}\cdot x^d

また、 1\leq i,j\leq N に対して形式的冪級数  g _ { i , j } を次のように定めます。

 \displaystyle
g_{i,j}=\sum_{d=0}^\infty \text{($j$ から $i$ への長さ $d$ キロメートルの道の本数)}\cdot x^d

このとき、 f g の間に次の関係式が成り立ちます。

 \displaystyle
f_i=\sum_{j=1}^N f_j\ast g_{i,j}+ \begin{cases}
1 & \text{if $\;i=1$} \\
0 & \text{otherwise}
\end{cases}

ここで、 \ast は畳み込み演算です。この関係式を行列を用いて表せば、 f に関する方程式が得られます。

 \displaystyle
\begin{pmatrix} f_1 \\ f_2\\ \vdots \\ f_N\end{pmatrix}=
\begin{pmatrix}
g_{1,1} & g_{1,2} & \cdots & g_{1,N} \\
g_{2,1} & g_{2,2} & \cdots & g_{2,N} \\
\vdots & \vdots & \ddots & \vdots \\
g_{N,1} & g_{N,2} & \cdots & g_{N,N}
\end{pmatrix}
\begin{pmatrix} f_1 \\ f_2\\ \vdots \\ f_N\end{pmatrix}
+
\begin{pmatrix} 1 \\ 0\\ \vdots \\ 0\end{pmatrix}

これを  f について解きます。

 \displaystyle
\begin{pmatrix}
f _ 1 \\ f _ 2\\ \vdots \\ f _ N
\end{pmatrix}=
\begin{pmatrix}
1 - g _ {1, 1} & g _ {1,2} & \cdots & - g _ {1,N} \\ - g _ {2,1} & 1 - g _ {2,2} & \cdots & - g _ {2,N} \\ \vdots & \vdots & \ddots & \vdots \\ - g _ {N,1} & - g _ {N,2} & \cdots & 1 - g _ {N,N}
\end{pmatrix} ^ {-1}
\begin{pmatrix}
1 \\ 0\\ \vdots \\ 0
\end{pmatrix}

ここで、右辺に現れる逆行列が存在するかが問題となります。これは、制約から  [x ^ 0]g _ { i , j } = 0 であり、対角成分にしか非零の定数項が存在しないことに注目するとその存在を示すことが出来ます*1

以上より、 T+1 次以上の項を切り捨てながら掃き出し法を行うことで、 O(N ^ 3 T\log T)逆行列 *2 が求まり、この問題を解くことが出来ました。

なお、逆行列のすべての成分を計算すると C++ で 4000 ms 程度かかってしまいますが、この問題では逆行列 (1,1) 成分しか必要ないためかなりの計算を省略することができ、1000 ms を切ることが出来ます。

実装

逆行列のすべての成分を計算する実装例 (C++, 4132 ms)

逆行列の (1,1) 成分だけを計算する実装例 (C++, 914 ms)

*1:対角成分にしか非零の定数項が存在しないことから、行列式

 \displaystyle
\sum _ { \sigma \in \mathfrak { S } _ N} \mathrm { sgn } ( \sigma ) \prod _ { i = 1 } ^ N (E-G)_{i,\sigma(i)}

において  \displaystyle [x ^ 0] \prod _ { i = 1 } ^ N (E-G)_{i,\sigma(i)} \neq 0 となり得る  \sigma は恒等置換  \sigma = \mathrm{Id} のみであり、このとき  \displaystyle [x ^ 0] \prod _ { i = 1 } ^ N ( 1 - g _ { i , i } ) = 1 です。従って、行列式が可逆なので正則であることが言えます。(可換環の元を係数に持つ形式的冪級数可換環を成し、可換環上の行列が正則であるための必要十分条件は、その行列式が可逆元であること (らしいです、私は代数に明るくないので、誤りがあれば指摘頂けるとありがたいです) )

あるいは、掃き出し法をシミュレートする過程で行列の対角成分の定数項が  1 のまま変化しない (つまり、可逆である) ことからも逆行列の存在が分かります。

*2:項を切り捨ててしまうと厳密な逆行列にはなりませんが、この問題を解く上では  T 次以下の項さえ正しければよいのでこのような操作をしても構いません。

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

こどふぉブログ 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


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

おまけ

紹介したテクニックを用いれば, \begin{bmatrix}N\\ N-i\end{bmatrix}\bmod P を全ての i=0,\ldots,N に対して  O(N\log N) で計算することができます.

 a _ {2n, i} を求める式をさらに変形します.

 \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}^\infty \dfrac{n^k}{k!}x^k\\
g_n:=&\sum_{k=0}^n a_{n,k}\cdot(n-k)!\cdot x^k
\end{aligned}\right)\\
&=[x^i]a_nh_n\quad\left(h_n:=\sum_{j=0}^n\dfrac{[x^j]f_ng_n}{(n-j)!}x^j\right)
\end{align}

上式より  a _ {2n} = a _ n h _ n が成り立ち,2 回の畳み込みで  a _ {2n} が求まることが分かります. f_n n 次までの項で打ち切ってよいことから, a_n から  a_{2n} を求める計算は  O(n\log n) でできることが分かりました.

 i 回目の反復で列の長さは  O(2 ^ i) であり,反復の回数は  O(\log N) 回であることから,全体の計算量は  O(N\log N) です.

実装例 (C++)

折り畳み

注意: 一部独自のクラスを用いているため,これをそのまま貼っても動きません

template <typename mint>
std::vector<mint> stirling_number1_reversed(int n) {
    std::vector<mint> fac(n + 1), fac_inv(n + 1);
    fac[0] = 1;
    for (int i = 1; i <= n; ++i) fac[i] = fac[i - 1] * i;
    fac_inv[n] = fac[n].inv();
    for (int i = n; i > 0; --i) fac_inv[i - 1] = fac_inv[i] * i;
    int l = 0;
    while ((n >> l) != 0) ++l;
    // 形式的冪級数のクラス
    FPS<mint> a {1};
    int m = 0;
    while (l --> 0) {
        // f, g は m 次の多項式
        FPS<mint> f(m + 1), g(m + 1);
        mint powm = 1;
        for (int i = 0; i <= m; ++i, powm *= m) {
            f[i] = powm * fac_inv[i];
            g[i] = a[i] * fac[m - i];
        }
        // f := f * g とし,(定数倍改善のために) f の m+1 次以降の係数を捨てる
        f *= g, f.pre_inplace(m);
        for (int i = 0; i <= m; ++i) f[i] *= fac_inv[m - i];
        // a = a * f とすることで a_{2n} が求められる
        // 次数 m を 2 倍し,m+1 次以降の係数を捨てる (必須)
        a *= f, m *= 2, a.pre_inplace(m);
        // a_{m} -> a_{m+1} の更新
        if ((n >> l) & 1) {
            a.push_back(0);
            for (int i = m; i > 0; --i) a[i] += m * a[i - 1];
            ++m;
        }
    }
    return a;
}


おまけ2

元の問題も同様に高速化することが出来ます. K+1 次以上の項を無視してよいので,次のように  a_{2n} を求めることが出来ます.

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

ここで, g _ n h _ n の係数を求める際に最大で  N/2 程度の階乗を計算する必要があります. N が小さい場合は前計算が出来るので問題とはなりませんでしたが, N が大きくなってくるとそうもいきません.

これは簡単には求まりませんが,1 つの階乗を  O(\sqrt{P}(\log P) ^ 2) あるいは  O(\sqrt{P}\log P) で求めるアルゴリズムが存在します 1.従って, 0!,\ldots,K! を前計算しておき,反復の度に都度  (n-\min(n,K))! を高速階乗計算により求めることで全体  O((K\log K + \sqrt{P}\log P)\log N) で問題を解くことが出来ます.

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


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