えびちゃんの日記

えびちゃん（競プロ）の日記です。

bitset 高速化が定数倍高速化という風潮に対するお気持ち表明

「bitset で 64 倍高速化できます」「64 は定数なのでオーダーは変わりません」のような説明は、競プロ界隈でたびたび見かけます。

AGC の解説 とかも。

「定数倍高速化を軽視する風潮に対するお気持ち表明」や「競プロのような入力に上限がある場合でもオーダーを信仰する風潮に対するお気持ち表明」などは、一旦おいておきます。

今回のメインは、「長さ \(n\) の bit 列の bitwise AND/OR や shift にかかる時間は \(\Theta(n)\) か？」（あるいは、そう仮定するのは妥当か？）という話です。

tl; dr

word RAM の下で、\(\Omega(\log(n))\) 倍の高速化と言えそうです。

補足

オーダー記法に慣れていない人が以下で怯えないように、さくっとまとめておきます。

  • \(\Theta(f(n))\)：ちょうど \(f(n)\) のオーダー。
    • 競プロ界隈で言うところの \(O(f(n))\) に近そう。
  • \(O(f(n))\)：高々 \(f(n)\) のオーダー。
    • \(O(n^3)\) と書いたとき、実際には \(\Theta(n^2)\) かもしれない。
  • \(\Omega(f(n))\)：少なくとも \(f(n)\) のオーダー。
  • \(o(f(n))\)：\(f(n)\) より真に小さいオーダー。
  • \(\omega(f(n))\)：\(f(n)\) より真に大きいオーダー。

小文字の方は今回は使いません。

導入

次のような問題を考えてみます。

\(0\) または \(1\) からなる長さ \(n\) の列 \(a = (a_0, a_1, \dots, a_{n-1})\) が与えられます。これが回文か判定してください。

  1. for \(i \in \{0, 1, \dots, n-1\}\) do
    1. if \(a_i \ne a_{n-i-1}\) then
      1. return no
  2. return yes

愚直には、このようなアルゴリズムが考えられるわけですが、これの計算量はどうなるでしょう？ 何の仮定も無しに \(O(n)\) 時間（+ worst \(\Theta(n)\) 時間）だと言い切っていいでしょうか？ \(i\) のインクリメントや \(n-i-1\) の計算、\(a\) へのアクセスが \(O(1)\) 時間でできるものと仮定しているはずです（もっと言えば、1 bit の等価判定を \(O(1)\) 時間と仮定しているはずです。さすがにこれくらいは仮定したいですが）。

読者の中には、「\(a\) を順方向と逆方向にシーケンシャルアクセスするだけでいいからそれらの仮定はいらないだろう」と言う人もおられるかもしれない*1ので、別の例も出します。

\(0\) 以上 \(2^n\) 未満からなる長さ \(n\) の整数列 \(a = (a_0, a_1, \dots, a_{n-1})\) が与えられます。その後、「与えられた \(i\) と \(j\) に対して \(a_i+a_j\) の値はなんですか？」という形式の \(m\) 個の質問に答えてください。

これは \(O(n+m)\) 時間と言えてほしいですが、何の仮定も無しには言えません。 そもそも各整数をどういう形式で表すのかについても触れていません。 「ひとつの整数を表すのに \(2^n -1\) 個の連続する bit を用い、そのうち（連続するとは限らない）\(x\) 個の bit が \(1\) で残りが \(0\) のとき、\(x\) を表すとする」などと言われては大変扱いにくいので、もう少しマシなものを要求したいです。 コンピュータで扱いにくいものの話をしてもしょうがないので、2 進法で表すのがいいかなとなります*2。

また、整数 \(i\) を用いて \(a_i\) に \(O(1)\) 時間でアクセスできることが要求されます。 \(a\) の長さは \(n\) であることから、\(\log(n)\) bits 程度の整数は \(O(1)\) 時間で扱えないとお話にならないことがわかります。 \(a_i+a_j\) の計算もしなきゃなので、\(\log(n)\) bits 程度の加算も \(O(1)\) 時間でしたいです。

word RAM の話

こうした要望に合った計算モデルとして、word RAM というのがあります*3。次のような仮定です。

  • \(U=2^u\) 個の word からなる配列をメモリと呼ぶ。
    • 入力はこのメモリに収まる必要がある。
  • word をポインタとして扱える。
    • word を添字として用いてメモリ（の全範囲）にアクセスできるということ。
  • 長さが \(O(1)\) の連続する word に対する読み書き、四則演算 (+, -, *, /, %)、ビット演算 (&, |, ^, !, <<, >>), 比較演算 (==, <, etc.) は \(O(1)\) 時間でできる。
    • ! はビット反転。C とかで言うところの ~。

word によってメモリにアクセスできることから、word は最低でも \(u\) bit 必要です。 word の bit 長を \(w\) とすると \(w\ge u\) です。

これらから、連続する \(u\) bit に対する各種演算が \(O(1)\) 時間でできることが従います。 ただし、「メモリから飛び飛びの \(u\) 個の bit を集めてきて、結合してひとつの整数と見なし、それらを用いた演算を行う」ということは \(O(1)\) 時間でできるとは言っていないことに注意してください。

本題

さて、元々の話に戻ります。

長さ \(n\) の bit 列に対するビット演算の計算量は？

連続する \(n\) 個の bit 列、すなわち \(n/w\) word に対する演算となります。 \(w\) は最低でも \(\log(n/w)\) bit 必要なことから \(w \ge \log(n/w)\)、すなわち \(w\cdot 2^w \ge n\) となります。

