AtCoder Beginner Contest 306 の E 問題で使ったふたつのヒープキューを使うと削除をシミュレートできるやつ

競技プログラミング（AtCoder Beginner Contest）

atcoder.jp

問題はこちら。

はじめに問題文を読んだとき、「平衡二分木か…？」とも思いましたが、競プロ界隈有名事実の「ヒープキューを２つ使えば、削除も対応可」を思い出しことなきを得ました。（一日置いて考え直してみると、座圧して Fenwick Tree ２本用意する感じでも解けなくもなさそう？　同じ値が複数あるときの取り扱いが面倒そうかも）

意外と python のサンプルスクリプトがネットに転がってなさそうだったので、次回はコピペで済ませたいなという思いから、この度パッと貼って使えるスニペットをご用意する運びとなりました。

以下、ふたつのヒープキューを使って、削除をシミュレートできるやつのスクリプトです。

import heapq

class DeletableHeapq:

    def __init__(self, q):
        # 一度シャローコピーしているので、計算量が心配な時はコピーを無くす
        # コピーをなくした場合は、引数に渡した q 自体もヒープキュー化されてしまうので注意
        self.q = [qi for qi in q] 
        heapq.heapify(self.q) 
        self.removed_q = []

    def pop(self):
        return heapq.heappop(self.q)

    def push(self, v):
        heapq.heappush(self.q, v) 

    def remove(self, v):
        heapq.heappush(self.removed_q, v) 
        for _ in range(len(self.removed_q)):
            if self.q[0] == self.removed_q[0]:
                heapq.heappop(self.q)
                heapq.heappop(self.removed_q)
            else:
                break

（追記）
とか言ってたら、記事を投稿した直後に他の方の Python での実装を発見しました。
[Python] heapqを応用して挿入/削除/最小値取得を効率的に行う - socha77’s blog
おそらくですが、こちらの方の実装の方が処理が少し早いはず。

AtCoder Beginner Contest 304 の F 問題に出てきた約数包除のサンプルコード

競技プログラミング（AtCoder Beginner Contest）

問題の元ネタはこちら。
ABC の F 問題で、約数系包除原理の問題が出るような時代になったみたいです。

atcoder.jp

約数系包除原理で検索すると、見つかる有名どころの解説記事はここらへんでしょうか。本記事ではあまり包除原理自体の説明には踏み込まないので、仕組みを知りたい方はこちらの記事に飛ぶのがおすすめです。

qiita.com

以下、ABC304F の AC コード を少しだけ書き換えたものを置いておきます。
約数包除系問題の解答コードの雛型として利用することを想定しています。

from itertools import combinations
import math

