任意 mod で二項係数を列挙する (2)


suisen-kyopro.hatenablog.com これを $\displaystyle O\left(\frac{N\log M}{\log\log M}\right)$ にします.

基本的なアイデア素因数分解であり, 前記事と同様です. 式 $(1)$ を用いるのも同じです. しかし, 今回セグメント木は使いません.

$$
C(N,i)=\begin{cases}\displaystyle 1& \text{(if $i=0$)} \\\displaystyle\frac{N-i+1}{i}\cdot C(N,i-1)& \text{(otherwise)}\end{cases}\tag{1}
$$

前回の記事では, \(M\) の素因数以外の素数についても指数を管理して累乗を計算していました. しかし, これには無駄な点が多いです. というのも, \(M^2\leq 2^{63}-1\) (符号付き 64 bit 整数の最大値) の範囲では, 実は \(M\) の素因数の個数は高々 $9$ つしかありません. *1

\(M\) と互いに素であれば \(\bmod M\) 上の乗算に関して逆元が存在するので, 除算が可能です. そこで, 「除算が可能な部分は逆元を使って楽に計算しよう」というのが今回のメインテーマです.

以下では, \(M\) の素因数は小さい順に \(p_1, \dots, p_K\) であるとします.

除算が可能な部分と不可能な部分を分けるため, \(X=C(N,i)\) を次のように分解します. $$X=Y\cdot {p_1}^{q_1}{p_2}^{q_2}\cdots {p_K}^{q_K}.$$ また, 上式において $Z={p_1}^{q_1}{p_2}^{q_2}\cdots {p_K}^{q_K}$ とします.

その他の整数 $A$ についても同様に $$A=(M \text{と素因数を共有しない部分}) \times (M \text{と素因数を共有する部分})$$ のように分割して考えます. 表記の簡単のため, $$\begin{align}S_M(A)&=(M \text{と素因数を共有しない部分}),\\ T_M(A)&=(M \text{と素因数を共有する部分})\end{align}$$ と定義しておきます *2. つまり, $S_M(X)=Y$, $T_M(X)=Z$ ということになります.

前置きはこれくらいにして, 本計算に入ります. $Z$ の指数部分 $q_i$ の更新方法は前記事と基本的に同様です. しかし, 今回は $p_1,\dots ,p_K$ 以外の素因数の指数は管理する必要がありません.

そこで, 式 $(2)$ を満たすような新しい配列 $d$ を用意します *3.

$$
d_i=\max \{p_j\mid p_j \text{は $i$ を割り切る} \} \cup\{0\}\quad (1\leq i\leq N) \tag{2}
$$

この $d_i$ は, 擬似的には以下のように求められます.

// primeFactors(M) は M の素因数を小さい順に列挙する
for (int prime : primeFactors(M)) {
    for (int i = prime; i <= N; i += prime) {
        d[i] = prime;
    }
}

計算回数は \(\displaystyle\sum_{p_i}\frac{N}{p_i}\) 回程度になりますが, これは次のように評価されます. $$\displaystyle\sum_{p_i}\frac{N}{p_i}\leq\sum_{p_i}N=NK=O\left(\frac{N\log M}{\log\log M}\right).\tag{3}$$

いま求めた $d$ を用いることで, 更新式に現れる $N-i+1$ および $i$ に対して, $S_M(N-i+1)$, $T_M(N-i+1)$, $S_M(i)$, $T_M(i)$ を効率的に求めることが出来ます. 以下にコードのイメージを示します.

int s1 = N - i + 1;
// div は p_1,...,p_K のいずれか
for (int div = d[s1]; d[s1] > 0; div = d[s1]) {
    int count = 0;
    do {
        s1 /= div;
        count++;
    } while (s1 % div == 0);
    // この時点で count は N-i+1 が div で割り切れる回数 (T_M(N-i+1) における div の指数) となっているので, Z の指数を更新する.
    // indexOf() は p_i -> i を表す仮の関数
    q[indexOf(div)] += count;
}
// この時点で s1=S_M(N-i+1) となっているので, Y を更新する.
Y *= s1;

int s2 = i;
for (int div = d[s2]; d[s2] > 0; div = d[s2]) {
    int count = 0;
    do {
        s2 /= div;
        count++;
    } while (s2 % div == 0);
    q[indexOf(div)] -= count;
}
// 実際は mod M での s2 の逆元を掛ける
Y /= s2;

計算量を考えましょう. Y/=s2 の部分で逆元が必要になっていますが, 合成数 $\bmod$ の場合でも $1,\dots,N$ に対する逆元の列挙は $O(N)$ で可能です*4.

割り算パートの解析は難しそうですが, 全体で $\displaystyle O\left(\frac{N\log M}{\log\log M}\right)$ であることが示せます. 各 $p_i$ について, $1,2,\dots,N$ の中に素因数として合計いくつ含まれるかが分かれば良いです. これは, 有名な 「$N!$ は素数 $P$ で何回割り切れるか?」と同じ問題です. そして, この問題の答えは $\displaystyle \sum_{i=1}^{\infty}\left\lfloor\frac{N}{P^i}\right\rfloor $ です *5. 従って, 計算量は $$\sum_{p_i}\sum_{j=1}^{\infty}\left\lfloor\frac{N}{{p_i}^j}\right\rfloor\leq\sum_{p_i}\sum_{j=1}^{\infty}\frac{N}{{p_i}^j}=\sum_{p_i}\frac{N}{p_i-1}=O\left(\frac{N\log M}{\log\log M}\right) \tag{4}$$ です (最後の変形は式 $(3)$ と同じ理屈です).

最後に素因数分解形 ${p_1}^{q_1}{p_2}^{q_2}\cdots {p_K}^{q_K}$ から $Z$ の値を復元するパートを考えます. 基本的には以下に示すコードのように愚直に計算します.

