しょぼんブログ

数学の色々とか様々とか

連載 しょぼんコーダー 番外編1 平方分割を実装してみた

連載 しょぼんコーダー

バグループ

聴いてください、バグループ
バグってない平方分割を知らない
バグってない平方分割が実装できない
バグってない平方分割を知らない
バグってない平方分割が実装できないよ

蟻本 pp169-170 を参考にして, 自分なりに平方分割を実装してみた.

Library Checker - Static RMQ Library Checker

区間最小値はセグメント木でも行けるが, 平方分割でも行ける.

コード

#include<bits/stdc++.h>
using namespace std;

int main(){
    int n,q;
    cin >> n >> q;
    int b = int(pow(n, 0.5) + 1);
    vector<int> a(n);
    for (int i=0; i<n; i++){
        cin >> a[i];
    }

    vector<int> dat(b, 1e9);
    for (int i=0; i<b; i++){
        for (int j=min(n, i*b); j<min(n, (i+1)*b); j++){
            dat[i] = min(dat[i], a[j]);
        }
    }

    for (int num=0; num<q; num++){
        int l, r;
        cin >> l >> r;
        int tl = l+b-(l%b), tr = r-(r%b), ans = 1e9;
        for (int i=l; i<min(r,tl); i++){
            ans = min(ans, a[i]);
        }
        for (int i=tl/b; i<tr/b; i++){
            ans = min(ans, dat[i]);
        }
        for (int i=max(tl,tr); i<r; i++){
            ans = min(ans, a[i]);
        }
        cout << ans << endl;
    }
}

説明

細かくみていく.

b

int b = int(pow(n, 0.5) + 1);

bはバケットのサイズでもあるが, バケットの数でもある. +1でうまい具合に調整してある.

前処理

vector<int> dat(b, 1e9);
for (int i=0; i<b; i++){
    for (int j=min(n, i*b); j<min(n, (i+1)*b); j++){
        dat[i] = min(dat[i], a[j]);
    }
}

前処理はこんな感じ. 各バケットについて計算している. ここは O(N)

（追記 2022/09/03）for文だと必要ないが, sortだとバケットの左端jを i*b としてしまうと配列外にアクセスしてしまう. これを防ぐために, ここの時点で for int j min(n, i*b) と書いている.

RMQ

for (int num=0; num<q; num++){
    int l, r;
    cin >> l >> r;
    int tl = l+b-(l%b), tr = r-(r%b), ans = 1e9;

    // (1)
    for (int i=l; i<min(r,tl); i++){
        ans = min(ans, a[i]);
    }

    // (2)
    for (int i=tl/b; i<tr/b; i++){
        ans = min(ans, dat[i]);
    }

    // (3)
    for (int i=max(tl,tr); i<r; i++){
        ans = min(ans, a[i]);
    }
    
    cout << ans << endl;
}

ここはクエリごとの処理. ここは O(Q\sqrt{N}). 蟻本と違うのは, 左から順番に見ていることと, while文を使っていないこと.

tl は, l から右に行って最も近いブロックを示し, tr は, r から左に行って最も近いブロックを示す.

蟻本でも tl, tr を使っているが, ちょっとだけ違う. (1) と (3) でそれぞれforの区間にmin, maxを使っている. これは下図のような状態を防ぐためである. 余計なものまで計算してしまい正しい結果が出ない.

左から順番に見ているのは, 左からやっていく場合に使いたいからである. 蟻本のような実装だと, tlは左から右へ移動するが, trは右から左に移動する. これでは操作を左からやっていきたいときに困る.

while文をfor文にしたのは, while文だとちょっと遅いからである（定数倍高速化）. while文で実装したときは1400msだったが, for文で実装したら858msだった.

while文
https://judge.yosupo.jp/submission/101976

for文
https://judge.yosupo.jp/submission/101979

ドドスコードの出力回数の期待値

このプログラムにおいて, 配列{"ドド","スコ"}からランダムに要素を標準出力する回数の期待値を求めました.
また, そのあと「ドド」の出る確率を pとしたときの出力回数の期待値を求め,
「ドド」の出る確率を何にしたら期待値を最小にできるかを調べました.

続きを読む

連載 しょぼんコーダー 4 先頭固定テクニック

連載 しょぼんコーダー

第4回では, 先頭固定テクニックとぼくが言う, 主に定数個入力数え上げで平等性が高い状況のときに使えるDPについて書きます.

問題1

例題4-1 N個の区別する玉を K個の区別しない箱に入れる方法の数はいくつですか. ただし, 各箱には必ず 1個以上の玉を入れるものとします.