式変形がんばるよ \(w\cdot 2^w = n\) の形の式を解くのは初等的にはつらいです。Lambert の \(W\) 関数 というのがあります。 一般には複素数で定義されますが、ここでは実数についてのみ考えます。 \[ ye^y = x \implies y = \begin{cases} W_0(x), & \text{if }x \ge 0; \\ W_{-1}(x), & \text{if } -1/e\le x\lt 0. \end{cases} \] ここでは \(x\ge 0\) のみ気にするので、\(W_0\) のみ考えます。 \[ \begin{aligned} w\cdot 2^w &\ge n \\ w\cdot 2^{w\cdot\log_2(e)\cdot\log_e(2)} &\ge n \\ w\cdot e^{w\cdot\log_e(2)} &\ge n \\ w\cdot\log_e(2)\cdot e^{w\cdot\log_e(2)} &\ge n\log_e(2) \\ w\cdot\log_e(2) &\ge W_0(n\log_e(2)) \\ w &\ge W_0(n\log_e(2))\cdot\log_2(e). \end{aligned} \]

また、\(W\) 関数について、以下のことが知られています*4。 \[ x\ge e \implies \ln(x)-\ln(\ln(x)) \le W_0(x). \] これらから、\(w\in \Omega(\log(n/\log(n)))\) が言えそうです。 がんばりおわり。

不等式を示すだけなら、\(W\) を使わなくてもできるかもです。

（追記）

まず、\(\log(n)\ge\log(n/\log(n))\) はいいとして、 \[ \begin{aligned} \log(n/\log(n)) &= \log(n) - \log(\log(n)) \\ &\ge \log(n) - \log(\sqrt{n}) \\ &= \log(\sqrt{n}) \\ &= \frac{1}{2}\log(n) \end{aligned} \] なので、\(\log(n/\log(n)) = \Theta(\log(n))\) でした。

やっぱり \(W\) を使わずにできましたね。

結局、\(w\in\Omega(\log(n))\) が言えたので、\(n/w\in O(n/\log(n))\) となります。というわけで、bitset 高速化は \(\Omega(\log(n))\) 倍の改善と言えそうです。

おまけ

別のモデルの話

たとえば、もっとやばい計算モデルを仮定して、「任意のサイズの整数の四則演算は \(O(1)\) でできます」というのがある程度の妥当性をもって言えれば、bitset 高速化は \(n\) をひとつぶん落としたと主張できそうです。

word RAM の元々の論文*6では、

これこれの論文では長さ \(n\) の整数配列 \(a\) のソートを \(O(n)\) 時間で達成したと述べているが、\(n^2 \log(\max a)\) の長さ（ビット？）の除算を \(O(1)\) 時間でできると仮定した上でのアルゴリズムであり、やばいだろ

といった旨のことが書かれていました。

（追記）上の式では \(w\) の下限にしか触れていないので、別に \(w=n\) のようにワードサイズがやばいモデルを仮定することもできそうです。が、ちょっと乱暴だと思います。また、メモリを無駄に使うことで「\(w\) がもっと必要でしょう」と主張することもできそうですが、無駄に使ったメモリのぶんの時間のせいで損をしそうです。

（追記）別にメモリを無駄に使わずとも「\(w\) がめちゃ大きいモデルです」というのを仮定するの自体は許される気がします。現実的にはそれが妥当ではなさそうなので word RAM を考えていそうではありますが。

www.cs.cmu.edu

今回の内容に関連する講義資料です。word RAM 以外にも、Turing machine や論理回路を用いた場合に計算量はどうなるか？といったことが書かれています。

許容する演算の話

word RAM においては word size の四則演算を \(O(1)\) と言っていますが、乗算や除算はこれに含めないとする立場もあるようです。 これは、特定の条件を満たす回路によって加算は実現できるが乗算は実現できないという事情によるようです。 n bit の乗算が M(n) 時間でできれば、Newton 法によって n bit の除算も O(M(n)) 時間でできるので、乗算を仮定するなら除算も仮定されそうです。 ここ嘘で、\(M(n)\in O(1)\) なので、\(O(\log(w))\) 時間とかになっちゃいます。無視してください。

また、立っている最上位 bit の位置は、乗算を許せば次の手法で \(O(1)\) 時間で得られます。

rsk0315.hatenablog.com

これにより、立っている最下位 bit が n & -n で得られることと合わせて、立っている最下位 bit の位置も \(O(1)\) 時間で得られます。

word RAM 上のアルゴリズムの話

こういう記事があります。

qiita.com

まとめ

\(O(n/w)\) などの \(w\) は定数ではなくて、入力サイズ \(n\) に依っていると見なす立場もあるよ、という話でした。

マシンが入力に依ることに納得がいかない人は、「大きいサイズの問題を、メモリの少ない昔の PC で解くのは不可能」みたいなことを考えてみるといいかもしれません？ 大きなメモリをちゃんと扱うことができるマシンであれば、それに伴っていくらか大きい値も扱えるようになる、といったような気持ちに近い気がします。

補足

「実際にはキャッシュとかの関係できっちり 64 倍高速化されるとは限らないから 定数 倍高速化ではない」とかそういう類の主張ではないです。

おわり

log は定数学派の人にとっては bitset 高速化は定数倍高速化ということがわかりました。いかがでしたか？

*1:「インクリメントはならし定数時間なので〜」とか思った人もいるかも？

*2:当然、文脈によってはコンピュータに縛られる必要はありませんが、ここでは一応競プロに近い文脈なので。

*3:入力に応じてサイズが変わるの自体は transdichotomous model と呼ばれるもので、単位時間でできる特定の演算を決めたものが word RAM かもしれません。

*4:Hoorfar, A., & Hassani, M. (2008). Inequalities on the Lambert W function and hyperpower function. J. Inequal. Pure and Appl. Math, 9(2), 5-9.

*5:このあたり、記法が怪しいので困りました。

*6:Fredman, Michael L., and Dan E. Willard. "Surpassing the information theoretic bound with fusion trees." Journal of computer and system sciences 47.3 (1993): 424-436.

素数判定するよ

\(n\) 以下の素数の個数を求める関数 \(\pi(n)\) を持っていたとします。 このとき、\(n\) の素数判定は \(\pi(n)-\pi(n-1) = 1\) かどうかでできます。

