Nyaanの日記

競プロ用(精進記録など) 乱文失礼します

CodeChef April Challenge 2020 Div.1

15位…と思ったら今見たら14位になっている(マラソンのrejudge中らしい、今後も変わりそう) 順位表
Mod 1e9+7と998...ばかりで個人的にかなり楽しめたセットでした

STRNO,SQRDSUB

  • 簡単枠 同値変形をすると解ける

ANSLEAK

  • ラソン枠 高速化+遷移を工夫しただけの山登りで94点
    • かなり雑な出来だと思ったが、やり込んだ人が少なかったからかDiv1では8番目の成績
  • 時間内に山を登り切れないと踏んで焼きなましに書き換えなかったが、高速化の甲斐あって最後の方は 10 ^ 6 回ループが回っていたので焼きなましにすべきだったね…

REBIT

PPDIV

  • 数学問 メビウス関数を使うと包除の部分は取り除けて   N + \sum _ { i = 2 } ^ {60} - \mu ( i ) \sum _ {2 \leq n \leq  \sqrt [ i ] {N}}  n ^ i \lfloor \frac{N}{n ^ i} \rfloor が答えになる
    • これは愚直に計算すると O(N^ {\frac{1}{2}} ) \  (N \leq 10 ^ {18})
  • 式の見た目がProject Eulerの知識を要求しそうでやばい
  • 冷静に考えてPE典型が800人に解かれるわけないだろ!となり式を睨むと、floorの部分の値で場合分けするいつものやつで \mathrm{O}( N ^ {\frac{1}{3} } ) に落ちて終了

FCTRE

  • 問題文を読むと、木上のMoをやるだけだとわかる
  • ところで木上のMoを書いたことがありますか?僕はありません
    • それどころか長さ2Nのオイラーツアーすら書いたことがなく…(本当に橙?笑)
  • うしさんの記事をガン見、parityを持つ実装を参考にしつつMoを書いてAC(うしさんに、感謝!)
  • 提出

LLLGRAPH

  • Line Graph何もわからん
  • Edmondのアルゴリズムを貼って部分点を通してあとは撤退…

TRIPWAYS

  • これ面白かった この記事の本題です
  • ↓問題へのリンク

www.codechef.com

  • 形式的冪級数で殴る(い つ も の)
    • 次数と日数を対応させた母関数を考えると、スタート地点が \frac{1}{1-L_1 x} 、遷移が \frac{x}{1-L_i x}になる
    • DPの部分は、分母を \Pi _i \left (1-L_i x  \right )で固定して分子の値を要素として持つと遷移が \mathrm{O} (N)、DP全体で\mathrm{O}(NM)となりTLに間に合う
  • 結局、「分母がN次以下の分数式を冪級数展開した時の x^{D_i} の係数を求めるクエリに Q個答えよ (N \leq 4000, Q \leq 500, D_i \leq 10 ^ {18} )」という問題に落ちる
    • →高速kitamasa法で解けた!勝ったな(確信)
      • → TLE… (kitamasa法を使った時の計算量は \mathrm{O} (QN \log N \log D)なので、当然落ちる)

  • 線形漸化式、kitamasa法じゃ間に合わない、分母が既に因数分解できている…といえば?
  • →部分分数分解!
  • 最初は変数変換の部分をFFTで実装したが、 \mathrm{O}( N ^ 2 \log N) でTLE…(Arbitrary Mod FFTは定数倍が重い…)
    • よくよく考えると変数変換はFFTせずとも部分分数分解に必要な次数だけ計算すればよく、そうするとlogが落ちて \mathrm{O} (N ^ 2)になる
  • 部分分数分解が出来れば、あとは各クエリに対して \mathrm{O} (N \log ( \min(D , Mod) ) ) で答えを求めればよい
  • 最終的に全体の計算量は \mathrm{O} (NM+N ^ 2+QN \log ( \min( D , Mod) ) ) となり、通る
    • 提出 結構勉強になった

Pythonで抽象化セグメントツリーを生やした