つまり, 第2種Stirling数 S(N, K)を求めてください.
ただし, 答えは 998244353で割った余りで出力してください.

時間制限：2秒、メモリ制限：1024 MB
制約

  • 1 \leq K \leq N \leq 300

入力

N K

考察

dp[i][j]を i個の区別する玉を j個の区別しない箱に各箱1個以上入れる方法とします.
次の O(NK)のDPが知られています.

mod = 998244353
n,k = map(int,input().split())
dp = [[0] * (k+1) for j in range(n+1)]
dp[0][0] = 1

for i in range(n):
	for j in range(k):
		dp[i+1][j+1] = (j+1) * dp[i][j+1] + dp[i][j]
		dp[i+1][j+1] %= mod

print(dp[n][k])

今回は, 先頭固定テクニックで解いてみます. dp[i][j]について考えます. ある1つの箱について考えます. その箱は m個の玉を含んでいるものとします. このとき, ある玉の組み合わせを固定すると, その組み合わせを含むものは, dp[i-m][j-1] 通りとなります.
それはそうなのですが, 実際にすべての組み合わせについてこれを足し合わせると（つまり \binom{i}{m}を掛けると）, 次の画像のように同じ結果になるものがダブルカウントしてしまうことがあります.

ダブルカウントしない方法として, 「今見ている配列のうち, 先頭を固定する」ということが考えられます. 先頭を固定すれば, 後続は \binom{i-1}{m-1}通りですが, 先頭が固定されていることにより遷移先で先頭が証人になっているので, 後々になってダブルカウントすることはありません. これが先頭固定テクニックです.
dp[i][j]は, 1\leq m\leq iのすべてについて dp[i-m][j-1] \times \binom{i-1}{m-1}の総和となります. x!, 1/y!を O(N + \log \text{MOD})で前計算することにより, このDPは O(N^2K)となります.

ACコード

mod = 998244353
N = 400
fact = [1]*(N+1)
factinv = [1]*(N+1)

for i in range(2, N+1):
	fact[i] = fact[i-1] * i % mod

factinv[-1] = pow(fact[-1], mod-2, mod)
for i in range(N-1, 1, -1):
	factinv[i] = factinv[i+1] * (i+1) % mod

def cmb(a, b):
	if (a < b) or (b < 0):
		return 0
	return fact[a] * factinv[b] % mod * factinv[a-b] % mod

n,k = map(int,input().split())
dp = [[0] * (k+1) for j in range(n+1)]
dp[0][0] = 1

for i in range(1,n+1):
	for j in range(1,k+1):
		for m in range(1,i+1):
			dp[i][j] += dp[i-m][j-1] * cmb(i-1, m-1) % mod
			dp[i][j] %= mod

print(dp[n][k])

時間計算量はあまりよくないけれど, 有名な遷移とはまた異なる遷移ができました.
余談ですが, \sum^{i}_{m=1} \binom{i-1}{m-1} dp[i-m][j-1] = \sum^{i-1}_{m=0} \binom{i-1}{m} dp[m][j-1]と変形できます. dp[i][j] = \sum^{i-1}_{m=0} \binom{i-1}{m} dp[m][j-1] = \sum^{i-2}_{m=0} \binom{i-1}{m} dp[m][j-1] + dp[i-1][j-1]に注目すると, 次の式が得られます.
\sum^{i}_{m=0} \left\{\binom{i+1}{m} - (j+1) \binom{i}{m}\right\} S(m, j) = 0 （ただし i, jは任意の非負整数）

問題2

例題4-2 ACM-ICPC Japan Alumni Group Contest 2019 2日目 H
問題文は長いので上のリンクを見てください.

時間制限：2秒、メモリ制限：256 MB
制約

  • 1 \leq N \leq 10^5

考察

長さ Nの順列 pのスコア f(p)の総和を A[N] = \sum f(p), スコアの2乗の総和を B[N] = \sum f(p)^2とします.
このとき, スコアの平均値は a=\frac{1}{N!} A[N], 分散は \frac{1}{N!} \sum (f(p)^2 - 2 f(p) a + a^2) = \frac{1}{N!} (B[N] - 2 A[N] a + a^2 N!) です.
よって, A[N], B[N]を求めればokです.