long Z = 1;
for (int j = 0; j < K; j++) {
    Z *= pow(p[j], q[j]);
}

しかし, この場合 pow(p[j], q[j]) の部分が $O(1)$ でなければ \(\displaystyle O\left(\frac{N\log M}{\log\log M}\right)\) は達成できません. そこで, $\displaystyle K=O\left(\frac{\log M}{\log\log M}\right)$ なので全て前計算してしまいましょう. $Z$ は $N!$ を割り切るので, 各 $p_i$ に対して $\displaystyle \sum_{j=1}^{\infty}\left\lfloor\frac{N}{{p_i}^j}\right\rfloor\leq\frac{N}{p_i-1}$ まで計算しておけば十分です. そして, この前計算に掛かる計算量は $(4)$ より \(\displaystyle O\left(\frac{N\log M}{\log\log M}\right)\) です.

以上より, 全体 \(\displaystyle O\left(\frac{N\log M}{\log\log M}\right)\) で $C(N,i)$ を列挙することが出来ました. 最後に, Java でのコードを掲載します. 篩パートで計算量に $O(\log N)$ を付けたくなかったので, 線形篩を用いています. また, \(M\) の素因数は $N$ 以下のものだけを列挙すれば十分なので, 篩の途中で $M$ の素因数を探しています.

verify 問題: https://atcoder.jp/contests/arc012/tasks/arc012_4
verify 提出: https://atcoder.jp/contests/arc012/submissions/18290575

*1:$2\cdot 3\cdot 5\cdot 7\cdot 11\cdot 13\cdot 17\cdot 19\cdot 23\cdot 29=6469693230$ で, \(6469693230^2\geq 2^{63}\) です.

*2:この辺りの分かりやすい/一般的な記号があれば教えて下さい...

*3:\(\max\) を付けているのは, 尺取りをするためです

*4: 参考 : 線形篩 – 37zigenのHP

*5: 参考 : ルジャンドルの定理(階乗が持つ素因数のべき数) | 高校数学の美しい物語

任意 mod で二項係数を列挙する


注意: 計算量解析が粗い可能性があります. より厳密に評価できるのであれば教えて頂けると幸いです.

$N=10^5$, $M=10^9$ くらいの場合に, 二項係数 ${}_NC_0, {}_NC_1, \dots, {}_NC_N \bmod M $ を $\displaystyle O\left(\frac{N(\log N)^2}{\log\log N}\right)$ で列挙します. 以降, 表記の簡単のため, 二項係数を $C(N,i)$ のように表すことにします.

