Complex FFT is (not) as bad as you think - えびちゃんの日記

これに触発されていろいろ調べたりがんばったりした記録です。

tl;dr: フ…上等だ…オレも一つ言っておくことがある
f64 の FFT でも誤差がいい感じに収まるような気がしていたが別にそんなことはなかったぜ！

ウオオオいくぞオオオ！

ご愛読ありがとうございました！*1

背景
本題
関連資料
あとがき
おわり

背景

多項式の積を求めるときは、離散 Fourier 変換 (discrete Fourier transform; DFT) を介して行うのが典型的です。この変換を高速に求めるアルゴリズムは、総称して 高速 Fourier 変換 (fast Fourier transform; FFT) と呼ばれています。都合のいい複素数を持ってきて行うのが典型的ですが、実部と虚部が（整数とは限らない一般の）実数になるため、コンピュータで高速に扱おうとする上では精度が心配になります。

多項式の積を求める際に、各係数をたとえば $998244353$ で割ったあまりを求めればよい場合*2は、$\Complex$ ではなく $\mathbb F_{998244353}$ 上で都合のいい元が取れるため、整数の演算のみで完結させることができます。特に 数論変換 (number-theoretic transform; NTT) と呼ばれますが、本質的なアルゴリズムは同じです。

ところで、$\mathbb F_{1000000007}$ を考えたい場合はしばしばありますが、$\mathbb F_{1000000007}$ では都合のいい元が取れないため数論変換ができません。代わりに、$\mathbb F_{754974721}$, $\mathbb F_{167772161}$, $\mathbb F_{469762049}$ などでそれぞれ数論変換を行い、中国剰余定理 (Chinese remainder theorem; CRT) などに基づいて $\mathbb F_{1000000007}$ での値を求める方法などが知られています。しかし、何度も数論変換をする必要があったり、並列化するのが大変だったり、高速化するにあたって Montgomery 乗算をしようとするとややこしくなってデバッグがしにくくなったりなど、うれしくない気持ちになることがしばしばあります。

そういった事情などから、古き良き $\Complex$ 上での FFT を好んでおられる方も多いです。多いのかな、たぶんある程度はいる気がします。冒頭に紹介した記事の筆者も同様のようです。

本題

素朴に f64 で FFT をしようとすると、誤差がとんでもないことになり、とても使えたものではありません。 NTT が想定されている問題のテストケース名に fft_killer と書かれていることもしばしばあります。

以下では、まず FFT の概要をおさらいした後、FFT を用いた畳み込みの（素朴な）方法について述べます。誤差を減らすための工夫や高速化のための工夫の例をいくつか挙げます。ここでは、FFT の高速化（キャッシュの考慮や非再帰など）ではなく、FFT をブラックボックスとして畳み込みをする場合の工夫を対象とします。また、誤差評価についても多少触れますが、（冒頭で匂わせた通り）結果的には反例がある形となったため、詳細には立ち入りません。

FFT について

下記では $\Complex$ 上の DFT について考えます。複素数列 $z = (z_0, z_1, \dots, z_{n-1})\in\Complex^n$ と、$1$ の原始 $n$ 乗根 $\omega$ に対して、DFT $F_{\omega}[z]\in\Complex^n$ は次のように定義されます。 $$ \left(\begin{matrix} F_{\omega}[z]_0 \\ F_{\omega}[z]_1 \\ \vdots \\ F_{\omega}[z]_{n-1} \end{matrix}\right) = \left(\begin{matrix} \omega^0 & \omega^0 & \cdots & \omega^0 \\ \omega^0 & \omega^1 & \cdots & \omega^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ \omega^0 & \omega^{n-1} & \cdots & \omega^{(n-1)^2} \end{matrix}\right) \left(\begin{matrix} z_0 \\ z_1 \\ \vdots \\ z_{n-1} \end{matrix}\right) $$ 左上成分を $(0, 0)$ 成分として、$(i, j)$ 成分が $\omega^{ij}$ で表されるような $n$ 次正方行列を $V_n(\omega)$ とすると、 $$ F_{\omega}[z]^{\top} = V_n(\omega) z^{\top} $$ と書けます。この $V_n(\omega)$ が、非常に疎な行列たちの（十分に少ない個数の）積で表せることに基づくアルゴリズムがあります。

$\gdef\krtimes{\stackrel{\text{Kr}}{\otimes}}$ $\gdef\diag{\operatorname{diag}}$

以降、ある非負整数 $t$ に対して $n = 2^t$ とします（必要なら、余分な $0$ で埋めることで帰着可能）。下記で定義される $n$ 次正方行列 $A_k$ ($1\le k\le t$) と $P_t$ を用いて、 $$ V_n(\omega) = A_t\cdots A_2 A_1 P_t $$ と書けます。ここで、$I_n$ を $n$ 次単位行列、$\krtimes$ を Kronecker 積、$\diag_n(a_j)$ を $\angled{a_0, \dots, a_{n-1}}$ を対角成分に持つ $n$ 次正方行列とし、 $$ \begin{aligned} A_k &= I_{n/2^k} \krtimes \left(\begin{matrix} I_{2^{k-1}} & \hphantom{-{}}\Omega_{2^{k-1}} \\ I_{2^{k-1}} & -\Omega_{2^{k-1}} \\ \end{matrix}\right), \\ \Omega_{2^{k-1}} &= \diag_{2^{k-1}}(\omega^{jn/2^k}) \end{aligned} $$ とします。$P_t$ は、各 $b_j\in\{0, 1\}$ に対して第 ${b_{2^t-1}\dots b_1b_0}_{(2)}$ 成分を第 ${b_0b_1\dots b_{2^t-1}}_{(2)}$ 成分に移す順列とします。

$n = 8 = 2^3$ の例