順列 pを対称群 i \mapsto p(i)の元だと見ると, 「サイクル森の1つの連結成分の数」は, 「巡回置換の長さ」（つまり G(p)は閉路からのみなる）と確認することができます.
その閉路に注目します. A[i] を求めてみましょう. mを, 「先頭に含まれるグループの連結成分の数」とします. このとき, そのようなグループの決め方は \binom{i-1}{m-1}個, さらにグループの中の閉路の作り方が (m-1)!個あります. （ m次対称群の長さ mの巡回置換は (m-1)!個あるので）
A[i-m]から乗算されるので, 閉路の大きさ mも掛けて, m A[i-m] \times \binom{i-1}{m-1} \times (m-1)! \times m が A[i]に加算されます.
m=1,2,\cdots, iを考慮して, A[i] = \sum^{i}_{m=1} A[i-m] \times \binom{i-1}{m-1} \times (m-1)! \times mとなります.

このように, 先頭固定テクニックは対称群（置換）と相性が良いです.

同様に B[i] = \sum^{i}_{m=1} B[i-m] \times \binom{i-1}{m-1} \times (m-1)! \times m^2です.
実はこれで O(N^2)コードができました. しかし間に合わないので工夫をします.

式変形をすると,
A[i] = \sum^{i}_{m=1} A[i-m] \times \binom{i-1}{m-1} \times (m-1)! \times m
= \sum^{i}_{m=1} A[i-m] \times (i-1)! / (m-1)! / (i-m)! \times (m-1)! \times m
= (i-1)! \sum^{i}_{m=1} A[i-m] / (i-m)! \times m
で, 順序を変えると
A[i] = (i-1)! \sum^{i-1}_{m=0} A[m] / m! \times (i-m)
= (i-1)! \left\{i \times \sum^{i-1}_{m=0} A[m] / m! - \sum^{i-1}_{m=0} A[m] / m! \times m\right\}

同様に
A[i] = (i-1)! \sum^{i-1}_{m=0} B[m] / m! \times (i-m)^2
= (i-1)! \left\{i^2 \times \sum^{i-1}_{m=0} B[m] / m! - 2i \sum^{i-1}_{m=0} B[m] / m! \times m + \sum^{i-1}_{m=0} B[m] / m! \times m^2 \right\}
です.

\sum^{i-1}_{m=0} A[m] / m!
\sum^{i-1}_{m=0} A[m] / m! \times m
\sum^{i-1}_{m=0} B[m] / m!
\sum^{i-1}_{m=0} B[m] / m! \times m
\sum^{i-1}_{m=0} B[m] / m! \times m^2
を持つ変数を管理すれば, この5つは, iを1つ増やすごとに O(1)で更新できます. よって累積和的に各 A[i], B[i]を O(1)で更新でき, 全体で O(N + \log \text{MOD})で解けます.

ACコード

mod = 10**9 + 7
N = 100256
fact = [1]*(N+1)
factinv = [1]*(N+1)

for i in range(2, N+1):
	fact[i] = fact[i-1] * i % mod

factinv[-1] = pow(fact[-1], mod-2, mod)
for i in range(N-1, 1, -1):
	factinv[i] = factinv[i+1] * (i+1) % mod

def cmb(a, b):
	if (a < b) or (b < 0):
		return 0
	return fact[a] * factinv[b] % mod * factinv[a-b] % mod

n = int(input())
a0 = 1
a1 = 0
b0 = 1
b1 = 0
b2 = 0
t = 0
s = 0
for i in range(1, n+1):
	t = fact[i-1] * (i*a0 - a1) % mod
	s = fact[i-1] * (i*i*b0 - 2*i*b1 + b2) % mod
	a0 = (a0 + t*factinv[i]) % mod
	a1 = (a1 + t*i*factinv[i]) % mod
	b0 = (b0 + s*factinv[i]) % mod
	b1 = (b1 + s*i*factinv[i]) % mod
	b2 = (b2 + s*i*i*factinv[i]) % mod

a = t*factinv[n]%mod
print((s - 2*t*a + a*a*fact[n])*factinv[n]%mod)

練習問題

ネタバレをあまりしたくない主義なので, リンクを隠してあります. わかりにくくてすみません.
1番目のリンク (難易度:AtCoder 橙)
2番目のリンク (難易度:yukicoder ★4)

ABC255のA~Fを解くまでの思考過程

AtCoder参加記

ぼくがABC255のA～F (+コンテスト後にEx)を解くまでの思考過程をまとめてみました.

atcoder.jp

A : 1分39秒
B : 5分04秒
C : 10分29秒 (1ペナ)
D : 13分18秒
E : 38分26秒
F : 64分47秒

で, 163位でした

本当に素の思考過程なので, 自分流の言葉とかふんわりとした気持ちを使ってます
あと, ネタバレを多く含むので展開式になってます