ABC332 勉強メモ

先日開催された、AtCoder Beginner Contest 332の勉強メモです。
問題への言及がありますので、今後解く予定の方はご注意ください。
使用言語はPython・提出はPyPy3です。

公式解説

A～C

解法
いずれも愚直に処理を行えばよいです。
Cは$O(N^2)$になりますが間に合います。
提出例: ABC332A, ABC332B, ABC332C O(N^2), ABC332C O(N)

D Swapping Puzzle

問題文はこちら

解法
バブルソート(隣り合う2要素を昇順交換するソート法)の逆操作を考えると、隣接する項の交換を繰り返すことで任意の順列を得ることができます。
なので、たとえば行の操作で各行を自由に並び替えることができます。
さらに行の操作と列の操作は互いに干渉しないので、結局各行各列を自由に並び替えることができます。
操作後の元々の行・元々の列の組み合わせを全探索し、各組に対してグリッドBと一致するか判定しましょう。
ここまでの計算量は$O(H! × W! × HW)$です。

最小操作回数を考えます。
例として行を考えますが、列も同様に考えてかまいません。

元々の行の位置は昇順にソートされた順列 $S = [1, 2, ･･･ , H]$ です。
これを隣接する項の交換で並び替えて、新しい順列 $P =[P_1, P_2, ･･･ , P_H]$を得たとします。
S→Pの並び替えの最小回数は、逆操作であるP→Sの並び替えの最小回数と一致します。従って、Pをバブルソートする最小回数を答えればよいです。
これは転倒数を考えれば$O(H^2)$で、BFSなら$O(H × H!)$で計算できます。
特に、BFSで行う場合は前計算しておいた方がよいでしょう。
提出例(BFS)

転倒数について

転倒数とは、$ i < j $なのに$ A_i > A_j $となる$(i,j)$の組数のことです。
たとえば $[5, 1, 2, 4, 3]$の転倒数は5、$[1,2,3,4,5]$の転倒数は0です。

転倒数はバブルソートの最小操作回数であることが知られています。(終)

転倒数とバブルソート 個人用メモ

証明が見つからなかったので、メモを残しておきます。

1. 隣接2項のswapで、転倒数は±1だけ変化する。
   順列$P =[P_1, P_2, ･･･ , P_i, P_j, ･･･ ,P_N]$を考える。
   隣接2項の$P_i,P_j$のswapにより、$P_i$と$P_j$の大小関係だけが変化し、それ以外の要素の大小関係は変化しない。
   従って、元々が$P_i < P_j$ならばswapで転倒数が1増加し、$P_i > P_j$ならばswapで転倒数が1減少する。
2. バブルソートの操作回数の下界は転倒数の値であり、これは達成可能である。
   一回の操作で転倒数は1しか減少しないので、最低でも転倒数の値だけswapが必要である。
   そして転倒数が1以上であれば、「隣接2項が転倒している」場所がかならず存在する。
   (帰納法で証明する。
   転倒数が1以上なのに、離れた2項しか転倒していない状況があると仮定する。
   この場合、隣接2項はすべて昇順ソートされているので、数列全体も昇順ソートされていることになる。
   すると転倒数は0となり、矛盾する)
   なので、転倒した隣接項を選びswapし続けることで、転倒数の回数でソートが完了する。(終)

BFSでの数え方

順列を頂点に見立てて、01-BFSで最短経路を数えればよいです。
これはダイクストラ法の応用で、辺のコストが常に1であることを利用してheapqではなくdequeやlistで行うものを指します。(heapqでもできます)
順列の要素数をNとすると、頂点数が$N!$、遷移の辺数が各頂点ごとに$N-1$個なので計算量は$O(N × N!)$となります。
詳細は提出例をご確認ください。(終)

E Lucky bag

問題文はこちら

解法
分散の計算式のうち、$\bar{x}$は実は定数です。
$V = \frac{1}{D} \sum_{i=1}^D (x_i - \bar{x})^2$ ･･･(1)