\(M\) が素数の場合は, 「1000000007 で割ったあまり」の求め方を総特集! 〜 逆元から離散対数まで 〜 - Qiita で紹介されているように, 任意の \(1\leq x\lt M\) に対して \(x\cdot x'\equiv 1 \pmod M\) を満たす $x'$ ($x$ の \(\bmod M\) での乗法逆元) が存在するので, 計算は容易です.

しかし, \(M\) が素数でない場合は, 逆元が存在するとは限りません. そこで, 別の方法を考えます.

まず, $C(N,i)$ は式 $(1)$ の漸化式で求まるので, 初期値を $X=1$ として, $i=1,2,\dots N$ の順に $X=C(N,i)$ を計算します.

$$
C(N,i)=\begin{cases}\displaystyle 1& \text{(if $i=0$)} \\\displaystyle\frac{N-i+1}{i}\cdot C(N,i-1)& \text{(otherwise)}\end{cases}\tag{1}
$$

ここで, 素因数分解の形で $X$ を管理することにします. つまり, 素数 $p_1, p_2, \dots$ を用いて, $X={p_1}^{q_1}{p_2}^{q_2}\cdots$ の形で保持します. 素数は無限個ありますが, $C(N,i)$ が持つ素因数は $N$ 以下なので, $N$ 以下の素数を列挙しておけば十分です. 以下では, $N$ 以下の素数の個数を $K$ とします.

まず, 素因数分解の指数 $q_1, q_2,\dots, q_K$ の更新を考えます. $N-i+1={p_1}^{a_1}{p_2}^{a_2}\cdots {p_K}^{a_K}$, $i={p_1}^{b_1}{p_2}^{b_2}\cdots {p_K}^{b_K}$とすると, $q_i\leftarrow q_i+a_i-b_i$ により更新すればよいです.

ただし, $\displaystyle K=O\left(\frac{N}{\log N}\right)$ なので, 全ての $i$ に対して毎回更新していると計算量が $\displaystyle O\left(\frac{N^2}{\log N}\right)$ となってしまいます. そこで, $a_i\neq 0$ または $b_i\neq 0$ であるような $i$ のみを選んで更新することにします. 定数倍高速化をしているだけに見えるかもしれませんが, 実はこのような $i$ の個数は, $\displaystyle O\left(\frac{\log N}{\log\log N}\right)$ 個です *1.

エラトステネスの篩において, 自身を割り切る素数を覚えておくようにすれば素因数分解は $O(\log N)$ で可能 *2 なので, 結局 $q$ の更新は全体 $O(N\log N)$ で出来ます.

素因数分解が求まっても, 元の値が高速に求まらなければ意味がありません. 例えば, 毎回愚直に $X={p_1}^{q_1}{p_2}^{q_2}\cdots {p_K}^{q_K}$ を計算していると, 累乗が $O(1)$ で計算できたとしても全体で $\displaystyle O\left(\frac{N^2}{\log N}\right)$ となってしまいます.

いま, 求めたいのは $X={p_1}^{q_1}{p_2}^{q_2}\cdots {p_K}^{q_K}$ です. $p_i$ は変化せず, $q_i$ は高速に更新できることが分かっています. そこで, $r_i={p_i}^{q_i}$ とし, これを \(\bmod M\) 上の乗算 $\times$ を二項演算とするセグメント木に載せます. すると, \(X=r_1\times r_2\times \cdots \times r_K\) は $O(\log K)=O(\log N)$ で計算できます. ただし, 更新箇所が $\displaystyle O\left(\frac{\log N}{\log\log N}\right)$ 個あるので, セグメント木パートの計算量は全体で $\displaystyle O\left(\frac{N(\log N)^2}{\log\log N}\right)$ です.

最後に, $r_i={p_i}^{q_i}$ の更新にかかる計算量を確認します. 二分累乗法により $O(\log q_i)$ で計算できますが, $q_i$ はどれくらい大きくなるでしょうか. $X$ は $N!$ を割り切るので, $N!$ を素因数分解した場合の $q_i$ の大きさを評価してみます. $S=\lceil \log_{p_i} N\rceil$ とすると, ${p_i}^S\geq N$ より $q_i\leq NS$ です. 従って, $\log_2 q_i\leq\log_2 N + \log_2 \lceil \log_{p_i} N\rceil=O(\log N)$ となります. 更新箇所は毎回 $\displaystyle O\left(\frac{\log N}{\log\log N}\right)$ 個なので, $r_i$ の更新パートは全体で $\displaystyle O\left(\frac{N(\log N)^2}{\log\log N}\right)$ となります.

以上より, 全体 $\displaystyle O\left(\frac{N(\log N)^2}{\log\log N}\right)$ で任意 mod での $C(N,i)$ を列挙することが出来ました.

最後に, Java での実装例を掲載して終わりにしたいと思います.

verify 問題: https://atcoder.jp/contests/arc012/tasks/arc012_4
verify 提出: https://atcoder.jp/contests/arc012/submissions/18243372

*1: 素因数 - Wikipedia

*2:全ての素因数は $2$ 以上なので, 割るごとに値は半分以下になります. 従って, 結局この割り算は高々 $\lceil \log_2{N} \rceil$ 回で終わります.

第四回 アルゴリズム実技検定 (バーチャル参加) [L-O]

AtCoder による 第四回 アルゴリズム実技検定 の問題を解きました.

f:id:suisen_kyopro:20201110204452p:plain

L - マンションの改築

問題文

難しい. 偶数番目を奇数番目に分けて, 差分を imos で管理する.

$A_1=H_1$, $A_i=H _ i-H _ {i-1}$ $(1\lt i\leq N)$ として, 下図のように要素を並べる.

f:id:suisen_kyopro:20201110223211p:plain

すると, クエリ 1 は上行の全要素を $+v$ して下行の全要素を $-v$ する操作である. 同様に, クエリ 2 は上行の $A_1$ 以外の要素を $-v$ して下行の全要素を $+v$ する操作である. クエリ 3 は, $A_i\leftarrow A_i+v$, $A_{i+1}\leftarrow A _ {i+1} -v$ とすればよい.

各クエリの後に 0 の個数が求まればよいので, 行ごとに値の頻度を Map で管理することを考える. クエリ 1 とクエリ 2 は Map の key を全て変化させる操作なので, 愚直に変更するわけにはいかない. そこで, 各行に対して下駄を表す変数 (プログラム中の g0, g1) を持ち, クエリ 1 とクエリ 2 ではこの下駄を $\pm v$ することにすればよい (クエリ 2 では  A_1 のみ別で処理する必要があるので注意).

    public static void solve(ExtendedScanner sc, FastPrintStream pw) {
        int n = sc.nextInt();
        int q = sc.nextInt();
        int n0 = n - n / 2;
        int n1 = n / 2;
        long[] h0 = new long[n0];
        long[] h1 = new long[n1];
        for (int i = 0; i < n; i++) {
            if ((i & 1) == 0) {
                h0[i >> 1] = sc.nextLong();
            } else {
                h1[i >> 1] = sc.nextLong();
            }
        }
        HashMap<Long, Integer> m0 = new HashMap<>();
        HashMap<Long, Integer> m1 = new HashMap<>();
        add1(m0, h0[0]);
        for (int i = n - 1; i > 0; i--) {
            if ((i & 1) == 0) {
                h0[i >> 1] -= h1[(i >> 1) - 1];
                add1(m0, h0[i >> 1]);
            } else {
                h1[i >> 1] -= h0[i >> 1];
                add1(m1, h1[i >> 1]);
            }
        }
        long g0 = 0, g1 = 0;
        while (q --> 0) {
            int t = sc.nextInt();
            if (t == 1) {
                int v = sc.nextInt();
                g0 += v;
                g1 -= v;
            } else if (t == 2) {
                int v = sc.nextInt();
                g0 -= v;
                g1 += v;
                sub1(m0, h0[0]);
                h0[0] += v;
                add1(m0, h0[0]);
            } else {
                int u = sc.nextInt() - 1;
                int v = sc.nextInt();
                if ((u & 1) == 0) {
                    sub1(m0, h0[u >> 1]);
                    h0[u >> 1] += v;
                    add1(m0, h0[u >> 1]);
                    if (u >> 1 < n1) {
                        sub1(m1, h1[u >> 1]);
                        h1[u >> 1] -= v;
                        add1(m1, h1[u >> 1]);
                    }
                } else {
                    sub1(m1, h1[u >> 1]);
                    h1[u >> 1] += v;
                    add1(m1, h1[u >> 1]);
                    if (u >> 1 < n0 - 1) {
                        sub1(m0, h0[(u >> 1) + 1]);
                        h0[(u >> 1) + 1] -= v;
                        add1(m0, h0[(u >> 1) + 1]);
                    }
                }
            }
            int ans0 = m0.getOrDefault(-g0, 0);
            int ans1 = m1.getOrDefault(-g1, 0);
            if (h0[0] == -g0) ans0--;
            pw.println(ans0 + ans1);
        }
    }

    static <Key> void add1(HashMap<Key, Integer> map, Key key) {
        map.putIfAbsent(key, 0);
        map.put(key, map.get(key) + 1);
    }

    static <Key> void sub1(HashMap<Key, Integer> map, Key key) {
        map.put(key, map.get(key) - 1);
    }

M - 筆塗り

問題文

HLD + 遅延セグ木をやるだけ. HLD は木上のパスに対するクエリを高速に処理できるので, 区間更新の遅延セグ木を載せれば OK.

実装例を貼ってもライブラリの中身がよく分からないと思うので省略. HLD や遅延セグ木の記事を見ると解けることが分かると思う.

N - マス目の穴埋め

問題文

二時間近くかかった, 大反省. 天才系の問題かと思って色々やったけど, 普通に DP だった...

まず, 入力の周りを 0 で囲っておく.

$$dp[i][S][T]=(i-1) 行目の 0,1 パターンが S で, i 行目の 0,1 パターンが T であるような場合の数$$

と DP テーブルを定義すればよい. 遷移は実装を見た方が早い気がする (希望があれば説明をちゃんと書きます).

    static final int H = 18;
    static final int W = 6;

    public static void solve(ExtendedScanner sc, FastPrintStream pw) {
        char[][] _b = sc.charArrays(H);
        int[][] b = new int[H + 2][W + 2];
        for (int i = 0; i < H; i++) {
            for (int j = 0; j < W; j++) {
                b[i + 1][j + 1] = _b[i][j] == '?' ? -1 : _b[i][j] - '0';
            }
        }
        long[][][] dp = new long[H + 2][1 << W][1 << W];
        dp[0][0][0] = 1;
        for (int i = 1; i <= H + 1; i++) {
            for (int s = 0; s < 1 << W; s++) {
                int[] z = toArray(s);
                for (int t = 0; t < 1 << W; t++) {
                    int[] y = toArray(t);
                    for (int u = 0; u < 1 << W; u++) {
                        int[] x = toArray(u);
                        boolean ok = true;
                        for (int j = 1; j <= W; j++) {
                            ok &= b[i][j] < 0 || b[i][j] == x[j];
                        }
                        if (!ok) continue;
                        if (valid(z, y, x)) {
                            dp[i][t][u] += dp[i - 1][s][t];
                        }
                    }
                }
            }
        }
        long ans = 0;
        for (int s = 0; s < 1 << W; s++) {
            ans += dp[H + 1][s][0];
        }
        pw.println(ans);
    }

    static int[] toArray(int s) {
        int[] x = new int[W + 2];
        for (int i = 0; i < W; i++) {
            x[i + 1] = (s >> i) & 1;
        }
        return x;
    }

    static boolean valid(int[] z, int[] y, int[] x) {
        int[][] xyz = new int[][]{z, y, x};
        for (int i = 1; i <= W; i++) {
            int[] c = new int[2];
            for (int d = 0; d < 4; d++) {
                int ni = i + Const.dx4[d];
                int nj = 1 + Const.dy4[d];
                c[xyz[nj][ni]]++;
            }
            if (++c[y[i]] < 3) {
                return false;
            }
        }
        return true;
    }

O - 宝箱

ボス問かと思いきや後半 4 問の中では一番簡単だった.

初めに $A_i$ は全てもらっておき, もらえなかった $A_i$ を損失として考えることにすると, 以下のようなグラフを構築して頂点 $0$ から頂点 $N$ までの最短経路問題を解けば OK.

f:id:suisen_kyopro:20201110222132p:plain

    public static void solve(ExtendedScanner sc, FastPrintStream pw) {
        int n = sc.nextInt();
        int m = sc.nextInt();
        long s = 0;
        long[] a = sc.longs(n);
        s = Arrays.stream(a).sum();
        Digraph<SimpleEdge> g = new Digraph<>(n + 1);
        for (int i = 0; i < m; i++) {
            int u = sc.nextInt();
            int v = sc.nextInt();
            long c = sc.nextLong();
            g.addEdge(new SimpleEdge(u - 1, v, c));
        }
        for (int i = 0; i < n; i++) {
            g.addEdge(new SimpleEdge(i, i + 1, a[i]));
            g.addEdge(new SimpleEdge(i + 1, i, 0));
        }
        long cost = new Dijkstra<>(g, 0).distance(n);
        pw.println(s - cost);
    }

第四回 アルゴリズム実技検定 (バーチャル参加) [G-K]

AtCoder による 第四回 アルゴリズム実技検定 の問題を解きました.

f:id:suisen_kyopro:20201110204452p:plain

G - 村整備

問題文

 N, M が小さいので, 全ての '#' に対して愚直にシミュレーションしても十分間に合う. BFS をしてもよかったが, UnionFind の方が個人的には楽.

    public static void solve(ExtendedScanner sc, FastPrintStream pw) {
        int n = sc.nextInt();
        int m = sc.nextInt();
        char[][] s = sc.charArrays(n);
        int sx = -1, sy = -1;
        int cnt = 1;
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                if (s[i][j] == '.') {
                    cnt++;
                    sx = i;
                    sy = j;
                }
            }
        }
        int ans = 0;
        for (int i = 0; i < n; i++) for (int j = 0; j < m; j++) {
            if (s[i][j] == '.') continue;
            boolean[][] g = GridUtil.build(s, '.');
            g[i][j] = true;
            DisjointSetUnion uf = new DisjointSetUnion(n * m);
            for (int r = 0; r < n; r++) for (int c = 0; c < m; c++) {
                if (!g[r][c]) continue;
                if (r < n - 1) {
                    if (g[r + 1][c]) {
                        uf.merge(r * m + c, (r + 1) * m + c);
                    }
                }
                if (c < m - 1) {
                    if (g[r][c + 1]) {
                        uf.merge(r * m + c, r * m + c + 1);
                    }
                }
            }
            if (cnt == uf.size(sx * m + sy)) {
                ans++;
            }
        }
        pw.println(ans);
    }