零成分は空白で表します。 $$ \begin{aligned} P_3 &= \left(\begin{matrix} 1 & & & & & & & \\ & & & & 1 & & & \\ & & 1 & & & & & \\ & & & & & & 1 & \\ & 1 & & & & & & \\ & & & & & 1 & & \\ & & & 1 & & & & \\ & & & & & & & 1 \\ \end{matrix}\right), \\ A_1 &= I_{8/2^1} \krtimes \left(\begin{matrix} I_{2^0} & \hphantom{-{}}\diag_{2^0}(\omega^{j\cdot 8/2^1}) \\ I_{2^0} & -\diag_{2^0}(\omega^{j\cdot 8/2^1}) \end{matrix}\right) \\ &= I_4 \krtimes \left(\begin{matrix} 1 & \hphantom{-{}}1 \\ 1 & -1 \end{matrix}\right) \\ &= \left(\begin{matrix} \omega^0 & \omega^0 & & & & & & \\ \omega^0 & \omega^4 & & & & & & \\ & & \omega^0 & \omega^0 & & & & \\ & & \omega^0 & \omega^4 & & & & \\ & & & & \omega^0 & \omega^0 & & \\ & & & & \omega^0 & \omega^4 & & \\ & & & & & & \omega^0 & \omega^0 \\ & & & & & & \omega^0 & \omega^4 \\ \end{matrix}\right), \\ A_2 &= I_{8/2^2} \krtimes \left(\begin{matrix} I_{2^1} & \hphantom{-{}}\diag_{2^1}(\omega^{j\cdot 8/2^2}) \\ I_{2^1} & -\diag_{2^1}(\omega^{j\cdot 8/2^2}) \end{matrix}\right) \\ &= I_2 \krtimes \left(\begin{matrix} 1 & & \hphantom{-{}}1 & \\ & 1 & & \hphantom{-{}}\omega^2 \\ 1 & & -1 & \\ & 1 & & -\omega^2 \\ \end{matrix}\right) \\ &= \left(\begin{matrix} \omega^0 & & \omega^0 & & & & & \\ & \omega^0 & & \omega^2 & & & & \\ \omega^0 & & \omega^4 & & & & & \\ & \omega^0 & & \omega^6 & & & & \\ & & & & \omega^0 & & \omega^0 & \\ & & & & & \omega^0 & & \omega^2 \\ & & & & \omega^0 & & \omega^4 & \\ & & & & & \omega^0 & & \omega^6 \\ \end{matrix}\right), \\ A_3 &= I_{8/2^3} \krtimes \left(\begin{matrix} I_{2^2} & \hphantom{-{}}\diag_{2^2}(\omega^{j\cdot 8/2^3}) \\ I_{2^2} & -\diag_{2^2}(\omega^{j\cdot 8/2^3}) \end{matrix}\right) \\ &= I_1 \krtimes \left(\begin{matrix} 1 & & & & \hphantom{-{}}1 & & & \\ & 1 & & & & \hphantom{-{}}\omega & & \\ & & 1 & & & & \hphantom{-{}}\omega^2 & \\ & & & 1 & & & & \hphantom{-{}}\omega^3 \\ 1 & & & & -1 & & & \\ & 1 & & & & -\omega & & \\ & & 1 & & & & -\omega^2 & \\ & & & 1 & & & & -\omega^3 \\ \end{matrix}\right) \\ &= \left(\begin{matrix} \omega^0 & & & & \omega^0 & & & \\ & \omega^0 & & & & \omega^1 & & \\ & & \omega^0 & & & & \omega^2 & \\ & & & \omega^0 & & & & \omega^3 \\ \omega^0 & & & & \omega^4 & & & \\ & \omega^0 & & & & \omega^5 & & \\ & & \omega^0 & & & & \omega^6 & \\ & & & \omega^0 & & & & \omega^7 \\ \end{matrix}\right). \quad\eod \end{aligned} $$

疎であることから、ベクトルとの積を高速に計算できます（右側から計算していけばよいということ）。実装の際には、陽に（たとえば $n\times n$ 要素の Vec などを使って）行列を表すのではなく、然るべき添字を考えて計算します。

また、$j$ 行目の非零成分が $j$ 列目と $j'$ 列目のとき、$j'$ 行目の非零成分も $j$ 列目と $j'$ 列目にあることがわかります。これを利用して、ベクトルを in-place に更新していくことができます。

実装例

fn fft(a: &mut [Complex]) {
    let n = f.len();

    let roots = {
        let mut roots = vec![(Complex::from(0.0)), Complex::from(1.0)];
        roots.reserve(n - 2);
        while roots.len() < n {
            let m = roots.len();
            for i in m / 2..m {
                roots.push(roots[i]);
                let k = (2 * i + 1 - m) * ((1 << 19) / m);
                // cos(k pi/2^{19}) + i sin(k pi/2^{19})
                roots.push(omega_2_20_k(k));
            }
        }
        roots
    };

    let rev = {
        let mut rev = vec![0; n];
        let shl = n.trailing_zeros() - 1;
        for i in 0..n {
            rev[i] = (i & 1) << shl | rev[i >> 1] >> 1;
        }
        rev
    };

    for i in 0..n {
        if i < rev[i] {
            f.swap(i, rev[i]);
        }
    }

    for k in (0..n.trailing_zeros()).map(|k| 1 << k) {
        for i in (0..n).step_by(2 * k) {
            for j in 0..k {
                // roots[j + k] = omega_{2k}^j
                let z = f[i + j + k] * roots[j + k];
                f[i + j + k] = f[i + j] - z;
                f[i + j] = f[i + j] + z;
            }
        }
    }
}

畳み込みについて

FFT は DFT をどう計算するかのレイヤーの話で、この節では DFT を使ってなにを計算するかに焦点を当てていきます。すなわち、FFT のことは一旦忘れても大丈夫です。