fn is_prime(n: usize) -> bool {
    n > 0 && prime_pi(n) - prime_pi(n - 1) == 1
}

たとえば \(O(n^{3/4}/\log(n))\) でできます*1。

お粗末さまでした。

おまけ

素数の数え上げが気になる人向けの記事：

rsk0315.hatenablog.com

*1:素直に \(O(n^{1/2})\) 時間の試し割り法をしましょう。

n 番目の素数を求めるよ

軽めの記事です。

\(n\) 番目の素数を \(p_n\)、\(n\) 以下の素数を \(\pi(n)\) と書くとします。 このとき、\(p_n\) は \(\pi(p_n) \ge n\) かつ \(\pi(p_n-1) \lt n\) を満たします。 よって、これは二分探索で求められます。

上限を決めるのが大変そうですが、指数探索っぽくすればいいです。 指数探索というのは、上限を \(2, 4, 8, \dots\) と倍々に増やしていき、true/false の境界を跨いだらその範囲で二分探索を行う方法です。 今回は、自明に \(p_n \gt n\) なので、上限の初期値はそれに基づいて決めてよいです（別に \(2\) から始めてもよいです）。

fn main() {
    assert_eq!(nth_prime(3), 5);
    assert_eq!(nth_prime(4118054813), 99999999977);
}

fn nth_prime(n: usize) -> usize {
    let mut lo = n;
    let mut hi = 2 * n;
    let small = |i| prime_pi(i) < n;
    while small(hi) {
        hi *= 2;
    }
    while hi - lo > 1 {
        let mid = lo + (hi - lo) / 2;
        if small(mid) {
            lo = mid;
        } else {
            hi = mid;
        }
    }
    hi
}