H - マス目のカット

問題文

正方形領域を固定すると, 一番多い数以外を変えるのが明らかに最適. 制約が小さいので,  n (プログラム中では t ) と正方形の角の座標を全て愚直に試しても間に合う.

    public static void solve(ExtendedScanner sc, FastPrintStream pw) {
        int n = sc.nextInt();
        int m = sc.nextInt();
        int k = sc.nextInt();
        int[][] s = new int[n][m];
        for (int i = 0; i < n; i++) {
            char[] row = sc.nextChars();
            for (int j = 0; j < m; j++) {
                s[i][j] = row[j] - '0';
            }
        }
        int ans = 1;
        for (int t = 1; t <= Math.min(n, m); t++) {
            for (int i = 0; i + t <= n; i++) for (int j = 0; j + t <= m; j++) {
                int[] cnt = new int[10];
                for (int r = i; r < i + t; r++) for (int c = j; c < j + t; c++) {
                    cnt[s[r][c]]++;
                }
                int max = IntArrays.max(cnt);
                int change = t * t - max;
                if (change <= k) {
                    ans = t;
                }
            }
        }
        pw.println(ans);
    }

I - ピザ

問題文

ここからは計算量削減を考えて解く必要がある.

 \displaystyle S=\sum _ {i=1}^N A _ i とする. 各  i に対して,  \displaystyle \sum _ {j=i}^k A _ {j \bmod n}\leq\frac{S}{2} を満たす最大の  k と, k+1 が求まればいい. 従って, 累積和+二分探索を用いると合計  O(N\log N) で解くことが出来る.

 A を二つ繋げた配列  B を用意して,  B に対して累積和を取ったり二分探索を行うと楽に実装できる.

    public static void solve(ExtendedScanner sc, FastPrintStream pw) {
        int n = sc.nextInt();
        long[] a = sc.longs(n);
        long sum = 0;
        for (long v : a) {
            sum += v;
        }
        long[] b = new long[2 * n];
        System.arraycopy(a, 0, b, 0, n);
        System.arraycopy(a, 0, b, n, n);
        long[] bsum = CumulativeSum.build(b);
        long ans = Const.LINF;
        for (int i = 0; i < n; i++) {
            int l = i - 1;
            int r = i + n;
            while (r - l > 1) {
                int m = (l + r) >> 1;
                long s = bsum[m + 1] - bsum[i];
                if (s << 1 <= sum) {
                    l = m;
                } else {
                    r = m;
                }
            }
            long x1 = bsum[l + 1] - bsum[i];
            long y1 = sum - x1;
            long x2 = bsum[r + 1] - bsum[i];
            long y2 = sum - x2;
            ans = Longs.min(ans, Math.abs(x1 - y1), Math.abs(x2 - y2));
        }
        pw.println(ans);
    }