列 $a = (a_0, a_1, \dots, a_{n-1})$ に対する DFT の第 $l$ 成分は次のように表せます。 $$ \begin{aligned} F_{\omega}[a]_l &= \sum_{j=0}^{n-1} a_j \omega^{jl}. \end{aligned} $$ 列 $b = (b_0, b_1, \dots, b_{n-1})$ に対しても同様です。これらの積は次のようになります。 $$ \begin{aligned} F_{\omega}[a]_l \cdot F_{\omega}[b]_l &= \left(\sum_{j=0}^{n-1} a_j \omega^{jl}\right) \left(\sum_{j=0}^{n-1} b_j \omega^{jl}\right) \\ &= a_0 b_0 \omega^{0l} \\ & \qquad\quad + (a_0 b_1 + a_1 b_0) \omega^{1l} \\ & \qquad\quad + (a_0 b_2 + a_1 b_1 + a_2 b_0) \omega^{2l} \\ & \qquad\quad + \cdots \\ & \qquad\quad + (a_0 b_{n-1} + \dots + a_{n-1} b_0) \omega^{(n-1)l} \\ & \qquad\quad + (a_1 b_{n-1} + \dots + a_{n-1} b_1) \omega^{nl} \\ & \qquad\quad + \cdots \\ & \qquad\quad + a_{n-1} b_{n-1} \omega^{2(n-1)l} \\ &= \sum_{k=0}^{n-1} \sum_{j=0}^k a_j b_{k-j} \omega^{kl} + \sum_{k=0}^{n-2} \sum_{j=k+1}^{n-1} a_j b_{k+n-j} \omega^{(k+n)l} \\ &= \sum_{k=0}^{n-1} \left( \sum_{j=0}^{n-1} a_j b_{(k-j)\bmod n} \right) \omega^{kl}. \end{aligned} $$ これは、$k$ 番目の要素が $c_k = \sum_{j=0}^{n-1} a_j b_{(k-j)\bmod n}$ で表される列 $c$ に対する $F_{\omega}[c]_l$ であると見ることができます。 $c_k = \sum_{j=0}^k a_j b_{k-j}$ でないのが気になりますが、たとえば、十分な個数の $0$ を列の末尾に追加することで帰着できます。

あとは、$F_{\omega}[c]$ から $c$ を復元できればよいです。これは、DFT の行列 $V_n(\omega)$ の逆行列が $V_n(\omega)^{-1} = \tfrac1n V_n(\omega^{-1})$ を満たすことから計算できます。

Claim 1: $V_n(\omega)^{-1} = \tfrac1n V_n(\omega^{-1})$ が成り立つ。

Proof

To-be-proved 1: $j=k \implies \sum_{l=0}^{n-1} \omega^{-jl}\cdot \omega^{lk} = n$.

$$ \begin{aligned} \sum_{l=0}^{n-1} \omega^{-jl}\cdot \omega^{lk} &= \sum_{l=0}^{n-1} \omega^{(k-j)l} \\ &= \sum_{l=0}^{n-1} \omega^0 \\ &= \sum_{l=0}^{n-1} 1 \\ &= n. \end{aligned} $$

To-be-proved 2: $j\ne k \implies \sum_{l=0}^{n-1} \omega^{-jl}\cdot \omega^{lk} = 0$.

$$ \begin{aligned} \sum_{l=0}^{n-1} \omega^{-jl}\cdot \omega^{lk} &= \sum_{l=0}^{n-1} \omega^{(k-j)l} \\ &= \frac{1-\omega^{(k-j)n}}{1-\omega^{k-j}} \\ &= \frac{1-(\omega^n)^{k-j}}{1-\omega^{k-j}} \\ &= \frac{1-1}{1-\omega^{k-j}} \\ &= 0. \quad\qed \end{aligned} $$

これらより、長さ $2n$ の列に対する FFT を $3$ 回行うことで畳み込みが求まります。

工夫について

FFT を用いて畳み込みをする際の工夫について考えていきます。

工夫①：虚部も使う

さて、たとえば $\mathbb F_{1000000007}$ 上での畳み込みを上記の枠組みで行うとき、$a_k$ や $b_k$ の虚部は $0$ ということになります。 DFT は複素数上での計算なのに、これはちょっともったいないです。ということで、$F_{\omega}[a+bi]$ を考えてみます。 $$ \begin{aligned} F_{\omega}[a+bi]_l &= \sum_{j=0}^{n-1} {(a_j+b_ji) \omega^{jl}} \\ &= \sum_{j=0}^{n-1} a_j \omega^{jl} + i\sum_{j=0}^{n-1} b_j \omega^{jl} \\ &= F_{\omega}[a]_l + F_{\omega}[b]_l \cdot i. \end{aligned} $$ ここから $F_{\omega}[a]_l$ と $F_{\omega}[b]_l$ を得るにはどうすればいいでしょうか？ $F_{\omega}[a]_l$ や $F_{\omega}[b]_l$ は実数ではないので、$F_{\omega}[a+bi]_l$ の実部や虚部を取り出してもうまくいきません。

ここで、$|\omega| = 1$ より $\omega\cdot\omega^{\ast} = 1$ で、$\omega^{-1} = \omega^{\ast}$ であることを考えます。$\bullet^{\ast}$ で複素共軛を表しています。これを踏まえて $F_{\omega}[a+bi]_{(n-l)\bmod n}^{\ast}$ を考えます。強引な流れです。 $$ \begin{aligned} F_{\omega}[a+bi]_{(n-l)\bmod n}^{\ast} &= \left( \sum_{j=0}^{n-1} {(a_j+b_ji) \omega^{j(n-l)}} \right)^{\ast} \\ &= \sum_{j=0}^{n-1} {(a_j+b_ji)^{\ast} (\omega^{-jl})^{\ast}} \\ &= \sum_{j=0}^{n-1} {(a_j-b_ji) \omega^{jl}} \\ &= F_{\omega}[a]_l - F_{\omega}[b]_l \cdot i. \end{aligned} $$ ここまで来ればどうにでもなりますね。 $$ \begin{aligned} F_{\omega}[a+bi]_l + F_{\omega}[a+bi]_{(n-l)\bmod n}^{\ast} &= 2 F_{\omega}[a]_l \\ F_{\omega}[a+bi]_l - F_{\omega}[a+bi]_{(n-l)\bmod n}^{\ast} &= 2i F_{\omega}[b]_l \end{aligned} $$ より $$ F_{\omega}[a]_l\cdot F_{\omega}[b]_l = \frac{(F_{\omega}[a+bi]_l)^2 - (F_{\omega}[a+bi]_{(n-l)\bmod n}^{\ast})^2}{4i} $$ などとして、所望の積が得られます*3。

これらより、長さ $2n$ の列に対する FFT を $2$ 回行うことで畳み込みが求まります。

