2023-11-07

u64 の平方根を f64 で計算するやつに関して

浮動小数点数

符号なしの 64-bit 整数型の値に関して平方根を計算したい局面はそこそこあります。しかし、多くの言語ではそのための関数が用意されておらず、浮動小数点数型用のメソッドを使う人が多くいます。

今回は、そうした方針が妥当か？（誤差が出うることを踏まえつつ、正しい値を復元するのは簡単か？）という話です。

記法・前提
考察
実装例
あとがき
おわり

記法・前提

浮動小数点数型に対する一般的な知識はあるとします。ない人向けの記事はそのうち書きます。

実数 $x$ を該当の浮動小数点数型の値に丸めたものを $\roundp x$ と書くことにします。

IEEE 754 の規格において、sqrt(x) は $\roundp{\sqrt x\,}$ を返すことを規定しているので、ここではそれを仮定します。

note: 他の関数の誤差に関して

GCC において、sin(x), log(x), cbrt(x) などは、$\roundp{\sin(x)}$, $\roundp{\ln(x)}$, $\roundp{\sqrt[3]{x}\,}$ を返すとは限りません*1。sqrt などは IEEE 754 の規格で required の枠に含まれるのに対し、sin などは recommended の枠に含まれるので、建前上は「recommended の方は提供していませんよ〜」ということなのでしょうか。今後、正確に丸めた (correctly-rounded) $\roundp{\sin(x)}$ を返す crsin(x) のような関数を提供するかも？とのことです。

www.gnu.org

任意の正規化数 $x$ に対し、$2^e \le x\lt 2^{e+1}$ が成り立つような $e$ が一意に存在します。この $e$ を用いて、$\hfloor x=2^e$ と定義します。 $\hfloor x\le x\lt 2\hfloor x$ や $\tfrac x2\lt\hfloor x\le x$ が成り立ちます。また、$x$ が正整数の下では $\hfloor x\le x\le 2\hfloor x-1$ や $\tfrac{x+1}2\le\hfloor x\le x$ が成り立ちます。

下記では、次のような関数を作ることを想定します。

/// $\floor{\sqrt n}$ を返す。
fn uint_floor_sqrt(u64 n) -> u64 {
    f64 tmp = (n as f64).sqrt();
    // tmp が持ちうる誤差に対して何かしらの補正をする
    todo!()
}

どのように補正を行うかについてはまだ与えないことにします（適当に考えた方針が「実はよかったよね」とするのではなく、考察ベースで自然なアルゴリズムを考えたいため）。考察次第では、.sqrt() をする前に何かしらをするかもしれません。

？？「仮数部が 64 bits 以上の long double がある環境なら long double でやればよくないか？」

それはそう... そうした型のない言語を想定しているとかの建前でやり過ごします。

考察

note: $n \le 2^{53}$ では $\roundp n=n$ ですが、それより大きい場合はそうとは限らないため、気をつける必要があります。

$\gdef\roundsqrt#1{\roundp{\sqrt{#1}\,}}$

Observation 1:

u32 に関してはすぐに全探索できるため、その範囲を f32（仮数部 24 bits）で計算してみます。

fn main() {
    for i in 0..=u16::MAX as u32 {
        for j in i * i..=i * i + 2 * i {
            let r = (j as f32).sqrt() as u32;
            assert_eq!(r, i, "failed on {}", j);
        }
    }
    println!("passed all tests.");
}

assertion `left == right` failed: failed on 16785407
  left: 4097
 right: 4096

やはり、これはどうやらうまくいかないようです。 $$ \floor{\roundsqrt{(2^{12}+1)^2-2}} = 2^{12}+1. $$

note: f64 においても、$\floor{\sqrt n}$ を固定して二分探索を使ったり、固定した際の最大値を試してから降順に見ていくなどで、同様の探索ができます。最小の反例は 4503599761588224 で、下記のようになっています。 $$ \floor{\roundsqrt{(2^{26}+1)^2-1}} = 2^{26}+1. $$ $(2^{\floor{p/2}}+1)^2$ より少し小さいあたりでだめになってしまう感じでしょうか。

Observation 2:

では、$+1$ の誤差は許容してみましょう。

fn main() {
    for i in 0..=u16::MAX as u32 {
        for j in i * i..=i * i + 2 * i {
            let r = (j as f32).sqrt() as u32;
            assert!([i, i + 1].contains(&r), "failed on {}", j);
        }
    }
    println!("passed all tests.");
}

passed all tests.

これはよいらしいです。

Observation 3:

ということで、u64 にして次のような関数を考えてみます。

fn u64_floor_sqrt(n: u64) -> u64 {
    let tmp = (n as f64).sqrt() as u64;
    let tmp_m1 = tmp.saturating.sub(1); // if tmp >= 1 { tmp - 1 } else { tmp };

    // ((tmp - 1) + 1)**2 <= n, in an overflow-safe way
    if tmp_m1 * (tmp_m1 + 2) < n { tmp } else { tmp_m1 }
}

$n = 0$ の場合や、$n = 2^{64}-1$ で $2^{32}$ が返ってきた場合などにオーバーフローしないように気をつけています。そうしたケースに気を払う必要がない状況では、もう少し雑にやっても大丈夫そうです。

fn main() {
    for i in 0..=u32::MAX as u64 {
        let actual_left = u64_floor_sqrt(i * i);
        let actual_right = u64_floor_sqrt(i * i + 2 * i);
        let expected = i;
        assert_eq!(actual_left, expected);
        assert_eq!(actual_right, expected);
    }
    println!("passed all tests.");
}

passed all tests.

これもよいらしいです。$\floor{\sqrt n}-1$ や $\floor{\sqrt n}+2$ になることはなさそうです。

$\roundp x$ や $\sqrt x$ の単調性から、u64_floor_sqrt は正当そうです。これで「u64 の範囲ではこれでいけます、めでたし」で終わってもいいにはいいんですが、たとえば C++ において unsigned __int128 の範囲を __float128 で計算するのを同じ方針で調べることはできないため、もう少し数学的に考えます。

もしかすると 128-bit では反例が見つかるかもしれないし、そうではないかもしれませんから。

ここまでの話をまとめると、次のようなことを示したいです。

Conjecture 4: それなりに大きい $n$ に対しても $\floor{\sqrt n} \le \roundsqrt{\roundp n} \lt \floor{\sqrt n}+2$ が成り立つ。

まず、浮動小数点数 $\hat x$ と実数 $x$ に対して、$\roundp x = \hat x$ となる十分条件について触れておきます。 $$ \roundp x=\hat x \impliedby \begin{cases} \hat x-\hat x\cdot 2^{-1-p} \le x\le \hat x+\hat x\cdot 2^{-p}, & \text{if }\hfloor{\hat x}=\hat x; \\ \hat x-\hfloor{\hat x}\cdot 2^{-p}\lt x\lt\hat x+\hfloor{\hat x}\cdot 2^{-p}, & \text{if }\hfloor{\hat x}\gt\hat x. \end{cases} $$ $\roundp x\ge \hat x$ などを示したいときは片側の不等式が成り立てばよいです。

$0\le k\le 2^{p-1}$ なる整数 $k$ に対して、整数 $n\in[k^2\ldots(k+1)^2)$ が conjecture 4 を満たすことを示します。

整数性から $n\in[k^2\ldots k(k+2)]$ であることと、単調性から $$\roundsqrt{\roundp{k^2}}\le\roundsqrt{\roundp n}\le \roundsqrt{\roundp{k(k+2)}}$$ であることに注意します。そのため、 $$ k\le \roundsqrt{\roundp{k^2}} \wedge k+2\gt \roundsqrt{\roundp{k(k+2)}} $$ を示せば十分です。考えている状況に合わせ、$\roundp{k^2}\lt\infty$ である前提とします。また、この範囲で $\roundp k=k$ が成り立つことに注意しておきます。加えて、$p\ge 3$ も仮定しておくことにします。

Lemma 5: $0\le k\le 2^p$ なる整数 $k$ に対して、$k\le \roundsqrt{\roundp{k^2}}$ が成り立つ。

disclaimer: もともと $k\le 2^p$ で示す予定だったのですが、Lemma 6 が $k\le 2^{p-1}$ でしか成り立たなさそうということが Lemma 5 を示した後にわかりました。Lemma 5 を直すのが面倒だったのでそのままになっています。

Proof

ある $e$ が存在して $k = 2^e$ のとき、$2^{2e} = 2^{p-1}\times 2^{2e+1-p}$ と表せるため、$\roundp{2^{2e}}=2^{2e}$ となる。 $\roundsqrt{\roundp{k^2}} = \roundsqrt{k^2} = \roundp k=k$ が成り立つ。すなわち、$k\le \roundsqrt{\roundp{k^2}}$ となる。

$k = 0$ のときも同様にして成り立つことが示せる。

次に、$k = 2^e$ なる $e$ が存在しないかつ $k\ge1$ とし、 $$ k-\hfloor{k}\cdot 2^{-p}\lt \sqrt{\roundp{k^2}} $$ を示す。これらを評価すると、 $$ \begin{aligned} k-\hfloor k\cdot 2^{-p} &\le k-\tfrac{k+1}2\cdot 2^{-p} \\ &= (1-2^{-1-p})k - 2^{-1-p}; \\ \sqrt{\roundp{k^2}} &\ge \sqrt{k^2-\hfloor{k^2}\cdot 2^{-p}} \\ &\ge \sqrt{k^2-k^2\cdot 2^{-p}} \\ &= k\cdot \sqrt{1-2^{-p}} \end{aligned} $$ となるので、$(1-2^{-1-p})k-2^{-1-p}\lt k\cdot\sqrt{1-2^{-p}}$ を示せばよい。非負であることは明らかであるから二乗の差を考え、 $$ \begin{aligned} &\phantom{{}={}} (k\cdot\sqrt{1-2^{-p}})^2 - ( (1-2^{-1-p})k-2^{-1-p})^2 \\ &= (1-2^{-p})k^2 - ( (1-2^{-1-p})k-2^{-1-p})^2 \\ %&= (1-2^{-p})k^2 - ( (1-2^{-1-p})^2k^2-2^{-p}\cdot(1-2^{-1-p})k+2^{-2-2p}) \\ %&= (1-2^{-p})k^2 - ( (1-2^{-p}+2^{-2-2p})k^2-2^{-p}\cdot (1-2^{-1-p})k+2^{-2-2p}) \\ %&= -(2^{-2-2p}\cdot k^2-2^{-p}\cdot (1-2^{-1-p})k+2^{-2-2p}) \\ %&= -2^{-2-2p}\cdot k^2+2^{-p}\cdot (1-2^{-1-p})k-2^{-2-2p} \\ %&= 2^{-p}\cdot (-2^{-2-p}\cdot k^2 + (1-2^{-1-p})k - 2^{-2-p}) \\ %&= 2^{-2p}\cdot (-2^{-2}\cdot k^2 + (2^p-2^{-1})k - 2^{-2}) \\ &= 2^{-2-2p}\cdot (-k^2 + (2^{2+p}-2)k - 1) \\ \end{aligned} $$ となる。

$-k^2+(2^{2+p}-2)k-1$ の部分が、主要項が負の $2$ 次関数になっていることに注意すると、これが $k=1$ および $k=2^p$ で正になっていることを示せば、$1\le k\le 2^p$ の範囲で正であることがわかる。 $$ \begin{aligned} -1^2 + (2^{2+p}-2)\cdot 1-1 &= 4\cdot(2^p-1); \\ -(2^p)^2 + (2^{2+p}-2)\cdot 2^p-1 %&= -2^{2p} + 2^{2+2p} - 2^{1+p} - 1 &= 2^p\cdot(3\cdot 2^p-2)-1 \end{aligned} $$ より、$p\ge 1$ でこれらは正となることがわかる。

以上より、$0\le k\le 2^p$ なる任意の整数 $k$ に対して $k\le\roundsqrt{\roundp{k^2}}$ が成り立つ。$\qed$

Lemma 6: $0\le k\le 2^{p-1}$ なる整数 $k$ に対して、$k+2\gt \roundsqrt{\roundp{k(k+2)}}$ が成り立つ。

Proof

任意の $0\le k\le 2^{p-1}$ に対して、 $$ \sqrt{\roundp{k(k+2)}} \lt (k+2)-\hfloor{k+2}\cdot 2^{-p} $$ を示せばよい。両辺を評価して $$ \begin{aligned} \sqrt{\roundp{k(k+2)}} &\le \sqrt{k(k+2)+\hfloor{k(k+2)}\cdot 2^{-p}} \\ &\le \sqrt{k(k+2)+k(k+2)\cdot 2^{-p}} \\ &= \sqrt{k(k+2)}\cdot \sqrt{1+2^{-p}}; \\ (k+2)-\hfloor{k+2}\cdot 2^{-p} &\ge (k+2)-(k+2)\cdot 2^{-p} \\ &= (k+2)\cdot (1-2^{-p}) \end{aligned} $$ となるので、 $$ \sqrt{k(k+2)}\cdot\sqrt{1+2^{-p}} \lt (k+2)\cdot(1-2^{-p}) $$ を示す。Lemma 5 同様にして、 $$ \begin{aligned} &\phantom{{}={}} ( (k+2)\cdot(1-2^{-p}))^2 - (\sqrt{k(k+2)}\cdot\sqrt{1+2^{-p}})^2 \\ &= (k+2)^2\cdot(1-2^{1-p}+2^{-2p}) - k(k+2)\cdot(1+2^{-p}) \\ &= 2^{-2p}(k+2)\cdot\left( (2^{2p}-2^{1+p}+1)(k+2)-(2^{2p}+2^p)k\right) \\ %&= 2^{-2p}(k+2)\cdot\left( (2^{2p}-2^{1+p}+1)k+2(2^{2p}-2^{1+p}+1)-(2^{2p}+2^p)k\right) \\ %&= 2^{-2p}(k+2)\cdot\left( (-2^{1+p}+1)k+2(2^{2p}-2^{1+p}+1)-2^p\cdot k\right) \\ &= 2^{-2p}(k+2)\cdot\left( (-2^{1+p}-2^p+1)k+2(2^{2p}-2^{1+p}+1)\right) \\ \end{aligned} $$ となる。$(-2^{1+p}-2^p+1)k+2(2^{2p}-2^{1+p}+1)$ の部分が、主要項が負の $1$ 次関数であることに注意すると、これが $k=0$ および $k=2^{p-1}$ で正であることを示せば、$0\le k\le 2^{p-1}$ の範囲でこれが正であることがわかる。$p\ge 3$ の前提で、 $$ \begin{aligned} &\phantom{{}={}} (-2^{1+p}-2^p+1)\cdot 0+2(2^{2p}-2^{1+p}+1) \\ &= 2^{1+2p}-2^{2+p}+1; \\ &\phantom{{}={}} (-2^{1+p}-2^p+1)\cdot 2^{p-1}+2(2^{2p}-2^{1+p}+1) \\ &= 2^{1+2p}+2^{p-1} - (2^{2p}+2^{2p-1}+2^{2+p}) + 2 \\ &\ge 2^{1+2p}+2^{p-1} - 2^{1+2p} + 2 \\ &\gt 2^{p-1}+2 \end{aligned} $$ より、これらは正であることがわかる。

すなわち、$0\le k\le 2^{p-1}$ なる整数 $k$ に対して、$k+2\gt \roundsqrt{\roundp{k(k+2)}}$ が成り立つ。$\qed$

よって、Lemmata 5–6 より、$0\le n\le 2^{2(p-1)}$ なる整数 $n$ に対して $$\floor{\sqrt{n}} \le \floor{\roundsqrt{\roundp n}} \le \floor{\sqrt{n}}+1$$ が成り立ちます。すなわち、上記の関数の正当性が示されたことになります。

note: (n as f64).sqrt() as u64 で計算しているのが $\floor{\roundsqrt{\roundp n}}$ です。n as f64 が $\roundp n$ に、x.sqrt() が $\roundsqrt x$ に、r as u64 が $\floor r$ に対応します。

実装例

Rust

fn u64_floor_sqrt(n: u64) -> u64 {
    let tmp = (n as f64).sqrt() as u64;
    let tmp_m1 = tmp.saturating.sub(1);
    if tmp_m1 * (tmp_m1 + 2) < n { tmp } else { tmp_m1 }
}

C++

uint64_t uint64_floor_sqrt(uint64_t n) {
    if (n == 0) return 0; 
    uint64_t tmp_m1 = std::sqrt(n) - 1;
    return tmp_m1 * (tmp_m1 + 2) < n ? tmp_m1 + 1 : tmp_m1;
}

Python

def int_floor_sqrt(n: int) -> int:
    tmp = int(n**0.5)
    return tmp if tmp**2 <= n else tmp - 1

下記のケースでテストしています。

$n$	$\floor{\sqrt n}$
$0$	$0$
$1$	$1$
$(2^{26}+1)^2-1$	$2^{26}$
$(2^{26}+1)^2$	$2^{26}+1$
$2^{64}-1$	$2^{32}-1$

あとがき

つかれました。おなかがすきました。

証明は一旦 iPad で下書きをしているのですが、組版する段階になって符号が間違っていることに気づいたり、指数部の符号も間違っていたり、示したい命題が違っていたりで散々でした。下書きの証明をそのまま写すのではなく、組版しながらまた別途考えているので、「あれ、こんな式書いた記憶なくない？」というところで気づいたりしています。

さておき、ある程度は浮動小数点数とも仲よくなってきたんじゃないかなという気持ちがあります。初心者の頃、「平方根を sqrt(n) で求めるのは誤差が出ちゃう」と聞いて、返り値の ±30 くらい前後を探すようなコードを書いていたのが懐かしいです。

また、今回の話を踏まえると、素数判定のときにループの上限を i <= sqrt(n) とするのは正当そうということもわかりますね。約数列挙についてはどうでしょう。あまり考えていませんが、少なくとも u64・f64 で全探索した限りでは、$n=i\cdot(i+1)$ のとき $\floor{\roundsqrt{\roundp n}}=i$ となるようだったので、重複して列挙されることはなさそうな気がします。

いかがでしたでしょうか。今回の話は、ここ数回の中でも競プロ文脈に比較的近そうな気がします。とはいえ、この記事を見て「なるほど！」となって、rated のコンテストでも f64 ベースの実装で平方根を求める人がいたらなかなか勇気があるなという気もします。そんなこともないかも。証明が嘘だったときに責任を取るのは嫌なので、各自で証明し直してほしい気はします。

当初はこの記事は「書く予定リスト」に含まれていなかったため、予告していた「浮動小数点数初心者向け」の記事は次の次になりそうです。ごめんね。

おわり

おわりです。

*1:「有限桁なので、数学的な厳密な値からは誤差が出て当然」という話ではなく、「無限精度で計算した値を正確に丸めた値が欲しいが、それとは異なる値が返ってきうる」という話です。

2023-11-04

(1.0 / 49.0) * 49.0 < 1.0 について

浮動小数点数

仮数部の精度が 53 bits の浮動小数点数型で計算すると下記のようになります。

>>> 1 / 49 * 49
0.9999999999999999
>>> f'{1/49*49:.100g}'
'0.99999999999999988897769753748434595763683319091796875'

今回はこれについて掘り下げます。

記法・前提知識
解説
疑問・興味
- 疑問 1 の解答
- 疑問 2 の解答
あとがき
おまけ
おわり

記法・前提知識

その文脈における浮動小数点数型と丸め方向において、実数 $x$ を丸めた値を $\roundp{x}$ と書くことにします。型としては double（IEEE 754 で言うところの binary64）を用いて、丸め方向に roundTiesToEven を採用している言語・処理系が多いので、ここではそれを考えます。