経緯

  • 先日、有志コンを開催した eeicpc #1 参加者のみなさんありがとうございました!
  • testerをやる際、C++だけでなくPypy3でも通したくなった(なんとなく)
  • その中に非可換モノイドを載せる問題があった(これ)
  • Python用のセグメント木(以下「セグ木」)を探したが、非可換モノイドの載るセグ木が見当たらない
  • (max,min,addなど可換でよければじゅっぴーさんのセグ木がよさそう)
  • 仕方がないので自分で生やそう!

ソースコードと使い方

  • 使用は自己責任でお願いします
  • (2020/03/27 18:22 指摘を受けて少し改良)
  • (2020/04/02 immutable版を作成、改良。2ベキセグ木の方が時間の定数倍が軽そうなので2ベキセグ木に。)

  • セグ木に載せるオブジェクトがmutableな場合

import copy
class SegmentTree:
    def __init__(self, N, func, I):
        self.sz = 2**(N-1).bit_length()
        self.func = copy.deepcopy(func)
        self.I = copy.deepcopy(I)
        self.seg = [copy.deepcopy(I) for i in range(self.sz * 2)]

    def assign(self, k, x):
        self.seg[k + self.sz] = copy.deepcopy(x)

    def build(self):
        for i in range(self.sz - 1, 0, -1):
            self.seg[i] = self.func(self.seg[2 * i], self.seg[2 * i + 1])

    def update(self, k, x):
        k += self.sz
        self.seg[k] = copy.deepcopy(x)
        while k > 1:
            k >>= 1
            self.seg[k] = self.func(self.seg[2 * k], self.seg[2 * k + 1])

    def query(self, a, b):
        L = copy.deepcopy(self.I)
        R = copy.deepcopy(self.I)
        a += self.sz
        b += self.sz
        while a < b:
            if a & 1:
                L = self.func(L, self.seg[a])
                a += 1
            if b & 1:
                b -= 1
                R = self.func(self.seg[b], R)
            a >>= 1
            b >>= 1
        return self.func(L, R)
  • immutableな場合
class SegmentTree:
    def __init__(self, N, func, I):
        self.sz = 2**(N-1).bit_length()
        self.func = func
        self.I = I
        self.seg = [I] * (self.sz * 2)
 
    def assign(self, k, x):
        self.seg[k + self.sz] = x
 
    def build(self):
        for i in range(self.sz - 1, 0, -1):
            self.seg[i] = self.func(self.seg[2 * i], self.seg[2 * i + 1])
 
    def update(self, k, x):
        k += self.sz
        self.seg[k] = x
        while k > 1:
            k >>= 1
            self.seg[k] = self.func(self.seg[2 * k], self.seg[2 * k + 1])
 
    def query(self, a, b):
        L = self.I
        R = self.I
        a += self.sz
        b += self.sz
        while a < b:
            if a & 1:
                L = self.func(L, self.seg[a])
                a += 1
            if b & 1:
                b -= 1
                R = self.func(self.seg[b], R)
            a >>= 1
            b >>= 1
        return self.func(L, R)
  • 使用例 : 区間最大クエリ
  • (抽象化してあるのでライブラリをいじらずに書ける)
# 要素数10のセグ木を宣言
seg = SegmentTree(10, max, -(2 ** 60))
# リストaで初期化
for i in range(10):
    seg.assign(i, a[i])
# 構築
seg.build()
# 0番目の要素を10に更新
seg.update(0, 10)
# 半開区間[3 , 6)に含まれる要素の最大値を取得
ma = seg.query(3, 6)

使ってみる

注意:行列+セグ木で有名な企業コンの問題のネタバレがあります
~~~~ 以下、ネタバレ防止のため改行 ~~~~



















  • 例1 : DISCO!
  • 「TL : 13s」という仰々しいTLのわりに解法はシンプルな問題
  • 行列をセグ木に乗せれば解ける(逆行列を用いた線形解法もあるようだ)
  • セグ木の要素数 10 ^ 6 と最大級、これが通れば大体通りそう
  • AtCoderなのでPython3+numpyで挑んでみる
# importとセグメント木は省略、適宜補ってください
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
S = readline().rstrip().decode()
Q = int(readline())
m = map(int, read().split())
LR = zip(m, m)
seg = SegmentTree(len(S), np.dot, np.identity(6, dtype=np.uint32))