工夫②：値を平方根付近で分ける

note: $10^9$ 程度の大きさの整数の積を $2^{19}$ 個足すと $5\cdot 10^{23}$ 程度になるが、これは f64 で誤差なく表せる整数の大きさを上回っており、計算途中の誤差を差し置いてもどうしようもないことがわかる。

$\mathbb F_p$ での畳み込みをしたい前提とします*4。ある $d$ を持ってきて、各 $k$ に対して $a_k \equiv a_k^{(0)} + a_k^{(1)} d \pmod p$ と $b_k \equiv b_k^{(0)} + b_k^{(1)} d \pmod p$ が成り立つような列 $a^{(0)}$, $a^{(1)}$, $b^{(0)}$, $b^{(1)}$ を考えます。このとき、 $$ \begin{aligned} &\phantom{{}={}} (a\times b)(x) \\ &= ( (a^{(0)}+a^{(1)}d)\times (b^{(0)}+b^{(1)}d) )(x) \\ &= (a^{(0)}\times b^{(0)})(x) + d\cdot (a^{(0)}\times b^{(1)} + a^{(1)}\times b^{(0)})(x) + d^2\cdot (a^{(1)}\times b^{(1)})(x) \end{aligned} $$ が成り立つため、下記のように（工夫①も適宜使いつつ）長さ $2n$ の列に対する FFT を $4$ 回行うことで畳み込みが求まります。$\odot$ は要素ごとの積です。

$F_{\omega}[a^{(0)}+b^{(0)}i]$
$F_{\omega}[a^{(1)}+b^{(1)}i]$
$F_{\omega^{-1}}[(F_{\omega}[a^{(0)}]\odot F_{\omega}[b^{(0)}])+(F_{\omega}[a^{(0)}]\odot F_{\omega}[b^{(1)}])\, i]$
$F_{\omega^{-1}}[(F_{\omega}[a^{(1)}]\odot F_{\omega}[b^{(0)}])+(F_{\omega}[a^{(1)}]\odot F_{\omega}[b^{(1)}])\, i]$

さて、$a_k = a_k^{(0)} + a_k^{(1)}d$ を決めるパートについて述べます。

まず、$0\le a_k\le p/2$ の場合を考えます。 $d = \ceil{\sqrt p}$ とし、 $$ (a_k^{(0)}, a_k^{(1)}) = \begin{cases} ( (a_k\bmod d), \floor{a_k/d}), & \text{if}\: (a_k\bmod d)\le d/2; \\ ( (a_k\bmod d)-d, \floor{a_k/d}+1), & \text{if}\: (a_k\bmod d)\gt d/2 \end{cases} $$ で定義します。$p/2\lt a_k\lt p$ の場合は、$p-a_k$ に対して上記を求め、符号をそれぞれ反転したものとします。これにより、 $$ \begin{aligned} |a_k^{(0)}| &\le \ceil{\sqrt p}/2, \\ |a_k^{(1)}| &\le \floor{\floor{p/2}/\ceil{\sqrt p}}+1, \\ % |a_k^{(0)} + a_k^{(1)}i| &\le \sqrt2\left(\ceil{\sqrt p}/2+1\right) |a_k^{(0)} + a_k^{(1)}i| &\le \ceil{\sqrt p}/\sqrt2+\sqrt2 \end{aligned} $$ が成り立ちます。$p$ が平方数でないときは $\ceil{\ceil{\sqrt p}/\sqrt{2}}$ で押さえられそうな気配がありましたが、あまり考えていません。

たとえば、$p = 1000000007$ と $a_k = 123456789$ に対しては、 $$ (a_k^{(0)}, a_k^{(1)}) = (597, 3904) $$ となります。また $p = 1000000007$ と $a_k = 499975442$ に対しては $$ (a_k^{(0)}, a_k^{(1)}) = (-15811, 15811) $$ となり、これが絶対値の最大ケースです。

$10^{4.5}$ 程度の整数の積の $2^{19}$ 個の和は $5\cdot 10^{14}$ 程度で、計算途中の誤差が十分に小さければ、f64 でもなんとかなりそうな気もしてきます。

工夫③：乱数を掛ける

え〜 😮‍💨 ここはあまり好きじゃない人もいると思います。えびちゃんはそうです。

法と互いに素な乱数 $\xi$ を取ってきます。法が素数 $p$ であれば、$\{1, 2, \dots, p-1\}$ から一様ランダムに取ってくればよいです。これに対して、下記が成り立ちます。 $$ \begin{aligned} (a\times b)(x) &= \sum_{k=0}^{n-1} \left(\sum_{j=0}^k a_j b_{k-j}\right) x^k \\ &= \sum_{k=0}^{n-1} \left(\sum_{j=0}^k a_j b_{k-j}\cdot \xi^k\right) \xi^{-k} x^k \\ &= \sum_{k=0}^{n-1} \left(\sum_{j=0}^k (a_j \xi^j) (b_{k-j} \xi^{k-j}) \right) (\xi^k)^{-1} x^k \\ \end{aligned} $$ すなわち、係数 $a_k$, $b_k$ をそれぞれ $(a_k\xi^k \bmod p)$, $(b_k\xi^k \bmod p)$ で置き換えた列の畳み込みを求めてから、$k$ 番目の項に $(\xi^k)^{-1}\bmod p$ を掛けて復元することができるということです。

ランダムケースにおいては誤差が大きくならなさそうっぽいという観測に基づいていますが、$\xi$ を固定すれば最悪ケースを容易に構成できるため、注意が必要です。ランダム性を入力側（入力全体からなる集合）からアルゴリズム側に移しているのに相当する気がします。平均計算量ではなく期待計算量を考えるときと似たような議論になる気がします。クイックソートをする前に一旦入力をシャッフルしておくみたいな話が数年前に話題になっていたような記憶もあります。

工夫④：虚部に巡回させる