知らない人向け：表せるもののうちで一番近いものに丸めるモード。この記事の内容では、ちょうど中間にあるケースは存在しないので、そういうことは気にしなくて大丈夫。詳しくは別の記事で解説します。

また、浮動小数点数 $x$, $y$ に対して、$x\otimes y \triangleq \roundp{x\times y}$ および $x\oslash y\triangleq \roundp{x\div y}$ とします。

たとえば、上記の例より $$ \begin{aligned} \roundp{\tfrac1{49}}\otimes 49 &= {\footnotesize 0.999999999999999888}{\scriptsize 97769753748434595763}{\tiny 683319091796875} \\ \end{aligned} $$ です。また、簡単な例として $\roundp1=1$ もあります。

$\gdef\hlog#1{\log_2\hfloor{#1}}$

正の実数 $x$ に対し、整数 $i$ を用いて $2^i$ と書ける有理数のうち、$x$ 以下で最大のものを $\hfloor x$ で表します。たとえば $\hfloor{3.5} = 2$、$\hfloor1 = 1$、$\hfloor{0.4} = 0.25$ などです。下記では「正の実数 $x$ に対し、$2^e \le x \lt 2^{e+1}$ なる唯一の整数 $e$ を考えたい」という局面が多く、これを $\hlog{x}$ と書けて便利です。

$\hlog{x} = i$ のとき、$2^{i+1-p}$ の位までが保持され、$2^{i-p}$ 以下の桁は丸められます。 double においては $p=53$ です。

解説

まず、1 / 49 の部分を考えます。

$$ \begin{aligned} \roundp{\tfrac1{49}} &= {\footnotesize 0.020408163265306}{\scriptsize 12082046367561360966}{\tiny 64704382419586181640625} \\ &= \texttt{14e5e0a72f0539}_{(16)} \times 2^{-58} \\ &= \tfrac1{49}\cdot(1-23\cdot 2^{-58}) \end{aligned} $$ となります。10 進の小数表記・16 進の指数表記と、誤差に着目した形式を挙げました。ここでは一番下の表記で考えるのがわかりやすいと思いますが、いくつかの表現方法があることを覚えておき、状況に応じて適切なものを選べるとよいでしょう。数学寄りの文脈で解析をするときは一番下の形が処理しやすい気がします。

さて、$\roundp{\tfrac1{49}}\otimes49$ を考えます。 $$ \roundp{\tfrac1{49}}\otimes 49 = \roundp{\roundp{\tfrac1{49}}\times 49} = \roundp{1-23\cdot 2^{-58}}. $$ $\tfrac12 \lt 1-23\cdot 2^{-58} \lt 1$ より、$\hlog{\roundp{\tfrac1{49}}\otimes 49}=-1$ とわかります。すなわち、$2^{-53}$ の位まで丸められます。

誤差に関して、$23\cdot 2^{-58} = \tfrac{23}{16}\cdot 2^{-54}$ より、$1\cdot 2^{-54} \lt 23\cdot 2^{-58} \lt 2\cdot 2^{-54}$ とわかります。よって、$\roundp{\tfrac1{49}}\times49$ は $1$ よりも $1-2^{-53}$ に近く、$1-2^{-53}$ に丸められます。

疑問・興味

ここで、いくつかの疑問が自然に出てくると思います。

$\roundp{\tfrac1x}\otimes x \gt 1$ となるような $x$ はある？
- ただし $\infty$ になるケースは除く
$x\in[1\ldots 2)$ なる $x$ のうちで $\roundp{\tfrac1x}\otimes x\lt 1$ となる最小のものは？
- すなわち $\hlog x=0$ の範囲
- $49$ の例から $x = \tfrac{49}{32}$ でも成り立つことはわかるので、存在性はわかる

ネタバレ

$\texttt{10000004e62386}_{(16)}\times 2^{970}$ など
- $\roundp{\tfrac1x}$ が正規化数になる範囲では存在しない
$\texttt{1000000f5cbf2a}_{(16)}\times 2^{-52}$

これについて考察していきます。

疑問 1 の解答

$\gdef\emin{e_{\text{min}}}$ $\gdef\emax{e_{\text{max}}}$

まず、正規化数の最小値を $1\times 2^{\emin}$、正規化数の最大値を $(2-2^{1-p})\times 2^{\emax}$ とおきます。 double においては、$1\times 2^{-1022}$ および $(2-2^{-52})\times 2^{1023}$ です。一般に、$|\emin|+1 = \emax$ が成り立ちます。

証明などは折りたたんでいます。

Claim 1: $x\in(1\ldots 2)$ のとき、$\roundp{\tfrac1x}\otimes x\in \{1-2^{-p}, 1\}$ が成り立つ。

Proof

このとき $\tfrac1x\in(\tfrac12\ldots1)$ であり、$\hlog{\tfrac1x}=-1$ である。そのため、ある $|\varepsilon|\le 2^{-1-p}$ を用いて $$ \roundp{\tfrac1x} = \tfrac1x + \varepsilon $$ と書ける。よって、 $$ \begin{aligned} \roundp{\tfrac1x}\times x &= (\tfrac1x + \varepsilon)\times x \\ &= 1 + \varepsilon x \end{aligned} $$ となる。$\varepsilon$ と $x$ の範囲から、$|\varepsilon x|\lt 2^{-p}$ とわかる。

$1$ と隣り合う二数は $1-2^{-p}$ と $1+2^{1-p}$ であり、$1$ に丸められる範囲は $$[1-2^{-1-p}\ldots 1+2^{-p}]$$ であるから、$1+\varepsilon x\lt 1+2^{-p}$ より $\roundp{\tfrac1x}\otimes x\le 1$ となる。同様にして $\roundp{\tfrac1x}\otimes x\ge 1-2^{-p}$ もわかる。$\qed$

Lemma 2: 浮動小数点数 $x\in(1\ldots2)$ に対し、$\roundp{\tfrac1x}\in(\tfrac12\ldots1)$ が成り立つ。

Proof

$1\oslash(1+2^{1-p})\lt 1$ および $1\oslash(2-2^{1-p})\gt\tfrac12$ を示せば十分。

$$ \begin{aligned} 1\div(1+2^{1-p}) &= 1 - 2^{1-p} + 2^{2(1-p)}\cdot (1+2^{1-p})^{-1} \\ &= 1 - 2^{-p} - 2^{-p} + 2^{-p}\cdot 2^{2-p}\cdot (1+2^{1-p})^{-1} \\ &= (1-2^{-p}) - 2^{-p}\cdot (1 - 2^{2-p}\cdot (1+2^{1-p})^{-1}) \\ &= (1-2^{-p}) - 2^{-p}\cdot (1 - 2\cdot (2^{p-1}+1)^{-1}) \\ &\lt 1-2^{-p}; \\ 1\div(2-2^{1-p}) % &= \tfrac12 \div (1-2^{-p}) \\ &= \tfrac12 \cdot (1+2^{-p}+2^{-2p}\cdot(1-2^{-p})^{-1}) \\ &\gt \tfrac12 \cdot (1+2^{-p}) \end{aligned} $$ より従う。$\qed$

Claim 3: $e \in[\emin\ldots\emax-2]$ に対して、$x\in(2^e\ldots2^{e+1})$ のとき、$\roundp{\tfrac1x}\otimes x\in \{1-2^{-p}, 1\}$ が成り立つ。

Proof

ある浮動小数点数 $\alpha\in(1\ldots2)$ を用いて $x = \alpha\times 2^e$ と書ける。また、Lemma 2 より $\roundp{\tfrac1\alpha}\in(\tfrac12\ldots1)$ が成り立つことに注意しつつ $$ \begin{aligned} \roundp{\tfrac1x} &= \roundp{\tfrac1\alpha\times 2^{-e}} \\ &= \roundp{\tfrac1\alpha}\times 2^{-e} \\ &= \roundp{\tfrac2\alpha}\times 2^{-1-e} \end{aligned} $$ が成り立つ。ここで $\emin\le -1-e\le \emax$ が成り立っていることに注意せよ。さて、 $$ \begin{aligned} \roundp{\tfrac1x}\times x &= (\roundp{\tfrac1\alpha}\times 2^{-e})\times (\alpha\times 2^e) \\ &= \roundp{\tfrac1\alpha}\times \alpha \end{aligned} $$ が成り立つので、Claim 1 から従う。$\qed$

Conjecture 4: $x\in(2^{\emax-1}\ldots 2^{\emax+1})$ のとき、$\roundp{\tfrac1x}\otimes x\in\{1-2^{-p}, 1\}$ が成り立つ。

Disproof

double において、$x = \texttt{10000004e62386}_{(16)}\times 2^{970}$ とすると $$\roundp{\tfrac1x}\otimes x = 1+2^{-52}\gt 1$$ が成り立つ。

$\tfrac1x\lt 2^{-(\emax-1)} = 2^{\emin}$ より、正規化数の最小値以下に丸められるので、ある $|\varepsilon|\le 2^{\emin-p}$ を用いて $$ \roundp{\tfrac1x} = \tfrac1x + \varepsilon $$ と書ける。よって、 $$ \begin{aligned} \roundp{\tfrac1x}\times x &= (\tfrac1x+\varepsilon)\times x \\ &= 1 + \varepsilon x \end{aligned} $$ となる。$|\varepsilon x| \lt 2^{\emax+1+\emin-p} = 2^{2-p}$ でしか抑えられず、適切な $x$ を選ぶことで反例を作れそうである。

Claim 5: $e\in[-(\emax+1)\ldots\emin)$ に対して、$x\in(2^e\ldots 2^{e+1})$ のとき、$\roundp{\tfrac1x}\otimes x\in\{1-2^{-p}, 1\}$ が成り立つ。

Proof

このとき $\tfrac1x\in(2^{-(e+1)}\ldots2^{-e})$ であり、$\hlog{\tfrac1x}=-(e+1)$ である。そのため、ある $|\varepsilon|\le 2^{-(e+1)-p}$ を用いて $$ \roundp{\tfrac1x} = \tfrac1x + \varepsilon $$ と書ける。よって、

$$ \begin{aligned} \roundp{\tfrac1x}\times x &= (\tfrac1x + \varepsilon)\times x \\ &= 1 + \varepsilon x \end{aligned} $$ となる。$\varepsilon$ と $x$ の範囲から、$|\varepsilon x|\lt 2^{-p}$ とわかる。 Claim 1 同様にして、成り立つことがわかる。$\qed$

note: Claim 3 では $x$ が正規化数の範囲をまとめて扱ったが、Claim 5 の $e$ の範囲を広げることでも示せそうである。

Claim 6: $e\in[\emin\ldots\emax]$ に対して、$x=1\times2^e$ のとき $\roundp{\tfrac1x}\otimes x=1$ が成り立つ。

Proof

$\roundp{\tfrac1x}=2^{-e}$ より従う。$\qed$

Claim 7: $\roundp{\tfrac10}\otimes0$ および $\roundp{\tfrac1{\infty}}\otimes\infty$ は NaN となる。また、$x\in(0\ldots 2^{-(\emax+2)}]$ では $\roundp{\tfrac1x}\otimes x = \infty$ となる。

Proof

前半は、$\roundp{\tfrac10}=\infty$ および $\roundp{\tfrac1{\infty}}=0$ であることと $0\otimes \infty$ が NaN であることから従う。後半は、該当の範囲で $\roundp{\tfrac1x}=\infty$ であることと $x\gt 0$ であることから従う。 $\qed$

よって、$\roundp{\tfrac1x}$ が正規化数となる範囲では $\roundp{\tfrac1x}\otimes x\in\{1-2^{-p}, 1\}$ であることがわかりましたが、そうでない範囲においては反例があることがわかりました。

たとえば、$x = \texttt{10000004e62386}_{(16)}\times 2^{970}$ に対して $\roundp{\tfrac1x}\otimes x = 1 + 2^{1-p}$ となります。また、最小性は示していませんが、$x = \texttt{18000003000000}_{(16)}\times 2^{970}$ に対して $\roundp{\tfrac1x}\otimes x = 1 - 2\cdot 2^{-p}$ となります。

疑問 2 の解答

$\roundp{\tfrac1x}\otimes x \lt 1$ なる $x\in(1\ldots2)$ の最小値について考えます。

Claim 8: $\varepsilon^2\lt 2^{-(p+1)}$ に対して $x=1+\varepsilon$ と書けるとき、$\roundp{\tfrac1x}\otimes x=1$ が成り立つ。

Proof

等比級数の和の公式から、 $$ \begin{aligned} \tfrac1x &= \tfrac1{1+\varepsilon} \\ &= 1 - \varepsilon + \tfrac{\varepsilon^2}{1+\varepsilon} \end{aligned} $$ である。ここで、$\tfrac{\varepsilon^2}{1+\varepsilon} \le \varepsilon^2 \lt 2^{-(p+1)}$ より、 $1\oslash x = 1-\varepsilon$ となる*1。この範囲では $$ \begin{aligned} \roundp{\tfrac1x}\times x &= (1-\varepsilon)\times(1+\varepsilon) \\ &= 1-\varepsilon^2 \end{aligned} $$ となり、$\varepsilon^2\lt 2^{-(p+1)}$ より $\roundp{\tfrac1x}\otimes x=1$ となる。 $\qed$

Claim 8 の証明から、$\tfrac{\varepsilon^2}{1+\varepsilon}$ による誤差が十分に大きくなれば $\roundp{\tfrac1x}\otimes x\lt 1$ が成り立ちそうです。 $1+2^{-(p+1)/2}$ 付近に解があることに期待して、仮数部を昇順に全探索してもよいかもしれません。

実際、仮数部が 53 bits である double や、long double のうち仮数部が 64 bits であるものについては比較的すぐ見つかりました。一方、仮数部が 113 bits である __float128（あるいはそういう処理系の long double）ではうまくいきませんでした。そこで、もう少し考察を進めます。

まず、整数 $m\in(2^{p-1}\ldots 2^p)$ に対して $x = m\times 2^{1-p}$ と表すことができ、 $$ \begin{aligned} \tfrac1x &= \tfrac{2^p}x\times 2^{-p} \\ &= \tfrac{2^{2p-1}}m\times 2^{-p} \end{aligned} $$ です。ここで、$m$ の範囲から $\tfrac{2^{2p-1}}m\in(2^{p-1}\ldots 2^p)$ であることがわかります。よって、 $$ \roundp{\tfrac1x}\in\left\{\floor{\tfrac{2^{2p-1}}m}\times 2^{-p}, \ceil{\tfrac{2^{2p-1}}m}\times 2^{-p}\right\} $$ となります。また、 $$ \Floor{\frac{2^{2p-1}}m} = \frac{2^{2p-1} - (2^{2p-1}\bmod m)}m $$ です。

丸めの際は近い方が採用されるため、$(2^{2p-1}\bmod m) \lt \tfrac m2$ のときは $\floor{\ldots}\times 2^{-p}$ の方が採用され、$(2^{2p-1}\bmod m)\gt \tfrac m2$ のときは $\ceil{\ldots}\times 2^{-p}$ の方が採用されます。$m$ の範囲から $\tfrac{2^{2p-1}}m$ が（2 進法で）無限小数となることがわかり、$(2^{2p-1}\bmod m)=\tfrac m2$ となるケースは存在しないことが従います。

よって、次のようになります。

Claim 9: $r = 2^{2p-1}\bmod m$ が $r\gt\tfrac m2$ を満たすとき、$\roundp{\tfrac1x}\otimes x=1$ が成り立つ。

Proof

該当の条件のとき、 $$ \begin{aligned} \roundp{\tfrac1x}\times x &= \ceil{\tfrac{2^p}x}\times 2^{-p} \\ &\gt \tfrac{2^p}x\times 2^{-p} \\ &= \tfrac1x \end{aligned} $$ であることと、Claim 1 より従う。$\qed$

すなわち、$\floor{\ldots}\times 2^{-p}$ が採用される方のみ考えればよいです。

Claim 10: $r = 2^{2p-1}\bmod m$ が $2^{p-2}\lt r\lt \tfrac m2$ を満たすとき、$\roundp{\tfrac1x}\otimes x=1-2^{-p}$ となる。

Proof

$r\lt \tfrac m2$ の範囲で考える。 $$ \begin{aligned} \roundp{\tfrac1x}\times x &= \left(\floor{\tfrac{2^{2p-1}}m}\times 2^{-p}\right)\cdot (m\times 2^{1-p}) \\ &= \left(\tfrac{2^{2p-1}-r}m\times 2^{-p}\right)\cdot (m\times 2^{1-p}) \\ % &= (2^{2p-1}-r)\cdot 2^{-(2p-1)} \\ &= 1 - r\cdot 2^{-(2p-1)} \end{aligned} $$ より、$r\cdot 2^{-(2p-1)} \gt 2^{-(p+1)}$ であれば $1-2^{-p}$ に丸められる。すなわち、$r\gt 2^{p-2}$ であることが必要十分である。$\qed$

Observation 11: $r = 2^{2p-1}\bmod m$ は次のようになっている。

$p=53$ とする。また、$m = 2^{p-1}+i$ とする。

$i$	$0$	$1$	$2$	$3$	$\cdots$	$47453132$	$47453133$	$47453134$	$47453135$
$r$	$0$	$2$	$8$	$18$	$\cdots$	${\tiny 4503599473218848}$	${\tiny 4503599663031378}$	${\scriptsize 178020282}$	${\scriptsize 367832819}$

ある $k$ に対し、$r = 2i^2 - k\cdot(2^{p-1}+i)$ となっている。

Proof

$$ \begin{aligned} 2^{2p-1} &= (2^p-2i)\times(2^{p-1}+i) + 2i^2 \\ &= (2^p-2i)\times(2^{p-1}+i) + \Floor{\tfrac{2i^2}{2^{p-1}+i}}\cdot (2^{p-1}+i) + (2i^2\bmod (2^{p-1}+i)) \\ &= \left(2^p-2i+\Floor{\tfrac{2i^2}{2^{p-1}+i}}\right)\times(2^{p-1}+i) + \left(2i^2-\Floor{\tfrac{2i^2}{2^{p-1}+i}}\cdot(2^{p-1}+i)\right) \end{aligned} $$ である。$k$ は $2i^2\div(2^{p-1}+i)$ の商に他ならない。$\qed$

そこで、次のような方針が考えられます。

$k$ を固定する。
その範囲で $r\gt 2^{p-2}$ を満たす $i$ の最小値 $i^{\star}$ を求める。
- ここは二分探索できる。
$i^{\star}$ で計算したあまり $r^{\star}$ が $r^{\star}\lt\tfrac m2$ を満たしていればそれが答え。
- 満たしていなければ次の $k$ を試す。

$k$ の上界はよくわからないものの、いくつか試したところでは十分小さい $k$ で停止しました。今回の目的は計算量改善ではないのでよしとします。

Example 12:

メジャーな処理系における _Float16, float, double, long double, __float128 を想定し、上記の計算をした。

type	minimum $x\in(1\ldots2)$ subject to $\roundp{\tfrac1x}\otimes x\lt 1$	$p$	$k$
`_Float16`	$\texttt{41c}_{(16)}\times2^{-10}$	$11$	$1$
`float`	$\texttt{801467}_{(16)}\times2^{-23}$	$24$	$6$
`double`	$\texttt{1000000f5cbf2a}_{(16)}\times2^{-52}$	$53$	$29$
`long double`	$\texttt{800000005a82799a}_{(16)}\times2^{-63}$	$64$	$0$
`__float128`	$\texttt{1000000000000022df0668c89bf13}_{(16)}\times2^{-112}$	$113$	$9$