ここで、prime_pi は 別の記事 で書いたやつを使うと、\(O(n^{3/4}/\log(n))\) time です。 また、\(p_n = O(n\log(n))\) が知られているので、全体で \(O( (n\log(n))^{3/4}/\log(n\log(n))\cdot \log(n\log(n)) = O( (n\log(n))^{3/4})\) time です。

お粗末さまでした。

眠れない夜は素数の個数でも数えましょう

別に、おふとんから出たくない朝とかに数えてもいいと思います。

\(n\) 以下の素数の個数を数えたいときどうしますか？ 各整数 \(2\le i\le n\) に対して素数判定して数えます？ 素数判定を試し割り法でやると \(\Theta(n\sqrt{n})\)、線形篩*1を使えば \(\Theta(n)\) とかですね。 これを \(o(n)\) でできたらうれしくないですか？という話です。

今回は \(O(\sqrt{n})\) space, \(O(n^{3/4}/\log(n))\) time の話がメインです。 定数倍があまり軽くない \(O(n^{2/3})\) space, \(O(n^{2/3})\) time の話も少し触れます。

（追記：2021/06/05）少し難しめかもですが、アイディア自体は水色もあれば理解できると思います。

（追記：2022/06/29）勉強会用に作ったスライドを公開しました：https://rsk0315.github.io/slides/prime-counting.pdf

前提知識

以下、\(n\) 以下の個数の素数を \(\pi(n)\) と書きます。

導入

素数の個数を数えるにあたって、素数を実際に列挙する必要はないというのが重要です。 実際に列挙していては、「与えられた素数の次の素数はなに？」というのに \(\Theta(1)\) 時間で答えられても、\(\Theta(n/\log(n))\) 時間かかってしまいます*2。

ここで、Eratosthenes の篩について考え直します。次のようなアルゴリズムです。

  1. \(1\) から \(n\) の bool 配列 \(a_i\) を用意する。\(a_1\) のみ false とする。
  2. \(i = 2, 3, \dots, \lfloor\sqrt{n}\rfloor\) に対して、次の処理をする。
    1. \(a_i\) が false なら次の \(i\) に移る。
    2. \(j = i, i+1, \dots, \lfloor n/i\rfloor\) に対して、\(a_{ij}\) を false にする。
let mut is_prime = vec![true; n + 1];
is_prime[1] = false;
for i in (2..=n).take_while(|&i| i * i <= n) {
    if !is_prime[i] { continue; }
    for j in i..=n / i {
        is_prime[i * j] = false;
    }
}

ここで、各 \(i\) に対して 2.2 のステップでいくつの \(a_{ij}\) が true から false になるかについて考えます。 この個数をうまく管理できれば、篩い終えた後にいくつ true になっているか、すなわち素数の個数がわかることになります。

この手法は、界隈では Lucy DP などと呼ばれることもあるようです*3。

解説

（追記 2021/05/19）画像とかをつけるとわかりやすくなる気がしたので、そのうちつけると思います。

\(S(v, p)\) を、\(i=p\) まで処理を終えた後に \(a_j\) が true である \(2\le j\le v\) の個数と定義します。 \(p\cdot p \gt v\) であれば \(i=p-1\) のときから配列は変化しないので、\(S(v, p) = S(v, p-1)\) とわかります。\(p\) が素数でない場合も同じです。

\(p\) が \(p\cdot p\le v\) なる素数の場合はむずかしいのでちょっと考えます。 \(j = p, p+1, \dots, \lfloor v/p\rfloor\) に対して \(a_{pj}\) が false になりますが、\(j\) が \(p\) 未満の素数で割り切れるものはすでに篩われています。よって、篩われる個数は「\(p, p+1, \dots, \lfloor v/p\rfloor\) のうちで篩われていないもの」であり、これは \(S(\lfloor v/p\rfloor, p-1)-S(p-1, p-1)\) です。すなわち、次が成り立ちます。 \[S(v, p) = S(v, p-1) - (S(\lfloor v/p\rfloor, p-1) - \pi(p-1) ).\]

さて、\(\pi(n)\) は \(S(n, \lfloor\sqrt{n}\rfloor)\) と等しいので、これを求めたいわけですが、\(S(v, \bullet)\) を求めるためには（その時点で見つかっている素数の個数と）\(S(\lfloor v/\bullet\rfloor, \bullet)\) の形で表せるものさえわかっていれば十分とわかります。 \(\lfloor \lfloor n / p\rfloor / q\rfloor = \lfloor n / pq\rfloor\) が重要です*4。

よって、\(v = 1, 2, \dots, \lfloor\sqrt{n}\rfloor\) と \(v = \lfloor n/1\rfloor, \lfloor n/2\rfloor, \dots, \lfloor n/\lfloor\sqrt{n}\rfloor\rfloor\) についてのみ管理します（サイズ \(O(\sqrt{n})\) の配列ふたつで管理できます）。 以下のコード中では \(\texttt{smalls}[i] = S(i, p)\)、\(\texttt{larges}[i] = S(\lfloor n/i\rfloor, p)\) とします。

fn main() {
    assert_eq!(prime_pi(10), 4);
    assert_eq!(prime_pi(100), 25);
    assert_eq!(prime_pi(100000000000), 4118054813);
}

fn prime_pi(n: usize) -> usize {
    let n_sqrt = floor_sqrt(n);

    let mut larges: Vec<_> = (0..=n_sqrt).collect();
    for i in 1..=n_sqrt { larges[i] = n / i - 1; }
    let mut smalls: Vec<_> =
        (0..n / n_sqrt).map(|x| x.saturating_sub(1)).collect();

    for p in 2..=n_sqrt {
        if p < smalls.len() {
            if smalls[p] <= smalls[p - 1] { continue; }
        } else {
            if larges[n / p] <= smalls[p - 1] { continue; }
        }
        let pc = smalls[p - 1];
        let q = p * p;
        for i in (1..=n_sqrt).take_while(|&i| n / i >= q) {
            // vi = n / i
            // dp[n / i] -= dp[n / i / p] - pc;
            let ip = i * p;
            let cur = *larges.get(ip).unwrap_or_else(|| &smalls[n / ip]) - pc;
            *larges.get_mut(i).unwrap_or_else(|| &mut smalls[n / i]) -= cur;
        }
        for i in (1..smalls.len()).rev().take_while(|&i| i >= q) {
            // vi = i
            // dp[i] -= dp[i / p] - pc;
            let cur = smalls[i / p] - pc;
            smalls[i] -= cur;
        }
    }
    larges[1]
}

fn floor_sqrt(n: usize) -> usize {
    if n <= 1 {
        return n;
    }
    let mut lo = 1;
    let mut hi = n;
    while hi - lo > 1 {
        let mid = lo + (hi - lo) / 2;
        match mid.overflowing_mul(mid) {
            (x, false) if x <= n => lo = mid,
            _ => hi = mid,
        }
    }
    lo
}

計算量解析

i に関するループがふたつありますが、ひとつめの計算量は次の値で抑えられます。以下の \(i\) は p に対応します。 \[ \begin{aligned} \int_2^{\sqrt{n}} \min\{\sqrt{n}, n/i^2\}\, di &= \int_2^{\sqrt[4]{n}} \sqrt{n}\,di + \int_{\sqrt[4]{n}}^{\sqrt{n}} n/i^2\,di\\ &= \sqrt{n}\cdot\sqrt[4]{n} + O(\sqrt{n}) + (\sqrt[4]{n}-1)\cdot\sqrt{n} \\ &= O(n^{3/4}). \end{aligned} \]

ふたつめも。 \[ \begin{aligned} \int_2^{\sqrt{n}} \max\{\sqrt{n}-i^2, 0\}\, di &= \int_2^{\sqrt[4]{n}} (\sqrt{n}-i^2)\, di + \int_{\sqrt[4]{n}}^{\sqrt{n}} 0\,di \\ &= 1/3\cdot (2n^{3/4} - 3\sqrt{n} + 1) \\ &= O(n^{3/4}). \end{aligned} \]

合わせて \(O(n^{3/4})\) です。実際には、素数の分布から \(O(n^{3/4}/\log(n))\) になっている気がします。

工夫

オーダーの改善ではないはずで、定数倍高速化だと思います。

上記のアルゴリズムで、篩う \(i\) を \(\lfloor\sqrt{n}\rfloor\) までではなく \(\lfloor \sqrt[4]{n}\rfloor\) までにしてみます。 このとき、篩われずに残っている合成数は、たかだか 3 つの素数の積で表せることがわかります*5。4 つかけると \(n\) を越えちゃうので。 このことから、全体をふたつの処理に分けて行うことが考えられます。

  • \(\lfloor\sqrt[4]{n}\rfloor\) 以下の素数で篩い、篩われた個数を引く。
  • たかだか 3 つの素数の積で表せる合成数の個数を引く。

合成数の個数を数える方について考えます。 篩われずに残っている整数（以前篩うのに使ったのは除く）を昇順に見ていきます。\(i_1, i_2, \dots\) とします。 各 \(k\) に対して、\(i_k\) の倍数である \(n\) 以下の整数の個数を足します。 また、\(k'\lt k\) に対して \(i_k\cdot i_{k'}\) の倍数である整数はすでに数えられているので、引きます。