J - ワープ

問題文

突然難しくなる.

実は AtCoder で上位互換 ARC061 C - すぬけ君の地下鉄旅行 が出題されており, 同じ考え方で解ける.

頂点  v_ {AB}, v _ {BC}, v _ {CA}, v _ {BA}, v _ {CB}, v _ {AC} を追加する. さらに, タイプ A の頂点からは頂点  v _ {AB} にコスト  X _ {AB} の有向辺, 頂点  v _ {AC} にコスト  X _ {AC} の有向辺, 頂点  v _ {BA} にコスト  0 の有向辺, 頂点  v _ {CA} にコスト  0 の有向辺を追加する. タイプ B, C についても同様に有向辺を追加する. 追加される辺の数は  4N 本である.

頂点および辺を追加したグラフ上で Dijkstra をすれば, 全体  O((N+M)\log N) でこの問題は解ける.

    public static void solve(ExtendedScanner sc, FastPrintStream pw) {
        int n = sc.nextInt();
        int m = sc.nextInt();
        long xab = sc.nextLong();
        long xca = sc.nextLong();
        long xbc = sc.nextLong();
        char[] s = sc.nextChars();
        // Digraph は有向グラフを表すクラス. SimpleEdge は始点, 終点, コストを持つ辺 のクラス
        Digraph<SimpleEdge> g = new Digraph<>(n + 6);
        for (int i = 0; i < m; i++) {
            int u = sc.nextInt() - 1;
            int v = sc.nextInt() - 1;
            long c = sc.nextLong();
            g.addEdge(new SimpleEdge(u, v, c));
            g.addEdge(new SimpleEdge(v, u, c));
        }
        int vab = n + 0, vbc = n + 1, vca = n + 2;
        int vba = n + 3, vcb = n + 4, vac = n + 5;
        for (int i = 0; i < n; i++) {
            if (s[i] == 'A') {
                g.addEdge(new SimpleEdge(i, vab, xab));
                g.addEdge(new SimpleEdge(vba, i, 0));

                g.addEdge(new SimpleEdge(vca, i, 0));
                g.addEdge(new SimpleEdge(i, vac, xca));
            } else if (s[i] == 'B') {
                g.addEdge(new SimpleEdge(i, vbc, xbc));
                g.addEdge(new SimpleEdge(vcb, i, 0));

                g.addEdge(new SimpleEdge(vab, i, 0));
                g.addEdge(new SimpleEdge(i, vba, xab));
            } else if (s[i] == 'C') {
                g.addEdge(new SimpleEdge(i, vcb, xbc));
                g.addEdge(new SimpleEdge(vbc, i, 0));

                g.addEdge(new SimpleEdge(vac, i, 0));
                g.addEdge(new SimpleEdge(i, vca, xca));
            }
        }
        var dij = new Dijkstra<>(g, 0);
        pw.println(dij.distance(n - 1));
    }