long double 以下の型については、$1$ から昇順に全探索したものと値が一致することを確認した。また、各型の値について、実際に $\roundp{\tfrac1x}\otimes x\lt 1$ となることを確認した。

あとがき

「浮動小数点数型」に対して、「謎の誤差が出て変な感じになる」というような感覚を持っている人はおそらく多いのではないかと思っていますが、こうした解析で遊んでいるとある程度「数学側のやり方で（手間は多少かかっても）ある程度まともに扱える」という気持ちになってくるのではないでしょうか。少なくともえびちゃんの中ではそういう感じになってきました。コンテスト中になにかできそうな気はまだしません。

とはいえまだ基礎パートの解説記事を書いていないので、取っかかりが掴めていない段階の人にはまだかもという気はしています。もう少し待っていてもらえるとうれしいです。

というのも、ある程度自分が親しんでいない段階では基礎用の記事になにを書くべきか判断しかねたからです。ボツになった下書きが二つか三つくらいあります。

その練習のために書いたのが、Sqrt Inequality の記事と floor(n * 0.1) の記事と今回の記事でした（最初にやるにしては重くないか？）。もう一つ練習に書こうとしているのがあるので、それを終えたら基礎用の記事を書こうと思っています。

こういう細々とした詰め詰めな解析も楽しいんですが、競プロ文脈で（コンテスト中に計算量のオーダーを見積もるくらいの手軽さで）誤差の見積もりができるようにもなれたらうれしいので、そういうところのお勉強もしていきたいです。

あまり多くの競プロ er が興味を持たないところに興味を持ってしまうところ、えびちゃんのよいところである気はしていますが、自分としてはさみしいところでもあります。

おまけ

#include <cstdio>
#include <string>
#include <quadmath.h>

using f16 = _Float16;
using f32 = float;
using f64 = double;
using f80 = long double;
using f128 = __float128;

std::string q(f128 x) {
    char buf[1000];
    quadmath_snprintf(buf, sizeof buf, "%.999Qg", x);
    return std::string(buf);
}

int main() {
    int x = 27355;
    printf("%.99g\n", f16(f16(f16(1.0) / f16(x)) * f16(x)));
    printf("%.99g\n", f32(f32(1.0F / f32(x)) * f32(x)));
    printf("%.99g\n", 1.0 / f64(x) * f64(x));
    printf("%.99Lg\n", 1.0L / f80(x) * f80(x));
    printf("%s\n", q(f128(1.0) / f128(x) * f128(x)).data());
}

0.99951171875
0.999999940395355224609375
0.99999999999999988897769753748434595763683319091796875
0.9999999999999999999457898913757247782996273599565029144287109375
0.99999999999999999999999999999999990370350278063820734720110287075363407309491758923059023800306022167205810546875

wandbox.org

5 つの型すべてで $\roundp{\tfrac1x}\otimes x\lt 1$ にさせる最小の非負整数です。実際には _Float16 では $\roundp{x}\ne x$ なのでちょっとずるいです。

自分で仮数部 $p$ bits の浮動小数点数型をエミュレートする型を作り、これら以外の型について遊ぶこともできます（遊びました）。

おわり

おわりです。

*1:$1+\varepsilon\lt 2$ が該当の浮動小数点数型で表せていることから、$\roundp{1-\varepsilon} = 1-\varepsilon$ です。

2023-10-30

floor(n / 10.0) と floor(n * 0.1) は等価？

浮動小数点数

これを double で計算したときに常に等しくなりますか？というクイズです。 n は double で表せる整数であることにします。

答え

>>> 6755399441055749 / 10
675539944105574.9
>>> 6755399441055749 * 0.1
675539944105575.0

より、前者の floor と後者の floor は異なることになります。すなわち、等価ではありません。

解説など

真であるなら証明をする必要があるんですが、偽だったので反例を挙げて終了です。それはそれとして、「これの範囲では成り立つ」という話をしても損はしないので、してみましょう。

浮動小数点数型についてある程度の知識は持っているものとします。そうでない人のための記事は近々（次か次の次くらいに）書く予定です。

$\gdef\fround#1{(\hspace{-.2em}[#1]\hspace{-.2em})}$

まず、0.1 と言っているのは、double（かつ、多くの処理系でデフォルトになっている丸め方向）の文脈においては $$ \begin{aligned} &\phantom{{}={}} {\footnotesize 0.100000000000000005}{\scriptsize 55111512312578270211}{\tiny 81583404541015625} \\ &= \texttt{ccccccccccccd}_{(16)} \times 2^{-55} \\ &= (\tfrac45(2^{52}-1)+1)\cdot 2^{-55} \\ % &= (\tfrac45(2^{-3}-2^{-55})+2^{-55}) \\ &= \tfrac1{10}+\tfrac15\cdot 2^{-55} \\ \end{aligned} $$ と等しい値になります*1。

実数 $x$ を浮動小数点数型（ここでは double）で表せる値に丸めたものを $\fround{x}$ と書くことにします*2。たとえば、$\fround{0.1} = 0.1 + \tfrac15\cdot 2^{-55}$ や $\fround{10} = 10$ です。

また、浮動小数点数 $x$, $y$ に対して $\fround{x\times y}$ および $\fround{x\div y}$ を、それぞれ $x\otimes y$ および $x\oslash y$ と書くことにします。浮動小数点数型の値同士の演算のみしか考えていないことに注意してください。

さて、今回やりたいのは $\floor{n\oslash 10}$ と $\floor{n\otimes\fround{0.1}}$ の比較です。

note: double の仮数部は $p = 53$ bits なので、$2^{53}$ 以上の（浮動小数点数型の）値 $x$ は整数となるため、$\floor{x}=x$ となります。

実数 $x\ge 0$ と整数 $k$ に対して $\floor{\fround{x}} = k$ となる条件は、$\fround{x}$ の整数部分が $k$ となることです。丸め誤差の関係で、$x$ の整数部分が $k$ であることとは同値ではありません。

さて、浮動小数点数で表せる整数 $n = 10q+r$ ($q\in\N, 0\le r\le9$) を考えます。$2^e\le 10q+r \lt 2^{e+1}$ とします。$0\le e\lt p$ としておきます。

$\fround{0.1} = \tfrac1{10} + \tfrac15\cdot 2^{-55} = \tfrac1{10}\cdot(1+2^{-54})$ なので、 $$ \begin{aligned} (10q+r)\cdot\fround{0.1} &= (q+\tfrac r{10})\cdot(1+2^{-54}) \\ &= q+(\tfrac r{10} + (q+\tfrac r{10})\cdot 2^{-54}) \end{aligned} $$ となります。一旦 $r$ が最大の $9$ である場合について考えます。 $$ \begin{aligned} (10q+r)\cdot\fround{0.1} &\le q+(\tfrac9{10} + \tfrac1{10}(10q+r)\cdot 2^{-54}) \\ &\lt q+\tfrac9{10} + \tfrac1{10}\cdot 2^{e+1-54} \\ &= q+\tfrac9{10} + \tfrac1{10}\cdot 2^{e-53}. \end{aligned} $$ $\tfrac1{10}(10q+r)\lt\tfrac18(10q+r)\lt 2^{e+1-3}$ より、$\tfrac1{10}(10q+r)$ の 0.5 ULP は $2^{e-56}$ 以下です。小数部分を考えて、$\tfrac9{10}+\tfrac1{10}\cdot 2^{e-53} \lt 1-2^{e-56}$ の範囲においては、$(10q+r)\cdot\fround{0.1}$ は $q$ の方に丸められます。これにより、$e\le51$ の範囲では $\floor{n\oslash 10} = \floor{n\otimes\fround{0.1}}$ となります。また、同様の考察により、$r\le 8$ では $e=52$ でも $\floor{n\oslash 10} = \floor{n\otimes\fround{0.1}}$ となることが示せます。

よって、$\floor{n\oslash 10} \ne \floor{n\otimes\fround{0.1}}$ となる最小の $n$ は、$e = 52$ かつ $r = 9$ においてこれを満たす最小の $q$ によって作られることがわかります（この時点では存在性は示していませんが、これから作るので大丈夫です）。 $(10q+r)\cdot\fround{0.1}$ の小数部分を考えて*3、 $$ \begin{aligned} \tfrac9{10} + (q+\tfrac9{10})\cdot 2^{-54} &\ge 1-2^{-4} = \tfrac{15}{16} \\ q+\tfrac9{10} &\ge \tfrac3{80} \cdot 2^{54} \\ &= \tfrac35\cdot 2^{50}. \end{aligned} $$ よって、条件を満たす最小の $q$ は $\ceil{\tfrac35\cdot 2^{50}-\tfrac9{10}} = \tfrac15\cdot(3\cdot2^{50}-2)$ となります。これにより、最小の $n$ は $10\cdot\tfrac15\cdot(3\cdot 2^{50}-2)+9 = 3\cdot2^{51}+5$ となります。これは、冒頭に挙げた例そのものです。

disclaimer: ここでは、該当の範囲で $(10q+r)\oslash 10 = q$ となることをちゃんと示していません。よしなに。

一般のケース

一般に、$\floor{n\oslash d}$ と $\floor{n\otimes\fround{\tfrac1d}}$ を考えてみます。 $d = (2k+1)\cdot 2^s$ とします。二冪のケースは指数部が変わるだけなので $k\ge 1$ とします。 $$ \tfrac1d = \tfrac{2^t}{2k+1}\times 2^{-(s+t)} $$ と表せます。ここで、$t$ は $2^{p-1}\le\tfrac{2^t}{2k+1}\lt 2^p$ となるように定めます。すなわち、$t = (p-1)+\ceil{\log_2(2k+1)}$ となります。 $$ \tfrac1d = \underbrace{\tfrac{2^{(p-1)+\ceil{\log_2(2k+1)}}}{2k+1}}_{\mu}\times 2^{-(s+t)}. $$ この $\mu$ の部分を $\floor{\mu}$ または $\ceil{\mu}$ に適切に丸めたものを $m$ として、 $$ \fround{\tfrac1d} = m\times 2^{-(s+t)} $$ となります。 $$ (2^{(p-1)+\ceil{\log_2(2k+1)}})\bmod (2k+1) \le k $$ であれば $m=\floor{\mu}$ に、そうでなければ $m=\ceil{\mu}$ になります。よって、$0\lt u\le k$ とおき、丸め方向に応じて適切な符号を用いることで $$ \fround{\tfrac1d} = \tfrac{2^t\pm u}{2k+1}\times 2^{-(s+t)} $$ と書くことができます。$2^t$ の持つ素因数が $2$ のみであることと、$2k+1\ge3$ が $2$ 以外の素因数のみを持つことから、$u=0$ にはなりません。

符号ごとに分けて考えます。

まず、$\fround{\tfrac1d} = \tfrac{2^t+u}{2k+1}\times2^{-(s+t)}$ のケースです。$\fround{\tfrac1d}\gt\tfrac1d$ に注意します。 $n = dq + r$ とし、$2^e\le d q+r\lt 2^{e+1}$ とします。

$$ \begin{aligned} (dq+r)\cdot\fround{\tfrac1d} &= (dq+r) \cdot \left(\tfrac{2^t+u}{2k+1}\times 2^{-(s+t)}\right) \\ &= (dq+r) \cdot \tfrac{2^t+u}d\cdot 2^{-t} \\ &= (q+\tfrac rd) \cdot (1+u\cdot 2^{-t}) \\ &= q + \tfrac rd + (q+\tfrac rd)\cdot u\cdot 2^{-t} \\ &= q + \tfrac rd + (q+\tfrac rd)\cdot u\cdot 2^{-(p-1)-\ceil{\log_2(2k+1)}}. \end{aligned} $$ あとは、先ほど同様に小数部分を考えていけばよいでしょう。疲れてきたので割愛です。

次に、$\fround{\tfrac1d} = \tfrac{2^t-u}{2k+1}\times2^{-(s+t)}$ のケースです。$\fround{\tfrac1d}\lt\tfrac1d$ に注意します。

いきなりですが、実はめちゃくちゃ反例があります。下記は Haskell の対話型ツール GHCi での出力です。

ghci> take 10 . map round $ filter (\x -> x * (1.0 / x) /= x / x) [1..]
[49,98,103,107,161,187,196,197,206,214]

$49\otimes\fround{\tfrac1{49}} \lt 1$ や $197\otimes\fround{\tfrac1{197}} \lt 1$ などが成り立ちます。先の $\fround{\tfrac1{10}}$ のケースでは身近な範囲（たとえば int に収まる程度）では大丈夫だったりしましたが、こちらに関しては身近な範囲でも全然事故りそうですね。注意として、浮動小数点数の精度によってこの反例は異なります。下記は Rust の対話型ツール evcxr での出力です*4。

>> (1..).filter(|&x| 1.0 / x as f32 * x as f32 != 1.0).take(10).collect::<Vec<_>>()
[41, 47, 55, 61, 82, 83, 94, 97, 107, 109]
>> (1..).filter(|&x| 1.0 / x as f64 * x as f64 != 1.0).take(10).collect::<Vec<_>>()
[49, 98, 103, 107, 161, 187, 196, 197, 206, 214]

この $49\otimes\fround{\tfrac1{49}}$ のようなものについては、書きたいことがいろいろあるので、別の記事で書くことにします。

また、興味のある読者は floor で満足せず、ceil について考えてみてもよいでしょう。ところで $\ceil{x} = -\floor{-x}$ や $\fround{-x} = \fround{x}$ などが成り立ちます。 f32 においては、仮数部が 24 bits しかないため、それなりに十分なテストを行うこともできます。C++ では floating-point promotion によって float が勝手に double になったりするため、Rust などを用いる方が楽そうです。

速度に関して

整数の除算が遅いというのは比較的有名な話かと思います。

そこで、$\sum_{i=1}^{10^{10}} \floor{\tfrac i{10}}$ を、整数型 (i64) と浮動小数点数型 (f64) の $i\oslash 10$ 版と $i\otimes \fround{\tfrac1{10}}$ 版で計測してみました。また、定数除算最適化の影響を考慮するため、除数を定数として用意するものと実行時に与えるものを比較しました。

もちろん値自体は $O(1)$ で（というか紙とペンで）求めることができて、 $$ \begin{aligned} \sum_{i=1}^{10^{10}} \floor{\tfrac i{10}} &= 9\sum_{i=0}^{10^9-1}i + \sum_{i=1}^{10^9}i \\ &= 10\sum_{i=1}^{10^9-1}i + 10^9 \\ &= 5\cdot (10^9-1)\cdot 10^9 + 10^9 \\ &= 5\cdot 10^{18} - 4\cdot10^9 \\ &= 4999999996000000000 \end{aligned} $$ です。これは $2^{63}$ よりも小さいことを指摘しておきます。

足し合わせる変数には i64 を用い、除算パートのみで上記の違いがある状態にします。なお、$10^{10} \lt 3\cdot 2^{51}+5$ なので、この範囲で正しい値を計算できることに注意してください。計測結果は下記の通りでした。

	`i64`	`f64` $\oslash\, {10}$	`f64` $\otimes\,\fround{0.1}$
定数	2.26 s	2.86 s	2.24 s
変数	5.75 s	2.90 s	2.20 s

どうやら、いずれも f64 で $\otimes$ を用いるものが最も速かったようです。あらあら。あと、やはり変数での除算は重いようですね。

今回は floor の方のみ実験しましたが、ceil の方においても適切に丸めた方を用いたりすることで高速に正しい値を得ることが可能そうな気がします。浮動小数点数型と仲よくなると、適切に証明した上で高速化をすることができそうでうれしいですね。

おきもち

コンテスタント目線では、この手の誤差評価などは “非本質” などと扱われがちな部分な気はしています。一方で、先ほど言及したように高速化に寄与する場合があり、知っていると役に立つ場面もあるのかなとも思います。

あるいは、テストケースを作成する側に回った際にも、「浮動小数点数型を用いたこうした解法は実は正当である」とか、「誤差がこうなることがわかるので、こういうケースを置いておけば落とせそう」とか、そうした考え方ができるとうれしそうです。AtCoder では（特に double を想定したケースは）考慮されていがちな気はしますよね（long double は見落とされがちな気もしますが）。

ソフトウェア開発の文脈だと、「数学的な証明をしました」で単体テストを書かないのは不安な感じがある気がします。一方、この手の浮動小数点数が絡むものにおいては、単体テストで適当にいくつかのケースを選んで「大丈夫そうですね」と言うのは不安な気持ちになります。テストしていないケースにおいて誤差でめちゃくちゃになることが容易に起こりうるからです。多少詳しい人でないと、先の例の $\fround{\tfrac1{49}}$ などの反例を挙げることは難しそうですよね。

とはいえ、浮動小数点数に明るくない人（世の中のマジョリティ）にそういう事情を説明してわかってもらうのは大変そうな気がしていて、うまいことやっていく必要がありそうだなという気持ちになったりしています。この辺の話は正規表現などのトピックでも当てはまりそうですね。

特になにかあったわけではないのですが、競プロ寄りの人々には興味のない話かなと思って薄めの文字で書きました。

あとがき

$\fround{x}$ の記法は、今回思いつきで導入したものになります。下記の前回の記事では $\hat x$ とか $\operatorname{round}(x)$ とか書いていたものですが、長めの式に対してもよい感じに書けるものがいいなと思ったので考えてみました。

rsk0315.hatenablog.com

経緯としては、最近（文脈自体はあまり関係ないものの）$[\![\texttt{a+b}]\!]$ のような括弧が使われているのを見て、「括弧をそうやって作るのも “アリ” なんだな」と思ったのがきっかけになります。操作としては round を表すので丸い括弧がいいなと思って、$(\hspace{-.2em}|x|\hspace{-.2em})$ や $\fround{x}$ などを試して決めたものです。 $(\hspace{-.2em}|x|\hspace{-.2em})$ もいいなと思ったのですが、別の文脈*5での使われを思い出したことと、$[\ldots]$ の方が Gauß 記号や $\floor{\ldots}$ $\ceil{\ldots}$ を連想させて、丸めに近いかな？という気持ちがあったことによります。定義は下記です*6。

$\gdef\fround#1{(\hspace{-.2em}[#1]\hspace{-.2em})}$

もしかしたらこの括弧に別の意味があったり、該当の概念に別の記法があったりするのかもしれませんが、あまりサーベイしていません。

おわり

最近は需要のなさそうな記事ばかり書いていますが、何らかの局面で誰かしらの役に立ったらいいなと思っています。あるいは他の人の役に立っていなくても、書いているときの自分が気持ちよくなっているので十分です。

*1:「誤差でぐちゃぐちゃになった意味不明な値」のようなイメージがある人もいるかもしれませんが、数学的に書ける値であることを意識するとよいかもです。

*2:仕様に応じて適切な丸め方向を採用するものとします。

*3:仮数部の偶奇に応じて等号を含むかが変わりますが、一旦等号を含むものから考えてみて、だめならその次のものを採用することにします。

*4:いろいろな対話型ツールがあってたのしいですね。

*5:The SATySFibook の record の部分です。こっちも別の出典があるのかも。分野がやや違いそうなので、同じ記号でも気にしてなくていいと思ってはいます。

*6:論文とかで何らかの独自の記法を使うときは、$\LaTeX$ での書き方を footnote にでも書いておいてくれと思っているタイプです。

2023-10-22

Sqrt Inequality の誤差評価

浮動小数点数