for i in range(len(S)):
    if S[i] == 'D':
        seg.seg[i + seg.sz][0][1] = 1
    elif S[i] == 'I':
        seg.seg[i + seg.sz][1][2] = 1
    elif S[i] == 'S':
        seg.seg[i + seg.sz][2][3] = 1
    elif S[i] == 'C':
        seg.seg[i + seg.sz][3][4] = 1
    elif S[i] == 'O':
        seg.seg[i + seg.sz][4][5] = 1
seg.build()

for L, R in LR:
    mat = seg.query(L - 1, R)
    print(mat[0][5])
  • ライブラリを貼ってちょろっと書くだけ、簡単だな!
  • →TLE どうして…

  • 仕方がないので重そうな部分を適当に書き替える 具体的には

# __init__関数内
self.seg = [copy.deepcopy(I) for i in range(self.sz * 2)]
  • が明らかにヤバそうなので、これを
self.seg = np.zeros(self.sz * 2 * 6 * 6, dtype=np.uint32).reshape(self.sz * 2, 6, 6)
for i in range(N):
  self.seg[i + self.sz] = I
  • に直す
  • →ギリギリ通りました( 12155ms ) 俺の勝ち! 提出



read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
def main():
 N = int(readline())
    m = map(int, read().split())
    seg = SegmentTree(N, gcd, 0)
    for i in range(N):
        seg.assign(i, next(m))
    seg.build()
    ans = 1
    for i in range(N):
        ans = max(ans, gcd(seg.query(0, i), seg.query(i+1, N)))
    print(ans)

main()
  • 今回はセグ木に載るオブジェクトがimmutableなのでimmutable版を使う
  • mutable提出 952ms
  • immutable提出 587ms
  • copy.deepcopy(I)が無くなった分だけはっきり軽くなっている

おわりに

  • C++最高!一番好きな競プロ用言語です!(は?)
  • 冗談はともかく、セグ木を使わざるを得ない問題は高速な言語を使った方がよさそう
  • (追記:なぜかというと、僕はpythonでDISCO!を通すのに2時間弱かかったので)
  • 高速に動かすのがかなり難しく、適切に書かないと \mathrm {O} ( N \log N ) , N \leq 2 \cdot 10 ^ 5がTLEする
  • Python使いならセグ木を使わない方向で強くなった方がよさそう
  • (例えばDISCO!は更新クエリがないので逆行列を使って線形で解くべき)
  • ただ、セグ木+Pythonがメチャクチャ弱いというわけでもなさそうなのが書いてみてわかったので記事にした
  • ライブラリの使用は自己責任でお願いします、問題点があったら教えてください

Educational Codeforces Round 78: F Cards

色々な解法があるようなので理解のために調べてまとめてみた。

問題概要

https://codeforces.com/contest/1278/problem/F

 JOKERが1枚含まれたm枚のカードから、カードを1枚引いて戻す操作をn回繰り返す。  JOKERを引いた回数をxと置くとき、x^kの期待値E\ [ x ^ k ]を求めよ。

TLを見たところ(1)~(4)の解法が考えられるようだ。

(1)自分の解法 (ゴチャゴチャやる)

p(x) := (JOKERをx枚引く確率)とおくとp(x)= {( \frac{1}{m} ) }^ n {( \frac{m-1}{m} ) }^ {n-x} \ _ n C _ x が成り立つ。 このときFPSP(T) = \sum _  { i = 0 } ^ {\infty} p(x) T ^ x のようにおくと、二項定理より P(T) = {( \frac{1}{m}T + \frac{m-1}{m} ) }^ n = m^{-n} \cdot (T + m - 1 ) ^ n が成り立つ。

ここで、A(T) = \sum _  { i = 0 } ^ {\infty} a(x) T ^ x なるA(T)について、「Tで微分してT倍する」という操作をk回行った関数をB_k(T)とおくと、
 B_0(T) = 1 \cdot a_0 + 1 \cdot a_1 T + 1 \cdot a_2 T ^ 2 + 1 \cdot a_3 T ^ 3 + \ldots
 B_1(T) = 0 \cdot a_0 + 1 \cdot a_1 T + 2 \cdot a_2 T ^ 2 + 3 \cdot a_3 T ^ 3 + \ldots
 B_2(T) = 0 \cdot a_0 + 1 \cdot a_1 T + 4 \cdot a_2 T ^ 2 + 9 \cdot a_3 T ^ 3 + \ldots
 B_3(T) = 0 \cdot a_0 + 1 \cdot a_1 T + 8 \cdot a_2 T ^ 2 + 27 \cdot a_3 T ^ 3 + \ldots