K - 転倒数

[問題文] (https://atcoder.jp/contests/past202010-open/tasks/past202010_k)

 A の転倒数を  I(A) と表すことにする. 列  X に対して列  Y を追加してできる列  Z=X+Y の転倒数  I(Z) がどうなるかを考えると,

$$ I(Z)=I(X)+I(Y)+\# \{ (i, j)\mid X_i\gt Y_j \} $$

が成り立つことが分かる. $I(X)$ は一つ前のクエリの答えであり, $I(Y)$ は前計算できる. 従って, $ \# \{ (i, j)\mid X_i\gt Y_j \} $ をどう高速に求めるかが問題となるが, ここで $1\leq X_i \leq 20$, $1\leq Y_j \leq 20$ の制約が活きる.

$$ \# \{ (i, j)\mid X_i\gt Y_j \}=\sum _ {u=1}^{20} \# \{ i\mid X_i=u \} \cdot \sum _ {v=1}^{u-1} \# \{ j\mid Y_j=u \}$$

と変形できるので, $C_u = \# \{ i\mid X_i=u \}$ を満たす配列 $C$ を持っておけばよい.

    public static void solve(ExtendedScanner sc, FastPrintStream pw) {
        int k = sc.nextInt();
        long[] inum = new long[k];
        int[][] a = new int[k][];
        int[][] cnta = new int[k][20];
        for (int i = 0; i < k; i++) {
            int n = sc.nextInt();
            a[i] = sc.ints(n, e -> e - 1);
            for (int v : a[i]) {
                cnta[i][v]++;
            }
            // 数列 A の転倒数は O(|A|log|A|) で求められる.
            inum[i] = InversionNumber.solveCompressed(a[i]);
        }
        int q = sc.nextInt();
        long[] cnt = new long[20];
        long ans = 0;
        while (q --> 0) {
            int i = sc.nextInt() - 1;
            ans += inum[i];
            for (int v = 0; v < 20; v++) {
                long s = cnta[i][v];
                long t = 0;
                for (int w = v + 1; w < 20; w++) {
                    t += cnt[w];
                }
                ans += s * (t % MOD) % MOD;
            }
            ans %= MOD;
            for (int v = 0; v < 20; v++) {
                cnt[v] += cnta[i][v];
            }
        }
        pw.println(ans);
    }

第四回 アルゴリズム実技検定 (バーチャル参加) [A-F]

AtCoder による 第四回 アルゴリズム実技検定 の問題を解きました.

f:id:suisen_kyopro:20201110204452p:plain

A - 中央値

問題文

ソートして  A,  B,  C と一致判定するのが楽?

    public static void solve(ExtendedScanner sc, FastPrintStream pw) {
        int a = sc.nextInt();
        int b = sc.nextInt();
        int c = sc.nextInt();
        int[] s = new int[]{a, b, c};
        Arrays.sort(s);
        if (s[1] == a) {
            pw.println('A');
        } else if (s[1] == b) {
            pw.println('B');
        } else {
            pw.println('C');
        }
    }

B - 電卓

問題文

 Y=0 の場合分けはいいとして,  100X/Y に小数点を挿入すると誤差が怖くない.

    public static void solve(ExtendedScanner sc, FastPrintStream pw) {
        int x = sc.nextInt();
        int y = sc.nextInt();
        if (y == 0) {
            pw.println("ERROR");
            return;
        } else {
            int res = x * 100 / y;
            int p0 = res % 10;
            int p1 = res / 10 % 10;
            int p2 = res / 100;
            pw.print(p2).print('.').print(p1).println(p0);
        }
    }

C - 隣接カウント

問題文

やるだけ. 配列外参照にだけ注意 (入力を受け取るときに周りを '.' で囲っても良い).

    public static void solve(ExtendedScanner sc, FastPrintStream pw) {
        int n = sc.nextInt();
        int m = sc.nextInt();
        char[][] s = sc.charArrays(n);
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                int c = s[i][j] == '#' ? 1 : 0;
                for (int d = 0; d < 8; d++) {
                    int ni = i + Const.dx8[d];
                    int nj = j + Const.dy8[d];
                    if (ni < 0 || ni >= n || nj < 0 || nj >= m) continue;
                    if (s[ni][nj] == '#') {
                        c++;
                    }
                }
                pw.print(c);
            }
            pw.println();
        }
    }

D - 分身

問題文

この辺りから問題のレベルが少し上がる.

まず,  x\geq (最左の忍者よりも左にあるマスの数),  y\geq (最右の忍者よりも右にあるマスの数) が必要. 加えて,  x+y\geq (隣り合う忍者の間の距離の最大値) が必要. 逆に, これで十分であることも言える.

    public static void solve(ExtendedScanner sc, FastPrintStream pw) {
        int n = sc.nextInt();
        char[] s = sc.nextChars();
        int max = 0;
        for (int i = 0; i < n;) {
            while (i < n && s[i] == '#') {
                i++;
            }
            if (i == n) break;
            int c = 0;
            while (i < n && s[i] == '.') {
                i++;
                c++;
            }
            max = Math.max(max, c);
        }
        int l = 0;
        while (s[l] == '.') {
            l++;
        }
        int r = n - 1;
        while (s[r] == '.') {
            r--;
        }
        r = n - 1 - r;
        if (l + r >= max) {
            pw.print(l).print(' ').println(r);
        } else {
            pw.print(l).print(' ').println(max - l);
        }
    }