DFT を用いた畳み込みの式を再掲します。 $$ \begin{aligned} F_{\omega}[a]_l \cdot F_{\omega}[b]_l &= \left(\sum_{j=0}^{n-1} a_j \omega^{jl}\right) \left(\sum_{j=0}^{n-1} b_j \omega^{jl}\right) \\ &= \sum_{k=0}^{n-1} \sum_{j=0}^k a_j b_{k-j} \omega^{kl} + \sum_{k=0}^{n-2} \sum_{j=k+1}^{n-1} a_j b_{k+n-j} \omega^{(k+n)l} \\ &= \sum_{k=0}^{n-1} \left( \sum_{j=0}^{n-1} a_j b_{(k-j)\bmod n} \right) \omega^{kl}. \end{aligned} $$ ここで、$\sum a_j b_{k-j}$ と $\sum a_j b_{k+n-j}$ が一つの項にまとめられている部分が不都合なのでした。これを $(\sum a_j b_{k-j}) + (\sum a_j b_{k+n-j})\,i$ のようにしてしまおうという試みです。工夫①とは別の虚部の使い方です。

$\psi$ を $i$ の原始 $n$ 乗根とします。$\omega$ の原始 $4$ 乗根という見方もできるかもしれません。工夫③で $\xi$ を使った部分に $\psi$ を使ってみます。 $$ \begin{aligned} &\phantom{{}={}} \left(\sum_{j=0}^{n-1} a_j (\psi\omega^l)^j\right) \left(\sum_{j=0}^{n-1} b_j (\omega^l)^j\right) \\ &= \sum_{k=0}^{n-1} \sum_{j=0}^k a_j b_{k-j} \psi^k\omega^{kl} + \sum_{k=0}^{n-2} \sum_{j=k+1}^{n-1} a_j b_{k+n-j} \psi^{k+n}\omega^{(k+n)l} \\ &= \sum_{k=0}^{n-1} \sum_{j=0}^k a_j b_{k-j} \psi^k\omega^{kl} + i\sum_{k=0}^{n-2} \sum_{j=k+1}^{n-1} a_j b_{k+n-j} \psi^k\omega^{kl} \\ &= \sum_{k=0}^{n-1} \left(\sum_{j=0}^k a_j b_{k-j} + i \sum_{j=k+1}^{n-1} a_j b_{k+n-j}\right) \psi^k\omega^{kl}. \end{aligned} $$ 工夫③同様にして $\psi^k$ の項を後から消せばよく、長さ $n$ の列に対しての DFT で済むことになりました。一方で工夫①と併用できないため、工夫②が必須である前提だと厳しいかもしれません。$1$ の原始 $4n$ 乗根が取れる $p$ に対する NTT であれば使い道がある気もします。

実装例

adamant の実装を載せます。いろいろ書いてあります。実はちゃんとわかっていません。

judge.yosupo.jp

誤差評価

三角関数について

一般に、table-maker’s dilemma により、三角関数の値を correct rounding で求めるのは簡単ではありません。すなわち、四則演算のように $0.5$ ULP で押さえられることを前提とするのは厳しいです。

rsk0315.hatenablog.com

また、計算したい $\cos(\theta)$ や $\sin(\theta)$ の引数 $\theta$ も無理数であり、これを求める段階でも誤差が生じます。 FFT したい列の長さを高々 $2^{20}$ 程度として、FFT のアルゴリズムを考慮すると、各整数 $0\le k\lt 2^{19}$ に対しての $\theta = k\pi/2^{19}$ のみ考えれば十分です。

記法としては、次のように書き分けます。

$\gdef\libmsin#1{\operatorname{\texttt{s\hspace{-0.04em}i\hspace{-0.04em}n}}(#1)}$ $\gdef\libmcos#1{\operatorname{\texttt{c\hspace{-0.04em}o\hspace{-0.04em}s}}(#1)}$ $\gdef\libmexp#1{\operatorname{\texttt{exp}}(#1)}$

$\theta$, $\cos(\theta)$, $\sin(\theta)$：各値の真の値
$\hat\theta$：$\theta$ を浮動小数点数で素朴に計算したもの
- $\hat\theta = 2\otimes k\otimes \roundcirc{\pi} \oslash 2^{20} = (k\otimes\roundcirc{\pi})\cdot 2^{-19}$
- 一般に $\hat\theta = \roundcirc{\theta}$ とは限らず、たとえば $k = 11$ が反例となる
$\roundcirc{\cos(\hat\theta)}$, $\roundcirc{\sin(\hat\theta)}$：$\hat\theta$ に対する三角関数の値を、該当の浮動小数点数型に正確に丸めたもの
- いわゆる「パソコンでふつうに計算したときに得られる値」であってほしい値
- table-maker’s dilemma により、実際に得るのはむずかしい
$\libmcos{\hat\theta}$, $\libmsin{\hat\theta}$：$\hat\theta$ に対する三角関数の値を数値計算ライブラリで計算したもの
- いわゆる「パソコンでふつうに計算したときに得られる値」の実際の値
- 用いる数値計算ライブラリに応じて結果が変わる

$\cos(\theta)$ と $\libmcos{\hat\theta}$ の誤差が実際の誤差評価で必要になるもので、$\roundcirc{\cos(\hat\theta)}$ と $\libmcos{\hat\theta}$ については「correctly rounded な値が得られないケースもあるんだねえ」というのを確かめる程度の気持ちで眺めます。

AtCoder の環境で実験をしました。libm-2.36 を使っていそうでした。他のジャッジ環境については知りません。

実験に関して

AtCoder 環境のコードテストでは 65535 bytes の出力を見られるため、欲しい三角関数の値たちを出力するのを繰り返しました（時間がかかりました）。別途それらの値たちの SHA-256 のハッシュ値を計算し、誤りがないことを確かめました。

rsk0315.hatenablog.com

その過程でこの記事を書き、手元環境で出力を再現できるようになったため、時間をかけたパートはほぼ無駄になりました。 __sin_fma や __cos_fma を使っていると思います。

実験したときは勘違いしていて $2^{20}$ まで試していました。

#include <math.h>
#include <stdio.h>

int main(void) {
  double pi = 0x1.921fb54442d18p+1;
  int n = 1 << 20;
  for (int k = 0; k < n; ++k) {
    double theta = 2 * k * pi / n;
    printf("%.13a\n", sin(theta));
  }
}