atcoder.jp

導入

これの long double 解法についての話です。下記のような解法です。

#include <cmath>
#include <cstdio>

constexpr long double eps = /* ??? */;

bool ok(long double a, long double b, long double c) {
  return std::sqrt(a) + std::sqrt(b) + eps < std::sqrt(c);
}

int main() {
  long double a, b, c;
  scanf("%Lf %Lf %Lf", &a, &b, &c);
  puts(ok(a, b, c) ? "Yes" : "No");
}

ただし、long double の仮数部は 64-bit であるとします（AtCoder 環境に準拠）。公式解説では 1.0E-14 に設定し、

ここで $\varepsilon = 10^{-14}$ の値を大きくしすぎたり小さくしすぎたりすると WA になってしまいます。なぜ $\varepsilon = 10^{-14}$ であるとうまくいくのかは説明できるのですが、より簡単で安全な方法は、全て整数でやってしまうことです。

と書かれています（1.0E-14 は double のリテラルであり、long double の 1.0E-14L より若干精度が劣りますが、ここでは問題になりません）。実際、

eps = 8.8817841970012523238e-16L;

では WA になり、

eps = 8.8817841970012523239e-16L;

では AC になります。また、

eps = 1.5099033134902128948e-14L;

では AC になり、

eps = 1.5099033134902128949e-14L;

では WA になります。

ざっくり 4 桁で見ると、$8.882\times 10^{-16} \le \text{\texttt{eps}} \le 1.509\times 10^{-14}$ の範囲であれば AC になります。

数値のざっくりとした出処は、$2^{\floor{\log_2(\sqrt{10^9})}-64} = 2^{-50}$ と $$ \underbrace{(\ceil{\tfrac12(10^9)^{-3/2}\cdot 2^{(64-1)-\floor{\log_2(\sqrt{10^9})}}}-\tfrac12)}_{8.5} \times \underbrace{2^{\floor{\log_2(\sqrt{10^9})}-(64-1)}\vphantom{\tfrac12}}_{2^{-49}} $$ です（これらに真に含まれる範囲）。

これをちゃんと考えていきましょうというのが今回のお話です。なお、後述しますが、上記の eps = 1.5099033134902128948e-14L による解法は hack することが可能です。

導入
- 定義・前提知識
方針
考察
- $\varepsilon_L'$ の上界
- $\varepsilon_U'$ の下界
hack について
ケースに関して
差の下界について
まとめ
参考
感情
おまけ
おわり

定義・前提知識

$\gdef\sqrtf{\text{\texttt{sqrt}}}$ $\gdef\roundf{\operatorname{round}}$

浮動小数点数型についての基本的な知識はあるものとします（表現の仕方、丸めによって誤差が生じることなど）。ない人のための記事はそのうち書く予定ではあります。

実数 $x$ に対し、（対象としている浮動小数点数型で）表せる最も近い数値*1を $\roundf(x)$ および $\hat x$ で表します（$x$ の部分が長いときには前者の記法を好んで使い、主に一文字の変数のときには後者の記法を好んで使うことにします）。浮動小数点数型の数値 $x$, $y$ に対して、$x\oplus y = \roundf(x+y)$ と $\sqrtf(x) = \roundf(\sqrt{x})$ と書くことにします。

和および平方根の計算では、正確に計算した値を正確に丸めた値を返すことが定められているので、誤差は 0.5 ULP です。すなわち、正確な値が $2^e$ 以上 $2^{e+1}$ 未満の範囲にあるとき、誤差の絶対値は $2^{e-p}$ 以下です。ここで、AtCoder 環境における long double においては $p = 64$ です。double や一部環境の long double では $p = 53$、別の環境の long double では $p=113$ だったりします。

また、正の実数 $x$, $y$, $x+y$ が正規化数の範囲に入るとき、$\roundf(x+y)$ と $(\hat x\oplus\hat y)$ の差の絶対値は 1 ULP 以下であることが示せます。すなわち、「丸めた値同士を足してから丸めた値」と「無限精度で足したものを一度だけ丸めた値」の差が 1 ULP で抑えられるということです*2。直感的には、$\roundf(x+y)=(\hat x\oplus\hat y)$ であるか、浮動小数点数型で表せるうちで隣り合う 2 数になります。

$\gdef\ulp{\operatorname{ulp}}$

実数 $x$ に対し、$\ulp(x) = 2^{\floor{\log_2(|x|)}+1-p}$ とします。$2^i\le x\lt2^{i+1}$ のとき、$\ulp(x) = 2^{i+1-p}$ となります。また、 $$ \ulp(x)\cdot 2^{p-1}\le x\lt \ulp(x)\cdot2^p $$ や $$ \ulp(x)\cdot 2^{p-1}\le \hat x\le \ulp(x)\cdot2^p $$ が成り立ちます。

方針

誤差のない世界で計算した結果と、誤差のある世界で計算した結果とで、大小関係が常に同じだと楽です（誤差を無視して計算して答えればよいので）。しかし、実際には次のようなケースが存在します。

$$ \sqrt{a} + \sqrt{b} = \sqrt{c} \wedge \sqrtf(a) \oplus \sqrtf(b) \lt \sqrtf(c). $$

そこで、適切な $\varepsilon$ を足すことで、上記の「大小関係が常に同じ」を成り立たせる方法が知られています（蟻本 p. 228 など）。 $$ \sqrt{a} + \sqrt{b} \lt \sqrt{c} \iff (\sqrtf(a) \oplus \sqrtf(b)) \oplus \varepsilon \lt \sqrtf(c). $$ ここで、平方根を求めた際やそれらの和を求めた際だけでなく、$\varepsilon$ を足した際にも丸め誤差が生じることに注意しておきましょう。

非常に嘆かわしいですが、この $\varepsilon$ をまともに考えて設定している人は、中級者以上でも非常に少ない印象があります*3。次のようなことを考える必要があるでしょう。

まず、$\sqrt{a} + \sqrt{b} \ge \sqrt{c}$ のケースにおいて、$\sqrtf(a) \oplus \sqrtf(b)$ が $\sqrtf(c)$ をどの程度下回るかを考える必要があります。「$\sqrt{a}+\sqrt{b}\ge\sqrt{c}$ のケースで No と答えるために、どの程度 $\varepsilon$ を大きくする必要があるか？」ということです。

また、$\sqrt{a} + \sqrt{b} \lt \sqrt{c}$ のケースにおいて、$\sqrtf(a) \oplus \sqrtf(b)$ が $\sqrtf(c)$ にどの程度近づくかを考える必要があります。「どの程度なら $\varepsilon$ を大きくしても、$\sqrt{a}+\sqrt{b}\lt\sqrt{c}$ のケースで Yes と答えるままでいられるか？」ということです。

未知の $\varepsilon_L$, $\varepsilon_U$ に対して $\varepsilon_L \le \varepsilon \le \varepsilon_U$ の範囲で AC できるとしたときに、$\varepsilon_L$ の上界 $\varepsilon_L'$ を求めるのが前者で、$\varepsilon_U$ の下界 $\varepsilon_U'$ を求めるのが後者です。$\varepsilon_L'\le\varepsilon_U'$ となるように求められれば、その範囲の $\varepsilon$ を使うことで十分となります。

逆に、そうした区間 $[\varepsilon_L\ldots\varepsilon_U]$ が空の場合はこの解法では正答できないことになります。

考察

一般性を失わず、$a\le b$ とします。また、制約の範囲において $a$, $b$, $c$ は注目している浮動小数点数型で正確に表せるので、そのことも仮定します。

$\varepsilon_L'$ の上界

$\sqrt{a}+\sqrt{b}\ge\sqrt{c}$ とし、$(\sqrtf(a)\oplus\sqrtf(b))-\sqrtf(c)$ を下から抑えます。丸めを一度しか行わない場合の値と隣り合った結果を得られることや、ULP の性質などに注意します。

まず、$\ulp(\sqrt{a}+\sqrt{b})\ge 4\ulp(\sqrt{c})$ の場合、 $$ \ulp(\sqrt{a}+\sqrt{b})\cdot 2^{p-1}\le \sqrt{a}+\sqrt{b} $$ より $$ \begin{aligned} \sqrtf(a)\oplus\sqrtf(b) &\ge \ulp(\sqrt{a}+\sqrt{b})\cdot 2^{p-1}-\tfrac12\ulp(\sqrt{a}+\sqrt{b}) \\ &\gt \ulp(\sqrt{a}+\sqrt{b})\cdot 2^{p-2} \\ &= \ulp(\sqrt{c})\cdot 2^{p} \\ &\ge \sqrtf(c) \end{aligned} $$ なので、$(\sqrtf(a)\oplus\sqrtf(b))-\sqrtf(c)\ge 0$ となります。

次に、$\ulp(\sqrt{a}+\sqrt{b})=2\ulp(\sqrt{c})$ の場合、 $$ \begin{aligned} \sqrtf(a)\oplus\sqrtf(b) &\ge \ulp(\sqrt{a}+\sqrt{b})\cdot 2^{p-1}-\tfrac12\ulp(\sqrt{a}+\sqrt{b}) \\ &= \ulp(\sqrt{c})\cdot 2^p - \ulp(\sqrt{c}) \\ &\ge \sqrtf(c) - \ulp(\sqrt{c}) \end{aligned} $$ なので、$(\sqrtf(a)\oplus\sqrtf(b))-\sqrtf(c)\ge -\ulp(\sqrt{c})$ となります。

最後に $\ulp(\sqrt{a}+\sqrt{b})=\ulp(\sqrt{c})$ の場合、

$$ \begin{aligned} \sqrtf(a)\oplus\sqrtf(b) &\ge \roundf(\sqrt{a}+\sqrt{b}) - \ulp(\sqrt{a}+\sqrt{b}) \\ &\ge \roundf(\sqrt{c}) - \ulp(\sqrt{a}+\sqrt{b}) \\ &= \sqrtf(c) - \ulp(\sqrt{c}) \end{aligned} $$ なので、こちらも $(\sqrtf(a)\oplus\sqrtf(b))-\sqrtf(c)\ge -\ulp(\sqrt{c})$ となります。

また、$\sqrtf(a)\oplus\sqrtf(b)\lt \sqrtf(c)$ の前提では $$\ulp(\sqrt{c})\cdot2^{p-1} \le \sqrtf(a)\oplus\sqrtf(b)\le \ulp(\sqrt{c})\cdot2^p$$ となります（刻み幅が $\ulp(\sqrt{c})$ の区間にあるということ）。

よって、差が $\ulp(\sqrt{c})$ 以下であることから、$(1-\tfrac12)\ulp(\sqrt{c})=\tfrac12\ulp(\sqrt{c})$ より真に大きい値を $\sqrtf(a)\oplus\sqrtf(b)$ に足すことで、$\sqrtf(c)$ 以上の値に丸められます。真に大きくする理由は、丸めの tiebreak で意図しない方に丸められるのを防ぐためです。

$c = 10^9$ かつ $p = 64$ の下では、これを満たす最小の値は $2^{\floor{\log_2(\sqrt{10^9})}-64}\cdot(1+2^{-63})$ であり、 $$ {\footnotesize 8.881784197001252324} {\scriptsize 3520183169} {\scriptsize 2018042652} {\tiny 79889712924636592690508241076940976199693977832794189453125} \times 10^{-16} $$ です。

なお、リテラルとして書いた際にも正確に丸められるため、$2^{\floor{\log_2(\sqrt{10^9})}-64}\cdot(1+2^{-64})$ より少しでも大きい値であれば上記の値に丸められます。大きすぎると別の値に丸められるので注意です。

constexpr long double eps = 8.8817841970012523238705358308233714632639944856462318296345254120538470488099846988916397094726562500000000000000000000000001e-16L;

/* 短めに書きたい人向け */
constexpr long double eps = 8.881784197001252324e-16L;
constexpr long double eps = 8881784197001252324e-34L;

もちろん、何らかの事情がない限りは、ぎりぎりの値を設定する必要はないため、9e-16L などと簡単な値を使いたいです。そのくらい大きくしても問題ないことを示すため、次の節が重要です。

$\varepsilon_U'$ の下界

$\sqrt{a}+\sqrt{b}\lt\sqrt{c}$ を考えます。まず、誤差のない世界での計算を考えます。

$$ (\sqrt a+\sqrt b+\sqrt c)(\sqrt a+\sqrt b-\sqrt c)(\sqrt a-\sqrt b+\sqrt c)(\sqrt a-\sqrt b-\sqrt c) $$ は整数となります。各因子の符号は $$ \underbrace{(\sqrt a+\sqrt b+\sqrt c)}_{\gt0}\underbrace{(\sqrt a+\sqrt b-\sqrt c)}_{\lt0}\underbrace{(\sqrt a-\sqrt b+\sqrt c)}_{\gt 0}\underbrace{(\sqrt a-\sqrt b-\sqrt c)}_{\lt0} $$ であり、整数性と合わせて $1$ 以上であることがわかります。符号が負のものは絶対値を考えるとして、 $$ |\sqrt a+\sqrt b-\sqrt c| \ge \frac1{(\sqrt a+\sqrt b+\sqrt c)(\sqrt a-\sqrt b+\sqrt c)\cdot|\sqrt a-\sqrt b-\sqrt c|} $$ であり、右辺の分母を上から評価したいです。

まず、$(\sqrt a+\sqrt b+\sqrt c)\lt 2\sqrt c$ です。次に、相加相乗平均の関係に気をつけつつ $$ \begin{aligned} (\sqrt a-\sqrt b+\sqrt c)\cdot|\sqrt a-\sqrt b-\sqrt c| &= -(\sqrt a-\sqrt b+\sqrt c)(\sqrt a-\sqrt b-\sqrt c) \\ &= -( (a+b-2\sqrt{ab})-c) \\ &= c - (a+b) + 2\sqrt{ab} \\ &\le c - (a+b) + (a+b) = c \end{aligned} $$ とわかります。よって、 $$ \sqrt c-(\sqrt a+\sqrt b) \gt \frac1{2c\sqrt c} $$ となります。すなわち、 $$ \sqrt a+\sqrt b\lt \sqrt c - \frac1{2c\sqrt c} $$ です。

$$ \begin{aligned} \sqrtf(a)\oplus\sqrtf(b) &\le \roundf(\sqrt{a}+\sqrt{b}) + \ulp(\sqrt{a}+\sqrt{b}) \\ &\le (\sqrt{a}+\sqrt{b})+\tfrac32\ulp(\sqrt{a}+\sqrt{b}) \\ &\lt \left(\sqrt c - \frac1{2c\sqrt c}\right) + \tfrac32\ulp(\sqrt{c}) \\ &= \sqrt c - \frac1{2c\sqrt c} + 3\cdot2^{\floor{\log_2(\sqrt{c})}-p}, \\ \sqrtf(c) &\ge \sqrt{c} - \tfrac12\ulp(\sqrt{c}) \\ &\ge \sqrt{c} - 2^{\floor{\log_2(\sqrt{c})}-p} \end{aligned} $$ なので、それらの差は $$ \begin{aligned} &\phantom{{}={}} \sqrtf(c) - (\sqrtf(a)\oplus\sqrtf(b)) \\ &\gt (\sqrt{c} - 2^{\floor{\log_2(\sqrt{c})}-p}) - \left(\sqrt c - \frac1{2c\sqrt c} + 3\cdot2^{\floor{\log_2(\sqrt{c})}-p}\right) \\ &= \frac1{2c\sqrt c} - 4\cdot2^{\floor{\log_2(\sqrt{c})}-p} \end{aligned} $$ となります。この差が $k\ulp(\sqrt{c})$ のとき、丸め誤差によって $$ (\sqrtf(a)+\sqrtf(b))\oplus\underbrace{(k-\tfrac12)\ulp(\sqrt{c})}_{\varepsilon} = \sqrtf(c) $$ になりえます。$\varepsilon$ を $(k-\tfrac12)\ulp(\sqrt{c})$ よりわずかに小さく設定すれば大丈夫そうです。

ということで、$c = 10^9$, $p = 64$ で計算してみると $$ \begin{aligned} \frac1{2c\sqrt c} - 4\cdot2^{\floor{\log_2(\sqrt{c})}-p} &= \tfrac12(10^9)^{-3/2} - 4\cdot 2^{\floor{\log_2(\sqrt{10^9})}-64} \\ &= \tfrac12\cdot 10^{-27/2} - 4\cdot 2^{-50} \\ &\approx 1.403 \times 10^{-14} \end{aligned} $$ となります。 $$ \underbrace{7.5\ulp(\sqrt{10^9})}_{\gt 1.332\times10^{-14}} \lt 1.403\times 10^{-14} \lt \underbrace{8\ulp(\sqrt{10^9})}_{\lt 1.422\times10^{-14}} $$ なので、$2k$ の整数性より $8\ulp(\sqrt{10^9})$ の差で下から抑えられそうです*4。 $7.5\ulp(\sqrt{10^9})$ 未満の最大の正規化数は $(1.875 - 2^{-63})\times 2^{16-63}$ であり、 $$ {\footnotesize 1.332267629550187848} {\scriptsize 43132080393349494087} {\tiny 7760882296602907258475934071384472190402448177337646484375} % {\footnotesize 1.332267629550187848} % {\scriptsize 4698394028} % {\scriptsize 2123965793} % {\tiny 88804411483014536292379670356922360952012240886688232421875} \times 10^{-14} $$ です。

こちらも、リテラルとして書いた際の丸めを考慮すると、$(1.875 - 2^{-64})\times 2^{16-63}$ より少しでも小さい値であれば上記の値に丸められます。

constexpr long double eps = 1.332267629550187848469839402821239657938880441148301453629237967035692236095201224088668823242187499999999999999999999999999e-14L;

/* 短めに書きたい人向け */
constexpr long double eps = 1.3322676295501878484e-14L;
constexpr long double eps = 13322676295501878484e-33L;

hack について

冒頭で、$8.5 \ulp(10^9)$ よりわずかに小さい値 1.5099033134902128948e-14L を使った解法は hack 可能と書きました。

hack できるケースの例は下記です。

239488342 239488343 957953370
247729794 247784320 991028225
248961062 248992620 995907363
249514080 249568802 998165761
249999995 249999996 999999982

これ以外にもたくさんあって、見つけた限りでも 649,152 個ありました。

ケースに関して

ローカル環境が M2 Mac なのですが、long double の精度が double と同じ $p=53$ でつらいです。数分かかるコードを AtCoder のコードテストに投げるわけにもいかないので、Docker でそういう環境を作って対処しました。$p=113$ の環境も作ることができたので満足です。

$\sqrt a+\sqrt b=\sqrt c$ かつ $\sqrtf(a)\oplus\sqrtf(b)\lt\sqrtf(c)$ の下で、誤差が大きくなるケースについても探索しました。 $\sqrt{c}\in\N$ となるケースでは（制約の範囲では）誤差が出ないので、square-free な整数 $m\ge 2$ に対して $i_a\sqrt m+i_b\sqrt m=(i_a+i_b)\sqrt m$ のケースを全探索しました。

8 268609842 268702562
8 799840008 800000000
1002528 236661768 268470792
99904500 467737920 999980820
100017228 400692747 901090683
100026368 400218632 900407048
107717610 430870440 969458490
109996887 439987548 989971983
249961152 250015923 999954147

など、235,172,054 個のケースで、$2^{-49} \approx 1.776\times 10^{-15}$ の誤差を出すことができました。