全福袋の重さの平均は、常に全グッズの重さの平均と等しいからです。
すなわち $\bar{x} = \frac{1}{D} \sum_{i=1}^D x_i = \frac{1}{D} \sum_{i=1}^N W_i$です。
なので結局$\sum_{i=1}^D (x_i - \bar{x})^2$を最小化すればよいことが分かります。

DP[d][S]: 既に決めた福袋がd袋、採用決定したグッズの集合がS(二進数表記)のときのΣの最小値
として、$O(3^N)$の部分集合列挙を行いましょう。
計算量は$O(D×3^N)$です。

なお、高速化すると$O(3^N)$になるようです。
詳細はNachiaさんのユーザ解説をご確認ください。

部分集合列挙の計算量解析や、「二乗の平均 - 平均の二乗」の別解の桁落ち、Pythonの高速化法などは長くなるので折りたたみ表示にしてあります。
提出例: O(D*3^N), O(D*3^N) 別解, O(3^N)

部分集合列挙

集合Sの部分集合Tは以下の手順で列挙できます。

部分集合列挙

T = S
while T >= 0:
    T = S & T
    #部分集合の処理
    T -= 1

計算量解析を行います。全体集合をUとします。
$\sum_{S⊂U} 2^{|S|}$ ･･･ Sの部分集合は$2^{|S|}$個ある
$ = \sum_{k=0}^N {} _ N \mathrm {C} _ k × 2^k$
$ = \sum_{k=0}^N {}_N \mathrm{C}_k × 2^k × 1^{N-k}$
$ = (2+1)^N = 3^N$
となるので、部分集合列挙の計算量は$O(3^N)$です。(終)

別解と桁落ち

元々の分散の定義を再掲します。
$V = \frac{1}{D} \sum_{i=1}^D (x_i - \bar{x})^2$ ･･･(1)

有名な公式ですが、分散には「二乗の平均 - 平均の二乗」という公式があります。
$V = \bar{x_i ^2} - (\bar{x})^2 = \frac{1}{D} \sum_{i=1}^D x_i^2 - (\bar{x})^2$ ･･･(2)

証明

(1)を展開すると
$V = \frac{1}{D} \sum_{i=1}^D (x_i^2 - 2 \bar{x} x_i + (\bar{x})^2 )$
$\bar{x_i^2}= \frac{1}{D} \sum_{i=1}^D x_i^2 $、平均の定義より$\frac{1}{D} \sum_{i=1}^D x_i = \bar{x}$を適応すると
$V = \bar{x_i^2} - 2\bar{x} × \bar{x} + (\bar{x})^2 = \bar{x_i^2} - (\bar{x})^2$
と変形できた。(終)

$(\bar{x})^2$は定数なので、結局$\sum_{i=1}^D x_i^2$を最小化すればよく、
DP[d][S]: d袋目まで決めたとき、採用済グッズ集合がSのときのΣの最小値
と定義すれば、$O(D × 3^N)$の部分集合列挙で解けます。
しかし これには罠があり、このまま計算を行うと桁落ちして巨大な誤差が発生してしまいます。

Pythonは整数に関しては桁数制限なく使用できます(厳密には$2^{63}$以上は多倍長整数となり、処理速度が著しく低下します)。
一方でfloat型は仮数部52bitしか保持できず、有効数字が16桁程度しかありません。
なのでfloat型で巨大な値を扱うと切り捨て処理が発生してしまいます。
特に 巨大数 - 巨大数 の場合、上位桁がなくなることでさらに有効数字が目減りしてしまいます。これを桁落ちといいます。

桁落ちの例

A = 1234567890123456789  #19桁の整数

#巨大整数の丸め誤差
B = float(A)
print(B)  #1.2345678901234568e+18  有効数字16桁程度に丸めてしまう

B = int(B)
print(B)  #1234567890123456768  intに戻しても誤差ぶんずれる
print(A-B)  #21

#桁落ち
A = 1234567890.0
B = A + 0.12345678901234567
print(B)  #1234567890.1234567  有効数字16桁程度に丸めてしまう
print(B-A)  #0.12345671653747559  有効数字7桁

本問の場合、巨大整数の有効数字の丸めが深刻です。
例として、分散が$1/4$となる2ケースを見てみます。