これの出力の SHA-256 のハッシュ値は下記でした。

1d20e0a9a38224f0bc155f370d83e00ddccb87fa95eb5853a6be53d4f2f22285

sin(theta) の部分を cos(theta) の部分にした際のハッシュ値は下記でした。

46718d7f62d74613e50b39b90baec3bd01d2388057d6f3ab86b9c93a4c298cab

おわりです。$\eod$

真の値に十分近い値は Sollya で計算しています。

$|{\cos(\theta) - \libmcos{\hat\theta}}|$ が最大となるのは $k = 350297$ のときで、$3.142\times 10^{-16}$ や $2.830\times 2^{-53}$ で上から押さえられる程度でした。また、$|{\sin(\theta) - \libmsin{\hat\theta}}|$ が最大となるのは $k = 501943$ のときで、$3.478\times 10^{-16}$ や $3.133\times 2^{-53}$ で上から押さえられる程度でした。

$\roundcirc{\cos(\hat\theta)}\ne\libmcos{\hat\theta}$ となる $k$ は $696$ 個あり、$k$ が最小なのは $k = 8434$ のときでした。誤差が最大なのは $k = 283410$ のときで $0.514$ ULP で押さえられる程度でした。また、$\roundcirc{\sin(\hat\theta)}\ne\libmsin{\hat\theta}$ となる $k$ は $722$ 個あり、$k$ が最小なのは $k = 13142$ のときでした。誤差が最大なのは $k = 496456$ のときで $0.515$ ULP で押さえられる程度でした。

また、点としての誤差は $k = 433455$ のとき最大で、各 $0\le k\lt 2^{19}$ について下記が成り立ちます。 $$ |(\libmcos{\hat\theta} + i\libmsin{\hat\theta}) - (\cos(\theta) + i\sin(\theta))| \lt 3.58277\times 2^{-53}. $$

複素数の四則演算

除算は今回使わない*5ので、三則演算と言うべきかもしれません？

実部と虚部を f64 で持っている複素数型 $\Complex_F$ を考え、それにおける +, -, * をそれぞれ $\oplus$, $\ominus$, $\otimes$ で書くことにします。また、$\alpha, \beta\in\Complex_F$ とします。 $$ \begin{aligned} |(\alpha\oplus\beta)-(\alpha+\beta)| &\le |\alpha+\beta| \cdot 2^{-53}, \\ |(\alpha\ominus\beta)-(\alpha-\beta)| &\le |\alpha-\beta| \cdot 2^{-53}, \\ |(\alpha\otimes\beta)-(\alpha\times\beta)| &\le |\alpha\times\beta| \cdot \sqrt8\cdot 2^{-53} \\ \end{aligned} $$ が成り立ちます。$\otimes$ の評価における係数 $\sqrt8\cdot 2^{-53}$ は、解析をがんばることで $\sqrt5\cdot 2^{-53}$ に落とせるらしいです。また、計算の際に FMA を用いることで $2\cdot 2^{-53}$ になるらしいです。詳細はあまり追っていません。

FFT のワンステップ

あるペア $(\alpha, \beta)$ と $1$ の冪根 $\omega$ に近い値 $\hat\omega$ に対して、$\alpha\oplus(\beta\otimes\hat\omega)$ と $\alpha\ominus(\beta\otimes\hat\omega)$ を計算することになります。各計算にかかる誤差項を $\delta_{\omega}$, $\delta_{\otimes}$, $\delta_{\oplus}$, $\delta_{\ominus}$ として $$ \begin{aligned} \alpha\oplus(\beta\otimes\hat\omega) &= (\alpha+(\beta\omega(1+\delta_{\omega}) )(1+\delta_{\otimes}) )(1+\delta_{\oplus}), \\ \alpha\ominus(\beta\otimes\hat\omega) &= (\alpha-(\beta\omega(1+\delta_{\omega}) )(1+\delta_{\otimes}) )(1+\delta_{\ominus}) \end{aligned} $$ と書けます。$|\omega| = 1$ に注意して $$ \begin{aligned} &\phantom{{}={}} | (\alpha\oplus (\beta\otimes\hat\omega) ) - (\alpha+\beta\omega) | \\ &= |\alpha\delta_{\oplus} + \beta\omega( (1+\delta_{\omega})(1+\delta_{\otimes})(1+\delta_{\oplus})-1)| \\ &\le |\alpha| |\delta_{\oplus}| + |\beta| |(1+\delta_{\omega})(1+\delta_{\otimes})(1+\delta_{\oplus})-1| \\ &\lt |\alpha| \cdot 2^{-53} + |\beta| |(1+3.58277\cdot 2^{-53})(1+\sqrt5\cdot 2^{-53})(1+2^{-53})-1| \\ &\lt 6.8184\cdot 2^{-53}\cdot (|\alpha| + |\beta|) \end{aligned} $$ となります。

また、任意の $z, w\in\Complex$ に対して $|z+w|^2\le 2(|z|^2+|w|^2)$ より $$ \begin{aligned} &\phantom{{}={}} | (\alpha\oplus (\beta\otimes\hat\omega) ) - (\alpha+\beta\omega) |^2 \\ &= |\alpha\delta_{\oplus} + \beta\omega( (1+\delta_{\omega})(1+\delta_{\otimes})(1+\delta_{\oplus})-1)|^2 \\ &\le 2\left( (|\alpha| |\delta_{\oplus}|)^2 + (|\beta| |(1+\delta_{\omega})(1+\delta_{\otimes})(1+\delta_{\oplus})-1|)^2 \right) \\ &\lt 2\left( (|\alpha| \cdot 2^{-53})^2 + (|\beta| |(1+3.58277\cdot 2^{-53})(1+\sqrt5\cdot 2^{-53})(1+2^{-53})-1|)^2 \right) \\ &\lt 2\cdot (6.8184\cdot 2^{-53})^2 \cdot(|\alpha|^2 + |\beta|^2) % &\lt (9.6433\cdot 2^{-53})^2 \cdot(|\alpha|^2 + |\beta|^2) \end{aligned} $$ となります。また、 $$ \begin{aligned} &\phantom{{}={}} |\alpha+\beta\omega|^2 + |\alpha-\beta\omega|^2 \\ &= (\alpha+\beta\omega)(\alpha+\beta\omega)^{\ast} + (\alpha-\beta\omega)(\alpha-\beta\omega)^{\ast} \\ &= (\alpha+\beta\omega)(\alpha^{\ast}+(\beta\omega)^{\ast}) + (\alpha-\beta\omega)(\alpha^{\ast}-(\beta\omega)^{\ast}) \\ &= (|\alpha|^2 + \alpha(\beta\omega)^{\ast} + \alpha^{\ast}\beta\omega + |\beta\omega|) + (|\alpha|^2 - \alpha(\beta\omega)^{\ast} - \alpha^{\ast}\beta\omega + |\beta\omega|) \\ &= 2(|\alpha|^2 + |\beta|^2) \end{aligned} $$ より、FFT のワンステップによりノルム（$2$-ノルム）は $\sqrt2$ 倍になることがわかります。