E - アナグラム

問題文

全ての並び替えを全探索する. Java には next_permutation が標準にないのでライブラリとして持っておくべき.

    public static void solve(ExtendedScanner sc, FastPrintStream pw) {
        int n = sc.nextInt();
        char[] s = sc.nextChars();
        for (int[] p : Permutation.ofAscending(n)) {
            char[] t = new char[n];
            for (int i = 0; i < n; i++) {
                t[i] = s[p[i]];
            }
            if (Arrays.equals(s, t)) continue;
            CharArrays.reverse(t);
            if (Arrays.equals(s, t)) continue;
            CharArrays.reverse(t);
            pw.println(t);
            return;
        }
        pw.println("None");
    }

F - 構文解析

問題文

問題名で身構えてしまうけど, 中身を読んで安心.

まず, HashMap<String, Integer> を用意して文字列の頻度を数えておく.  c_i=頻度が i である文字列の数 となる  c を用意すると,  K 番目の文字列の頻度  x を求めることができる.

このとき,  c_x\gt 1 なら "AMBIGUOUS". そうでなければ頻度  x の文字列は一意に定まるので, Map を走査して見つければ OK.

    public static void solve(ExtendedScanner sc, FastPrintStream pw) {
        int n = sc.nextInt();
        int k = sc.nextInt();
        String[] s = sc.strings(n);
        HashMap<String, Integer> map = new HashMap<>();
        for (String t : s) {
            map.putIfAbsent(t, 0);
            map.put(t, map.get(t) + 1);
        }
        int[] c = new int[100001];
        for (var e : map.values()) {
            c[e]++;
        }
        int ans = 0;
        int cnt = 0;
        for (int i = 100000; i >= 1; i--) {
            cnt += c[i];
            if (cnt >= k) {
                if (c[i] > 1) {
                    pw.println("AMBIGUOUS");
                    return;
                } else {
                    ans = i;
                    break;
                }
            }
        }
        for (var e : map.entrySet()) {
            if (e.getValue() == ans) {
                pw.println(e.getKey());
                return;
            }
        }
    }

ACL Begginer Contest F - Heights and Pairs 解説

先日,ACLC 2 の問題が爆破された影響により代替コンテストとして ACL Begginer Contest が開催されました.私は 62:49 で全完し,過去最高の順位を獲得することが出来ました. *1

このコンテストは急に開催が決まった影響もあってか,公式の解説が AC コードのみとなっています.中身を読んでも何が何やらといった感じなので,とりあえずボス問の解説を書くことにしました.


ここから本題です.まず,解説する問題を以下に引用しておきます ( 引用元 : Heights and Pairs )

問題文

$2N$ 人の人 (  1 番から  2N 番まで ) がいます。 人 $i$ の身長は $h_i$ です。 以下の条件を満たすように、$N$ 個の人のペアを作る方法は何通りありますか? 答えを modulo $998,244,353$ で求めてください。

  • どの人もちょうど一つのペアに含まれる。
  • どのペアも、そのペアに属する二人の人の身長が異なる。

ある $p$ と $q$ に対し、人 $p$ と人 $q$ がペアになったかどうかが異なる場合、異なる方法であるとみなします。

制約

  • $1 \leq N \leq 50,000$
  • $1 \leq h_i \leq 100,000$
  • 入力は全て整数である。

まず,そのまま数えるのは難しそうなので包除原理の適用を考えます.$P_i$ を「 $2N$ 人の中から身長が同じであるペアを $i$ 組選ぶ場合の数」,$Q_i$ を「$2i$ 人を $i$ 個のペアに分ける場合の数」と定義すると,求める値は次のようになります.

\[ \sum _ {i=0}^ {N} (-1) ^ i \cdot P _ i \cdot Q _ {N-i} \]

一般に  \displaystyle Q _ M=\frac{{} _ {2M}C _ {2} \cdot {} _ {2M-2}C _ {2} \cdot\ \cdots\ \cdot {} _ {2}C _ {2}}{M!}=\frac{(2M)!}{2 ^ M \cdot M! } なので,あとは $P_i$ が高速に求まればこの問題は解けます.

身長の大小に意味は無く,等しいかのみが重要なので,身長に対するヒストグラムを考えます.つまり,身長 $j$ の人が何人いるかをカウントしておきます.これを $S_j$ とします.

例えば ${S_j}=[3, 1, 2, 4, 1, 0, 5]$ として,$i=4$ 組の身長が同じペアを選ぶとすると,以下のような選び方が考えられます.

f:id:suisen_kyopro:20200930030045p:plain
図 1 : ペアの数が $[1,0,0,1,0,0,2]$ であるような選び方の例

f:id:suisen_kyopro:20200930030048p:plain
図 2 : ペアの数が $[1,0,1,0,0,0,2]$ であるような選び方の例

ここで,身長が $j$ の人たちから $k$ 個のペアを作る場合の数を $A(j, k)$ とすると,$P_i$ は以下のように表せることが分かります.ただし,式中の $a _ j$ は身長が $j$ の人から選ぶペアの数で,図 1 における $[1,0,0,1,0,0,2]$ や図 2 における $[1,0,1,0,0,0,2]$ を指します.

\[ P _ i = \sum _ {\sum a _ j = i}\prod _ {j} A(j, a _ j) \]

$A(j, k)$ は $S_j$ 人の中から $2k$ 人を選んで $k$ 組のペアを作る場合の数なので  \displaystyle A(j,k)={} _ {S _ j}C _ {2k}\cdot Q _ k=\frac{S _ j!}{2 ^ k\cdot(S _ j-2k)!\cdot k!} です.