sample imput 1 (sample output 1: 0.25)

2 2
10 11

$W_i$が小さい場合、「二乗の平均 - 平均の二乗」でもACできます。
しかし平均値が巨大になった次のケースではどうでしょう。

sample input 2 (sample output 2: 0.25)

2 2
90000010 90000011

$\frac{DP[D][{1,2,･･･,N}]}{D} - (\frac{ΣW}{D})^2$だと0.25ではなく0.0を出力するはずです。
DP[D][all]の値が約$1.6×10^{16}$と巨大になってしまい、有効数字16桁程度のfloat型にすると小数点以下どころか整数部分すら精度が保てません。
当然ジャッジの許容誤差の$10^{-6}$を全く満たせずWAとなってしまいます。

より顕著な例はこちらです。「二乗の平均 - 平均の二乗」だと22.22... ではなく0.0を出力します。

sample imput 3 (sample output 3: 22.2222222 ... )

15 3
100000000 100000000 100000000 100000000 100000000 100000000 100000000 100000000 100000000 100000000 100000000 100000000 100000000 100000000 99999990

本問は$ΣW_i$が整数なので、式変形を行うことで最後の除算以外を整数型で扱えるようになります。
あるいはFractionやDecimalといった高精度演算を一時的に併用することでもACが可能です。
import時間を含めて+50ms程度で処理できます。(終)

Python 高速化法

本問はPythonだと実行時間制限がかなり厳しいです。
hirayuu_Atさんのユーザ解説とも重複しますが、高速化の工夫を示します。

1. DPの一次元化
DP[d][S]のような二次元配列は遅いので、一次元化しましょう。
DP[d * (1<<N) + S] ないし DP[(d << N) + S] のように対応すればよいです。
趣旨が異なりますが、処理順を工夫すればDP[S]のように添字を減らす高速化も可能です。詳細は省略します。

2. すべてfloat型で計算する
int→floatの変換コストが損なので、可能な限りfloat型にしましょう。

3. 枝狩りを行う
直感的に、空の福袋ができることは損だとわかります。
なので空の福袋が生じうる組み合わせの処理を省くことで高速化できます。

証明

要素数Nの集合Uを最小分散を満たすようD個に分割し、要素数の昇順に並び替える。
$|S_1| \leq |S_2| \leq ･･･ \leq |S_D|, sum(|S_i|) = |U| = N \geq D$
となる。また、平均を$\bar{x}$で表す。

このなかに空集合があると仮定し、帰納法で矛盾を示す。
昇順ソートしたので$|S_1| = 0, |S_D| \geq 2$である。
$S_D$の総和を$X$、要素のひとつを$d$として、これを$S_1$に移動させる。
移動前の分散は
$D × V_{before} = (X- \bar{x})^2 + (0- \bar{x})^2 + ･･･$
移動後の分散は
$D × V_{after} = (X-d- \bar{x})^2 + (0+d- \bar{x})^2 + ･･･$
$ = (X- \bar{x})^2 - 2d(X - \bar{x}) + d^2 + (0- \bar{x})^2 + 2d(0- \bar{x}) + d^2 + ･･･$
となる(･･･の部分は$S_2$から$S_{D-1}$の部分であり、同一である)。
この差を考えると
$D × (V_{after} - V_{before})$
$ = - 2d(X - \bar{x}) + d^2 + 2d(0- \bar{x}) + d^2$
$ = - 2d(X - d)$
となる。$D,d,X-d > 0$なので、$V_{after} - V_{before} < 0$である。
これは最小分散の仮定と矛盾する。(終)

具体的には、集合Sの要素数と現在の分割数を確認します。
分割しすぎだったり、要素を決めすぎて今後空袋が生じるようであれば処理を行わないことで高速化が可能です。
この枝狩りは特にDが大きい場合に効果的です。

計算回数の確認コード

処理回数の確認