これらより、誤差のベクトルのノルムは、真の値のノルムの $6.8184\cdot 2^{-53}$ 倍で押さえられることがわかります。

note: 誤差のノルムが十分小さいんだから個々の要素の誤差も小さいだろ〜というアプローチをしようとしましたが、要素ごとの積を求めた後のノルムってどんなことになるんだ？そもそも別に十分小さくなくない？などになり、参っています。

畳み込み全体

しんどいよ〜

$|\alpha+\beta\omega|\le|\alpha|+|\beta|$ に注意します。 $\alpha' \gets \alpha\oplus(\beta\otimes\omega)$ の更新をすると $$ |\alpha'| \le (1 + 6.8184\cdot 2^{-53})(|\alpha|+|\beta|) $$ となります。工夫②での $a_k^{(0)}$, $a_k^{(1)}$ の上界に注意し、 $$ \begin{aligned} z_0 &= 15811\sqrt2, \\ z_{k+1} &= 2(1+6.8184\cdot 2^{-53})\cdot z_k \end{aligned} $$ で定義される $z_k$ は、更新を $k$ 回行った後の $|\alpha|$, $|\beta|$ の上界となります。また、$k$ 回目の更新で発生する誤差の上界は $z_k - 2^k z_0$ とできます。

よって、$20$ 回の更新で累積する誤差は $$ \begin{aligned} \sum_{k=1}^{20} z_k - 2^k z_0 &= \sum_{k=1}^{20} {( (2(1+6.8184\cdot 2^{-53}) )^k - 2^k) z_0} \\ % &= z_0 \sum_{k=1}^{20} {\left( (2(1+6.8184\cdot 2^{-53}) )^k - 2^k\right)} \\ % &= z_0 \left( 2(1+6.8184\cdot 2^{-53})\frac{1-(2(1 + 6.8184\cdot 2^{-53}) )^{20}}{1-2(1+6.8184\cdot 2^{-53})} - 2^{21}+2 \right) % &\lt 3.0164\cdot 10^{-8} z_0 \\ % &\lt 6.7446\cdot 10^{-4} &\lt 1.3813 \cdot 2^{-11} \end{aligned} $$ です。よって、誤差項 $|\delta|\lt 1.3813\cdot 2^{-11}$ に対して $\hat F_{\omega}[a]_l = F_{\omega}[a]_l + \delta$ が得られることになります。cancellation がたくさん起きるはずなのと、最終的には誤差が $0.5$ 未満であることを示したいので、相対誤差ではなく絶対誤差で評価しようとしています。

記号をよしなに導入しつつ、$F_{\omega}[a]_l\le 2^n z_0$ で評価することにして $$ \begin{aligned} &\phantom{{}={}} |(F_{\omega}[a]_l + \delta_a)(F_{\omega}[b]_l + \delta_b)(1+\delta_{\otimes}) - F_{\omega}[a]_l F_{\omega}[b]_l| \\ &\le |F_{\omega}[a]_l| |F_{\omega}[b]_l| |\delta_{\otimes}| + (|F_{\omega}[a]_l| |\delta_b| + |F_{\omega}[b]_l| |\delta_a| + |\delta_a| |\delta_b|)(1+|\delta_{\otimes}|) \\ &\lt 2\cdot 2^{20} z_0 \cdot 1.3813\cdot 2^{-11} + 1.9080\cdot 2^{-22} \\ &\lt 1.8888 \cdot 2^{24} \end{aligned} $$ を得ます。これを $n$ で割ったものに対して逆変換をすることになりますが、逆変換をする前の段階で、誤差の上界が $1.8888\cdot 2^4\approx 30.221$ くらいになっています。困りました。

道は二つあって、「かしこい評価をがんばって、実は逆変換をした後も誤差が $0.5$ より小さいことを示せる」か「実は反例があるのでがんばって探す」かです。

証明ができないので、反例見つかれーッって言いながらランダムケースを投げまくってる
ダメ人間の末路
— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2025年6月4日

twitter.com

その他

$\sum_{k=0}^{n-1} \sin(\frac{k\pi}n) \lt \frac{2n}{\pi}$ が成り立ちます。区分求積法などで示せます。評価に役立つかと思ったのですが気のせいでした。

反例

簡単のため、乱数の部分を固定します。乱数がリークできたら落とせるというような感じです。

工夫②の出力が $z = 15811-15811i$ となるような整数を構成して、それだけからなる列を考えたり、乱数で少しずつ散らしたりすると、反例がそれなりに見つかりました。

adamant の実装や、自前のもう少し素朴な実装についていくつか試しましたが、工夫①と②による実装では $(\mathbb F_{1000000007})^{2^{19}}$ の列の畳み込みには耐えなさそうでした。

展望など

double-double を使って $2^{-100}$ 程度の精度で行うことでなんとかする
- 工夫①を使って DFT $2$ 回で済ませる
  - 実測がめちゃくちゃ遅かった（それはそうか）
値を 3 つに分ける
- 工夫②の要領で $d_1$, $d_2$ を使って分ける
  - 計 $8$ 回の FFT を行うことになる
  - double-double 解法ほどは遅くなかった
