# 再現できる数値実験：乱数種・環境・設定の記録と決定性 > 数値実験をコード・入力・設定・乱数種・環境からの写像とみて、再現の 3 つの水準を区別する。和の順序による差の上界 2γΣ|x|、内容ハッシュの衝突確率、種による決定性、種を変えた平均の標準誤差を証明し、実行記録の書き方を示す。 > https://rikai.mugen-giken.com/computer-science/research-computing/reproducible-experiments ## 0. この記事の要点 同じプログラムを同じように動かしたのに、あとで数値が少し違って出る、ということがあります。この記事は、その違いが「気にしなくてよい誤差」なのか「どこかの間違い」なのかを見分ける考え方と、あとから同じ結果を作り直すための記録の取り方を扱います。読み終えると、結果に何を添えて残せばよいか、末尾の桁がずれたときにどう判断するかを、理由つきで説明できるようになります。 - 数値実験は「コード・入力・設定・乱数種（乱数を作り始める最初の値）・環境」の組から出力への写像です。この 5 つを記録していない結果は、あとから同じものを作れません。 - 再現には 3 つの水準があります。**ビット単位**の再現、**許容誤差つき**の再現、**統計的**な再現です。どの水準を約束するかを先に決めます。 - 浮動小数の加法は結合的ではありません。ただし、$n$ 個の数をどんな順序・どんな括弧づけで足しても、2 つの結果の差は $2\gamma_{n-1}\sum|x_i|$ 以下です（$\gamma_k=ku/(1-ku)$、$u$ は 1 回の丸めで生じる相対誤差の上限）。順序を変えて結果が変わったとき、差がこの範囲なら丸めの差、超えていれば別の誤りです。 - 擬似乱数生成器（決まった規則で乱数らしい数を順に出す仕組み）は状態機械なので、乱数を消費する回数と順序が同じなら、結果は種の関数です。1 つの生成器を並列に取り合うと、種を固定しても再現しません。実験ごとに生成器を分けます。 - $K$ 個の独立な種の平均の標準誤差は $s/\sqrt{K}$ で推定できますが、種が 1 個では推定できません。種の集合は実行前に決め、都合のよい種だけを選びません。 - 出力には必ず**実行記録**（設定・種・版・入力と出力の SHA-256。内容から決まる固定長の指紋のような値）を添えます。$b$ ビットのハッシュで $N$ 個の入力が衝突する確率は $N(N-1)/2^{b+1}$ 以下です。 ## 1. 動機 半年前に自分が出した図を、査読者の求めで作り直すことになったとします。スクリプトは残っています。ところが実行すると、数値が小数第 6 位から合いません。考えられる原因はいくつもあります。乱数の種を記録していなかった、NumPy を上げた、設定の一部をコードに直接書いていて書き換えた、計算機が変わってスレッド数が違う、入力ファイルを誰かが更新した、などです。 困るのは、「合わない」ことそのものより、**その差が無害かどうかを判定できない**ことです。末尾の数桁の違いは、和の順序が変わっただけかもしれませんし、入力が 1 行ずれているのかもしれません。前者は放置してよく、後者は結論を変えます。 この記事は、この判定に根拠を与えます。作法（種を記録する、版を記録する、ハッシュを添える）を並べるだけでなく、それぞれの作法が**何を保証し、何を保証しないか**を命題として述べ、証明します。浮動小数点数の誤差モデルは[浮動小数点数と丸め誤差](/mathematics/numerical-analysis/floating-point-arithmetic)の章のものを使い、ここでは二重に定義しません。 ## 2. 実験を写像とみる ### 2.1. 何が結果を決めるか **数値実験**とは、次の 5 つの組 $(c,x,\theta,\sigma,\eta)$ に出力 $y$ を対応させる手続きである。 - $c$: **コード**（実行するプログラムの全体。版管理の 1 つのコミットで特定する） - $x$: **入力**（データファイルの内容） - $\theta$: **設定**（刻み幅・反復回数・モデルの係数など、実行時に選ぶ値） - $\sigma$: **乱数種** - $\eta$: **環境**（言語処理系と数値ライブラリの版、線形代数ライブラリ BLAS の実装、スレッド数、CPU の種類） 出力が組 $(c,x,\theta,\sigma,\eta)$ だけで決まる、すなわち写像 $R:(c,x,\theta,\sigma,\eta)\mapsto y$ が定まるとき、実験は**決定的**であるという。同じ組に対して実行ごとに出力が変わりうるとき、**非決定的**であるという。 「同じスクリプトを実行したのに結果が違う」のは、5 つのどれかが変わったか、実験が非決定的かのどちらかです。非決定性の源は限られています。 | 源 | 何が起きるか | 対処 | |---|---|---| | スレッドごとの部分和 | 和をとる順序がスレッド数や割り当てで変わる（3 節） | 分割を固定する。差を 3 節の上界と比べる | | 集合・辞書の走査順 | Python の `set` の走査順は文字列のハッシュの乱数化で実行ごとに変わりうる（`dict` は Python 3.7 以降は挿入した順に走査するが、挿入の順序が `set` の走査などで決まっていれば同じことが起きる） | 走査の前に整列する | | 時刻・プロセス番号から作った種 | 実行ごとに種が変わる | 種を設定に書く（5 節） | | 1 つの生成器の取り合い | 乱数を受け取る順序が実行ごとに変わる（5 節） | 実験ごとに生成器を分ける | | GPU の一部の演算 | 同じ入力でも実行ごとに加算の順序が変わる実装がある | 決定的な実装を選ぶか、許容誤差つきで比べる | ### 2.2. 再現の 3 つの水準 2 回の実行の出力を $y$, $y'$ とする（実数の配列とし、成分を $y_j$, $y'_j$ と書く）。 1. **ビット単位の再現**: $y$ と $y'$ がバイト列として一致する。 2. **許容誤差つきの再現**: あらかじめ決めた $\tau_{\mathrm{abs}},\tau_{\mathrm{rel}}\ge0$ について、すべての $j$ で $|y_j-y'_j|\le\tau_{\mathrm{abs}}+\tau_{\mathrm{rel}}|y_j|$ が成り立つ。 3. **統計的な再現**: 出力が乱数に依存する量の推定値であるとき、別の種の集合で求めた推定値の信頼区間が、元の信頼区間と重なる。 水準 1 は、同じ環境で同じ組を与えたときに求めるものです。環境 $\eta$ が違えば、同じコードでも和の順序や初等関数の末尾のビットが変わるので、水準 1 は期待できません。そのときに求めるのが水準 2 で、許容誤差 $\tau$ の根拠を与えるのが 3 節の定理です。水準 3 は、乱数を使う実験の結論が「たまたまその種だったから」ではないことを確かめるためのもので、6 節で扱います。 ### 2.3. 実行記録 1 回の実行の**実行記録**（manifest）とは、次の項目を 1 つのファイルにまとめたものである。 - 設定 $\theta$ の全体と乱数種 $\sigma$ - コードの版（コミットの識別子） - 環境 $\eta$（処理系とライブラリの版、機種） - 入力 $x$ の内容ハッシュと、出力 $y$ の内容ハッシュ - 開始と終了の時刻 入力と出力は大きいので、内容そのものではなく**内容ハッシュ**（ の暗号学的ハッシュ関数 $H$ の値）を記録します。Git がコミットをハッシュで同定するのと同じ考え方で（）、ハッシュが一致すれば内容が一致するとみなします。この「みなす」がどれだけ安全かは 4 節で評価します。環境を丸ごと固定したいときは、[コンテナ](/computer-science/software-engineering/containers-and-kubernetes#def-container)（アプリを必要なライブラリごと隔離して動かす仕組み）のイメージ（）の識別子を環境の欄に書きます。イメージも内容で名前が付くので、同じ識別子なら同じライブラリの組です。 ## 3. 和の順序はどこまで結果を変えるか ### 3.1. 加法は結合的でない 以下、$\mathbb{F}$ は精度 $p$ の浮動小数点数の集合（）、$u=2^{-p}$ は単位丸め、$\mathrm{fl}$ は最近接偶数丸め（）です。binary64 では $p=53$ です。 $x=M\cdot2^{e-p+1}\in\mathbb{F}$ を正規化数（$2^{p-1}\le M\lt2^p$）で、$M$ が偶数、$e_{\min}+1\le e\le e_{\max}$ とし、$h=2^{e-p}$（$x$ の右隣との間隔の半分）とおく。このとき $h\in\mathbb{F}$ で、 $$ \mathrm{fl}\bigl(\mathrm{fl}(x+h)+h\bigr)=x,\qquad\mathrm{fl}\bigl(x+\mathrm{fl}(h+h)\bigr)=x+2h $$ である。2 つの括弧づけの結果は、$x$ の位置での浮動小数点数の間隔 $2h$（1 ulp）だけ異なる。とくに binary64 で $x=1$, $h=2^{-53}$ とすると、$(1+2^{-53})+2^{-53}$ は $1$、$1+(2^{-53}+2^{-53})$ は $1+2^{-52}$ と計算される。 ここで示したいのは「括弧の位置を変えるだけで結果が 1 ulp ずれる例が、いつでも作れる」ことです。証明は、小さい数 $h$ が $x$ の右隣までの距離のちょうど半分であることを使い、丸めが偶数側に倒れる向きを追います。指数の範囲の場合分けは、$h$ や $2h$ が浮動小数点数として表せること（極端に小さい数を含めて）を確かめるためのもので、本筋は「左の括弧づけでは $x$ に戻り、右の括弧づけでは $x+2h$ に届く」という一点です。 **$h$ と $2h$ が $\mathbb{F}$ に属すること.** $t\in\{e-p,\ e-p+1\}$ について $2^t\in\mathbb{F}$ を示します。$t\ge e_{\min}$ なら、$2^t=2^{p-1}\cdot2^{t-p+1}$ は指数 $t$（$e_{\min}\le t\le e\le e_{\max}$）の正規化数です。$t\lt e_{\min}$ なら、$2^t=2^{t-e_{\min}+p-1}\cdot2^{e_{\min}-p+1}$ と書けます。整数 $2^{t-e_{\min}+p-1}$ は、$t\ge e-p\ge e_{\min}+1-p$ から指数が $0$ 以上なので $1$ 以上、$t\lt e_{\min}$ から $2^{p-1}$ 未満なので、これは非正規化数です。 **$x$ の右隣.** $M$ は偶数で $M\lt2^p$ なので $M\le2^p-2$、したがって $M+1\le2^p-1$ で、$x+2h=(M+1)\cdot2^{e-p+1}$ は指数 $e$ の正規化数です。$\mathbb{F}\cap[2^e,2^{e+1}]$ は間隔 $2^{e-p+1}=2h$ の等間隔な点列（）なので、$x$ と $x+2h$ の間に $\mathbb{F}$ の元はありません。 **左の括弧づけ.** 実数 $x+h$ は $x$ と $x+2h$ のちょうど中点なので、$\nabla(x+h)=x$, $\Delta(x+h)=x+2h$ で、両者は $x+h$ から等距離です。最近接偶数丸めは仮数 $M$ が偶数の方を選びます。$x$ の仮数は $M$（偶数）、$x+2h$ の仮数は $M+1$（奇数）なので、$\mathrm{fl}(x+h)=x$ です。もう一度 $h$ を足しても同じ計算で、$\mathrm{fl}(\mathrm{fl}(x+h)+h)=\mathrm{fl}(x+h)=x$ です。 **右の括弧づけ.** $h+h=2h\in\mathbb{F}$ なので $\mathrm{fl}(h+h)=2h$（$\mathbb{F}$ の元は丸めで動きません）。$x+2h\in\mathbb{F}$ なので $\mathrm{fl}(x+2h)=x+2h$ です。 binary64 の $x=1$ は $M=2^{52}$（偶数）, $e=0$ で、$h=2^{-53}$ です。
x（仮数は偶数） x+h（真ん中） x+2h（仮数は奇数） 間隔 2h fl(x+h)=x
この反例の仕組みは、「小さい数を 1 つずつ足すと毎回切り捨てられるが、先にまとめてから足すと届く」というものです。交換法則 $\mathrm{fl}(a+b)=\mathrm{fl}(b+a)$ は成り立つ（厳密な和 $a+b$ が同じなので）ので、和の結果を変えるのは**括弧づけ**と、どの数をどの位置に置くかです。 ### 3.2. どんな順序でも成り立つ上界 「任意の括弧づけ・任意の順序」を木で表します。 $x_1,\ldots,x_n\in\mathbb{F}$（$n\ge1$）の**和の木**とは、葉がちょうど $n$ 枚あり、葉以外の頂点（**内部頂点**）がちょうど 2 つの子をもつ根つき木で、葉に $x_1,\ldots,x_n$ を 1 つずつ割り当てたものである。各頂点の**値**を、葉では割り当てた $x_i$、内部頂点では 2 つの子の値 $a,b$ から $\mathrm{fl}(a+b)$ と定め、根の値を木の**計算結果**という。葉 $i$ の**深さ** $d_i$ を、根から葉 $i$ までの道の上にある内部頂点の個数とし、$\max_id_i$ を木の**高さ**という。 前から順に足す逐次和は、高さ $n-1$ の一本道の木です。隣どうしを足していく対ごとの和は、高さ $\lceil\log_2n\rceil$ の釣り合った木です。スレッドごとに部分和をとってから合わせる計算も、SIMD 命令（1 つの命令で複数の数を同時に処理する命令）で 8 本の部分和を並行にとる計算も、それぞれ 1 つの和の木です。
+++ +++ x₁x₂x₃x₄ x₁x₂x₃x₄ 深さ 3, 3, 2, 1 深さ 2, 2, 2, 2
$x_1,\ldots,x_n\in\mathbb{F}$ の和の木の高さを $h$ とし、$hu\lt1$ とする。すべての内部頂点で、2 つの子の値の厳密な和が $0$ であるか、絶対値が $x_{\min}$ 以上 $x_{\max}$ 以下である（オーバーフローもアンダーフローも起きない）とする。このとき計算結果 $\hat s$ は $$ \hat s=\sum_{i=1}^nx_i(1+\theta_i),\qquad|\theta_i|\le\gamma_{d_i}\le\gamma_h $$ を満たす（$\gamma_k=ku/(1-ku)$、$\gamma_0=0$）。また $h\le n-1$ である。 この補題は、どんな順序・括弧づけの和でも、各項 $x_i$ が受ける丸めの回数は「根から葉までに通る加算の個数」だけで、その分だけ $(1+\delta)$ の因子が掛かる、と言っています。証明は、木の下から上へ、部分木ごとに「各項に何個の因子が掛かっているか」を数える帰納法です。 頂点 $v$ を根とする部分木の葉の集合を $\Lambda(v)$、葉 $i\in\Lambda(v)$ から $v$ までの道の上にある内部頂点の個数を $d_i(v)$ と書きます。次を、部分木の葉の枚数についての帰納法で示します。 > 頂点 $v$ の値は $\sum_{i\in\Lambda(v)}x_i\,\pi_i(v)$ と書ける。ここで $\pi_i(v)$ は $d_i(v)$ 個の因子 $(1+\delta)$, $|\delta|\le u$ の積である。 $v$ が葉 $i$ なら値は $x_i$ で、$d_i(v)=0$、$\pi_i(v)=1$（空の積）です。$v$ が内部頂点で、子を $v_1,v_2$、その値を $a,b$ とします。仮定から四則の標準モデル（）が使えて、$v$ の値は $(a+b)(1+\delta_v)$, $|\delta_v|\le u$ です。帰納法の仮定を $a$ と $b$ に使うと、$v$ の値は $$ \sum_{i\in\Lambda(v_1)}x_i\,\pi_i(v_1)(1+\delta_v)+\sum_{i\in\Lambda(v_2)}x_i\,\pi_i(v_2)(1+\delta_v) $$ です。$\Lambda(v)$ は $\Lambda(v_1)$ と $\Lambda(v_2)$ の交わらない和集合で、$i\in\Lambda(v_k)$ について $d_i(v)=d_i(v_k)+1$ なので、$\pi_i(v)=\pi_i(v_k)(1+\delta_v)$ とおけば主張の形です。 根 $r$ に適用すると $\hat s=\sum_ix_i\,\pi_i(r)$ で、$\pi_i(r)$ は $d_i$ 個の因子の積です。$d_iu\le hu\lt1$ なので、誤差因子の積の補題（）から $\pi_i(r)=1+\theta_i$, $|\theta_i|\le\gamma_{d_i}$ です（$d_i=0$ なら $\theta_i=0$）。$\gamma_k$ は $k$ について増加（ の 1）なので $\gamma_{d_i}\le\gamma_h$ です。 **$h\le n-1$.** 内部頂点の個数を $m$ とします。根以外の頂点はちょうど 1 つの親をもつので、辺の本数は（頂点の総数）$-1=m+n-1$ です。一方、辺は内部頂点から子へ 2 本ずつ出るので $2m$ 本です。よって $2m=m+n-1$、$m=n-1$ です。根から葉への道は同じ頂点を 2 度通らないので、その上の内部頂点は $n-1$ 個以下です。 逐次和の誤差（）と対ごとの和の誤差（）は、この補題で木を一本道・釣り合った木にとった場合です。 $x_1,\ldots,x_n\in\mathbb{F}$（$n\ge2$）、$(n-1)u\lt1$ とする。同じ $x_1,\ldots,x_n$ の 2 つの和の木の計算結果を $\hat s_1,\hat s_2$、高さを $h_1,h_2$ とし、どちらの木でも の仮定（オーバーフローもアンダーフローも起きない）が成り立つとする。$s=\sum_ix_i$ とおくと $$ |\hat s_k-s|\le\gamma_{h_k}\sum_{i=1}^n|x_i|\quad(k=1,2),\qquad |\hat s_1-\hat s_2|\le(\gamma_{h_1}+\gamma_{h_2})\sum_{i=1}^n|x_i|\le2\gamma_{n-1}\sum_{i=1}^n|x_i| $$ である。 から $h_k\le n-1$ なので $h_ku\lt1$ で、$\hat s_k=\sum_ix_i(1+\theta^{(k)}_i)$, $|\theta^{(k)}_i|\le\gamma_{h_k}$ です。したがって $\hat s_k-s=\sum_ix_i\theta^{(k)}_i$ で、三角不等式から $|\hat s_k-s|\le\gamma_{h_k}\sum_i|x_i|$ です。もう一度三角不等式を使って $$ |\hat s_1-\hat s_2|\le|\hat s_1-s|+|s-\hat s_2|\le(\gamma_{h_1}+\gamma_{h_2})\sum_i|x_i| $$ で、$\gamma_{h_k}\le\gamma_{n-1}$（ の 1）から最後の不等式が出ます。 この定理から、順序を変えて結果が変わったときの**判定規則**が得られます。 1. 2 つの結果の差 $|\hat s_1-\hat s_2|$ と、$2\gamma_{n-1}\sum|x_i|$（binary64 ではほぼ $2.2\times10^{-16}\,(n-1)\sum|x_i|$）を比べる。 2. 差が上界の内側なら、順序の違いによる丸めの差として説明がつく。水準 2（許容誤差つき）の再現として受け入れてよい。 3. 差が上界を超えていたら、和の順序では説明できない。入力・設定・コードのどれかが違うので、そちらを調べる。 注意が 2 つあります。第一に、上界は最悪の場合の評価で、実際の差はふつうずっと小さくなります（下の例）。上界の内側にあることは「丸めの差で説明できる」ことを言うだけで、「ほかに誤りがない」ことまでは言いません。第二に、定理が扱うのは**同じ浮動小数点数**を足す順序の違いです。足す項そのものの計算が違う場合（積和を 1 回の丸めで済ませる融合積和が使われた、初等関数の実装が違う）は、項ごとの差を別に見積もって足します。 乱数種を固定した標準正規乱数 $10^6$ 個（binary64）を、20 通りの順列で並べ替えて足します。足し方は NumPy の `np.sum`（対ごとの和の変形）と、前から順に足す逐次和（`np.cumsum` の最後の成分）の 2 通りです。正しく丸めた和は `math.fsum` で求めます。 ```python # Python 3.9.6 / NumPy 2.0.2（arm64）／乱数種 20261008 import math import numpy as np n = 10**6 u = 2.0 ** -53 # binary64 の単位丸め gamma = (n - 1) * u / (1 - (n - 1) * u) rng = np.random.default_rng(20261008) x = rng.standard_normal(n) # 正負が混じる 10^6 個 exact = math.fsum(x) # 正しく丸めた和 bound = 2 * gamma * math.fsum(np.abs(x)) pair, seq = [], [] for _ in range(20): # 20 通りの順列 y = x[rng.permutation(n)] pair.append(float(np.sum(y))) # 対ごとの和（NumPy） seq.append(float(np.cumsum(y)[-1])) # 逐次和 assert math.fsum(y) == exact # fsum は順序によらない both = pair + seq print(f"fsum {exact!r}") print(f"np.sum の幅 {max(pair) - min(pair):.3e} 異なる値 {len(set(pair))}") print(f"逐次和の幅 {max(seq) - min(seq):.3e} 異なる値 {len(set(seq))}") print(f"40 個全体の幅 {max(both) - min(both):.3e}") print(f"fsum からの最大のずれ {max(abs(s - exact) for s in both):.3e}") print(f"上界 2γΣ|x| {bound:.3e}") ``` ```text fsum 833.5802334626452 np.sum の幅 1.478e-12 異なる値 10 逐次和の幅 2.068e-10 異なる値 19 40 個全体の幅 2.068e-10 fsum からの最大のずれ 1.439e-10 上界 2γΣ|x| 1.771e-04 ``` **検算.** (a) 既知値: `math.fsum` は正しく丸めた和を返すので順序によりません。20 通りの順列すべてで同じ値になることを `assert` で確かめています。(b) 別手法: 同じ数の和が、`np.sum` では 10 通り、逐次和では 19 通りの異なる値になりました。 のとおり、順序で結果が変わります。(c) 上界: 40 個の結果の最大と最小の差 $2.1\times10^{-10}$ は、 の上界 $1.8\times10^{-4}$ の内側です。`fsum` は正しく丸めた和で、$s$ との差は $u|s|\approx9\times10^{-14}$ 以下なので、「`fsum` からのずれ」は $|\hat s_k-s|$ をほぼそのまま表しており、これも $\gamma_{n-1}\sum|x_i|$ の内側です。 上界は実際の幅より 6 桁大きくなっています。上界はすべての丸め誤差が同じ向きにそろう最悪の場合を数えるのに対し、実際には符号が打ち消し合うからです。この例から「実際の差は上界の $10^{-6}$ 倍」という一般的な主張は導けません。定理が保証するのは上界だけです。 なお、この出力のうち `np.sum` の行は**環境によって変わります**。`np.sum` が内部でどの順に足すか（和の木の形）は NumPy の版と CPU の命令セットで違うからです。同じ機械で Python 3.12.14 / NumPy 2.5.2 に替えて同じコードを実行すると、`np.sum の幅` は $9.095\times10^{-13}$、異なる値は 8 個になりました。`fsum`・逐次和・上界の行は同じでした（乱数の列と、足す順序が決まっている計算は変わらないため）。版を変えて末尾の桁が変わる、という 1 節の状況の実例で、どちらの環境の値も上界の内側にあります。 ## 4. 内容ハッシュと数値の書き出し ### 4.1. ハッシュが一致すれば同じ内容か 実行記録は、入力と出力を SHA-256（$b=256$ ビット）のハッシュで同定します。異なる内容が同じハッシュになる（衝突する）と同定を誤るので、その確率を評価します。 ハッシュ関数 $H$ の出力が $\{0,1\}^b$ 上の一様乱数として振る舞うと仮定する。相異なる $N$ 個の内容のうち、どれか 2 つのハッシュが一致する確率 $P$ は $$ P\le\frac{N(N-1)}{2^{b+1}} $$ を満たす。 同じ仮定のもとで、誕生日境界（）から $P\le\binom{N}{2}2^{-b}$ です。$\binom{N}{2}=N(N-1)/2$ なので、右辺は $N(N-1)/2^{b+1}$ です。 SHA-256 で $N=10^{12}$ 個のファイルを同定しても、上界は $10^{24}/2^{257}\approx4\times10^{-54}$ です。一方、ハッシュの先頭 8 桁の 16 進数（$b=32$）だけを鍵にすると、$N=10^5$ で上界は $10^5(10^5-1)/2^{33}\approx1.16$ となり、何も保証しません。実行記録には**ハッシュの全桁**を書きます。表示のために短縮するのは構いませんが、短縮した値を同一性の判定に使ってはいけません。 なお、この評価は偶然の衝突についてのものです。悪意のある相手が衝突を作ろうとする場合の安全性は、ハッシュ関数の衝突困難性（）の問題で、ここでは扱いません。 設定・種・版・SHA-256 を集めて JSON にする関数です。設定は辞書 1 つにまとめ、コードの中に数値を埋めません。設定のハッシュは、鍵の順序と区切りを固定した JSON 文字列からとります（辞書の書き順が違うだけで別のハッシュにならないようにするため）。 ```python # Python 3.9.6 / NumPy 2.0.2（arm64） import hashlib, json, platform, sys from datetime import datetime, timezone import numpy as np def sha256(data: bytes) -> str: return hashlib.sha256(data).hexdigest() def canonical(obj) -> bytes: # 鍵の順序と区切りを固定した JSON return json.dumps(obj, sort_keys=True, separators=(",", ":")).encode() def experiment(config, x): # 実験本体: 設定と入力と種だけで決まる rng = np.random.default_rng(config["seed"]) noise = rng.standard_normal(x.shape) return np.cumsum(x + config["sigma"] * noise) def run(config, x): started = datetime.now(timezone.utc).isoformat() out = experiment(config, x) manifest = { "config": config, "config_sha256": sha256(canonical(config)), "input_sha256": sha256(x.tobytes()), "output_sha256": sha256(out.tobytes()), "code_version": "(git rev-parse HEAD の出力を入れる)", "env": {"python": platform.python_version(), "numpy": np.__version__, "machine": platform.machine(), "byteorder": sys.byteorder}, "started": started, "finished": datetime.now(timezone.utc).isoformat(), } return out, manifest x = np.arange(1000, dtype=np.float64) / 7 out, m = run({"seed": 20261008, "sigma": 0.5}, x) text = json.dumps(m, sort_keys=True, indent=1) # ファイルに書く内容 m2 = json.loads(text) # 読み戻して再実行 out2, m3 = run(m2["config"], x) print(m["config_sha256"][:16], m["input_sha256"][:16], m["output_sha256"][:16]) print("入力の一致", m3["input_sha256"] == m2["input_sha256"]) print("出力の一致", m3["output_sha256"] == m2["output_sha256"]) out3, m4 = run({"seed": 20261009, "sigma": 0.5}, x) print("種を変えた出力", m4["output_sha256"][:16]) ``` ```text 7ce2658066add4e2 2acb4f9dc64dd69f efe3e9d0349d858f 入力の一致 True 出力の一致 True 種を変えた出力 2b78acec23b32161 ``` **検算.** (a) 再実行: 読み戻した設定で実行した出力のハッシュが、記録したものと一致しました（同じ環境でのビット単位の再現）。(b) 感度: 種を 1 だけ変えると出力のハッシュは別の値になり、ハッシュが設定の違いを拾うことがわかります。(c) 範囲: 一致を確かめたのは同じ機械・同じ版での再実行です。別の環境でハッシュが一致しなくても、すぐに誤りとは言えません。そのときは出力の数値を比べ、差を 3 節の判定規則にかけます。表示した 16 桁は見やすさのための短縮で、比較は全桁で行っています。 ### 4.2. 浮動小数を 10 進で書いて読み戻す 設定や結果を JSON に書くと、binary64 の数は 10 進の文字列になります。読み戻したときに元の数に戻らなければ、設定のハッシュは同じでも計算は別物になります。 $\mathbb{F}$ を binary64、$x\in\mathbb{F}$ を $x_{\min}\le x\lt x_{\max}$ の正規化数とし、$10^k\le x\lt10^{k+1}$ となる整数 $k$ をとる。$10^{k-16}$ の整数倍（有効数字 17 桁の 10 進数）のうち $x$ に最も近いものの 1 つを $d$ とすると、$\mathrm{fl}(d)=x$ である。 $d$ のとり方から $|d-x|\le\tfrac12\cdot10^{k-16}$ です。 $2^e\le x\lt2^{e+1}$ とし、$y\in\mathbb{F}$, $y\ne x$ を任意にとります。$|y-x|\ge ux$ を示します。$\mathbb{F}\cap[2^e,2^{e+1}]$ は間隔 $2u\cdot2^e$ の等間隔な点列で両端を含みます（。$e=e_{\max}$ のときは右端を除きますが、$x\lt x_{\max}$ なので $x$ の右隣はこの点列の中にあります）。$x\gt2^e$ なら、$x$ の両隣はこの点列の点で、$|y-x|\ge2u\cdot2^e\gt ux$ です（$x\lt2^{e+1}$）。$x=2^e$ なら、右隣との距離は $2u\cdot2^e$、左隣との距離は、$e\gt e_{\min}$ のとき $[2^{e-1},2^e]$ の間隔 $2u\cdot2^{e-1}=u\cdot2^e$、$e=e_{\min}$ のとき $[0,2^{e_{\min}}]$ の間隔 $2^{e_{\min}-p+1}=2u\cdot2^e$ で、いずれにしても $|y-x|\ge u\cdot2^e=ux$ です。 $x\ge10^k$ と $2^{53}=9007199254740992\lt10^{16}$ から $$ |d-x|\le\tfrac12\cdot10^{k-16}\lt\tfrac12\cdot2^{-53}\cdot10^k\le\tfrac12ux\le\tfrac12|y-x| $$ です。したがって $|d-y|\ge|y-x|-|d-x|\gt|d-x|$ で、$x$ は $d$ に最も近いただ 1 つの $\mathbb{F}$ の元です。$x\lt x_{\max}$ なので $x$ の右隣 $y_+\in\mathbb{F}$ があり、$d\lt x+\tfrac12|y_+-x|\lt x_{\max}$ から $\mathrm{fl}(d)$ は定義されて、$\mathrm{fl}(d)=x$ です。 16 桁では足りません。$0.1+0.2$ の計算結果を 16 桁で書くと `0.3` になり、読み戻すと別の数です（17 桁なら `0.30000000000000004` で元に戻ります）。Python での使い分けは次のとおりです。 | 書き方 | 往復するか | 向く用途 | |---|---|---| | `repr(x)`（`json.dumps` も同じ） | する。往復する最短の 10 進（17 桁以下）を出す | 設定ファイル、人が読む記録 | | `f"{x:.17g}"` | する（） | 桁数をそろえた表 | | `x.hex()` と `float.fromhex` | する。2 進の値をそのまま書くので丸めが入らない | ビット単位で比べたい値、別言語との受け渡し | | `f"{x:.6g}"` など短い桁 | しない | 表示だけ。再計算の入力には使わない | ## 5. 乱数種と決定性 ### 5.1. 生成器は状態機械 **擬似乱数生成器**とは、有限集合 $S$（状態の集合）、写像 $T:S\to S$（遷移）、$g:S\to V$（出力）、種の集合 $\Sigma$ から $S$ への写像 $\iota$（初期化）の組である。種 $\sigma\in\Sigma$ に対し $$ s_0=\iota(\sigma),\qquad v_k=g(s_k),\qquad s_{k+1}=T(s_k)\quad(k=0,1,2,\ldots) $$ で定まる列 $(v_k)_{k\ge0}$ を、種 $\sigma$ の**出力列**という。 **乱数を消費する決定的なプログラム**とは、集合 $Q$（プログラムの状態）、初期状態 $q_0\in Q$、停止状態の集合 $Q_{\mathrm{halt}}\subset Q$、写像 $F:Q\times V\to Q$、$\mathrm{out}:Q_{\mathrm{halt}}\to Y$ の組である。生成器から値を 1 つ受け取るたびに $q_{k+1}=F(q_k,v_k)$ と進み、初めて $q_k\in Q_{\mathrm{halt}}$ となったところで $\mathrm{out}(q_k)$ を出力して止まる。 要するに、生成器は「種から決まる数の列を順に取り出す機械」で、乱数を使うプログラムは「取り出した値に応じて自分の状態を更新していく機械」です。たとえば 1 つの生成器を共有する 2 つの消費者 A, B のうち、A が最初に受け取る値が $v_0$ なのか $v_1$ なのかで結果が変わる、という状況をこの形で扱います。 $q_0$ は のコード・入力・設定で決まる部分です。生成器の中身（周期の長さ、出力列の統計的な性質）はこの記事では扱わず、状態機械であることだけを使います。 の記号を使う。 1. 出力列 $(v_k)$ は種 $\sigma$ だけで決まる。 2. 乱数を消費する決定的なプログラムが止まるとき、その出力は $(q_0,\sigma)$ だけで決まる。 3. 2 つの消費者 A, B が 1 つの生成器を共有し、それぞれ値を 1 つ受け取ってそのまま出力するとする。A が先に受け取る実行では出力の組は $(v_0,v_1)$、B が先に受け取る実行では $(v_1,v_0)$ で、$v_0\ne v_1$ ならば種が同じでも 2 つの実行の結果は異なる。 4. A, B がそれぞれ自分専用の生成器（種 $\sigma_{\mathrm{A}},\sigma_{\mathrm{B}}$）をもち、互いの状態を読み書きしないならば、A の出力は $(q^{\mathrm{A}}_0,\sigma_{\mathrm{A}})$ だけ、B の出力は $(q^{\mathrm{B}}_0,\sigma_{\mathrm{B}})$ だけで決まり、A と B の手続きがどんな順序で交互に進んでも結果は変わらない。 1. $k$ についての帰納法で、$s_k$ が $\sigma$ だけで決まることを示します。$s_0=\iota(\sigma)$ です。$s_k$ が $\sigma$ で決まれば $s_{k+1}=T(s_k)$ も決まります。したがって $v_k=g(s_k)$ も $\sigma$ だけで決まります。 2. $k$ についての帰納法で、$q_k$ が $(q_0,\sigma)$ だけで決まることを示します。$k=0$ は自明です。$q_k$ が決まれば、1 から $v_k$ は $\sigma$ で決まるので、$q_{k+1}=F(q_k,v_k)$ も決まります。止まる番号（初めて $q_k\in Q_{\mathrm{halt}}$ となる $k$）も列 $(q_k)$ で決まるので、出力 $\mathrm{out}(q_k)$ も $(q_0,\sigma)$ だけで決まります。 3. 生成器は呼ばれた順に $v_0,v_1$ を返します。A が先なら A が $v_0$、B が $v_1$ を受け取り、B が先ならその逆です。出力の組（A の出力, B の出力）は $(v_0,v_1)$ と $(v_1,v_0)$ で、$v_0\ne v_1$ なら異なります。 4. A の状態の列は、A の初期状態と A の生成器が返す値の列だけから $q^{\mathrm{A}}_{k+1}=F^{\mathrm{A}}(q^{\mathrm{A}}_k,v^{\mathrm{A}}_k)$ で決まります。A の生成器を呼ぶのは A だけなので、A が $k$ 番目に受け取る値は種 $\sigma_{\mathrm{A}}$ の出力列の第 $k$ 項 $v^{\mathrm{A}}_k$ で、B がいつ何回進んだかによりません。したがって 2 を A に適用でき、A の出力は $(q^{\mathrm{A}}_0,\sigma_{\mathrm{A}})$ だけで決まります。B も同じです。 3 は、スレッドやプロセスが 1 つの生成器を取り合う並列計算で実際に起こります。どちらが先に呼ぶかは OS のスケジューラが決めるので、実行ごとに変わります。4 が対処で、**実験（あるいはワーカー）ごとに生成器を分けます**。分けた生成器の出力列が互いに重ならず独立とみなせるかは生成器の設計の問題で、確率的計算のテーマの擬似乱数の章で扱います。ここでは、NumPy がそのために用意している `SeedSequence.spawn` の使い方だけを示します。 4 つの実験（それぞれ標準正規乱数 $10^4$ 個の標本平均）を、順序 $0,1,2,3$ と順序 $3,1,0,2$ で実行します。(a) は親の種から `spawn` で作った子を実験ごとに渡し、(b) は 1 つの生成器を全実験で使い回します。 ```python # Python 3.9.6 / NumPy 2.0.2（arm64）／親の種 20261008 import numpy as np def trial(rng, m=10**4): # 1 回の実験: 標本平均 return rng.standard_normal(m).mean() K = 4 # (a) 実験ごとに独立な生成器（SeedSequence.spawn） def separate(order): children = np.random.SeedSequence(20261008).spawn(K) res = [None] * K for k in order: res[k] = trial(np.random.default_rng(children[k])) return res # (b) 1 つの生成器を全実験で共有 def shared(order): rng = np.random.default_rng(20261008) res = [None] * K for k in order: res[k] = trial(rng) return res for name, f in (("spawn", separate), ("共有", shared)): r1, r2 = f([0, 1, 2, 3]), f([3, 1, 0, 2]) print(name, [f"{v:+.6f}" for v in r1]) print(name, [f"{v:+.6f}" for v in r2], "一致" if r1 == r2 else "不一致") vals = np.array(separate(range(K))) print(f"平均 {vals.mean():+.6f} s {vals.std(ddof=1):.6f}" f" 標準誤差 {vals.std(ddof=1) / np.sqrt(K):.6f}") ``` ```text spawn ['+0.004517', '-0.005054', '+0.000175', '+0.012945'] spawn ['+0.004517', '-0.005054', '+0.000175', '+0.012945'] 一致 共有 ['+0.002960', '+0.002897', '-0.004041', '-0.013365'] 共有 ['-0.004041', '+0.002897', '-0.013365', '+0.002960'] 不一致 平均 +0.003146 s 0.007615 標準誤差 0.003807 ``` **検算.** (a) 命題との対応: 生成器を分けた場合は実行の順序を変えても各実験の結果がビット単位で一致し（ の 4）、共有した場合は同じ 4 つの値が別の実験に割り当てられます（同 3。実験 1 だけは両方の順序で 2 番目に実行されるので一致しています）。(b) 既知値: 標準正規乱数 $10^4$ 個の標本平均の標準偏差は $1/\sqrt{10^4}=0.01$ で、4 個の値から求めた $s=0.0076$ はこれと同じ大きさです。(c) 最終行は次の節の標準誤差の計算で、4 個の平均 $0.0031$ は真の値 $0$ から標準誤差 $0.0038$ の 1 倍以内にあります。 ### 5.2. 種の決め方と、版による出力列の違い 種は、実行のたびに思いつきで選ぶのではなく、設定から決まるようにします。よく使うのは次の 2 つです。 - **親の種と実験番号から派生させる。** 親の種を設定に 1 つ書き、$k$ 番目の実験には `SeedSequence(親の種).spawn(K)[k]` を渡す。上の例の方法です。 - **設定のハッシュから導く。** 設定の正規化した JSON の SHA-256 の先頭 8 バイトなどを整数にして種にする。設定を変えると種も変わり、同じ設定なら同じ種になります。 時刻やプロセス番号から種を作ると、実行ごとに種が変わります。そうしたい場合でも、**使った種を実行記録に書き出せば**再現はできます。 環境 $\eta$ に版を含める理由の 1 つは、同じ種でも**生成器が違えば出力列が違う**ことです。 の 1 が言うのは「生成器 $(S,T,g,\iota)$ を固定すれば」出力列が種で決まることでした。NumPy の古い `RandomState(0)` は Mersenne Twister（MT19937）を、`default_rng(0)` は現在は PCG64 を使うので、同じ種 `0` から最初に出る $[0,1)$ の一様乱数は、手元の NumPy 2.0.2 でそれぞれ $0.5488135\ldots$ と $0.6369616\ldots$ です。NumPy は `default_rng` が選ぶ生成器と、分布に従う乱数を作る算法を将来の版で変えうると明言しています（NEP 19）。長く残す結果では、`np.random.Generator(np.random.PCG64(seed))` のように生成器の名前を明示し、NumPy の版を実行記録に書きます。 ## 6. 種を変えた繰り返し 乱数を使う実験の結果は、種を変えれば変わります。1 つの種の結果だけを報告すると、読み手はその値がどれくらい揺らぐのかを知ることができません。 $X_1,\ldots,X_K$ を独立で同じ分布に従う確率変数（$K$ 個の種の結果のモデル）とし、$E[X_i]=\mu$, $\mathrm{Var}(X_i)=\sigma^2\lt\infty$ とする。$\bar X=\frac1K\sum_{i=1}^KX_i$ とおく。 1. $E[\bar X]=\mu$, $\mathrm{Var}(\bar X)=\sigma^2/K$ である。 2. $K\ge2$ のとき、$S^2=\frac{1}{K-1}\sum_{i=1}^K(X_i-\bar X)^2$ は $E[S^2]=\sigma^2$ を満たす。したがって $S^2/K$ は $\mathrm{Var}(\bar X)$ の不偏推定量である。 3. $K=1$ のとき、分散の不偏推定量は存在しない。正確には、分布の族 $\mathcal{P}$ が 1 点に集中する分布（$P(X=c)=1$）をすべての $c\in\mathbb{R}$ について含み、さらに分散が正の分布を 1 つ以上含むならば、「$\mathcal{P}$ のどの分布に従う $X_1$ についても $E[h(X_1)]=\mathrm{Var}(X_1)$」となる関数 $h:\mathbb{R}\to\mathbb{R}$ はない。 1 と 2 は、種を変えた平均の揺らぎが $\sigma^2/K$ で、それがデータから推定できるという標準的な計算です。3 は、種が 1 個だと揺らぎを測れないことを、「どんな関数 $h$ も、分散 0 の分布と分散が正の分布の両方では合わない」という背理法で示します。 1. 期待値の線形性（ から従う基本性質）から $E[\bar X]=\frac1K\sum_iE[X_i]=\mu$ です。$X_i$ と $X_j$（$i\ne j$）は独立なので $\mathrm{Cov}(X_i,X_j)=0$ です（）。 の 2 と 4 から $$ \mathrm{Var}(\bar X)=\frac{1}{K^2}\mathrm{Var}\Bigl(\sum_{i=1}^KX_i\Bigr)=\frac{1}{K^2}\sum_{i=1}^K\mathrm{Var}(X_i)=\frac{\sigma^2}{K}. $$ 2. $X_i-\bar X=(X_i-\mu)-(\bar X-\mu)$ を 2 乗して $i$ について足します。$\sum_i(X_i-\mu)=K(\bar X-\mu)$ なので $$ \sum_{i=1}^K(X_i-\bar X)^2=\sum_{i=1}^K(X_i-\mu)^2-2(\bar X-\mu)\cdot K(\bar X-\mu)+K(\bar X-\mu)^2=\sum_{i=1}^K(X_i-\mu)^2-K(\bar X-\mu)^2 $$ です。期待値をとると、$E[(X_i-\mu)^2]=\sigma^2$ と、1 から $E[(\bar X-\mu)^2]=\mathrm{Var}(\bar X)=\sigma^2/K$ なので、右辺の期待値は $K\sigma^2-\sigma^2=(K-1)\sigma^2$ です。$K-1\gt0$ で割って $E[S^2]=\sigma^2$ です。$E[S^2/K]=\sigma^2/K=\mathrm{Var}(\bar X)$ も従います。 3. そのような $h$ があったとします。$c\in\mathbb{R}$ を任意にとり、$X_1$ が $c$ に集中する分布に従うとすると、$E[h(X_1)]=h(c)$、$\mathrm{Var}(X_1)=0$ なので $h(c)=0$ です。$c$ は任意だったので $h$ は恒等的に $0$ です。すると $\mathcal{P}$ の中の分散が正の分布について $E[h(X_1)]=0\ne\mathrm{Var}(X_1)$ となり、仮定に反します。 2 から、$K$ 個の種の結果の平均 $\bar x$ には、標準誤差の推定値 $s/\sqrt{K}$（$s$ は $S$ の実現値）を添えて報告できます。3 は、種が 1 個では揺らぎの大きさを**データから知る方法がない**ことを言っています。1 個の値は、分散 $0$ の分布から出たのか、分散の大きい分布から出たのかを区別できないからです。 の水準 3 を確かめるには、少なくとも 2 個、実際には区間が意味をもつ程度の個数の種が必要です。 この命題は $X_1,\ldots,X_K$ が**あらかじめ決めた種**の結果であることを前提にしています。20 個の種で実行して結果のよい 5 個だけを報告すると、報告された 5 個はもう同じ分布からの独立な標本ではなく、平均は偏り、$s/\sqrt{K}$ は揺らぎを過小に見せます。種の集合は実行の前に決めて設定に書き、実行した種はすべて報告します。途中で失敗した種も、失敗したという記録を残します。 ## 7. 実務のまとめ ここまでの命題を、作業の約束事に直します。このテーマの以降の章と、数値計算の各章のコード例は、この約束事に従っています。 | 約束事 | 根拠 | |---|---| | 設定は 1 つのファイル（辞書）にまとめ、コードに数値を埋めない | : 設定 $\theta$ を特定できなければ写像 $R$ の引数が決まらない | | 出力には実行記録を必ず添え、ハッシュは全桁を書く | , | | 浮動小数は `repr`（17 桁以下で往復）か `float.hex` で書く | | | 環境が変わって末尾の桁が変わったら、差を $2\gamma_{n-1}\sum\lvert x_i\rvert$ と比べる | | | 実験ごと・ワーカーごとに生成器を分ける | の 3 と 4 | | 生成器の名前とライブラリの版を記録する | の 1 は生成器を固定したときの主張 | | 種の集合を先に決め、全部の結果を標準誤差つきで報告する | | コード例の先頭の 1 行（処理系とライブラリの版、機種、乱数種）は、実行記録を 1 行に縮めたものです。書いてあるのはその例を実際に実行した環境で、このテーマでは前半の 3 章が Python 3.9.6 / NumPy 2.0.2、後半の 3 章が Python 3.12.14 / NumPy 2.5.2 です（python-flint を使う例は後者）。載せた出力は、先頭の行の環境での実行結果です。別の版では、`np.sum` のように足す順序が実装に任されている計算の末尾の桁や、数値の表示の丸め方が変わることがあります。実行時間は環境によって大きく変わるので、このテーマでは時間の計測値を載せず、演算量の式や反復回数で比べます。配列計算の書き方は[NumPy で書く配列計算](/computer-science/research-computing/numpy-array-computing)、誤差の見積もりと桁数の選び方は[精度の選び方](/computer-science/research-computing/precision-and-error-budget)で扱います。 ## 8. 演習 binary32（$p=24$）で $(a+b)+c\ne a+(b+c)$ となる $a,b,c$ を から 1 組作り、2 つの結果とその差を求めよ。 $x=1=2^{23}\cdot2^{-23}$ は $M=2^{23}$（偶数）, $e=0$ の正規化数なので、$h=2^{-24}$ として $a=1$, $b=c=2^{-24}$ ととります。命題から $\mathrm{fl}(\mathrm{fl}(1+2^{-24})+2^{-24})=1$、$\mathrm{fl}(1+\mathrm{fl}(2^{-24}+2^{-24}))=1+2^{-23}$ で、差は $2^{-23}\approx1.19\times10^{-7}$（$1$ の右隣との間隔）です。NumPy では `np.float32(1) + np.float32(2**-24) + np.float32(2**-24)` と `np.float32(1) + (np.float32(2**-24) + np.float32(2**-24))` で確かめられます。 binary64 の正の数 $n=10^8$ 個の和を 2 台の計算機で求めたところ、結果は約 $5.0\times10^{7}$ で、2 台の差は $3\times10^{-3}$ だった。(1) この差は和の順序の違いで説明できるか。(2) 一方の計算機が逐次和、他方が高さ $27$ の釣り合った木で足していたとわかった場合、上界はどう変わるか。 (1) 正の数なので $\sum|x_i|=s\approx5.0\times10^7$ です。$(n-1)u\approx10^8\times1.11\times10^{-16}=1.11\times10^{-8}$ で、$\gamma_{n-1}\approx1.11\times10^{-8}$。 の上界は $2\gamma_{n-1}\sum|x_i|\approx2\times1.11\times10^{-8}\times5.0\times10^7\approx1.1$ です。差 $3\times10^{-3}$ はこの内側なので、和の順序の違いで説明できます。 (2) 定理の精密な形 $(\gamma_{h_1}+\gamma_{h_2})\sum|x_i|$ を使います。$h_1=n-1$, $h_2=27$ で、$\gamma_{27}\approx27u\approx3.0\times10^{-15}$ は $\gamma_{n-1}$ に比べて無視でき、上界は約 $1.11\times10^{-8}\times5.0\times10^7\approx0.56$ です。差はこの内側にあります。なお、釣り合った木の側の誤差は $\gamma_{27}\sum|x_i|\approx1.5\times10^{-7}$ 以下なので、差 $3\times10^{-3}$ のほとんどは逐次和の側の誤差だとわかります。 実行記録のファイル名に、出力の SHA-256 の先頭 $t$ 桁の 16 進数を使いたい。$N=10^6$ 個の相異なる出力について、どれか 2 つのファイル名が衝突する確率を $10^{-6}$ 以下に抑えるには、$t$ をいくつ以上にすればよいか。 を使って答えよ。 16 進 $t$ 桁は $b=4t$ ビットです。$N(N-1)/2^{b+1}\le10^{-6}$ となればよいので、$2^{4t+1}\ge10^6\cdot N(N-1)$ が十分条件です。右辺は $10^{18}$ 未満で、$10^{18}\lt2^{60}$（$2^{10}=1024\gt10^3$ を 6 乗して $2^{60}\gt10^{18}$）なので、$4t+1\ge60$、すなわち $t\ge15$ で足ります。$t=14$ では $2^{57}\approx1.44\times10^{17}\lt10^{6}\cdot N(N-1)\approx10^{18}$ で、上界は約 $6.9\times10^{-6}$ となり $10^{-6}$ を超えます。したがって $t\ge15$ とします。これは上界にもとづく十分条件です。 2 つのワーカー A, B が 1 つの生成器を共有し、A は値を 2 つ、B は値を 1 つ受け取るとする。A の 2 回の受け取りと B の 1 回の受け取りが交互に進む順序をすべて挙げ、それぞれで A, B が受け取る値を、種 $\sigma$ の出力列 $v_0,v_1,v_2$ で表せ。$v_0,v_1,v_2$ が相異なるとき、A が受け取る値の組は何通りあるか。 受け取りの順序は、3 回のうち B が何番目かで決まり、BAA, ABA, AAB の 3 通りです。生成器は呼ばれた順に $v_0,v_1,v_2$ を返すので、 - BAA: B が $v_0$、A が $(v_1,v_2)$ - ABA: A が $v_0$、B が $v_1$、A が $v_2$。A は $(v_0,v_2)$ - AAB: A が $(v_0,v_1)$、B が $v_2$ です。$v_0,v_1,v_2$ が相異なれば、A が受け取る組は $(v_1,v_2)$, $(v_0,v_2)$, $(v_0,v_1)$ の 3 通りで、すべて異なります。種が同じでも、スケジューラが選ぶ順序によって A の結果は 3 通りに分かれます。A と B に別々の生成器を渡せば、 の 4 により順序によらなくなります。 ## 参考文献 - N. J. Higham, *Accuracy and Stability of Numerical Algorithms*, 2nd ed., SIAM, 2002 — 第 4 章（和の誤差解析。和の木による統一的な扱い）。 - D. Goldberg, "What every computer scientist should know about floating-point arithmetic", *ACM Computing Surveys* 23 (1991), 5–48. [doi:10.1145/103162.103163](https://doi.org/10.1145/103162.103163) — 2 進と 10 進の往復に必要な桁数。 - G. K. Sandve, A. Nekrutenko, J. Taylor, E. Hovig, "Ten simple rules for reproducible computational research", *PLoS Computational Biology* 9 (2013), e1003285. [doi:10.1371/journal.pcbi.1003285](https://doi.org/10.1371/journal.pcbi.1003285) - NumPy 公式ドキュメント, [Parallel random number generation](https://numpy.org/doc/stable/reference/random/parallel.html) — `SeedSequence.spawn` による生成器の分け方。 - R. Kern, [NEP 19 — Random number generator policy](https://numpy.org/neps/nep-0019-rng-policy.html) — 版をまたいだ出力列の互換性についての NumPy の方針。 - Python 公式ドキュメント, [Floating-point arithmetic: issues and limitations](https://docs.python.org/3/tutorial/floatingpoint.html) — `repr` が往復する最短の 10 進を返すことと `float.hex`。