ここで、\(i_k\) の倍数である \(n\) 以下の整数は、\(S(\lfloor n/i_k\rfloor, \lfloor\sqrt[4]{n}\rfloor) - \pi(i_k -1)\) と表せます。\(i_k i_{k'}\) の方は \(S(\lfloor n/i_k i_{k'} \rfloor, \lfloor\sqrt[4]{n}\rfloor) - \pi(i_{k'} -1)\) です。

実装

アイディアは上記の通りですが、実装上の工夫が少しあるので下で説明します。

変数の概要は次の通りです。

  • \(\texttt{smalls}[i] = S(i, p)\)
  • \(\texttt{roughs}[i]\)：消されていない \(i\) 番目の整数
    • \(1\) は残る、篩うのに使った素数は消える
    • サイズは \(\texttt{s}\) で管理される
  • \(\texttt{larges}[i] = S(\lfloor n/\texttt{roughs}[i]\rfloor, p)\)
  • \(\texttt{pc} = \pi(p-1)\)
  • \(\texttt{skip}[i]\)：\(i\) が篩われていれば true
    • ただし偶数については更新しない

\(p\) の変化に伴って更新されるので場所によって ±1 がありますが、各ループの最後で上記の定義に合うように管理します。 最初のループの開始前においては、\(p = 2\) での値です。

fn main() {
    assert_eq!(prime_pi_fast(100000000000), 4118054813);
}

fn prime_pi_fast(n: usize) -> usize {
    if n <= 3 {
        return n.saturating_sub(1);
    }
    let v = floor_sqrt(n);
    let mut smalls: Vec<_> = (0..=v).map(|i| (i + 1) / 2).collect();
    let mut s = (v + 1) / 2;
    let mut roughs: Vec<_> = (0..s).map(|i| 2 * i + 1).collect();
    let mut larges: Vec<_> =
        (0..s).map(|i| (n / (2 * i + 1) + 1) / 2).collect();
    let mut skip = vec![false; v + 1];

    let mut pc = 0;
    for p in (3..=v).step_by(2) {
        if skip[p] { continue; }
        let q = p * p;
        pc += 1;
        if q * q > n { break; }
        skip[p] = true;
        for i in (q..=v).step_by(2 * p) {
            skip[i] = true;
        }
        let mut ns = 0;
        for k in 0..s {
            let i = roughs[k];
            if skip[i] { continue; }
            let d = i * p;
            let x = if d <= v {
                larges[smalls[d] - pc]
            } else {
                smalls[n / d]
            };
            larges[ns] = larges[k] + pc - x;
            roughs[ns] = i;
            ns += 1;
        }
        s = ns;
        let mut i = v;
        for j in (p..=v / p).rev() {
            let c = smalls[j] - pc;
            let e = j * p;
            while i >= e {
                smalls[i] -= c;
                i -= 1;
            }
        }
    }

    let roughs = roughs;
    let pc = pc;

    let mut res: usize = larges[0] + (s + 2 * (pc - 1)) * (s - 1) / 2
        - larges[1..s].iter().sum::<usize>();

    for l in 1..s {
        let q = roughs[l];
        let m = n / q;
        let e = smalls[m / q] - pc;
        if e <= l { break; }
        let t: usize = roughs[l + 1..=e].iter().map(|&r| smalls[m / r]).sum();
        res += t - (e - l) * (pc + l - 1);
    }
    res
}

fn floor_sqrt(n: usize) -> usize { /* 上と同じなので省略 */ }

実装の説明

冷静になると自明かも。わからなければ、下記の説明を読みつつ、各箇所で値を適宜出力したりしながら確かめるとよいと思います。

roughs[ns] の更新

roughs[..s] が篩われていない整数です。 ns が次のループにおける s です。 ns = 0 で初期化しておき、roughs[ns] に値が入ったら ns += 1 していきます。前から詰め直しているということですね。

larges[smalls[d] - pc] についての説明

smalls[d] - pc は「篩われずに残っている整数のうち d は何番目か？」を意味します*6。

smalls[d] - pc は「d 以下にある篩われていない整数」から「篩うのに使った整数」を除いたものの個数です。 ただし、この段階ではまだ pc は \(\pi(p -1)-1\) であることに注意です。これに、\(1\) のぶんの個数を足して、所望の値であることがわかります。 roughs[それ] が篩われていない d 番目の整数であることと larges の定義から、結局 \(n/d\) 以下で篩われていない整数の個数になります。

smalls[i] の更新

j = i / p を保って更新されるようにします。 i >= j * p である間 smalls[i] -= smalls[j] - pc で更新され、smalls[i] -= smalls[i / p] - pc になっています。

res の初期化

初期化の際に、（上で説明した）\(i_k\) の倍数を先に引いています。等差数列の公式などを思い出すと吉です。

res の更新

  1. q (\(i_{k'}\)) を固定
  2. q より大きい r (\(i_k\)) を固定
  3. n / (q * r) として表せる篩われていない整数の個数を求める

\(\sqrt[4]{n} \lt i_{k'} \lt i_k\) から、\(\lfloor n / i_k i_{k'}\rfloor \le \lfloor\sqrt{n}\rfloor\) であり、m / r = n / (q * r) で smalls にアクセスしても問題ありません。

3. で求めた個数から、（前半パートですでに数えた）\(\sqrt[4]{n}\) 以下の素数を含む個数を除く処理が -= (e - l) * (pc + l - 1) です。 \(\lfloor n / i_k i_{k'}\rfloor\) として表せるうちで数えたい最小値は \(i_k'\) なので、2. で固定するのは \(\lfloor n / i_{k'}^2\rfloor\) 以下の範囲です。

計算量解析

前半パートでは、各 \(i \le n^{1/4}\) の素数に対して \(O(\sqrt{n})\) 程度の処理をするので \(O(n^{3/4}/\log(n))\) time です。篩の部分も \(O(\sqrt{n}\log(\log(n))) \subset o(n^{3/4}/\log(n))\) です。 後半パートでは、雑に見積もって \(n^{1/4} \lt i \le n^{1/2}\) の各素数に対して \(O(n/i^2)\) 程度の処理をしています。 \[ \int_{n^{1/4}}^{n^{1/2}} \frac{n}{i^2}\, di = n^{3/4} - n^{1/2}. \] であり、素数の分布から \(O(n^{3/4}/\log(n))\) になりそうです。

\(n = 10^{13}\) でも 1 秒前後くらいで終わりました。 \(\log(n)\) で割れる部分のちゃんとした（お気持ちでない）証明ができたら追記します。

その他

\(\sqrt[3]{n}\) とかで区切って \(O(n^{2/3})\) 時間にできるのかな？ 区切り方を変えれば \(O(n^{3/4-\varepsilon})\) にはできそうな気もするけど、定数倍とか諸々のことはわかりません。

Dirichlet 積などを使った方法 とかもあります。軽めの定数倍では実装できませんでした。

参考文献

おわり

おやすみなさいにゃ。

（記事を書いていたら眠れなくなりました）

*1:実装はこれとか。素因数分解の一意性をうまく使いつつ配る DP をします。

*2:素数定理から。

*3:Project Euler での Lucy_Hedgehog 氏の投稿 に由来するようです。

*4:ABC 141 D の解説 とかに証明があります。

*5:\(\lfloor\sqrt[3]{n}\rfloor\) までにすればたかだか 2 つにできます。

*6:rank/select 演算を知っている人は、それに似たことを考えるとよいと思います。

intXX_t に関して

どうしてこんなことに...

まえおき

C++ には次の型があります：

  • signed char
  • unsigned char
  • char
  • short signed int (= short)
  • short unsigned int (= unsigned short)
  • signed int (= int)
  • unsigned int (= unsigned)
  • signed long int (= long)
  • unsigned long int (= unsigned long)
  • signed long long int (= long long)
  • unsigned long long int (= unsigned long long)

char の signedness（符号つきかどうか）は処理系定義なのですが、signed char とも unsigned char とも異なる型であることが規定されています。

また、signed long long int などについては、順番はどうでもよくて、たとえば long int signed long なども同じ型として扱われます。

ともかく、signedness などを無視すると、次の 5 つがあることになります。

  • char
  • short
  • int
  • long
  • long long

ところで、扱えるビット長は 8, 16, 32, 64 は 4 種類です*1。 鳩ノ巣原理から、上記の型について、少なくとも 1 組はビット長が同じものが存在しますね。

よくある環境では、次のようになっているはずです。

  • char → 8 bit
  • short → 16 bit
  • int → 32 bit
  • long → 32 bit OR 64 bit
  • long long → 64 bit

（char が 64 bit の処理系とか、short が 32 bit の処理系とかも規格上は許されてはいますね。これらの型同士の大小関係とか、最低限のビット長くらいしか規定されていません。）

long が 64 bit の環境をみんなが使っていれば「64 bit 整数を使いたいなら long を使いましょう。おわり」にできるんですが...

多くのジャッジ環境（Codeforces は別）は long が 64 bit なので、手元環境が Windows でない人は long と書けばいい気がします。

本題

「そもそも 64 bit 整数を使っている」みたいなことを明示したいのに「long」とかの名前を使いたくないという感情はあるので、ビット長で名付けられた型を使いたい気分になります。

そこで、std::int64_t みたいなのを使うようになるのですが、これは罠が多いので、それに関する記事です。

まず、intXX_t は上記の型のどれかの typedef である（あれらと別の型として存在しているわけではない）ことに注意です。

前述の通り、long が 64 bit の処理系とそうでない処理系があるので、int64_t は long だったり long long だったりします。 このことは、意外と厄介です。

リテラルのつくりかた

では、int64_t の定数を自分で定義するときはどうすればいいでしょう。

unsigned の 0u とか、long の 0L のような suffix は用意されていません。

auto x = INT64_C(0);

INT64_C というマクロが定義されており、これを使うとよいです（めんどくさいね）。 このマクロは処理系が定義しておいてくれているもので、適宜 0L とか 0LL とかに置き換えてくれます。

int64_t(0) や (int64_t)0 のようなキャストでも別にいいです。

多用するなら、リテラルを定義しちゃいましょう。 ユーザ定義リテラルは _ で始めないと未定義動作（予約されている識別子を勝手に使ったことになる）なので、そういう感じにします*2。

int64_t operator ""_i64(unsigned long long) { return n; }

これを定義しておけば、0_i64 のようにして作ることができます。

次のように適当に決め打ちしていいじゃんと侮っていると、面倒なことになります。

int64_t x = 0LL;

std::max

もちろんこれに限るわけではないですが、ありがちなので。

int64_t x = /* 計算で出てきた値 */;
int64_t y = std::max(0L, x);  // ???

int64_t が long の処理系では 0L と x の型が同じなのでコンパイルが通りますが、long long の処理系では異なるので通りません。

Note: std::max は 2 つの引数の型が同じである必要があります。

scanf / printf

std::cin を何らかの理由で使いたがらない人向けです。

int64_t x;
scanf("%ld", &x);

これは、int64_t が long long である環境では未定義動作になるはずです（scanf はフォーマット文字列で指定される型と実際の型が異なると未定義動作になるので）。

そこで、これもやっぱりマクロが用意されていて、次のようにすると安全です。 受け取った 2 整数を順番を入れ替えてスペース区切りで出力する例です。

#include <cinttypes>
int main() {
    int64_t x, y;
    scanf("%" SCNd64 "%" SCNd64, &x, &y);
    printf("%" PRId64 " %" PRId64 "\n", y, x);
}

SCNd64 などが "ld" とかに展開されるように処理系が定義してくれています。 "%ld" となっていないのは、"%02ld" のようなのを作りやすいようにするためだと思います。

Note: C++ の文字列リテラルはくっつけると文字列の意味で連結することができ、"%" "ld" "%" "ld" と書けば "%ld%ld" と書いたのと同じことになります。

どうせ競プロで使うぶんにはジャッジ環境であってさえいればいいので、"%ld" とか書いてしまえばいいんですが、気持ち悪さを抱えながら書くことになります。 そもそも long とかを考えたくないからこういう型を使っているのに... という気分になります。

余談：intmax_t という型もあって、これについては %jd というのが使えます。

int8_t

これは std::cin / std::cout で入出力する人向けです。

前述のように、これらの型は typedef であって、「専用に用意された整数型」ではないです。 つまるところ、int8_t は char だったりします。

int8_t x = 65;
std::cout << x << '\n';

のように書くと...？ A とかが出ちゃいます。あーあ。

ACL の話

これ忘れてました。追記です。

ACL が int64_t をサポートしていないというか、long をサポートしていないというのが正しい気はします（ジャッジ環境で int64_t == long なので）。

そもそも定義されない処理系もある

こんな処理系を触ることはないでしょうので、これは気にしません。

そのほか

int_fast32_t とか int_least64_t とかについても同じような問題があります。 これらはそもそもどのビット長なのかすらわからないので、参りがちです。

int_fast32_t は別に「普段の 32 bit 整数とは異なるすごい整数型」ではなく、「処理系が適当に決めた最低 32 bit ある整数型」なので、必ずしも常に速いとすら限りません。

どうせ人々は自作ライブラリをぺたりしてるのでしょうから、自分で定義する方がいいんじゃないかなと思います。好きな名前をつけられますし。

using i32 = int;
using u64 = unsigned long long;
...

むかしの

*1:すごい処理系であれば別にこれに縛られはしないと思うのですが、とりあえず現在よくある環境たちについての話ということにします。

*2:グローバルでアンダースコアで始めるのは、リテラルなら許されるんですかね？

ICPC 2020 Asia Yokohama Regional 参加記

あー、おわりです。

えびちゃんの他の参加記は ここから*1。

えびちゃんは B3 の頃から ICPC に出ていました。 four-t で 3 回と、今回 tsutaj で出ました。チームメイトに大変恵まれていたなぁというのをとても思います。 今回で老害化です。

tsutaj（チーム）はえびちゃんと、たぶきゅんと、monkukui からなるチームです。tsutaj（個人名）は個人です。

いろいろ思い出はありますが、それは別の機会に書くことにします*2。

Day -10 くらい

システムテストコンテスト

先々週にあったやつです。

オタクなので、verdict 集めをしました*3。

f:id:rsk0315:20210317221111p:plain
verdict 集め

オタクなので、OUTPUT-LIMIT の境界を探る遊びをしました。 こういうのはどうせ 2 べき bytes だと思うので、指数部をにぶたんしました。

雑に出力すると TIMELIMIT になったり、前に作った大量出力用のコードがバグっていて期待通りでないバイト数を出力したりで少しわたわたしました。

223 bytes ではセーフで、224 bytes ではアウトでした。 224-1 ならセーフなのか、223+1 からアウトなのか、中間（10 MB ± 1 とか？）なのかで無駄な探索をしたりしましたが、結局、8 MB がセーフの上限だと突き止めました。 ジャッジに負荷がかかるかな？とか、攻撃と見なされたら嫌だな〜とか思いつつも、提出数自体が少なかったので、いいかなと思ってやっちゃいました。

f:id:rsk0315:20210317221752p:plain
一応チーム slack に報告する人

そんなことより static_assert(__cplusplus >= 201703L); とかを調べるべきなんですが、完全に忘れていました（1 敗）。

Day 1

受付

任意文字列発声脆弱性のコーナー。

文脈がないと伝わらないね（The Atama チームの受付が終わったときの様子）。

開会式

コンテストシステムのチュートリアルのコーナー。

これ言われなかった（例年、ログインの際に P-A-S-S-W-O-R-D と言いながらパスワードを打つのをやっていた気がします）。

WA を投げて clar で不満などを言うコーナーです。\(ax^2+bx+c=0\) を解いてねという問題です。

a=b=c=0 について触れられていなかったり、それにも関わらず CORRECT になったり、それで人々がわくわくしていたりして、面白かったです。

プラクティス

この前の日くらいに Judging Details が公開されて、OLE が > 4 MB と書かれていて、ほんとかな〜と思っていました。 念のため OLE 境界のコードふたつを投げて確認しました。

clar を投げたら対応してもらえてうれしかったです。

f:id:rsk0315:20210317222356p:plain
なけなしの英語力さん

f:id:rsk0315:20210317222439p:plain
A new announcement about this issue

別に深刻な問題ではないと思うんですが、不備を報告して対応してもらうとうれしくなりますね。

うれしくなっていたら static_assert(__cplusplus >= 201703L); のことを忘れていました（合わせて 1 勝 2 敗）。

Day 2

本番

A 困っていましたね。

B 貪欲？ いや、できなくない？ DP しますか？ 一生バグらせた後たぶきゅんに丸投げ。たぶきゅんも一生バグらせた後、えびちゃんが貪欲っぽいのを書いて CORRECT。う〜〜ん、ごめん。

monkukui が J を通していてえらい。

たぶきゅんが、「（概要の説明）〜〜（考察）〜〜（ここまで数秒）...ので、G を書きます。」とシームレスに書き始めていたのがよかったです。かっこよし。

この次たぶきゅんが H を書いてたかも？ なんかペースが速い。

たぶきゅんが書いている間に slack で E の概要を伝えると、たぶきゅんが「2 年前の JAG 夏で見覚えある」と言っていたので、探して貼ったりしました。E も通していました。

monkukui が I の DP をがんばっていて、デバッグとかをいろいろがんばったけど、そこで終わり。

後々聞くところによれば、C に手をつけておく方がよかったかも。う〜〜ん、大くやしい。

YES_No

f:id:rsk0315:20210317223818p:plain
YES_No 表記の元ネタ（？）

eat_ice、盛り上げ界の tourist。

任意文字列発声の脆弱性（失敗しがち）イベント。つたじ（かなしい）。

Yes を言わせられなかったのくやしいね。

お顔大公開のコーナーで、お顔が大公開されて恥ずかしかったです。 んわーとなってタブを閉じてしまいました。 後で懇親会のときに、かわいいと言ってもらえたのでうれしかったです。

えらい先生が「Yes/No おじさん」「___ KING ___ 様」などと呼んでいたのが面白かったです。 優勝コメントも面白かったですね。

オタクなので気になるんですが、___KING___ ではなく ___ KING ___ な気がします。

結局 6 完 20 位で、中央値よりは真に上でしたが、そういう意味では去年と同じなので、うむむむという感じです。

Closing Party

___ KING ___ 様がスタッフ様と話しているのを聞いていると、（たぶん F 問題？）「あれが難しいつもりというのがかなりびっくりで〜」のようなことを言っていて、ひえーとなりました。

いろいろな人と話したり、いつもの勢と話したり、びーとくんや olphe をつんつんしたりしました。

みさわさんにつんつんしてもあまり話してくれなかったりした*4のですが（マイクが不調だったりしたのかな？）、なんやかんやで、フローの講義をしてくれました。

climpet さんがえびちゃんについてきて「あなたと JAG」とだけ言ってきたのがちょっと怖かったです。

この辺からねむくなったり、おなかが空いたりしていました。

「朝早かったからね〜」と言われたのでつい「バカがつくったスケジュール」と口走ったら、なんかきゃっきゃされました。

近くにえらい先生がいて危ない気がしたので、念のためフォロー（？）のつもりで元ネタの紹介をしました。 これが最古じゃなかったらすみません。

楽しい時間はすぐ終わり。

おわり

あーー。じわじわ実感してきています。

えびちゃんと組んでくれた人たちありがとうです。 えびちゃんの力不足やオタク言動などで負担をかけたかと思いますが、長いこと組んでくれてうれしいです。 えびちゃん的には、このチーム（four-t / tsutaj）に身を置けたのが幸せです。

忙しさから、チーム練が疎かになったりしがちだったので申し訳ないです。

tsutaj は全員今年で引退なのですが、来年以降も弊学の若い人たちにがんばってほしいです。 チーム名に困ったら rsk0315 を使ってもいいですよ（？？）。 と思ったんですが、Yes/No おじさんが読むのに苦労しそうなので、TAB や monkukui などを使う方がいい気がします。

ステマ？

JAG に入るなどして、ちゃんと恩返しをしなきゃですね。

*1:他のと言いつつたぶんこの記事も含まれるんですが。

*2:書いているうちに今日の記憶が曖昧になると困るので。

*3:TOO-LATE は別で集めました。

*4:みさわさんが好きだけど, 好かれてはいない. か？

自分で未定義動作を書いて逆ギレする人かわいそう

お店に入るときにドアが開いているか確認せずに歩いていって、ドアにぶつかって怒っている人がいたらおかしい話じゃないですか*1。

未定義動作を起こして怒るのはこれと似ている気がします。 競プロでよくあるコーナーケースなども同じです。

「いま書いている処理は、ありえる入力のどれに対しても適切に動くか？」というのをちゃんと意識するといい気がします。

昔の記事 では、未定義動作を書いたらどうなりうるかがメインだった気がします。 今回は、競プロの人々がよくやっている or やりそうな未定義動作の紹介集みたいなのに近いつもりです。

「複雑な書き方をするとこわれます」というのはまた今度書くとして、「簡単にこわれます」というのがメインです。

以下では、int は 32 bit の 2 の補数表現ということで書きます。AtCoder など多くの環境と同じです。

演算子

多くの演算子、入力によってはすぐ未定義動作を起こしてしまう...（かなしい）

+

#include <climits>
int main() {
  int bad = INT_MAX + 1;
}

オーバーフローはだめ、それはそう。

-

#include <climits>
int main() {
  int bad = INT_MIN - 1;
}

負のオーバーフローも気をつけましょう。

*

#include <climits>
int main() {
  int bad = INT_MAX * 2;
  int bad2 = INT_MIN * (-1);
}

オーバーフロー定期。

/

#include <climits>
int main() {
  int bad = 1 / 0;
  int bad2 = INT_MIN / (-1);
}

ゼロ除算は当然気をつけましょう。 オーバーフローにも気をつけましょう。

%

#include <climits>
int main() {
  int bad = 1 % 0;
  int bad2 = INT_MIN % (-1);
}

INT_MIN % (-1) は 0 でしょうと思うんですが、x / y がオーバーフローするときは x % y も未定義です。

<<

int main() {
  int bad = 1 << 100;
  int bad2 = 1 << (-5);
}

x << y は、y が負の場合や、x の（汎整数拡張後の）ビット幅以上の場合は未定義です。

[]

int main() {
  int a[5] = {31, 41, 59, 26, 53};
  int bad = a[5];
}

配列外参照もだめ。

標準関数

<cmath> の関数とかは予想に反しがちなときがあるかも？ 適宜調べておきましょう。

en.cppreference.com

<algorithm> の関数、特に sort など順序が絡むものは誤る人が多いので気をつけましょう。 「順序としてこわれているものを第三引数に渡しちゃお〜」とすると、無限ループになったり Segfault になったりします*2。

GCC 拡張の関数

gcc.gnu.org

__builtin_clz __builtin_ctz

x のうち先頭・末尾の 0 のビットの個数を数えるものです。 x == 0 のとき（すなわち 1 が現れないとき）の動作は未定義です。

ビット幅を返してくれるとは言っていません。

飽きたよね

えびちゃんも飽きました。

何をやるにしても未定義動作が絡んできて、覚えておくのが大変すぎてつらいというのは、えびちゃんもわかります。 「こんな処理が未定義じゃなかったら、やばいだろ」というような感覚を養うといいかもしれません。

「未定義動作なんて非本質すぎて考えたくない」と思ってしまうタイプの人には、C++ 以外の言語が向いているかもしれません？

ねこ

典型コーナーケースに対する「自分で場合分けを放棄して WA に逆ギレする人かわいそう」という記事や、「自分で定数倍の遅いコードを書いて TLE に逆ギレする人かわいそう」のような記事も書こうとしていました。「誤読して逆ギレする人かわいそう」という気持ちはありますが、記事に書くことはあまりないような気もします？

人々が何らかの対象に逆ギレしているのを見て、お気持ち表明したくなったらまた書きます。

*1:え、もしかして実在しますか？

*2:未定義なので、もちろん、こうなることは保証されていませんが。