以上より,各 $j$ に対して形式的冪級数 $f_j$ を

\[ f _ j=\sum _ {k=0} ^ {\lfloor S _ j/2\rfloor}A(j, k)\cdot x ^ k=\sum _ {k=0} ^ {\lfloor S _ j/2\rfloor}\frac{S _ j!}{2 ^ k\cdot(S _ j-2k)!\cdot k!}\cdot x ^ k \]

で定義すると, \displaystyle \prod _ {j} f _ j=\sum _ {i=0} ^ {N}P _ i\cdot x ^ i となるので,左辺の積が高速に求まればよいです.

 \displaystyle \sum _ {j} S _ j=2N なので,$f_j$ の次数の和は高々 $N$ です.従って,FFT とマージテクを用いれば左辺は  O(N(\log N) ^ 2) で求めることが出来ます.

あとは,求めた $P _ i$ と $Q _ i$ を用いて冒頭に示した和を計算すれば良いです.階乗の逆元などを適切に前処理しておくことで,全体  O(N(\log N) ^ 2) でこの問題を解くことができます.

*1: 参加人数が少なかったというのもありますが...

ABC 169 F - Knapsack for All Subsets を愚直 DP からの変形のみで解く




先日の AtCoder Begginer Contest 169 に出場し,何とか全完したものの,F 問題で \(O(N^2S)\) の愚直 DP からの高速化に失敗して時間を溶かしてしまいました.そこで,復習を兼ねて \(O(N^2S)\) の DP を元にこの問題を解くことが出来ないかを考えてみました.問題文へはF - Knapsack for All Subsetsからどうぞ.



以降,数列 \(A\) は 1-indexed で考えます.



愚直 DP についてかなり簡単に説明します.\(A_{x_1}+A_{x_2}+\cdots+A_{x_j}=S\) の寄与を考えることで

$$ dp[i][j][k] := i \text{番目まで見て,}j \text{個使って和が}k \text{であるような選び方の個数} $$

で定義される DP テーブルにおける \(\displaystyle \sum_{j=0}^{N}dp[N][j][S]*2^{N-j}\) が答えとなることが分かります (\(S\gt 0\) なので \(dp[N][0][S]=0\) ですが,後に扱いやすいよう \(j=0\) を含めています).このテーブルは,\(dp[0][0][0]=1\) で初期化を行い,

$$ dp[i][j][k]=dp[i-1][j][k]+dp[i-1][j-1][k-A_i] \tag{1}$$

で更新することが出来ます.しかし,これは計算量が全体で \(O(N^2S)\) なので到底間に合いません.


高速化の手掛かりとして,必要な答えの形を思い出してみましょう.\(\displaystyle \sum_{j=0}^{N}dp[N][j][S]*2^{N-j}\) でした.この形を元に,次のように新しい DP テーブルの定義をしてみます.

$$ dp'[i][k]=\sum_{j=0}^{i}dp[i][j][k]*2^{i-j} \tag{2}$$

答えの形から \(N\rightarrow i, S\rightarrow k\) と書き換えただけです.次元が1つ減っていることにお気づきでしょうか.これは,新しいテーブルの更新が十分高速であればこの問題が解けることを表しています.

式 \((1)\) を用いて式 \((2)\) を計算してみると,

$$ \begin{eqnarray}dp'[i][k]&=&\sum_{j=0}^{i}dp[i][j][k]*2^{i-j} \\&=&\sum_{j=0}^{i}(dp[i-1][j][k])*2^{i-j}+\sum_{j=1}^{i}dp[i-1][j-1][k-A_i]*2^{i-j}\\&=&2*\sum_{j=0}^{i}(dp[i-1][j][k])*2^{(i-1)-j}+\sum_{j=0}^{i-1}dp[i-1][j][k-A_i]*2^{(i-1)-j}\\&=&2*dp'[i-1][k]+dp'[i-1][k-A_i]\end{eqnarray}$$

 となり,新しい DP テーブルは \(O(1)\) で更新できることが分かり,全体 \(O(NS)\) の計算量となるのでこの問題は解けました.

 

以下,Java で記述したコードを掲載します.メモリ消費量を抑えるため,in-place な更新を行っています.この場合,\(k\) のイテレーション方向に気を付けてください.

import java.util.Scanner;

public class Main {
    public static void main(String[] args) {
        final int MOD = 998244353;
        final var sc = new Scanner(System.in);
        final int N = Integer.parseInt(sc.next());
        final int S = Integer.parseInt(sc.next());
        final int[] A = new int[N + 1];
        for (int i = 1; i <= N; i++) { // 1-indexed
            A[i] = Integer.parseInt(sc.next());
        }
        sc.close();
        final int[] dp = new int[S + 1];
        dp[0] = 1;
        for (int i = 1; i<= N; i++) {
            for (int k = S; k >= 0; k--) { // descending
                dp[k] += dp[k];
                if (dp[k] >= MOD) dp[k] -= MOD;
                if (k - A[i] >= 0) {
                    dp[k] += dp[k - A[i]];
                    if (dp[k] >= MOD) dp[k] -= MOD;
                }
            }
        }
        System.out.println(dp[S]);
    }
}

ちなみに,私は純粋な dp 解法は思い浮かばなかったので形式的冪級数を用いて解きました.\(\left[x^S\right]\prod_{i=1}^N\left(2+x^{A_i}\right)\) が答えとなるのですが,こちらでも同様の漸化式 (dp の更新式) が立ちます.この解法は maspy 様が[AtCoder 参加感想] 2020/05/31:ABC 169 | maspyのHPにて大変分かりやすく解説されていますので,是非読んでみてください.