FFT の各ステップを radix-2 ではなく radix-4 や radix-8 などにする
- 改善はするはずだが反例はあるはず
  - 実際、adamant の解法は radix-4 のはず

double-double の解法においては、$32$ 以下の各非負整数 $k$ に対して $\sin(\frac{k\pi}{2^{19}})$, $\sin(\frac{k\pi}{2^{13}})$, $\sin(\frac{k\pi}{2^{7}})$, $\cos(\frac{k\pi}{2^{19}})$, $\cos(\frac{k\pi}{2^{13}})$, $\cos(\frac{k\pi}{2^{7}})$ の近似値を計算しておき、加法定理などを用いて復元しました。

3 つに分ける解法においては、$(d_1, d_2) = (10^3, 10^6)$ として $a_k = a_k^{(0)} + a_k^{(1)}d_1 + a_k^{(2)}d_2$ とすることで、各値の絶対値を $500$ 以下にできます。上記と同様の評価をして、$20$ 回の順変換で $1.3978\cdot 2^{-16}$、逆変換の前までで $1.935\cdot 2^{-6}$ 程度に収まりました。逆変換をした後の最終的な誤差の上界は $1.9311\cdot 2^{13}$ とかになってしまいました。さすがにこれは、もっとかしこくやればどうにかなるような気もします。

素朴にやると「初期値の上界が $2^{20}\cdot (500\sqrt 2)^2\approx 1.908\cdot 2^{38}$ の列に対しての FFT」と同程度の誤差が出ることになるので、逆変換であることを使って評価しないといけないような気がしています。

下記の問題は長さが高々 $2^{18}$、値が高々 $100$ なので、ざっくり似た範囲な気がします。この問題の FFT 解法の正当性がどう示されているのかは知らないですし、$(2^{18}, 100)$ と $(2^{20}, 500)$ で反例の有無が異なるかも知りません。ざっくり評価だと $100$ 倍くらい違う気がするので、やっぱり似ていると見なせないかもしれません。

atcoder.jp

あるいは、乱数 $\xi$ に対して $a(\xi x) b(\xi x)$ を求めるのではなく、各係数を小さくするような $u, v\in \mathbb F_p$ をうまいこと持ってきて $(u a(x) )(v b(x) )\cdot (uv)^{-1}$ を計算するみたいなアプローチはどうでしょう。あまり考えていません。

あとがき

定数倍をよくするために、誤り率を $\varepsilon\gt 0$ にしようというアプローチとあまり相容れません。「オーダー改善の文脈なら許せるのか？」「許せないのなら SHA-2 を使う資格がないし、よくある署名の方法にも納得しないべきでは？」などの声もあります。茂野吾郎と申します。本文中で触れたように、確率パートが入力にあるかアルゴリズムにあるか、というような違いはある気もします。ふわふわした立場で話している気もします。

応用の文脈においては、「$\varepsilon\gt 0$ だから嫌い！」とするのではなく、「その応用先に対して、$\varepsilon\gt 0$ でも許容できるか？」を考えて判断するのが妥当だと思います。それはそれとして、一般論での嗜好の問題として「速度よりも正当性が先」という価値観がある気がします。もちろん、逆に「確率 $0.1^{10^{100}}$ での誤りを許容しないがために時間を $100$ 倍かけます」のような態度を理解できないであろう人もいる気はします。

競プロで想定解を用意するという文脈においては、そういう誤りをあまり許容できない側な気がします。正しい値のファイルを用意するというだけなら 2 sec より多い時間をかけられる気もしますし、それはフェアではないような気もします。

最初から「今から誤り率 $\varepsilon\gt 0$ の話をしますよ〜」という前置きがあって話されるのならいいですが、さも「こういう問題を（正確に）解くアルゴリズムの話をしますよ〜」のような顔をして「さて！ここでは誤答も許すことにしちゃいます！」みたいなことを言われると萎えるというのもありますね。クエリ定数時間の簡潔データ構造の話だと思ったらクエリ $\Theta(\log(n) )$ 時間のコンパクトデータ構造をお出しされるのと似ています*6。簡潔データ構造が実際にうれしいかはさておき、期待の裏切られ感で嫌な気持ちになるというやつです。

誤差評価に関しては、たいへんだなぁという気持ちになりました。ランダムケースで数値計算をして「毎回同じ方向に上界と同程度の大きさの誤差が出ることはほぼないだろうから、なんかいい感じになるだろ〜」みたいなノリで解析されるとやっぱりびっくりしちゃいます。各状況での最悪ケースの積み重ねで悲観的すぎる評価をするのも悪手な気もしつつ、どうしたらいいのかなという感じです。

各 $\mathbb F_{u 2^e+1}$ で何度も FFT をするのが嫌だな〜という話だった気がするのに、$\Complex_F$ でやるにしても精度の問題で結局何回も FFT をやらされる羽目になる上、求めた値の正当性もちゃんとわからず、ん〜〜という気持ちです。どうしたものですかね。

気が向いたらサーベイをちゃんとやりつつ、気が向かなかったら別のトピックのお勉強をしたり気ままに遊んだりしようかなという気持ちです。しばらくは浮動小数点数はおなかいっぱいかもという気がしていますが、明日には浮動小数点数関連の別のトピックのお勉強を始めているかもしれませんし、なんとも言えません。

おわり

おわりです。お粗末さまでした。

*1:元ネタ。同じ元ネタのパロディの宣伝。

*2:あるいは、係数体を $\mathbb F_{998244353}$ とする多項式を考えている場合？

*3:catastrophic cancellation 的な意味で、$(u+v)(u-v)$ として計算するのがいいか $u^2-v^2$ として計算するのがいいかは場合によります。ここでは、相対誤差ではなく絶対誤差が重要になるはずで、どちらでもいいような気がしています。たぶん。

*4:一般の整数を法とした畳み込みでも使える手法ではありますが、あまり区別しないで話すことにします。

*5:$i$ での形式的な除算などをしようとしている際は浮動小数点演算を介さないのでノーカンです。

*6:無限回擦っているが、競プロ界隈で簡潔データ構造の話をしているのは老人しかいないため、若い人には「無限回シリーズ」だと認識されていなさそう。

背景

本題