上記の $\sqrt{a}+\sqrt{b}=\sqrt{c}$ となるケース (2,809,012,244 個) と、後述の $\sqrt{a}+\sqrt{b}\lt\sqrt{c}$ の差が極小となるケース (250,056,424 個) すべてについて、今回設定した eps が正しく判定できていることを確認しました。前者は 11 分、後者は 1 分程度で済みました。

差の下界について

$\sqrt{a}+\sqrt{b}\lt\sqrt{c}$ のケースにおける $\sqrt{c}-(\sqrt{a}+\sqrt{b})$ の下界についてです。

この節の内容は eps の設定には直接関係ない（範囲を絞るだけなら前の節だけで十分）ですが、自分で解決できなかったので載せておきます。数学に自信のあるフォロヮをお待ちしております。

$$ f(c) = \min_{\substack{\sqrt{a}+\sqrt{b}\lt\sqrt{c'}\\c'\le c\\a, b, c\in\N}} \sqrt{c}-(\sqrt{a}+\sqrt{b}) $$ を考えます。

Conjecture 1: $c\le 10^5$ くらいで実験したところ、$c$ を昇順に見ていったときに $f(c)$ の値が更新されるのは、$c$ が下記の形で書けるとき (if and only if) になりそうです。

$c$	$a$	$b$	$c\bmod 4$	$f(c)$
$2(2k+1)$	$k$	$k+1$	$2$	$\varepsilon_0(c)$
$(2k+1)(2k+3)$	$k(k+1)$	$(k+1)(k+2)$	$3$	$\varepsilon_1(c)$
$2(k+1)(k+2)$	$\tfrac12 k(k+1)$	$\tfrac12(k+2)(k+3)$	$0$	$\varepsilon_2(c)$
$(2k+1)(6k+3\pm2)$	$k(3k\pm1)$	$(k+1)(3(k+1)\pm1)$	$1$	$\varepsilon_3(c)$

ただし、$\varepsilon_i(x)$ は $x\ge i$ で定義される関数 $$ \varepsilon_i(x) = \sqrt{x} - \frac12\left(\sqrt{x+i-2\sqrt{ix+1}} + \sqrt{x+i+2\sqrt{ix+1}}\right) $$ です。

Conjecture 2: 任意の $i$ について $\varepsilon_i(x) \sim \tfrac12x^{-3/2}$ が成り立ちそうです。

Wolfram|Alpha 的には合っていそうでしたが、数学力が足りず、自分では示せませんでした。

Conjecture 3: また、任意の $i$ について $x\ge i+1$ ならば $\varepsilon_i(x) \lt \varepsilon_{i+1}(x)$ が成り立ちそうです。

グラフをプロットした感じは合っていそうでしたが、手計算はうまくいかず、投げ出してしまいました。

まとめ

Sqrt Inequality に関して、$c\le 10^9$ かつ $p=64$ であれば $$(1+2^{-63})\times 2^{-50} \le \text{\texttt{eps}} \le (1.875 - 2^{-63})\times 2^{-47}$$ の範囲で設定すればよさそうということがわかりました。10 進 4 桁だと $$8.882\times 10^{-16} \le \text{\texttt{eps}} \le 1.332\times 10^{-14}$$ となります。

また、ジャッジ上で AC となるコードを hack できうるケースをたくさん見つけることができました。

浮動小数点数の誤差の評価が大変ということもわかりました（実際にはもっとゆるっと評価していい部分をがんばって詰めようとしたため、必要以上に大変になったという部分はあります）。これからもやっていきたいと思いました。

参考

twitter.com

こういうの、下からの評価の方が上からの評価よりも難しくなりやすく。逆数を上から評価しましょう。

x_0 = \sqrt{a}+\sqrt{b}+\sqrt{c}
x_1 = \sqrt{a}+\sqrt{b}-\sqrt{c}
x_2 = \sqrt{a}-\sqrt{b}+\sqrt{c}
x_3 = \sqrt{a}-\sqrt{b}-\sqrt{c}

とすると、x_0x_1x_2x_3 は整数になります。
— maspy (@maspy_stars) 2020年3月14日

このツイートでは $$\sqrt{c}-(\sqrt{a}+\sqrt{b})\ge \tfrac1{27}\cdot (10^9)^{-3/2}$$ だけを示していますが、方針で詰まっているときに読んでとても参考になりました。やっぱり数学に明るくて経験もある方は違うなと思いました。

ぽかいこちゃんリスペクト。

感情

浮動小数点数の誤差については未証明でもおっけーみたいな風潮ありますよね。未証明貪欲についてあれこれ言っている人が double で適当に計算して AC することに対して平気なのはダブルスタンダードなのではないでしょうか。

実用的にやばい例をまだ思いつけていないのですが、ある述語 $P(x)$ が、無限精度で計算したときは「$x\in(-\infty\ldots \alpha]$ では true となり、$x\in(\alpha\ldots{+\infty})$ では false となる」という条件を満たすとしても、$P(x)$ を浮動小数点数型で計算した場合にその性質が満たされるとは限りません。

たとえば $x \le 0 \vee (\tfrac 1x)\cdot x = 1$ という条件は上記の性質を満たします ($\alpha=+\infty$) が、浮動小数点数型では多くの反例があります。

$p$	型	反例
$24$	`float` (32-bit)	$41$
$53$	`double` (64-bit)	$49$
$64$	`long double` (80-bit)	$41$
$113$	`long double` (128-bit)	$43$

また、$[1\ldots 2)$ の範囲の最小の反例は、上記の $p$ の順に次のようになります。~~$p=113$ については探索中です*5。~~ $$ \begin{aligned} \texttt{801467}_{(16)}\times2^{-23} &= {\footnotesize 1.000622630119323730}{\scriptsize 46875}, \\ \texttt{1000000f5cbf2a}_{(16)}\times2^{-52} &= {\footnotesize 1.000000057228997096}{\scriptsize 81428248586598783731}{\tiny 4605712890625}, \\ \texttt{800000005a82799a}_{(16)}\times 2^{-63} &= {\footnotesize 1.000000000164636126}{\scriptsize 99697816044164255799}{\tiny 842067062854766845703125}. \end{aligned} $$

$p = 113$ は、全探索だと無理そうだったので、考えて求めました。別途記事にします。$1+(2^{57}+\varepsilon)\cdot 2^{-112}$ くらいだったので、$1$ から昇順に求めるのをやめて正解でした。 $$ \begin{aligned} &\phantom{{}={}} \texttt{1000000000000022df0668c89bf13}_{(16)}\times 2^{-112} \\ &= {\footnotesize 1.000000000000000030}{\scriptsize 24593730708203819677}{\tiny 89344748946071241091121526548576685378133532822175766341388225555419921875}. \end{aligned} $$

C/C++ で計算すると、float が勝手に double に変換される (floating-point promotion) ので誤差が生じないように見えましたが、$p=24$ で計算していないので不正ですね。なお、各値が正規化数となる範囲においては、上記のように $(1.0\oslash x)\otimes x = 1-2^{-p}$ になることはあるものの、上方向にズレて $1+2^{1-p}$ とはならないことを示せます（別の記事に書きそう）（非正規化数についてはまだ考えていません）。

↓ 追記：上記二つのトピックに対して書きました ↓

rsk0315.hatenablog.com

二分探索で何かしらの計算をしていて、的外れな値で意図しない true/false が返ってきて変なことになる可能性は全然あると思っていて、そうしたことに無頓着なのはどうなの？という気が最近しています。こわいので、最近は、可能な限り有理数で二分探索するようにしています。

atcoder.jp

とはいえ、計算途中でオーバーフローしないことの担保すら “非本質な細部” のように思っている人も多数いそうな感じがしますし、要求しすぎでしょうか。実際、これらの評価をコンテスト中にやれと言われたら、それが楽しいとはあまり思わなさそうです（というか時間が足りなさそう）（別のジャンルの競技としてはありかもしれませんが）。

それはそれとして、元の問題の解説で「整数でやりましょう」で済ませるのではなく「$\varepsilon=10^{-14}$ だとうまくいくことを説明できる」と書いてあるのがやっぱりすごいなと思いました。

ABC 169 C なども誤差に関する教育的な問題です。コンテスト当時は long double で解くだけで AC になってしまうようでしたが、現在はケースが追加されているようです。こちらについても誤差評価を考えてみてはいかがでしょうか。

おまけ

Python の Decimal で精度を上げて計算した場合も、同様にたくさんのケースで hack が可能です。

atcoder.jp

2 2 8
8 8 32
2 8 18

この子たちは、小さいケースながら、多くの精度・丸め設定において誤判定させることができました。

おわり

お疲れさまでした。

*1:複数ある場合は仮数部が偶数であるものとします。いわゆる roundTiesToEven です。

*2:$\roundf(x+y)$ と $\hat x\oplus\hat y$ の ULP は異なるケースがありますが、ここでは小さい方を採用するとします。

*3:多くの問題は、このような “細部” に気を払わずに解けるようになっていますが、そのせいで「これは職人がノリで決める値だ」「よくわからないが適当に書いておけばよい」と思考停止している人が多いのはあまり好ましいことではない気がします。

*4:$2\ulp(\sqrt{a}+\sqrt{b})=\ulp(\sqrt{c})$ となるケースで $k$ は半整数となることがありそうです。

*5:$1.9\times 10^{12}$ 個調べましたがまだ出ていません。見つかり次第追記します。もしかしてこれ $1+2^{-56}$ くらい未満には存在しなかったりしますか？

2023-09-20

s-factorization を求めたり run enumerate を解いたり

文字列アルゴリズム ABC-EF くらい

judge.yosupo.jp

こちらです。これを解くアルゴリズムの一つとして下記があります。

Kolpakov, Roman, and Gregory Kucherov. "Finding maximal repetitions in a word in linear time." In 40th Annual Symposium on Foundations of Computer Science (Cat. No. 99CB37039), pp. 596–604. IEEE, 1999.

これをする過程で s-factorization というのが出てくるのですが、ABC-EF くらいの文字列アルゴリズムの例題としてよさそうな気がしたので紹介です。

導入
s-factorization
- 解法
s-factorization 以外
- 解法
提出
所感
おわり

導入

run というのは、ランレングス圧縮 (RLE) のランと同じで、繰り返しを意味する用語です。

$\gdef\ttstr#1{\text{\texttt{#1}}}$

上記の論文中では run という言い回しはされていなくて、repetition と呼ばれています。長さ $n$ の文字列 $w = a_0a_1\dots a_{n-1}$ の周期とは、任意の $p\le i\lt n$ に対して $a_{i-p} = a_i$ を満たすような最小の正整数 $p$ を指します。たとえば、$\ttstr{aaaa}$ の周期は $1$、$\ttstr{abcabca}$ の周期は $3$、$\ttstr{abcbc}$ の周期は $5$ です。$n$ が $p$ の倍数になっている必要はないことに注意してください。

文字列 $w = a_0a_1\dots a_{n-1}$ と、非負整数 $k, r$ ($0\le r\lt n$) に対して、 $$w^{k+r/n} = \underbrace{ww\ldots w}_{k\text{ times}}\, a_0a_1\dots a_{r-1}$$ と定義します。$|w^{k+r/n}| = |w|\cdot(k+r/n) = nk+r$ であることに気づいておきます。

一般に、$w^{e/n} = w^{\floor{e/n} + (e\bmod n)/n}$ とします。指数部は、$n$ を分母とする非負の分数に対して定義されるものとします。

長さ $n$ の文字列 $w$ が、ある文字列 $s$ と非負整数 $k$ を用いて $w = s^{2+k/n}$ と書けるとき、$w$ を repetition と呼びます。たとえば、$\ttstr{abcabca} = \ttstr{abc}^{7/3}$ と書けるので $\ttstr{abcabca}$ は repetition です。$\ttstr{abcab}$ は repetition ではありません。

与えられた文字列 $s$ の部分文字列のうち、repetition であるものを全て列挙することを考えます。ただし、部分文字列は文字列としてではなく添字の選び方として区別することにします。

しかし、$s = \ttstr{a}^n$ を考えると、$s$ の長さ $2$ 以上の部分文字列全てが repetition であり、最悪 $\Theta(n^2)$ 個の repetition が存在することになります。これは「列挙してください」という問題に対して相性が悪いです（愚直にやってもよいことになってしまうため）。

そこで、maximal という概念を導入します。 $s = a_0a_1\dots a_{n-1}$ の部分文字列 $s_{l, r} = a_l\dots a_{r-1}$ が maximal repetition であるとは、$s_{l, r}$ の周期 $p\ge 2$ に対して、$s_{l-1, r}$ または $s_{l, r+1}$ が $p$ より大きい周期を持つ（または $l-1\lt 0$ または $r\ge n$ が成り立つ）ことを言います。 maximal repetition は $O(n)$ 個しかないことが示されている（これ難しいです）ので、これであれば $O(n)$ 時間で列挙できる可能性があるわけです。

mississippi における maximal repetition の列挙は次のようになります。

mississippi
 ~~~~~~~     # (iss)^{7/3}
  ~~         # (s)^2
     ~~      # (s)^2
        ~~   # (p)^2

以下において、文字列や文字同士の連結を $\concat$ で表します。$\ttstr{ab}\concat\ttstr{c} = \ttstr{abc}$ などです。

s-factorization

文字列 $w = a_0a_1\dots a_{n-1}$ の s-factorization とは、$w = u_0u_1\dots u_{k-1}$ なる文字列の列 $(u_0, u_1, \dots, u_{k-1})$ であって、次のように構成されるものです。ただし、各 $u_i$ を factor と呼びます。便宜上、先頭 $i$ 個の factor の長さを $k_i = |u_0\dots u_{i-1}|$ と表します ($k_0 = 0$)。

手順：先頭 $i$ 個の factor $a_0\dots a_{k_i-1}$ に文字 $a_{k_i}$ が含まれないとき、$u_i = w_{k_i}$ とします（1 文字のみからなる文字列）。そうでないとき、$u_i$ は次の条件を満たす最大の $m$ を用いて $u_i = a_{k_i}\dots a_{k_i+m-1}$ とします ($|u_i| = m$)。

条件：$a_{k_i}\dots a_{k_i+m-1}$ は、$a_0\dots a_{k_i+m-2}$ の部分文字列である。

直感的に言えば、「$u_i$ は $u_0\dots u_i$ に（overlap を許容して）2 回以上現れる」とか「$u_i$ は $a_0\dots a_{k_i-2}$ に 1 回以上現れる」のような感じです。たとえば、$\ttstr{abaabbabaaba} = (\ttstr{a}, \ttstr{b}, \ttstr{a}, \ttstr{ab}, \ttstr{ba}, \ttstr{baab}, \ttstr{a})$ や $\ttstr{aaaaa} = (\ttstr{a}, \ttstr{aaaa})$ です。

解法

さて、これを求めていきます。

$\gdef\prefix#1#2{{#1[\ldots#2]}}$ $\gdef\suffix#1#2{{#1[#2\ldots]}}$ $\gdef\lcp{\operatorname{lcp}}$

文字列 $w = w_0w_1\dots w_{n-1}$ に対して、接頭辞 $j$ を $\prefix wj = w_0\dots w_{j-1}$、接尾辞 $j$ を $\suffix wj = w_j\dots w_{n-1}$ で定義します。$0\le j\le n$ に対して $w = \prefix wj\concat\suffix wj$ が成り立ちます。

また、文字列 $s$ と $t$ の最長共通接頭辞 (LCP) の長さを $\lcp(s, t)$ と書きます。$m = \lcp(s, t)$ としたとき、$\prefix sm = \prefix tm$ かつ、$\prefix s{m+1}\ne \prefix t{m+1}$（または $m+1\gt |s|$ または $m+1\gt |t|$）を満たします。

$u_i$ の構成の手順のうち後者の方（それより前の factor に $w_{k_i}$ が含まれる方）の条件を言い換えます。 $j\lt k_i$ に対して $\suffix wj$ と $\suffix w{k_i}$ の最長共通接頭辞の長さを $m_j$ とするとき、$|u_i| = \max_j m_j$ となります。

というわけで、以下のクエリを処理できればよいことになります。

長さ $n$ の文字列 $w$ が与えられるので、各 $0\lt i\lt n$ に対して下記を求めてください。 $$ \max_{0\le j\lt i} \lcp(\suffix wj, \suffix wi). $$

これは suffix array と LCP array を用いて答えられます。 suffix array において、接尾辞 $i$ に近い方が LCP が長いことから、接尾辞 $i$ から最も近くにある接尾辞 $j$ ($j\lt i$) をなんとかして探し、それとの LCP を計算すればいいことになります。

これは、適切にスタックを管理しながら suffix array を昇順（・降順）に走査することで、$i$ から左側（・右側）に最も近い $j$ を $O(n)$ 時間で探すことができます（そうした $j$ が存在しないケースには注意）。

(0, ..., n-1) の順列 (p_0, ..., p_{n-1}) が与えられます。各 0 ≦ i < n に対して、「j < i かつ p_j < p_i を満たす最大の j（存在しなければその旨）」を求めてください（全体で O(n) 時間で）
— えびちゃん🍑🍝🦃 (@rsk0315_h4x) September 15, 2023

i昇順にstackでp_iより小さいものが出てくるまでpopしてp_iをpushですか
— tute (@tute7627) September 15, 2023

アルファベットサイズが $n^{O(1)}$ であれば、suffix array は SA-IS を用いて $O(n)$ 時間で構築できます。 LCP array や RMQ も $O(n)$ 時間で構築できるので、全体で $O(n)$ 時間でクエリ処理できます。

クエリの答えが $0$ だったとき、s-factorization の構成手順で言うところの前者の方になることに注意すると、s-factorization を $O(n)$ 時間で求められたことになります。

s-factorization 以外

s-factorization を求めた後は、次のような問題を解かされます。

文字列 $t = t_0t_1\dots t_{m-1}$ と $u = u_0u_1 \dots u_{n-1}$ が与えられます。各 $1\le i\le n$ に対して、$t$ と $t\concat \prefix ui$ の最長共通接尾辞の長さを求めてください。

全体で $O(m+n)$ 時間でお願いします。

この $t$, $u$ は s-factorization に基づいて構成される文字列ですが、これを解く上では s-factorization の性質を使うわけではないので、一旦忘れても問題ありません。

解法

$\gdef\rev{\operatorname{rev}}$

求めたいものは、 $$ \begin{aligned} &\phantom{{}={}} {\lcp}(\rev(t), \rev(t\concat\prefix ui)) \\ &= {\lcp}(\rev(t), \rev(\prefix ui)\concat \rev(t)) \end{aligned} $$ と言い換えられます。ここで、 $$ \rev(\prefix u{i+1})\concat \rev(t) = u_i\concat\rev(\prefix ui)\concat \rev(t) $$ であることに注意すると、$\rev(u)\concat\rev(t)$ の各接尾辞に対して、$\rev(t)$ との LCP を求めればよいことがわかります。

よって、$t$ にも $u$ にも現れない文字 $\bot$ を用いて $\rev(t) \concat \bot \concat \rev(u) \concat \rev(t)$ を構成し、それに対して Z algorithm を適用すればよいです。

提出

LC #161740

所感

文字列界隈の線形時間アルゴリズム天才がち。

競プロ頻出のやつは「ふつう」みたいな感想になってきますが、あまり出てこないやつを調べるとびっくりしがちです。

s-factorization を求めてどうするのかとか、所望のものが線形個しかないことの証明とかは論文を読んでほしいですが、記事中で紹介した部分は ABC-E から F くらい（勝手なイメージで 450–475 くらい？）な気がするので、文字列アルゴリズムに苦手意識のある人は練習として解いてみるといい気がします。

おわり

おわりです。

2023-09-19

unsigned の exact division などをうまくやる

おつまみ記事

背景としては、こちらの記事です。

rsk0315.hatenablog.com

1431655766, 1717986920, 1431655782, 1431655768, 1840700272, 1431655966, -1431655056 などさまざまなマジックナンバーが登場します。「$x$ が $3$ の倍数であることがわかっているとき、$1431655766x \equiv \tfrac23 x\pmod{2^{32}}$ となる」のような性質を持ちます。

一般に除算は重い処理なので、結果が同じになるような別の命令（乗算など）に置き換える最適化が典型的に行われるのですが、その過程で出てきた値になります。

今回は、そうした値を自分で求められるようになってみようという回です（コンパイル時定数であれば、多くの場合はコンパイラが勝手に最適化してくれるとは思います）。

exact division というのは、あまりの出ない除算のことを指しています。

本題
例
一般化
- 不等式評価
- 具体例
おわり

本題

まず、上記の記事でも書いていますが、$1431655766 = \tfrac13(2^{32}+2)$ です。$2^{32}+2\equiv 0\pmod{3}$ であることに注意しておきましょう。 $$ \begin{aligned} 3x\times 1431655766 &= 3x\times\tfrac13(2^{32}+2) \\ &= x\times (2^{32}+2) \\ &\equiv x\times(0+2) = 2x \pmod{2^{32}} \end{aligned} $$ となるわけです。

以降、$ax\times n\equiv bx \pmod{2^{32}}$ となることを $ax\xmapsto{n} bx$ と表すことにします。上記で言えば $3x\xmapsto{1431655766}2x$ です。

例

では、手始めに $3x\xmapsto{n}1x=x$ となるような $n$ について考えてみましょう。合同式で $2^{32}$ が消えてくれることだけに目が行って、$\tfrac13(2^{32}+1)$ か？と思ってしまうと危ないです。$2^{32}+1\not\equiv 0\pmod{3}$ なので、$\tfrac13(2^{32}+1)$ が整数になってくれないんですね。

というわけで、$q\cdot 2^{32}+1 \equiv 0\pmod{3}$ なる $q$ があってくれるとうれしいわけです。これは競プロ er ならもう全員わかっていて*1、 $$q \equiv (-1)\cdot (2^{32})^{-1} \pmod{3}$$ ですね。逆元を使いたくなったときは $\gcd(2^{32}, 3) = 1$ であることをちゃんと意識した上でやりましょう。

$3^{-1}\bmod 2^{32}$ であれば Newton 法を使ってうまくやる方法があるのですが、$(2^{32})^{-1} \bmod 3$ だとうまい方法はあるんでしょうか。わかりません。$x^{-1} \bmod y$ がわかっても $y^{-1}\bmod x$ を求めるのには役立たなさそうな気がします。

というわけで Euclid の互除法なりなんなりで求めると、$(2^{32})^{-1} \bmod 3 = 1$ とわかります（そもそも $2^{32}\equiv 1\pmod 3$ なので）。よって、$q \equiv -1 \equiv 2$ とわかります。

すなわち、$3x\xmapsto{\tfrac13(2\cdot 2^{32}+1)}x$ というわけです。${\tfrac13(2\cdot 2^{32}+1)}=2863311531$ です。

一般化

$a=1$ の場合は単に $x\xmapsto{b} bx$ とできるので、以降は除外します。

奇数 $a\gt 1$ に対して、$2^k\cdot ax\xmapsto{n} bx$ なる $n$ を求めたいです。そもそもシフト演算は高速ですから、$(2^k\cdot ax)\gg k = ax$ を用いて、$ax\xmapsto{n}bx$ に帰着させることができます。 $a$ は奇数なので、$\gcd(2^{32}, a) = 1$ が保証されて、うれしいです。

やりたいことは $\tfrac 1a(q\cdot 2^{32}+b)$ の形の整数を探すことで、それは $q\cdot 2^{32}+b \equiv 0 \pmod a$ を解くことによって得られます。 $$q \equiv (-b)\cdot (2^{32})^{-1} \pmod a$$ ですね。

気づいてしまったのですが、 $$ q \equiv (-b) \cdot (2^{-1})^{32} \pmod a $$ であり、$a$ が奇数であることから $2^{-1}\bmod a = \tfrac{a+1}2$ なので、 $$ q \equiv (-b) \cdot (\tfrac{a+1}2)^{32} \pmod a $$ ですね。32 乗の部分は繰り返し二乗法っぽくやりましょう。

warning: $a = 2^{32}-1$ のときは $\tfrac{a+1}2$ の計算過程でオーバーフローしうるので注意しましょう。$\floor{\tfrac a2}+1$ とするのがよいでしょうか。

不等式評価

より注意深くなるために、$\tfrac1a(q\cdot 2^{32}+b)$ の大きさについて見積もっておきましょうか。 $1\le q\le a-1$ であることと、$a$, $b$ が 32-bit 符号なし整数型であることから、 $$ \begin{aligned} \tfrac1a(q\cdot 2^{32} + b) &\le \tfrac1a( (a-1)\cdot 2^{32} + (2^{32}-1)) \\ &= \tfrac{a-1}a\cdot 2^{32} + \tfrac1a\cdot2^{32}-\tfrac1a \\ &\lt 2^{32} \end{aligned} $$ です。整数性から、$\tfrac1a(q\cdot 2^{32}+b)\le 2^{32}-1$ ですね。よって、予め前計算を 64-bit 整数などで行った後は、32-bit に収まる値になっていることがわかります*2。

なお、$a + b = 2^{32}$ となるケースでは $\tfrac1a(q\cdot 2^{32}+b) = 2^{32}-1$ となるようです。 $$ \begin{aligned} q &\equiv (-b) \cdot (2^{32})^{-1} \\ &= -(2^{32}-a)\cdot (2^{32})^{-1} \\ &\equiv -(2^{32})\cdot (2^{32})^{-1} \equiv -1 \pmod{a}. \end{aligned} $$ よって、 $$ \begin{aligned} \tfrac1a(q\cdot 2^{32}+b) &\equiv \tfrac1a(-2^{32}+b) \\ &= \tfrac1a(-a) \\ &= -1 \equiv 2^{32}-1 \pmod {2^{32}} \end{aligned} $$ となります。

具体例

というわけで、もう好きなようにできてしまうわけです。試しに $271x\xmapsto{n} 314x$ となるような $n$ でも求めてみましょうか。 $$ \begin{aligned} q &\equiv (-314)\cdot (\tfrac{271+1}2)^{32} \\ &\equiv 228 \cdot 136^{32} \\ &\equiv 228\cdot 99 \equiv 79 \pmod {271} \end{aligned} $$ と出ました。よって、所望のマジックナンバー $n$ は $\tfrac1{271}(79\cdot 2^{32}+314) = 1252038438$ になりそうです。

>>> 123456760 % 271
0
>>> 123456760 * 1252038438 % 2**32
143045840
>>> 123456760 // 271 * 314
143045840

あっていそうですね。

おわり

こういう類のことを知っていると、コンパイラの出したアセンブリに不思議な数が入っていても驚かずに済んだりしそうです。

コンパイラの出したアセンブリを眺めて気になったところについて考えてみると、新しい発見があって面白いですね。

*1:わからなかった人はこれから競プロ er になりましょう。

*2:とはいえ、仮にオーバーフローしたとしても $2^{32}$ を法として乗算しているだけなので、$n\bmod 2^{32}$ を考えればいいだけのような気もしてきます？（たぶん）。

2023-09-19

Clang の k 乗和の最適化を眺める

ボリューム大きめ個人的新発見あり

Clang が $\sum_{i=0}^n i$ を $n(n+1)/2$ にしてくれることは有名です*1。

また、$\sum_{i=0}^n i^2$ も $n(n+1)(2n+1)/6$ にしてくれます。その過程では、

    unsigned v1 = n * (n - 1) * (n - 2) / 2 * 1431655766u;
    unsigned v2 = n * (n - 1) / 2;
    return 3 * v2 + v1 + n;

のような計算をしていました。ここで、$1431655766 = (2^{32}+2)/3$ であり、 $$ x\equiv 0\pmod{3} \implies (x\times 1431655766)\bmod 2^{32} = \tfrac23 x $$ が成り立ちます。

今回は、$\sum_{i=0}^n i^k$ の最適化を、より大きい $k$ についてやってもらったら別の発見があるのではないか？という記事です。発見がなかったらお蔵入りになる予定だったのですが、発見があったのでよかったです。

追記
見てみる
- 0 乗和
- 1 乗和
- 2 乗和
- 3 乗和
  - ちょっと整理
- 4 乗和
- 5 乗和
- 6 乗和
- 7 乗和
- 8 乗和
- 9 乗和・10 乗和
- 11 乗和
- それ以降
数学
関連資料
所感
- 追記
おわり

追記

lpha-z.hatenablog.com

↑ こっちの記事を読んだ方がいいかも。

下記では、吐かれたアセンブリを観察してエスパーする流れになっていますが、上記の記事では LLVM のソースコードを見ていてえらいです。下記からも学びがあるとは思いますが、Clang・LLVM の最適化手法自体を学ぶというよりは、その最適化手法を見て学びを得るスタイルに近い気がします。

見てみる

godbolt.org

いつもお世話になっております。x86-64 clang 16.0.0 で -O3 に設定します。

アセンブリに関しては下記の記事でざっくり説明したので、ある程度は読めると思って詳しい解説はしません。

rsk0315.hatenablog.com

0 乗和

こういうのは 0 からやりましょう。$0^0 = 1$ ということにします。

unsigned sum(unsigned n) {
    unsigned res = 0;
    for (unsigned i = 0; i <= n; ++i) res += 1;
    return res;
}

sum(unsigned int):
        lea     eax, [rdi + 1]
        ret

$\sum_{i=0}^n i^0 = n + 1$ ということですね。よしよしという感じです。

1 乗和

unsigned sum(unsigned n) {
    unsigned res = 0;
    for (unsigned i = 0; i <= n; ++i) res += i;
    return res;
}

sum(unsigned int):
        mov     ecx, edi
        lea     eax, [rdi - 1]
        imul    rax, rcx
        shr     rax
        add     eax, edi
        ret

$\sum_{i=0}^n i = \tfrac12 n(n-1) + n$ としていそうです。

2 乗和

.pow(2) 的なものが整数にないので、職人が気持ちを込めて一つずつ * i を書いていきます。

unsigned sum(unsigned n) {
    unsigned res = 0;
    for (unsigned i = 0; i <= n; ++i) res += i * i;
    return res;
}

sum(unsigned int):
        mov     eax, edi
        lea     ecx, [rdi - 1]
        imul    rcx, rax
        lea     eax, [rdi - 2]
        imul    rax, rcx
        shr     rax
        imul    edx, eax, 1431655766
        shr     rcx
        lea     eax, [rcx + 2*rcx]
        add     eax, edi
        add     eax, edx
        ret

最後の lea 以降で eax に答えを足していっていますが、次のような感じです。整理パートはただの手の運動です。

$$ \begin{aligned} \sum_{i=0}^n i^2 &= \underbrace{\tfrac32 n(n-1)}_{\text{\texttt{3*rcx}}} + \underbrace{\vphantom{\tfrac12} n}_{\text{\texttt{edi}}} + \underbrace{\tfrac12 n(n-1)(n-2) \times 1431655766}_{\text{\texttt{edx}}} \\ &\equiv \tfrac32 n(n-1) + n + \tfrac13 n(n-1)(n-2) \pmod{2^{32}} \\ &= \tfrac12 n(n-1) + n + n(n-1) + \tfrac13 n(n-1)(n-2) \\ &= \tfrac12 n(n+1) + \tfrac13 n(n-1)(3+(n-2)) \\ &= \tfrac12 n(n+1) + \tfrac13 n(n-1)(n+1) \\ &= \tfrac16 n(n+1)(3+2(n-1)) \\ &= \tfrac16 n(n+1)(2n+1). \end{aligned} $$

3 乗和

unsigned sum(unsigned n) {
    unsigned res = 0;
    for (unsigned i = 0; i <= n; ++i) res += i * i * i;
    return res;
}

sum(unsigned int):
        mov     eax, edi
        lea     ecx, [rdi - 1]
        imul    rcx, rax
        lea     eax, [rdi - 2]
        mov     edx, ecx
        lea     esi, [rdi - 3]
        imul    rsi, rax
        imul    rsi, rcx
        shr     rcx
        lea     r8d, [8*rcx]
        sub     r8d, ecx
        imul    edx, eax
        and     edx, -2
        shr     rsi, 2
        and     esi, -2
        add     r8d, edi
        lea     eax, [r8 + 2*rdx]
        add     eax, esi
        ret

長くなってきました。レジスタも rsi や r8d などが登場してきました。r8d は r8 の下位 32 bits (dword) です。目新しそうなポイントとしては and edx, -2 のあたりでしょうか。

今から人間向けに解釈するので少々お待ちください。各レジスタの最終的な値に基づいて高級言語っぽく書くと、次のようになりました。

using ul = unsigned long;
unsigned sum(unsigned n) {
  unsigned edx = (ul(n - 2) * ul(n - 1) * ul(n)) & -2;
  unsigned long rsi = ul(n - 3) * ul(n - 2) * ul(n - 1) * ul(n) / 4;
  unsigned esi = rsi & -2;
  unsigned long r8d = ul(n - 1) * ul(n) / 2 * 7 + n;
  return ul(r8d) + 2 * ul(edx) + esi;
}

-2 == 0xfffffffe、すなわち ~1（最下位 bit 以外が立っている）です。つまり k & -2 というのは以下を意味します。数式中では & は $\wedge$ で表します。

$$ k \wedge (-2) = \begin{cases} k, & \text{if }k\equiv 0\pmod{2}; \\ k-1, & \text{if }k\equiv 1\pmod{2}. \end{cases} $$

さて、$n(n-1)(n-2)$ は $2$ の倍数ですし、$n(n-1)(n-2)(n-3)/4$ も $2$ の倍数ですから、& -2 は行わなくても値は変わらなさそうです*2。ということでもう少し書き換えます。

using ul = unsigned long;
unsigned sum(unsigned n) {
  unsigned edx = ul(n - 2) * ul(n - 1) * ul(n);
  unsigned long rsi = ul(n - 3) * ul(n - 2) * ul(n - 1) * ul(n) / 4;
  unsigned esi = rsi;
  unsigned long r8d = ul(n - 1) * ul(n) / 2 * 7 + n;
  return ul(r8d) + 2 * ul(edx) + esi;
}

うまいこと整理する方法が思いつかなかったので端折りますが、結果は所望のものになっています。

$$ \begin{aligned} \sum_{i=0}^n i^3 &= \underbrace{\tfrac72 n(n-1)+n}_{\text{\texttt{r8d}}} + \underbrace{\vphantom{\tfrac12} 2n(n-1)(n-2)}_{\text{\texttt{2*edx}}} + \underbrace{\tfrac14 n(n-1)(n-2)(n-3)}_{\text{\texttt{esi}}} \\ &= \tfrac14 n^2(n+1)^2. \end{aligned} $$

ちょっと整理

アセンブリで計算しているものを見るに、 $\gdef\perm#1#2{{{}_{#1}\mathrm{P}_{#2}}}$ $$\sum_{i=0}^n i^k = c_{k, 0}\cdot \perm{n}{k+1} + c_{k, 1}\cdot \perm{n}{k} + \dots + c_{k, k}\cdot \perm{n}{1}$$ のような形式で書けるような $c_{\ast, \ast}$ を求めている感じなのでしょうか。$\perm{n}{j}$ は $j$ 次式であることに注意すると、所望の多項式に最高次から順に定めていくことができて、一意に定まりそうです。少し考えると定数項が $0$ であることもわかります。

ここまでの $c_{k, j}$ の表を書いてみましょう。

$k$ \ $j$	$0$	$1$	$2$	$3$
$1$	$\tfrac12$	$1$	-	-
$2$	$\tfrac13$	$\tfrac32$	$1$	-
$3$	$\tfrac14$	$2$	$\tfrac72$	$1$

ところで、これはもう多項式補間をしてくださいという問題に見えますね。$k$ を固定したとき、$k+1$ 次の多項式になって、定数項を含めて $k+2$ 個の係数を求めたいので、先頭 $k+2$ 個の $k$ 乗和をを求めれば定めることができるわけです。定数項が $0$ なのはわかるので、以下ではそれを除いて考えます。

$(i, j)$ 成分が $\perm{i}{k+1-j}$ であるような $k\times k$ 行列 $A_k$、$i$ 成分が $\sum_{u=0}^i u^k$ であるようなベクトル $b$ に対して、$Ax=b$ なる $x$ が $x=(c_{k, 0}, c_{k, 1}, \dots, c_{k, k})^{\top}$ を満たしそうです。$4$ 乗和の係数を先読みしちゃいましょう。

$$ \begin{pmatrix} 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 2 & 2 \\ 0 & 0 & 6 & 6 & 3 \\ 0 & 24 & 24 & 12 & 4 \\ 120 & 120 & 60 & 20 & 5 \end{pmatrix} \cdot \begin{pmatrix} c_{4, 0} \\ c_{4, 1} \\ c_{4, 2} \\ c_{4, 3} \\ c_{4, 4} \end{pmatrix} = \begin{pmatrix} % 1^4 \\ 1 \\ % 1^4 + 2^4 \\ 17 \\ % 1^4 + 2^4 + 3^4 \\ 98 \\ % 1^4 + 2^4 + 3^4 + 4^4 \\ 354 \\ % 1^4 + 2^4 + 3^4 + 4^4 + 5^4 979 \end{pmatrix} $$

www.wolframalpha.com

$c_4 = (\tfrac15, \tfrac52, \tfrac{25}3, \tfrac{15}2, 1)^{\top}$ とのことです。

$5$ 乗和もやっちゃいましょう。

www.wolframalpha.com

$c_5 = (\tfrac16, 3, \tfrac{65}4, 30, \tfrac{31}2, 1)^{\top}$ とのことです。

表も見直します。

$k$ \ $j$	$0$	$1$	$2$	$3$	$4$	$5$
$1$	$\tfrac12$	$1$	-	-	-	-
$2$	$\tfrac13$	$\tfrac32$	$1$	-	-	-
$3$	$\tfrac14$	$\tfrac63$	$\tfrac72$	$1$	-	-
$4$	$\tfrac15$	$\tfrac{10}4$	$\tfrac{25}3$	$\tfrac{15}2$	$1$	-
$5$	$\tfrac16$	$\tfrac{15}5$	$\tfrac{65}4$	$\tfrac{90}3$	$\tfrac{31}2$	$1$

え、あ、なんで？こわい、これ $c_{k, j} = \frac{{k+1 \brace k+1-j}}{k+1-j}$ でしょうか。正当性もそうなりそうな理由もなにもわかりません、どうしてでしょう。なにかしらをがちゃがちゃすると出てくるのでしょうか。

びっくりして取り乱しましたが、ここで ${n\brace k}$ は第二種 Stirling 数です*3。

正当性については記事の後半で示します。

en.wikipedia.org

n 次多項式 f(x) を x^i の線形結合の形から P(x, i) の線形結合の形に書き換えるのって高速にできたりする？無理かね（必要なら f(0), ..., f(n) の値も得られているのを仮定してよい）
— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2023年9月17日

https://t.co/V04MYjYJIj
O(N(logN)^2)です！！！
— hotman (@hotmanww) 2023年9月17日

$k$ 乗和の多項式自体は多項式補間で $O(n\log(n)^2)$ 時間でできるので、上記を組み合わせることで、${n\brace 0}, {n\brace 1}, \dots, {n\brace n}\pmod{998244353}$ をまとめて $O(n\log(n)^2)$ 時間で求められそうです。びっくりしました（あってますよね？）。

↑ そもそも $n$ を固定した際の第二種 Stirling 数自体は $O(n\log(n))$ 時間で求められました。← それを利用する方法で $k$ 乗和を $O(k\log(k))$ 時間で求めることもできそうです*4。

ともかく、Clang 様が吐くアセンブリがどうなるか予想できるようになったと思います。

4 乗和

気を取り直して、元のコーナーを進めましょう。気づいたんですが、これ += i * i * i * i とか書いてあるだけの C++ のコードは貼る必要がないですね。Clang 様のお出ししたアセンブリがこちらになります。

sum(unsigned int):
        mov     ecx, edi
        lea     eax, [rdi - 1]
        imul    rax, rcx
        lea     ecx, [rdi - 2]
        imul    rcx, rax
        lea     edx, [rdi - 3]
        imul    rdx, rcx
        lea     esi, [rdi - 4]
        imul    rsi, rdx
        shr     rsi, 3
        imul    esi, esi, 1717986920
        shr     rcx
        imul    ecx, ecx, 1431655782
        shr     rdx, 3
        lea     edx, [rdx + 4*rdx]
        shr     rax
        lea     eax, [rax + 4*rax]
        lea     r8d, [rax + 2*rax]
        add     ecx, edi
        add     ecx, esi
        lea     eax, [rcx + 4*rdx]
        add     eax, r8d
        ret

immutable に書くと次のようになります。ecx がすごいことになっています。

unsigned sum(unsigned n) {
  unsigned ecx =
      unsigned(ul(n - 2) * ul(n - 1) * n / 2) * 1431655782u + n +
      unsigned((ul(n - 4) * (n - 3) * (n - 2) * (n - 1) * n) / 8) * 1717986920u;
  unsigned edx = 5 * (ul(n - 3) * (n - 2) * (n - 1) * n / 8);
  unsigned r8d = 3 * 5 *(ul(n - 1) * n / 2);
  return ecx + 4 * edx + r8d;
}

まずは * 1431655782u と * 1717986920u について考えましょう。 $1431655782 = \tfrac13\,(2^{32}+50)$, $1717986920 = \tfrac25\,(2^{32}+4)$ です。つまり、次のようになります。 $$ \begin{aligned} 3x\times 1431655782 &= 3x\times \tfrac13\,(2^{32}+50) \\ &= x\times(2^{32}+50) \\ &\equiv 50x \pmod{2^{32}}, \\ 5x\times 1717986920 &= 5x\times \tfrac25\,(2^{32}+4) \\ &= 2x\times(2^{32}+4) \\ &\equiv 8x \pmod{2^{32}}. \end{aligned} $$

$\tfrac12 n(n-1)(n-2)$ は $3$ の倍数、$\tfrac18 n(n-1)(n-2)(n-3)(n-4)$ は $5$ の倍数であることに注意すると、 $$ \begin{aligned} \sum_{i=0}^n i^4 &= \underbrace{\tfrac{50}3\,\tfrac12\,n(n-1)(n-2) + n + \tfrac85\,\tfrac18\, n(n-1)(n-2)(n-3)(n-4)}_{\text{\texttt{ecx}}} + {} \\ &\phantom{{}={}} \qquad \underbrace{4\cdot\tfrac58 n(n-1)(n-2)(n-3)}_{\text{\texttt{4*edx}}} + \underbrace{3\cdot \tfrac52 n(n-1)}_{\text{\texttt{r8d}}} \\ &= \tfrac15\, \perm{n}{5} + \tfrac52\, \perm{n}{4} + \tfrac{25}3\, \perm{n}{3} + \tfrac{15}2\, \perm{n}{2} + n \end{aligned} $$ となります。

これは先ほど求めた係数と一致しています。展開して $n^i$ の線形結合で表すことはもうしません。

5 乗和

Clang ちゃんはまだ音を上げないみたいです。

sum(unsigned int):
        mov     ecx, edi
        lea     eax, [rdi - 1]
        imul    rax, rcx
        lea     edx, [rdi - 2]
        imul    rdx, rax
        lea     r8d, [rdi - 3]
        imul    r8, rdx
        lea     ecx, [rdi - 4]
        imul    rcx, r8
        lea     esi, [rdi - 5]
        imul    rsi, rcx
        shr     rsi, 4
        imul    esi, esi, 1431655768
        shr     r8, 3
        mov     r9d, r8d
        shl     r9d, 7
        lea     r8d, [r9 + 2*r8]
        shr     rdx
        imul    edx, edx, 60
        shr     rax
        mov     r9d, eax
        shl     r9d, 5
        sub     r9d, eax
        shr     rcx, 3
        lea     eax, [rcx + 2*rcx]
        add     r8d, edi
        add     r8d, edx
        add     r8d, r9d
        add     r8d, esi
        lea     eax, [r8 + 8*rax]
        ret

immutable に直すのは職人が手作業でやっていて、大変です。

unsigned sum(unsigned n) {
  unsigned edi = n;
  unsigned esi =
      ((ul(n - 5) * (n - 4) * (n - 3) * (n - 2) * (n - 1) * n) / 16) *
      1431655768u;
  unsigned r8d = ul(n - 3) * (n - 2) * (n - 1) * n / 8 * 128 +
                 2 * ul(n - 3) * (n - 2) * (n - 1) * n / 8;
  unsigned edx = (ul(n - 2) * (n - 1) * n / 2) * 60;
  unsigned r9d = ul(n - 1) * n / 2 * 31;
  unsigned eax = 3 * ul(n - 4) * (n - 3) * (n - 2) * (n - 1) * n / 8;
  return edi + edx + r9d + esi + r8d + 8 * eax;
}

* 1431655768u について考えます。$1431655768 = \tfrac13\,(2^{32}+8)$ なので、 $$3x\times1431655768 \equiv 8x\pmod{2^{32}}$$ です。慣れたものですね。

$$ \begin{aligned} \sum_{i=0}^n i^5 &= \underbrace{\vphantom{\tfrac12}n}_{\text{\texttt{edi}}} + \underbrace{60\cdot\tfrac12 n(n-1)(n-2)}_{\text{\texttt{edx}}} + \underbrace{31\cdot \tfrac12 n(n-1)}_{\text{\texttt{r9d}}} + {} \\ &\phantom{{}={}}\qquad \underbrace{\tfrac83\tfrac1{16} n(n-1)(n-2)(n-3)(n-4)(n-5)}_{\text{\texttt{esi}}} + {} \\ &\phantom{{}={}}\qquad \underbrace{128\cdot\tfrac18 n(n-1)(n-2)(n-3) + 2\cdot\tfrac18n(n-1)(n-2)(n-3)}_{\text{\texttt{r8d}}} + {} \\ &\phantom{{}={}}\qquad \underbrace{8\cdot 3\cdot\tfrac18 n(n-1)(n-2)(n-3)(n-4)}_{\text{\texttt{8*eax}}} \\ &= \tfrac16\, \perm{n}{6} + 3\, \perm{n}{5} + \tfrac{65}4\, \perm{n}{4} + 30\, \perm{n}{3} + \tfrac{31}2\, \perm{n}{2} + n \end{aligned} $$

先の表と同じになっています。

すごい最適化をしているはずなのに新鮮味がなくなってきましたね。流れがわかってきた証拠です。

6 乗和

もう少し続けます。もう少しで流れが変わるので。

sum(unsigned int):
        mov     eax, edi
        lea     ecx, [rdi - 1]
        imul    rcx, rax
        lea     eax, [rdi - 2]
        imul    rax, rcx
        lea     r8d, [rdi - 3]
        imul    r8, rax
        lea     r9d, [rdi - 4]
        imul    r9, r8
        lea     esi, [rdi - 5]
        imul    rsi, r9
        lea     edx, [rdi - 6]
        imul    rdx, rsi
        shr     rdx, 4
        imul    edx, edx, 1840700272
        shr     rax
        imul    eax, eax, 1431655966
        shr     r8, 3
        imul    r8d, r8d, 700
        shr     r9, 3
        imul    r9d, r9d, 224
        shr     rcx
        mov     r10d, ecx
        shl     r10d, 6
        sub     r10d, ecx
        shr     rsi, 4
        imul    ecx, esi, 56
        add     eax, edi
        add     eax, r8d
        add     eax, r9d
        add     eax, r10d
        add     eax, ecx
        add     eax, edx
        ret

職人も慣れてきたので作業が早くなってきました。最初に各レジスタに $\perm{n}{i}$ を詰めて、あとは賢く係数合わせをするだけですね。

unsigned sum(unsigned n) {
  unsigned edi = n;
  unsigned edx = ul(n - 6) * (n - 5) * (n - 4) * (n - 3) * (n - 2) * (n - 1) *
                 n / 16 * 1840700272u;
  unsigned eax = ul(n - 2) * (n - 1) * n / 2 * 1431655966u;
  unsigned r8d = ul(n - 3) * (n - 2) * (n - 1) * n / 8 * 700;
  unsigned r9d = ul(n - 4) * (n - 3) * (n - 2) * (n - 1) * n / 8 * 224;
  unsigned r10d = ul(n - 1) * n / 2 * 63;
  unsigned ecx =
      ul(n - 5) * (n - 4) * (n - 3) * (n - 2) * (n - 1) * n / 16 * 56;
  return eax + edi + r8d + r9d + r10d + ecx + edx;
}

* 1840700272u と * 1431655966u を考えます。

職人さんは次のようなことをして求めています。

>>> 2**32 / 1840700272  # とりあえず割る
2.3333333304358854
>>> 3 * 2**32 / 1840700272  # 分母に 3 くらいの値がありそうなので 3 を掛ける
6.999999991307656
>>> 7 * 1840700272 % 2**32  # 7 に近いので、7 に掛けたときの挙動を見る
16
>>> divmod((3*2**32 + 16), 7)  # 検算
(1840700272, 0)

$1840700272 = \tfrac17\,(3\cdot 2^{32}+16)$, $1431655966 = \tfrac13\,(2^{32}+602)$ で、 $$ \begin{aligned} 7x\times 1840700272 &= x\times(3\cdot 2^{32}+16) \\ &\equiv 16x \pmod{2^{32}}, \\ 3x\times 1431655966 &= x\times(2^{32} + 602) \\ &\equiv 602x \pmod{2^{32}} \end{aligned} $$ です。

edx を見るに、$\perm{n}{7}/16$ は $7$ の倍数なので、$1840700272$ は $\tfrac{16}7$ と読み替えてよさそうです。同様に eax の $\perm{n}{3}/2$ は $3$ の倍数なので、$1431655966$ は $\tfrac{602}3$ と読み替えられます。

あとは、所望の係数だけ持ってくれば十分でしょう。

$$ \sum_{i=0}^n i^6 = \tfrac17\,\perm n7 + \tfrac72\,\perm n6 + 28\,\perm n5 + \tfrac{175}2\,\perm n4 + \tfrac{301}3\,\perm n3 + \tfrac{63}2\,\perm n2 + n. $$

係数列を低次の方から並べると $(\tfrac11, \tfrac{63}2, \tfrac{301}3, \tfrac{350}4, \tfrac{140}5, \tfrac{21}6, \tfrac17)$ です。あっていそうですね。適宜 Stirling 数の表を調べてください。

7 乗和

残念ですが、まだ流れは変わりません。しかも長いです。

sum(unsigned int):
        mov     eax, edi
        lea     ecx, [rdi - 1]
        imul    rcx, rax
        lea     edx, [rdi - 2]
        imul    rdx, rcx
        lea     eax, [rdi - 3]
        imul    rax, rdx
        lea     esi, [rdi - 4]
        imul    rsi, rax
        shr     rax, 3
        imul    eax, eax, 3402
        lea     r8d, [rdi - 5]
        imul    r8, rsi
        shr     rsi, 3
        imul    esi, esi, 1680
        shr     rdx
        imul    edx, edx, 644
        shr     rcx
        mov     r9d, ecx
        shl     r9d, 7
        sub     r9d, ecx
        lea     ecx, [rdi - 6]
        mov     r10d, r8d
        imul    r10d, ecx
        and     r10d, -16
        lea     r11d, [rdi - 7]
        imul    r11, rcx
        imul    r11, r8
        shr     r11, 3
        and     r11d, -16
        shr     r8, 4
        imul    ecx, r8d, -1431655056
        add     eax, edi
        add     eax, edx
        add     eax, r9d
        add     eax, esi
        lea     eax, [rax + 4*r10]
        add     eax, r11d
        add     eax, ecx
        ret

職人さんはかなり慣れてきましたが、次から流れが変わります。かなしいね。

unsigned sum(unsigned n) {
  unsigned edi = n;
  unsigned eax = ul(n - 3) * (n - 2) * (n - 1) * n / 8 * 3402;
  unsigned esi = ul(n - 4) * (n - 3) * (n - 2) * (n - 1) * n / 8 * 1680;
  unsigned edx = ul(n - 2) * (n - 1) * n / 2 * 644;
  unsigned r9d = ul(n - 1) * n / 2 * 127;
  unsigned r10d =
      ul(n - 6) * (n - 5) * (n - 4) * (n - 3) * (n - 2) * (n - 1) * n;
  // r10d &= -16;
  unsigned long r11 = ul(n - 7) * (n - 6) * (n - 5) * (n - 4) * (n - 3) *
                      (n - 2) * (n - 1) * n / 8;
  // r11 &= -16;
  unsigned ecx =
      ul(n - 5) * (n - 4) * (n - 3) * (n - 2) * (n - 1) * n / 16 * -1431655056u;
  eax += edi;
  eax += edx;
  eax += r9d;
  eax += esi;
  eax += 4 * r10d;
  eax += r11;
  eax += ecx;
  return eax;
}

and r10d, -16 とand r11d, -16 があります。-16 == 0xfffffff0 で、$k \wedge (-16) = \floor{k/16}\cdot 16$ です。しかし、前回同様、r10d は $2\cdot 4\cdot 2=16$ の倍数、r11d も $2\cdot 4\cdot 2\cdot 8/8 = 16$ の倍数なので、何もしないのと同様に見えます。

あとは -1431655056u です。32-bit 符号なしなので 2863312240u と同じです。 $2863312240 = \tfrac23\,(2^{32}+1064)$ なので、$3x\times 2863312240\equiv 2128 \pmod{2^{32}}$ です。

r10d に関しては eax への寄与が 4 * であることに気をつけつつ、求めているものは次のようになります。

$$ \sum_{i=0}^n i^7 = \tfrac18\,\perm n8 + 4\,\perm n7 + \tfrac{133}3\,\perm n6 + 210\,\perm n5 + \tfrac{1701}4\,\perm n4 + 322\,\perm n3 + \tfrac{127}2\,\perm n2 + n. $$

係数列を低次の方から並べると $(\tfrac11, \tfrac{127}2, \tfrac{966}3, \tfrac{1701}4, \tfrac{1050}5, \tfrac{266}6, \tfrac{28}7, \tfrac18)$ です。あっていそうですね。

8 乗和

来ました。流れが変わります。変わった結果、SIMD を使った泥臭い最適化になってしまったので、（がんばって書いたのですが）この節は読まなくてもいいです。

.LCPI0_0:
        .long   0                               # 0x0
        .long   1                               # 0x1
        .long   2                               # 0x2
        .long   3                               # 0x3
.LCPI0_1:
        .long   4                               # 0x4
        .long   4                               # 0x4
        .long   4                               # 0x4
        .long   4                               # 0x4
.LCPI0_2:
        .long   8                               # 0x8
        .long   8                               # 0x8
        .long   8                               # 0x8
        .long   8                               # 0x8
sum(unsigned int):
        inc     edi
        xor     ecx, ecx
        mov     eax, 0
        cmp     edi, 8
        jb      .LBB0_4
        mov     ecx, edi
        and     ecx, -8
        pxor    xmm0, xmm0
        movdqa  xmm1, xmmword ptr [rip + .LCPI0_0] # xmm1 = [0,1,2,3]
        movdqa  xmm3, xmmword ptr [rip + .LCPI0_1] # xmm3 = [4,4,4,4]
        movdqa  xmm4, xmmword ptr [rip + .LCPI0_2] # xmm4 = [8,8,8,8]
        mov     eax, ecx
        pxor    xmm2, xmm2
.LBB0_2:                                # =>This Inner Loop Header: Depth=1
        movdqa  xmm5, xmm1
        paddd   xmm5, xmm3
        movdqa  xmm6, xmm1
        pmuludq xmm6, xmm1
        pshufd  xmm7, xmm1, 245                 # xmm7 = xmm1[1,1,3,3]
        pmuludq xmm7, xmm7
        pshufd  xmm8, xmm5, 245                 # xmm8 = xmm5[1,1,3,3]
        pmuludq xmm5, xmm5
        pmuludq xmm8, xmm8
        pmuludq xmm7, xmm7
        pmuludq xmm6, xmm6
        pmuludq xmm8, xmm8
        pmuludq xmm5, xmm5
        pmuludq xmm6, xmm6
        pshufd  xmm6, xmm6, 232                 # xmm6 = xmm6[0,2,2,3]
        pmuludq xmm7, xmm7
        pshufd  xmm7, xmm7, 232                 # xmm7 = xmm7[0,2,2,3]
        punpckldq       xmm6, xmm7              # xmm6 = xmm6[0],xmm7[0],xmm6[1],xmm7[1]
        paddd   xmm0, xmm6
        pmuludq xmm5, xmm5
        pshufd  xmm5, xmm5, 232                 # xmm5 = xmm5[0,2,2,3]
        pmuludq xmm8, xmm8
        pshufd  xmm6, xmm8, 232                 # xmm6 = xmm8[0,2,2,3]
        punpckldq       xmm5, xmm6              # xmm5 = xmm5[0],xmm6[0],xmm5[1],xmm6[1]
        paddd   xmm2, xmm5
        paddd   xmm1, xmm4
        add     eax, -8
        jne     .LBB0_2
        paddd   xmm2, xmm0
        pshufd  xmm0, xmm2, 238                 # xmm0 = xmm2[2,3,2,3]
        paddd   xmm0, xmm2
        pshufd  xmm1, xmm0, 85                  # xmm1 = xmm0[1,1,1,1]
        paddd   xmm1, xmm0
        movd    eax, xmm1
        jmp     .LBB0_5
.LBB0_4:
        mov     edx, ecx
        imul    edx, ecx
        imul    edx, edx
        imul    edx, edx
        add     eax, edx
        inc     ecx
.LBB0_5:
        cmp     edi, ecx
        jne     .LBB0_4
        ret

$O(1)$ 時間じゃない気配を感じますが、とりあえず解読しましょう。命令もレジスタも新しいものがいろいろ見えます。dq とついている命令がたくさんあります。レジスタも、xmm に番号がついたものが登場しています。

inc 命令は、引数を increment する命令です。cmp 命令は、引数を compare する命令です。

xmm0, xmm1, ... たちは、それぞれ 128-bit のレジスタです。本来は、8-bit 整数 16 個を並列に処理したり、64-bit 浮動小数点数 2 個を並列に処理したりなどいろいろな命令に対応しているのですが、ここでは 32-bit 整数 4 つを並列に処理するのにしか使っていないので、そういう前提で話すとして、xmm0 = [a, b, c, d] のような表記でそれらを表すことにします。

命令の意味合いは次のようになります。

        movdqa  xmm, [e, f, g, h]               # xmm = [e, f, g, h]
        paddd   [a, b, c, d], [e, f, g, h]      # a += e; b += f; c += g; d += h;
        pmuludq [a, b, c, d], [e, f, g, h]      # [a, b] = a * e; [c, d] = c * g;
        pxor    xmm, xmm                        # xmm = [0, 0, 0, 0]
        pshufd  xmm, [e, f, g, h], 245          # xmm = [f, f, h, h]
        punpckldq   [a, b, c, d], [e, f, g, h]  # xmm = [a, e, b, f]

補足が必要そうなのは、pmuludq pshufd punpckldq でしょうか。

pmuludq は、偶数番目の要素同士の積を 64-bit で計算し、結果を格納します。ここでは上位 32-bit ずつは使っていないため、a *= e; c *= g だと解釈しても大丈夫でしょう。

pshufd は、第三引数で指定された添字に従って [e, f, g, h] の要素をコピーします。 $245 = 3311_{(4)}$ なので、[[1], [1], [3], [3]] 番目を取得しています（下位桁が先頭側に来ます）。ここでは、「[0] と [2] に所望の値を入れたい。[1] と [3] はどうでもよい」のような使われ方が多いです。

punpckldq は、先頭側二つの要素を交互に詰めています。

さて、このあたりがわかっていれば概ね読めるでしょう。処理の内容は次のようになっています。

範囲	処理内容
`.LBB0_2` より前	定数やループ回数の初期化
`.LBB0_2` から `.LBB0_4` まで	値 8 つごとにまとめて 8 乗を計算
`.LBB0_4` から `.LBB0_5` まで	端数を 1 つずつ計算
`.LBB0_5` より後	return

$\gdef\register#1{r_{\text{\texttt{#1}}}}$ 数式中では、たとえば ecx の値は $\register{ecx}$ のように表すことにします。

初期化の段階では $\register{ecx} = \floor{\tfrac{n+1}8}\cdot 8$ とし、$8k\lt \register{ecx}$ なる $k$ について $(8k+0)^8 + (8k+1)^8 + \dots + (8k+7)^8$ をまとめて計算する準備をします。

$k$ 回目 ($0\le k\lt \floor{\tfrac{n+1}8}$) に .LBB0_2 に到達した時点では、各レジスタは次のようになっています。 xmm0, xmm2 が出力、xmm1 がループ変数、xmm3, xmm4 がループ変数の増分（ステップ？ stride？）を持っている定数です。 $$ \begin{aligned} \register{xmm0} &= \left[\sum_{i=0}^{k-1} (8i+0)^8, \sum_{i=0}^{k-1} (8i+1)^8, \sum_{i=0}^{k-1} (8i+2)^8, \sum_{i=0}^{k-1} (8i+3)^8\right], \\ \register{xmm2} &= \left[\sum_{i=0}^{k-1} (8i+4)^8, \sum_{i=0}^{k-1} (8i+5)^8, \sum_{i=0}^{k-1} (8i+6)^8, \sum_{i=0}^{k-1} (8i+7)^8\right], \\ \register{xmm1} &= [8k+0, 8k+1, 8k+2, 8k+3], \\ \register{xmm3} &= [4, 4, 4, 4], \\ \register{xmm4} &= [8, 8, 8, 8]. \end{aligned} $$

ループ内では、xmm5 から xmm8 を用いて繰り返し二乗法をしつつ、$\register{xmm6} = [8k+0, 8k+1, 8k+2, 8k+3]$ や $\register{xmm5} = [8k+4, 8k+5, 8k+6, 8k+7]$ を計算しています。

.LBB0_2 のループが終了すると、pshufd を駆使しつつ $\register{xmm0}+\register{xmm2}$ を計算したりして、 $$\register{eax} = \sum_{k=0}^{\floor{\tfrac{n+1}8}\cdot 8-1} k^8$$ として .LBB0_4 のループに向かいます。

.LBB0_4 では、繰り返し二乗法を使いつつ $\register{edx} = k^8$ を求め、$\register{eax}$ に足していきます。 .LBB0_4 は端数に関する処理のため、高々 7 回しか行われません。ecx がループ変数です。

最終的に $\register{eax} = \sum_{i=0}^n i^8$ になり、これを返して終了です。

9 乗和・10 乗和

8 乗和と同じように xmm を使っていました。大変なのでもう解説はしません。目新しい部分はなさそうです。興味のある読者は自分でやってみるとよいでしょう。

11 乗和

xmm が使われなくなりました。やる気がなくなったのでしょうか。あるいはその方が効率がよいと判断したのでしょうか。

sum(unsigned int):
        lea     edx, [rdi + 1]
        test    edi, edi
        je      .LBB0_1
        mov     esi, edx
        and     esi, -2
        xor     ecx, ecx
        xor     eax, eax
.LBB0_6:                                # =>This Inner Loop Header: Depth=1
        mov     edi, ecx
        imul    edi, ecx
        imul    edi, edi
        imul    edi, ecx
        mov     r8d, edi
        imul    r8d, ecx
        imul    r8d, edi
        add     r8d, eax
        lea     eax, [rcx + 1]
        mov     edi, eax
        imul    edi, eax
        imul    edi, edi
        imul    edi, eax
        imul    eax, edi
        imul    eax, edi
        add     eax, r8d
        add     ecx, 2
        cmp     esi, ecx
        jne     .LBB0_6
        test    dl, 1
        je      .LBB0_4
.LBB0_3:
        mov     edx, ecx
        imul    edx, ecx
        imul    edx, edx
        imul    edx, ecx
        imul    ecx, edx
        imul    ecx, edx
        add     ecx, eax
        mov     eax, ecx
.LBB0_4:
        ret
.LBB0_1:
        xor     ecx, ecx
        xor     eax, eax
        test    dl, 1
        jne     .LBB0_3
        jmp     .LBB0_4

皆さんもうある程度読めるようになっていると思われるので、詳細な解説はしません。大まかには、ループ変数 ecx を 2 ずつ増やしていき、各ループでは $\register{ecx}^{11} + (\register{ecx}+1)^{11}$ を計算しています。上限は $\register{esi} = \floor{\tfrac{n+1}{2}}\cdot 2$ です。

ここはあまり特筆すべき点はないかな？と思っていたのですが、そんなことはありませんでした。「累乗なんて繰り返し二乗法などを駆使してオーダーを落とすのは当然でしょう」という感覚が競プロ er 的にはある気がして、コンパイラがオーダーを落としてくれるのを当然がっていましたが、冷静になるとこれも総和同様に賢くやってくれているものの一つですね。

ところで、繰り返し二乗法は乗算回数の観点では最適とは限らないんですよね。 addition-chain exponentiation とか Knuth's power tree とかで調べると楽しそうなものが出てきます。一般に最適な回数を求めるのは NP-complete らしいです。コンパイラが出したコードは繰り返し二乗法に基づいているように見えました。乗算回数を減らそうとすると一時利用のレジスタが増えがちなので不都合なのでしょうか。

↓ 遊んでいた様子

unsigned pow(unsigned n) {
    return 
      n * n * n * n * n * n * n * n * n * n *
      n * n * n * n * n * n * n * n * n * n *
      n * n * n * n * n * n * n * n * n * n *
      n
    ;
}

pow(unsigned int):
        mov     eax, edi
        imul    eax, edi
        imul    eax, edi
        imul    eax, eax
        imul    eax, edi
        imul    eax, eax
        imul    eax, edi
        imul    edi, eax
        imul    eax, edi
        ret

$\gdef\mulgets{\xleftarrow{\times}}$

$\register{eax} \gets \register{edi}$ ($\register{eax} = n^1$)
$\register{eax} \mulgets \register{edi}$ ($\register{eax} = n^2$)
$\register{eax} \mulgets \register{edi}$ ($\register{eax} = n^3$)
$\register{eax} \mulgets \register{eax}$ ($\register{eax} = n^6$)
$\register{eax} \mulgets \register{edi}$ ($\register{eax} = n^7$)
$\register{eax} \mulgets \register{eax}$ ($\register{eax} = n^{14}$)
$\register{eax} \mulgets \register{edi}$ ($\register{eax} = n^{15}$)
$\register{edi} \mulgets \register{eax}$ ($\register{edi} = n^{16}$)
$\register{eax} \mulgets \register{edi}$ ($\register{eax} = n^{31}$)

最後に返す前に edi の方に掛けているのがよくわかりませんが、そうした方が都合がよいのでしょうか。

乗算回数で言えば、たとえば次のようにすれば一回減らせそうです。

$\register{eax} \gets \register{edi}$ ($\register{eax} = n^1$)
$\register{eax} \mulgets \register{edi}$ ($\register{eax} = n^2$)
$\register{eax} \mulgets \register{edi}$ ($\register{eax} = n^3$)
$\register{eax} \mulgets \register{eax}$ ($\register{eax} = n^6$)
$\register{ecx} \gets \register{eax}$ ($\register{ecx} = n^6$)
$\register{eax} \mulgets \register{eax}$ ($\register{eax} = n^{12}$)
$\register{eax} \mulgets \register{eax}$ ($\register{eax} = n^{24}$)
$\register{eax} \mulgets \register{ecx}$ ($\register{eax} = n^{30}$)
$\register{eax} \mulgets \register{edi}$ ($\register{eax} = n^{31}$)

それ以降

12–16 乗和は xmm あり、その後 17–40 乗和までは確認しましたが xmm なしでした。

数学

いろいろ示します。

Lemma 1: $$ \sum_{k=0}^n \textstyle{n \brace k}\, \perm{x}{k} = x^n $$

Proof

組合せ的に解釈します。

$1$ 番から $n$ 番までのボールがあって、各々を $x$ 色のうちのいずれかに塗ることを考えます。これは当然 $x^n$ 通りの塗り方があります。

別の数え方をします。同じ色で塗るボールごとに分けることを考えると、分け方は ${n\brace k}$ 通りあります（${n\brace k}$ は $n$ 要素の集合を $k$ 個の部分集合に分割する通り数なので）。この $k$ 個のグループにどの色を割り当てるかは $\perm{x}{k}$ 通りあります。すなわち、$\sum_{k=0}^n {n\brace k}\, \perm{x}{k}$ 通りです。

よって、$\sum_{k=0}^n {n\brace k}\, \perm{x}{k} = x^n$ となります。$\qed$

Theorem 2: $$ \sum_{i=1}^{k+1} \tfrac{1}{i} {\textstyle{k+1\brace i}}\,\perm{n}{i} = \sum_{i=1}^n i^k. $$

Proof

まず $k = 0$ のとき、 $$ \begin{aligned} \sum_{i=1}^{1} \tfrac{1}{i} {\textstyle{1\brace i}}\,\perm{n}{i} &= \tfrac11 {\textstyle{1\brace 1}}\,\perm{n}{1} \\ &= n = \sum_{i=1}^n 1. \end{aligned} $$

以下、$k\ge 1$ を固定する。左辺が $n$ に関する $k+1$ 次式であることに注意する。

$n = 0$ のとき、 $$ \sum_{i=1}^{k+1} \tfrac{1}{i} {\textstyle{k+1\brace i}}\,\perm{0}{i} = \sum_{i=1}^0 i^k = 0 $$ が成り立つ。次に、 $$ \sum_{j=1}^{k+1} a_j\,\perm nj = \sum_{j=1}^n j^k $$ なる $(a_1, \dots, a_{k+1})$ を考える。

$n = 1$ のとき、$n\lt j\implies \perm nj = 0$ に注意して $$ \sum_{j=1}^{k+1} a_j\,\perm 1j = a_1 = \sum_{j=1}^1 j^k = 1 = \tfrac11 {\textstyle {k+1\brace 1}}. $$

$n = i$ で固定すると、$\perm nj$ が $0$ でない範囲に注意して $$ \sum_{j=1}^{k+1} a_j\,\perm{n}{j} = \sum_{j=1}^i a_j\,\perm{i}{j} $$ となる。すなわち、$a_1, \dots, a_i$ のみの式となる。

$j\lt i$ に対して $a_j = \tfrac 1j {\textstyle{k+1\brace j}}$ が成り立つとき、$a_i = \tfrac 1i {\textstyle{k+1\brace i}}$ が成り立つことを示す。すなわち、 $$ \sum_{j=1}^{i-1} \tfrac 1j {\textstyle{k+1\brace j}}\,\perm ij + a_i\,\perm ii = \sum_{j=1}^i j^k $$ を $a_i$ について解く。

$$ \begin{aligned} a_i &= \frac1{i!} \left(\sum_{j=1}^i j^k - \sum_{j=1}^{i-1} \tfrac 1j {\textstyle{k+1\brace j}}\, \perm ij\right) \\ &= \frac1{i!} \left(i^k + \sum_{j=1}^{i-1} j^k - \sum_{j=1}^{i-1} \tfrac 1j {\textstyle{k+1\brace j}}\, \perm ij\right). \end{aligned} $$ 両辺に $i\cdot i!$ を掛けたり、帰納法の仮定を用いたりして変形を進める。 $$ \begin{aligned} i\cdot i!\cdot a_i &= i^{k+1} + i\sum_{j=1}^{i-1} j^k - i\sum_{j=1}^{i-1} \tfrac 1j {\textstyle{k+1\brace j}}\, \perm ij \\ &= i^{k+1} + i\sum_{j=1}^{i-1} \tfrac 1j{\textstyle{k+1\brace j}}\,\perm {i-1}j - i\sum_{j=1}^{i-1} \tfrac 1j {\textstyle{k+1\brace j}}\, \perm ij \\ &= i^{k+1} + \sum_{j=1}^{i-1} \tfrac 1j{\textstyle{k+1\brace j}}\,(\perm i{j+1} - i\cdot\perm ij). \end{aligned} $$ Lemma 1 を使いつつ、$\perm ij \gt 0$ の範囲や ${k+1\brace 0} = 0$ であることなどに注意して、 $$ \begin{aligned} i\cdot i!\cdot a_i &= \sum_{j=0}^{k+1} {\textstyle{k+1 \brace j}}\, \perm ij + \sum_{j=1}^{i-1} \tfrac 1j{\textstyle{k+1\brace j}}\,(\perm i{j+1} - i\cdot\perm ij) \\ &= {\textstyle {k+1\brace i}}\,\perm ii + \sum_{j=1}^{i-1} {\textstyle{k+1 \brace j}}\, \perm ij + \sum_{j=1}^{i-1} \tfrac 1j{\textstyle{k+1\brace j}}\,(\perm i{j+1} - i\cdot\perm ij) \\ &= {\textstyle {k+1\brace i}}\,i! + \sum_{j=1}^{i-1} \tfrac 1j{\textstyle{k+1\brace j}}\,\underbrace{(j\cdot\perm ij + (i-j)\cdot\perm ij - i\cdot\perm ij)}_0. \\ \end{aligned} $$ これにより、$a_i = \tfrac1i {\textstyle {k+1\brace i}}$ を得る。

$1\le i\le k+1$ に対して、$n=i$ の際に等式が成り立つように $a_i$ を定めたので、$n=0$ のケースと合わせて $k+2$ 個の点で等式が成り立っている。左辺は $k+1$ 次の多項式なので、任意の $n$ について等式が成り立つ。

よって、任意の $k\ge 1$ についても示されたので、任意の $k\ge 0$ に対して $$ \sum_{i=1}^{k+1} \tfrac{1}{i} {\textstyle{k+1\brace i}}\,\perm{n}{i} = \sum_{i=1}^n i^k $$ が成り立つ。$\qed$

Lemma 1 によって ${\textstyle{k+1\brace j}}$ の形を作るために、両辺に $i$ を掛けたところがおもしろかったです。

記事の冒頭では $i=0$ から足していましたが、$0$ 乗和において上式から自然に出てくる値を考慮すると、$i=1$ から足す方がよいかなとなってそうしました。

所感

元々、Clang が直接 $\tfrac12 n(n+1)$ を計算しているのではなく何かしら特殊なことをしていそうなことは知っていたのですが、実際にやってみるとたしかに便利そうな形でやっていそうで納得でした。第二種 Stirling 数であれば DP で手軽に求められますし、$\tfrac1a(b\cdot 2^{32}+c)$ などの形の定数を用意してコンパイラが最適化するのもよくあることだと思うので、そういう最適化をしてくれるの自体はなるほどなぁという感じです。そうまでして $k$ 乗和を最適化しようとしてくれるのはすごいなと思います。

内部実装については知らないので実際に DP で求めているかどうかなどは知りません。単なる $k$ 乗和でなく $k$ 次式の場合でもうまくやっているはずですし、別のうまい方法などがあるのかもしれません。

また、この記事では unsigned での最適化を検証しましたが、$k\ge 31$ のとき $1^k+2^k\ge 2^{31}$ ですから、signed であれば $n\ge 2$ のときオーバーフローして未定義になるはずです。$n\le 0$ であれば $0$、$n = 1$ のとき $1$ なので、次のような最適化も可能なはずです。

signed sum(signed n) {
    // 1 以上 n 以下の 31 乗和を返す
    return n > 0;
}

試した限りでは、このような最適化は行われていないようでした。

また、同様の概念として Bernoulli 数というものもあると記憶しているのですが、最適化をやる上ではあまり相性がよくないのかな？と思いつつ、実はちゃんと調べていません。

今回、アセンブリを元にして 3 乗和から 7 乗和で C++ 風のコードを書きましたが、それに関しては 0 から -1u まで任意の値を入力して、愚直な方法と値が同じになることは検証いたしました。

追記

lpha-z.hatenablog.com

↑ こっちの記事では、（たぶん、）より広く $k$ 次式和に関して自然に扱えそうな気がします。

おわり

ここ最近は重めの記事をぽんぽん書いていて、読者の人々が大変そうですね。

*1:私は知っていますという程度の意味。

*2:オーバーフローした場合でもそれは変わらないですし、なぜやっているのかはよくわかりませんでした。何らかの事情があるか、見落としがあるかだと思います。

*3:行き当たりばったりで書いていたので本当に取り乱しました。右端が $1$ になっていることとその隣が $\tfrac12(2^k-1)$ になっていること、あと左端が $\tfrac1{k+1}$ になっていることはすぐ気づくと思いますが、分母を $k+1-j$ に揃えたらどうかというのは $k=5$ を計算したあたりまで気づきませんでした。3, 6, 7 の並びに既視感があって Stirling 数の表を見に行ってびっくりしました。

*4:とはいえ線形篩を使えば $1\le i\le k+1$ に対する $i^k$ たちを $O(k)$ 時間で求められるので、等差数列になっている場合の多項式補間をちゃんとやったりすることで $O(k)$ 時間を達成できそうです。

記法・前提

考察

実装例

あとがき

おわり

記法・前提知識

解説

疑問・興味

疑問 1 の解答

疑問 2 の解答

あとがき

おまけ

おわり

答え

解説など

一般のケース

速度に関して

おきもち

あとがき

おわり

導入

定義・前提知識

方針

考察

$\varepsilon_L'$ の上界

$\varepsilon_U'$ の下界

hack について

ケースに関して

差の下界について

まとめ

参考

感情

おまけ

おわり

導入

s-factorization

解法

s-factorization 以外

解法

提出

所感

おわり

本題

例

一般化

不等式評価

具体例

おわり

追記

見てみる

0 乗和

1 乗和

2 乗和

3 乗和

ちょっと整理

4 乗和

5 乗和

6 乗和

7 乗和

8 乗和

9 乗和・10 乗和

11 乗和

それ以降

数学

関連資料

所感

追記

おわり