#高速bit_count
def popcnt(N):
    N = (N & 0x5555555555555555) + ((N >>  1) & 0x5555555555555555)
    N = (N & 0x3333333333333333) + ((N >>  2) & 0x3333333333333333)
    N = (N & 0x0F0F0F0F0F0F0F0F) + ((N >>  4) & 0x0F0F0F0F0F0F0F0F)
    N = (N & 0x00FF00FF00FF00FF) + ((N >>  8) & 0x00FF00FF00FF00FF)
    N = (N & 0x0000FFFF0000FFFF) + ((N >> 16) & 0x0000FFFF0000FFFF)
    N = (N & 0x00000000FFFFFFFF) + ((N >> 32) & 0x00000000FFFFFFFF)
    return N

N = 15
for D in range(1,N+1):
    cnt1 = cnt2 = 0
    for d in range(D):
        for S in range(1<<N):
            cnt1 += 2 ** popcnt(S)
            #枝狩りの処理
            if popcnt(S) >= d+1 and N - popcnt(S) >= D - (d+1):
                cnt2 += 2 ** popcnt(S)
    print(D,cnt1,cnt2,cnt2/cnt1)

枝狩りなしは$D×3^N$で、処理回数の最大値はD = 15の時で$2.1×10^8$です。
一方、枝狩りありだと最大でもD = 8の$8.2×10^7$で抑えられると分かります。(終)

参考: 実行時間比較
枝狩りなし
二次元DP float型 枝狩りなし: 2103ms(TLE)
一次元DP int型 枝狩りなし: (TLE3)
一次元DP float型 枝狩りなし: 1188ms

枝狩りあり
二次元DP int型 枝狩りあり: 1191ms
二次元DP float型 枝狩りあり: 629ms
一次元DP int型 枝狩りあり: 1081ms
一次元DP float型 枝狩りあり: 485ms

DP一次元化・float型の統一・枝狩り いずれも3割以上の高速化ができます。(終)

F Random Update Query

問題文はこちら

解法
遅延評価セグメントツリーを使用します。
これは区間作用値の反映を遅らせることで、区間更新クエリを$O(logN)$で行えるようにしたセグメントツリーです。
実装の詳細はPAST上級本をご確認ください。

なお、本問のように区間更新・一点取得の場合、N個の遅延とN個のノードがあれば十分です。
なので双対セグメントツリーを改造したような実装をすることで、空間計算量を半減できます。

双対セグメントツリーについてはtatyamさんの解説が詳細です。
本問は作用素の交換法則が成立しないので、遅延反映ができる実装にする必要があります。(双対セグ木と遅延セグ木の中間の実装が必要です)

クエリ順読み
公式解説と同様です。
ARC008D タコヤキオイシクナールを区間更新で行えばよいです。
$d = L_i - R_i + 1$とすると、要素の変更により期待値は$E← \frac{d-1}{d}E + \frac{1}{d}X_i$になります。
また、$x←a_1 x+b_1$の変換後に$x←a_2 x+b_2$の変換を行う操作は
$x←a_1 a_2 x + ( a_2 b_1 + b_2 )$と表すことができます。
なので遅延に一次関数を乗せ、遅延合成を上記の一次関数の合成と見立てる実装を行えばよいです。
なお、$998244353 \leq A_i \leq 10^9$がコーナーケースになりやすいので注意してください。
提出例(双対セグ木寄り, 810ms)

クエリ逆読み
各$A_i$に対して、答えへの寄与を考えます。
$A_i$の値は総変更回数を問わず最後に変更された値になります。
具体的には$d = L_i - R_i + 1$とすると、$1/d$の確率で値が$X_i$になります。
これをクエリの逆順で見ることで、値ごとの変更確率が順に求まります。

ノードと遅延を(一回も変更されない確率, 答えへの寄与)のペアで定義します。
$A_i$の期待値は これまで一回も変更されなかった確率 × 今回変更される確率 × $X_i$ の総和で求めることができます。
提出例(1977ms)

$A_i$順に処理
各$A_i$に対して、ARC008D タコヤキオイシクナールを解きます。
セグメントツリーにクエリごとの作用値を乗せ、各$A_i$ごとに機能しているクエリの全幅作用値を取得すればよいです。
計算量は$O((N+M)logM)$です。
提出例(1120ms)

おわりに

おわりです。
コーナーケースがおおくてむずかしいなとおもいました。