def dividers_pie(n):
    # 前処理のエラトステネスの篩
    prime_flg = [True] * (n + 1)
    prime_flg[0] = False
    prime_flg[1] = False
    for i in range(2, math.ceil(n ** 0.5) + 1):
        if not prime_flg[i]:
            continue
        for j in range(2, n // i + 1):
            prime_flg[i * j] = False
    primes = [i for i, flg in enumerate(prime_flg) if flg]

    # 素因数と、素因数の個数を求めておく
    prime_factors = [p for p in primes if n % p == 0]
    pf_count = len(prime_factors)

    ans = 0
    # 場合によっては、このループを 1 始まりなどにする必要あり。
    # 元ネタの ABC304F では 1 始まりにしてた
    for i in range(0, pf_count + 1):
        for choosed_pfs in combinations(prime_factors, i):
            factor = 1
            for pf in choosed_pfs:
                factor *= pf
            if factor > n:
                continue
            # ここも問題によりけり。ABC 304F では符号は逆
            sign = 1 if i % 2 == 0 else -1
            # あとは hogehoge を問題に合わせて求める
            ans += sign * hogehoge

    return ans

計算量につきましては、エラトステネスの篩の部分が O(n \times log(log(n))) ほど。

    for i in range(0, pf_count + 1):
        for choosed_pfs in combinations(prime_factors, i):

ここのループが、 O(2^{x}) （ただし、 x は n の素因数の個数）となりますが、 n \leq 10^{5} のとき x \leq 6、 n \leq 10^{6} のとき x \leq 7 なので大したループにはなりません。

「ネストすごろく」のゴール確率を計算してみた

SNS 関連

震源地はこちら。

検索した感じ、この記事ができた時点で一番端的な解答に見えた方のツイートがこちら。（とか言ってたら、ぎりぎり公式解答がこの記事より先に出てて泣いた）

というわけで、やっていきましょう。

確率空間 (S, \Sigma, \mu) および、可算無限個の確率変数 X_{1}, X_{2}, \dots を用意します。
任意の i \in \mathbb{N} について、 X_{i} \colon S \rightarrow \lbrace 1,2,3,4,5,6 \rbrace です。
また X_{1}, X_{2}, \dots の組は独立であり、任意の i \in \mathbb{N}, j \in \lbrace 1, 2, 3, 4, 5, 6 \rbrace について \mu(\lbrace s \in S | X_{i}(s) = j \rbrace) = \frac{1}{6} とします。

有限回の移動でゴールに到達したことを表す論理式 \psi を置きます。
このとき A = \lbrace s \in S | \psi(s) \rbrace とすると、求めたい値は以下のようになります。

\displaystyle{\mu(A)}

つぎに、最も深いネストが第 n 層（ n \in \mathbb{N}）以内という条件下で、ゴールに到達したことを表す論理式 \psi_{n} を置きます。
このとき、以下の３つが成り立ちます。

  • 任意の n \in \mathbb{N}, s \in S において、 \psi_{n}(s) \rightarrow \psi(s)
  • 任意の s \in S について、ある n \in \mathbb{N} が存在し \psi(s) \rightarrow \psi_{n}(s)
  • n \leq m ならば、任意の s \in S について \psi_{n}(s) \rightarrow \psi_{m}(s)

A_{n} = \lbrace s \in S | \psi_{n}(s) \rbrace とおくと、以下が成り立ちます。

  • 任意の n \in \mathbb{N} において、 A_{n} \subset A_{n + 1}
  • さらに A_{n} \subseteq S より、 \lbrace A_{n} \rbrace が収束する
  • \displaystyle{ \lim_{n \rightarrow \infty} A_{n} = A}

一番最後の主張は \displaystyle{ \lim_{n \rightarrow \infty} A_{n} = B} としたとき、先述の \psi と \psi_{n} の性質より B \backslash A = \emptyset かつ A \backslash B = \emptyset から示せます。

以上より、 \displaystyle{ \lim_{n \rightarrow \infty} \mu(A_{n}) = \mu(A)} が成り立ちます。

次に、ゴールからのマス目の距離が k \in \lbrace 1, 2, 3, 4, 5, 6 \rbrace のマスからスタートして、最も深いネストが第 n 層（ n \in \mathbb{N}）以内という条件下でゴールに到達できる確率を p_{n, k} と置きます。 p_{n, 6} = \mu(A_{n}) である点に注意してください。

このとき、以下が成り立ちます。

\displaystyle{ p_{n+1, 1} = 1 \\ p_{n+1, 2} = \frac{p_{n, 6}}{6} p_{n+1, 1} + \frac{5}{6} \\ p_{n+1, 3} = \frac{p_{n, 6}}{6} (p_{n+1, 1} + p_{n+1, 2}) + \frac{4}{6} \\ p_{n+1, 4} = \frac{p_{n, 6}}{6} (p_{n+1, 1} + p_{n+1, 2} + p_{n+1, 3}) + \frac{3}{6} \\ p_{n+1, 5} = \frac{p_{n, 6}}{6} (p_{n+1, 1} + p_{n+1, 2} + p_{n+1, 3} + p_{n+1, 4}) + \frac{2}{6} \\ p_{n+1, 6} = \frac{p_{n, 6}}{6} (p_{n+1, 1} + p_{n+1, 2} + p_{n+1, 3} + p_{n+1, 4} + p_{n+1, 5}) + \frac{1}{6} }

これをうまいこと式変形すると、以下を得ます。

\displaystyle{ p_{n+1, 1} = 1 \\ p_{n+1, 2} = (\frac{p_{n, 6}}{6} + 1) p_{n+1, 1} - \frac{1}{6} \\ p_{n+1, 3} = (\frac{p_{n, 6}}{6} + 1) p_{n+1, 2} - \frac{1}{6} \\ p_{n+1, 4} = (\frac{p_{n, 6}}{6} + 1) p_{n+1, 3} - \frac{1}{6} \\ p_{n+1, 5} = (\frac{p_{n, 6}}{6} + 1) p_{n+1, 4} - \frac{1}{6} \\ p_{n+1, 6} = (\frac{p_{n, 6}}{6} + 1) p_{n+1, 5} - \frac{1}{6} }

それぞれ代入することで、

\displaystyle{ p_{n+1, 6} = \frac{1}{7776} p_{n, 6}^{5} + \frac{29}{7776}p_{n, 6}^{4} + \frac{55}{1296}p_{n, 6}^{3} + \frac{25}{108}p_{n, 6}^{2} + \frac{5}{9}p_{n, 6} + \frac{1}{6} }

を得ます。
手計算で求めるのは大変なので、以下の python スクリプトで算出しています。

from fractions import Fraction

# 係数をリストで保持
ans = [0, 0, 0, 0, 0, 1]

for _ in range(5):
    # (p_(n, 6) / 6 + 1) 倍
    ans = [
        (
            ans[i] + ans[i + 1] * Fraction(1, 6)
            if i < 5
            else ans[i]
        )
        for i
        in range(6)
    ]
    # 1/6 を引く
    ans[-1] = ans[-1] - Fraction(1, 6)
    

print(ans)

ちなみに、

\displaystyle{ f(p) = \frac{1}{7776} p^{5} + \frac{29}{7776}p^{4} + \frac{55}{1296}p^{3} + \frac{25}{108}p^{2} + \frac{5}{9}p + \frac{1}{6} }

とおくと、 f(p) について以下が成り立ちます。

  • 連続である
  • 狭義単調増加である
  • 下に凸である
  • 方程式 p = f(p) は区間 \lbrack 0, 1\rbrack において、解 p = \alpha, 1 を持つ。ただし、 \frac{1}{6} \lt \alpha \lt 1 を満たす
  • f'(\alpha) \lt 1 である

以上より、 0 \lt p \lt \alpha ならば、

\displaystyle{ p - \alpha \lt f'(\alpha)(p - \alpha) \lt f(p) - f(\alpha) \lt 0 }

が成り立ちます。

以上から、

\displaystyle{ p_{1, 6} = \frac{1}{6} \\ p_{n+1, 6} = f(p_{n, 6}) }

に対し、

\displaystyle{ q_{1} - \alpha= \frac{1}{6} - \alpha \\ q_{n+1} - \alpha = f'(\alpha) (q_{n} - \alpha) }

を用いて下から評価すると、

\displaystyle{ q_{n} - \alpha = f'(\alpha)^{n-1} (\frac{1}{6} - \alpha) \leq p_{n, 6} - \alpha \lt 0 }

となるため、 \displaystyle{ \lim_{n \rightarrow \infty} p_{n, 6} = \lim_{n \rightarrow \infty} \mu(A_{n}) = \alpha} となります。

ちなみに \alpha がいくらくらいになるかというと

\displaystyle{ \alpha \fallingdotseq 0.5476832304276367 }

となるようです。

以下、算出に使用した python スクリプトです。

from fractions import Fraction

l = [
    Fraction(1, 7776), 
    Fraction(29, 7776), 
    Fraction(55, 1296), 
    Fraction(25, 108), 
    Fraction(5, 9), 
    Fraction(1, 6)
]

# ここも Fraction 型にすると計算時間がとてもかかるので、ans は float 型に。
ans = 1/6

while True:
    prev_ans = ans
    ans = sum(
        li * (prev_ans ** (5 - i))
        for i, li
        in enumerate(l)
    )
    if abs(ans - prev_ans) < 0.0000001:
        print(ans)
        break

ウィルソンの定理

競プロのアルゴリズムまとめ

素数に関する定理で、面白いものを見つけたのでメモ。 自然数 p について、素数全体の集合を \mathbb{P} と置くとき、

(p-1)! \equiv \begin{cases} {-1 \quad (p \in \mathbb{P})} \\ {2 \quad (p = 4}) \\ {0 \quad (otherwise)} \end{cases} (mod \ p)

表題の定理の場合は p が素数 \Leftrightarrow (p-1)! = -1 \ (mod \ p) となっています。
p が素数でないときの、 p で割った余りがいくつになるかも、ひと手間加えれば示せるので、先のような書き方にしています。

以下、証明です。

p = 1 のとき、
(1 - 1)! = 1 \equiv 0 \ (mod\ 1) より成立します。

p=2 のとき、
(2-1)! = 1 \equiv -1 \ (mod \ 2) より成立します。

p = 4 のとき、
(4 - 1)! = 6 \equiv 2 \ (mod\ 4) より成立します。

p が 4 以外の合成数のとき、
p = mn となる m, n \in \mathbb{N} が存在します。 ただし m, n はどちらも 2 以上 p-1 以下かつ、 \neg (m=2 \land n=2) です。
m \neq n のとき、適当な l \in \mathbb{N} を用いて (p-1)! = lmn と表せるため、 (p-1)! \equiv 0 \ (mod \ p) です。
m = n のとき、 m = n \gt 2 となることから、 2m \leq (p-1) です。ここから、適当な l \in \mathbb{N} を用いて (p-1)! = lm2m = 2lmn と表せるため、 (p-1)! \equiv 0 \ (mod \ p) です。

p が 2 以外の素数のとき、
法 p における原始根 e が存在します。ここから、 e^{1} e^{2} \dots e^{p-1} \equiv (p-1)! \ (mod \ p) となります。 左辺を e^{1} e^{2} \dots e^{p-1} = e^{\frac{p(p-1)}{2}} と式変換し、これを b とおくと、 b^{2} = e^{p(p-1)} \equiv 1^{p} \equiv 1 \ (mod \ p) となり、 p が素数なことから、 b \equiv \pm 1 \ (mod \ p) となります。
b \equiv 1 \ (mod \ p) とすると、 \frac{p(p-1)}{2} = k(p-1) を満たす k \in \mathbb{N} が存在することになりますが、 p が奇数であることより矛盾。よって b \equiv -1 \ (mod \ p) となります。

競プロのアルゴリズムまとめのカテゴリを着けましたが、たぶん競プロで使うことはないです。素数判定は O(\sqrt{p}) で可能ですし。どちらかというと、定理を用いて計算飛ばしたりするのに使う可能性の方が高いかも。
もしくは証明の途中に登場した原始根が、NNT（数論変換）を理解する上での前提知識となっているので、そちらの方が競プロでは有用かもしれません。

AtCoder Heuristic Contest 015 参加しました

競技プログラミング（AtCoder Heuristic Contest）

atcoder.jp

ご参加された皆々様、お疲れさまでした。

得点：121,539,661 点
順位：184位（1447人中）

でした。ヒューリスティックで青パフォは初だったのでうれしい。

解法はこちらです。
いわゆる貪欲 + モンテカルロ系解法（のはず）です。

atcoder.jp

一番最初に出した解答はこちら。
問題文を読み込んで必要そうな関数を実装し終えたあと、毎ステップごとに（連結成分の大きさの二乗和）/（各キャンディーの個数の二乗和）が最大になるように傾けていきます。ただの貪欲。

atcoder.jp

評価関数を調整するのが正しいのかなぁと考えながらも、いいアイディアが思いつかなかったので、処理時間が許す限り、１手先の盤面をシュミレーションする方針に。
サンプリングする数を増やすごとに点数が伸びたのでこれはこれでよかったのかな。

Twitter でも流れてきましたが、きわめてシンプルなルールベースでも、もっと高得点が出せるらしいです。悲しい。

def main():
  f = list(map(int, input().split()))
  prev = None
  for fi in f[1:]:
    _ = input()
    move = None
    if fi == 1:
      move = "B"
    elif prev == "B":
      move = "F"
    elif fi == 2:
      move = "R"
    elif fi == 3:
      move = "L"
      
    print(move)
    prev = move
     
      
if __name__ == "__main__":
  main()

これで、133,637,155 点出ます。

AtCoder Beginner Contest 267 の G 問題の解説に出てきた Eulerian number

競技プログラミング（AtCoder Beginner Contest）

元ネタその１

atcoder.jp

元ネタその２

en.wikipedia.org

早速ですが python スクリプトをドン。

def eulerian_number(n, m, MOD=998244353):
    ans = 0
    sign = 1
    comb = 1    
    for i in range(m+1):
        ans += sign * comb * pow((m + 1 - i), n, MOD) 
        ans %= MOD
        sign *= -1
        comb *= (n + 1 - i)
        comb *= pow(i + 1, -1, MOD)
        comb %= MOD
    return ans

Wikipedia の Explicit formula セクションの式を、愚直に実装したものになります。
\displaystyle{\Sigma} の部分については、Wikipedia の通りに i = m+1 まで足し上げると、 (m + 1 - i)^{n} = 0 となってしまい無駄になるため、python の for ループの上限を i = m までにしています。

で、スクリプトを見ていただければわかる通り、この実装の Eulerian number の計算量は O(m (log(MOD) + log(n) )) となっています。
解説記事は O(n) のアルゴリズムがあるとのことなので、他の計算方法があるということなんでしょうか...？

AtCoder Beginner Contest 266 の E 問題の解析解

競技プログラミング（AtCoder Beginner Contest）

問題はこちら

atcoder.jp

で、AC コードはこちら

atcoder.jp

第一引数が float の場合の python の pow 関数の計算量について詳しくないのですが、もし int の場合と同じであれば計算量 O(log(N)) の解答となります。

やっていることは、公式解説をもとに、簡単なさんすうの知識で計算量を改善しているだけです。以下解説。

まずは解説の式を見ていただきつつ、有理数は全順序集合であることから、解説の式は以下の７パターンに書き下せることがわかります。

\begin{align} f(N) = \begin{cases} \frac{7}{2} & (f(N-1) \lt 1) \cr \frac{1}{6} \times f(N-1) + \frac{10}{3} & (1 \leq f(N-1) \lt 2) \cr \frac{1}{3} \times f(N-1) + 3 & (2 \leq f(N-1) \lt 3) \cr \frac{1}{2} \times f(N-1) + \frac{5}{2} & (3 \leq f(N-1) \lt 4) \cr \frac{2}{3} \times f(N-1) + \frac{11}{6} & (4 \leq f(N-1) \lt 5) \cr \frac{5}{6} \times f(N-1) + 1 & (5 \leq f(N-1) \lt 6) \cr f(N-1) & (6 \leq f(N-1) \end{cases} \end{align}

また、漸化式の形から f(N) が広義単調増加関数であることと、特に f(N) \lt 6 の範囲においては狭義単調増加関数であること、そして f(1) = \frac{7}{2} かつ \displaystyle{\lim_{N \rightarrow \infty} f(N) = 6} から考慮が必要なのは 3 \leq f(N-1) \lt 6 の条件を満たす３つの式だけであることがわかります。

よって f(N) は以下のように書き下せることがわかります。
漸化式が高々一次方程式である場合簡便な解析式が与えられることは周知の事実であり、かつパラメータの取得手法も有名であり、その手法を使用することで以下の式を得ます。

\begin{align} f(N) = \begin{cases} - \frac{3}{2} \times (\frac{1}{2}) ^ {N - 1} + 5 & (N \leq 2) \cr - \frac{5}{4} \times (\frac{2}{3}) ^ {N - 2} + \frac{11}{2} & (N \leq 5) \cr - \frac{47}{54} \times (\frac{5}{6}) ^ {N - 5} + 6 & (5 \lt N) \cr \end{cases} \end{align}

これを、コードに書き下したものが先述の AC コードとなります。