のように係数がTのべき数倍されていき、一般にB _ k ( T ) = \sum _ { i = 0 } ^ {\infty} { i ^ k \cdot a _  x T ^ x } となる。

求めたい値は E [ T  ^  k  ] = \sum _ { i = 0 } ^ {\infty} {  p(i) i ^ k  } であるため、P(T)に上記の操作をk回行ってT=1を代入すればよいとわかり、あとは適切な式変形をすると O ( k ^ 2 ) のDPに持ち込める。

(2)モーメント母関数

モーメント母関数の意味と具体例 | 高校数学の美しい物語
積率母関数 - Wikipedia
詳しくはよく知らないが、リンク先を読んでやることをまとめると
(a) P(T)T = e ^ t を代入する
(b) k次の係数を求めると、それが答え
となる。すごい(こなみかん) よってM ( e ^ t ) = m ^ {-n} ( e ^ t + m - 1 ) ^ n  k次の項をFPSのpowを使ってO(k \log k)で求めればよい。

(3) 多項式補間

P(T)のa微分(0 \leq a \leq k )にT=1を代入していくと以下の式が導かれる。
(  {} ^ {(k)}は下降階乗冪)
 E [ 1 ]= 1 \  , \ E [ x ]= \frac{n}{m} \ , \ E [ x ( x - 1 ) ] = \frac { n ( n - 1 ) } { m ^ 2 } \ \ldots \  , \  E \ [ x ^ {(k)} ] = \frac{n ^ {(-n)} } { m ^ k }
これが求まれば、あとは
 x ^ k = a_0 \cdot 1 + a_1 \cdot x + \ldots + a _ k  \cdot x ^ {(k)}
を満たす係数をニュートン補間で求めればO(k ^ 2 )  E [ x ^ k ] が求まる。
(a_nは性質の良い数列なので式変形や包除などにより係数を直接導出することも可能。後述。)

(4) スターリング数

editorialはこの解き方で解説されているようだ。
確率変数x _ i を「 i番目にジョーカーを引いたら1、引かなかったら0」となるようにとる。この時、
 E [ x ^ k ] = E [ ( x _ 1 + x _ 2 + \ldots + x _ n ) ^ k ]
であり、さらに期待値の線形性を用いると、
 E [ x ^ k ] = \sum _ { 0 \leq i \leq k } { (i種類の確率変数の積) \cdot (その確率)  }
 = \sum _ {0 \leq i  \leq k} { (i種類のものを全て使うようにkコ選ぶ ) {} _ n C _ i \cdot \frac { 1 } {m ^ {i} }}
となる。(i種類のものを全て使うようにkコ選ぶ)というのはいわゆる第二種スターリング数なので、包除や漸化式を用いてO(k ^ 2)で計算できる。

(1)と(2)は同値

P ( T ) = M ( {e} ^ t ) について、以下の操作は変数変換すれば等しいとわかる。
(a) P(T)x微分してx倍する操作をk回繰り返してx=1を代入する
(b)M ( e ^ t ) t微分する操作をk回繰り返してt=0を代入する
また、(b)の操作は M ( e ^ t ) k次の項を求めるのに等しい。

(3)と(4)は同値

途中で補間により求めた数列a_nは実は第二種スターリング数の定義である。
よってやってることはだいたい同じ(論理の飛躍さん…)

(2)と(4)は同値

(4)で得られた式を第二種スターリング数  S ( k , i ) を用いて書き直す。  E [ x ^ k ] = \sum _ {0 \leq i  \leq k}{ S(k , i) \cdot\frac {  {} _ n C _ i   } {m ^ {i} } }
ところで第二種スターリング数 S ( k  ,  i ) は、実は { ( e ^ t - 1) } ^ i  t ^ k 次の係数である。
(形式的)べき級数と数え上げの写像12相との関係性 前編 - Senの競技プログラミング備忘録
また、(2)の式を変形すると
 E [ x ^ k ] = [ t ^ k ] M ( e ^ t ) = [ t ^ k ] ( m ^ {-n} ( ( e ^ t - 1 ) + m ) ^ n ) = [ t ^ k ] \sum _ { 0 \leq i \leq k } (e ^ t - 1 ) ^ i \cdot \frac { {} _ n C _ i  } { m ^ i }
となる。(  e ^ t - 1 は1次以上であることからsumの範囲がk以下に限定されることを用いた)
よってスターリング数のFPS表記と合わせると確かに(2)と(4)は同じ式である。

おわりに

気になったので色々調べてみてまとめたが、正直よくわかってないところも多い。間違いがあったらご指摘いただけると助かります。

Beginners Typhoon After Contest#01 EX - No Bingo!

初めてブログ記事を書いた
フォントやレイアウトが壊れまくっており厳しい…

www.hackerrank.com
昨日ようやく解説の意味がわかったので自分なりの理解をまとめておく。

題意の問題は以下のように読み替えられる。

順列p _ iについて、以下の条件を満たすp _ iの場合の数を求めよ。
条件 \  \ldots \ p _ i = i ,p _ j = N + 1 - jを満たすi,jがともに一つ以上存在する。


場合分けをする。
A \ldots 全事象の場合の数
B \ldots 全てのiに対してp _ i \neq iとなる場合の数
C \ldots 全てのiに対してp _ i \neq N + 1 - iとなる場合の数
D \ldots 全てのiに対してp _ i \neq iかつp _ i \neq N + 1 - iとなる場合の数

とすると、求める場合の数はA - (B + C) + Dとなる。
A,B,Cは容易に求まるのでDについて考える。

まずNが偶数の時について考える。(このときN = 2Mとおく。)

a _ n := (少なくともn個のiについてp _ i = iまたはp _ i = N + 1 - iとなる順列の場合の数)
とおくと、包除原理よりD = \sum _ {i=0} ^ N ( a _ i (-1) ^ i ) が成り立つ。
よってa_nを求めればDの値が求まることになるが、直接求めるのは難しい。

ここで1N2N-1\ldotsMM+1をペアにして考える。
これらのM個の組のうちいくつの組が条件を満たすかで場合分けすることを考える。
条件( \ast ) を以下のように定める。

 条件( \ast ) \  \ldots \  とある組(i , N + 1 - i)で{ \bf "斜めのラインを防いだ穴"}が1つ以上ある。
 すなわち、 p _ i = i , p _ i = N + 1 - i , p _ {N + 1 - i } = i , p _ {N + 1 - i} = N + 1 - i
 のうち1つ以上を満たす。

次に数列f _ kの第n項を以下のように定める。

( f _ k ) _ n := (少なくともk個の組について( \ast )  を満たすとき、
k個の組の中で{\bf"斜めのラインを防いだ穴"}の個数がn個となる順列の場合の数)


この時、1対1の性質から a _ n = \sum _ {k=0} ^ M ( f _ k ) _ n が成り立つため、
f _ kの和を高速に求めれば題意の問題が解けることがわかる。

さて、組(i , N + 1 - i)において斜めのラインを防いだ穴がどのように配置され得るかを考えると、
2個配置される \ldots (i , i)(N+1-i,N+1-i) \ , \ (i,N+1-i)(N+1-i,i)の2通り
1個配置される \ldots (i , i)\ ,\ (N+1-i,N+1-i) \ , \ (i,N+1-i)\ , \ (N+1-i,i)の4通り
となる。
この事実をもとに形式的冪級数を導入すると f _ k = \ _ n C _ k  ( 2 x ^ 2 + 4 x ) ^ k となるため、
二項定理より a = \sum _ {k=0} ^ M ( f _ k ) = (2 x ^ 2 + 4 x + 1 ) ^ M が従う。
(この辺りはmaspy氏の形式的冪級数の解説記事を読むと理解しやすいと思う。)

N = 2 M + 1と表される場合も、同様の考察を経ることで a = (2 x ^ 2 + 4 x + 1 ) ^ M (x + 1)を導ける。
あとは適切な実装によりこの問題を高速に解くことが出来る。
提出
Programming Problems and Competitions :: HackerRank

感想
とても難しい…こういう考察を素早く進められるようになりたい。