§E20.33擬似乱数生成器

最終更新

計算機上の擬似乱数生成器は、有限個の値をとる状態を決まった規則で更新し、各状態から値を読み出して値の列を作る。空でない有限集合SS、写像T ⁣:S→ST\colon S\to S、集合ZZへの写像O ⁣:S→ZO\colon S\to Zの組(S,T,O)(S,T,O)としてこの仕組みを定めたものを有限状態生成器といい、TTを遷移、OOを出力写像という。初期状態s0s_0を定めると、状態列sj=Tj(s0)s_j=T^j(s_0)と出力列zj=O(sj)z_j=O(s_j)はただ一つに決まるので、遷移、出力写像、初期状態が等しい二つの生成器の出力列は一致し、同じ計算を再現することができる。他方、状態は有限個しかないから、状態列はある項以降で周期的になり、その最小周期は∣S∣|S|以下である。出力列の最小周期は状態列の最小周期を割り、それより小さいことがある。たとえばs∈{0,…,15}s\in\{0,\ldots,15\}を5s+15s+1を1616で割った余りへ写す遷移では、初期状態00からの状態列は1616個の状態をちょうど一度ずつ通って周期1616で戻るが、各状態を22で割った余りの列は0,1,0,1,…0,1,0,1,\ldotsである。本記事では、有限状態生成器の周期を調べ、その出力が相互独立で一様な値の列とどこで異なるかを明らかにする。

1 有限状態生成器

定義 1.1.

  1. SSを空でない有限集合、ZZを集合とし、T ⁣:S→ST\colon S\to SとO ⁣:S→ZO\colon S\to Zを写像とする。三つ組(S,T,O)(S,T,O)を 有限状態生成器 (finite-state generator) といい、SSを状態集合、TTを遷移、OOを出力写像という。s0∈Ss_0\in Sに対してsj:=Tj(s0)s_j:=T^j(s_0)、zj:=O(sj)z_j:=O(s_j)(j∈N≥0)(j\in\N)と置き、s0s_0を 初期状態 (initial state)、(sj)j≥0(s_j)_{j\ge0}を 状態列 (state sequence)、(zj)j≥0(z_j)_{j\ge0}を 出力列 (output sequence) という。
  2. AAを集合とし、(xj)j≥0(x_j)_{j\ge0}をAAの元の列とする。整数μ≥0\mu\ge0とλ≥1\lambda\ge1が、任意の整数j≥μj\ge\muに対してxj+λ=xjx_{j+\lambda}=x_jを満たすとき、(xj)(x_j)はμ\mu以降で周期λ\lambdaをもつという。あるμ\mu以降である周期をもつ列を 最終周期的 (eventually periodic) といい、00以降で周期をもつ列を 純周期的 (purely periodic) という。最終周期的な列に対し、あるμ\mu以降の周期となる整数λ≥1\lambda\ge1の最小値を、その列の 最小周期 (minimal period) という。

命題 1.2.(S,T,O)(S,T,O)を有限状態生成器、OOの終域をZZとし、s0∈Ss_0\in Sからの状態列を(sj)(s_j)、出力列を(zj)(z_j)とする。

  1. 整数μ≥0\mu\ge0とλ≥1\lambda\ge1の組であって、s0,…,sμ+λ−1s_0,\ldots,s_{\mu+\lambda-1}が相異なり、かつsμ+λ=sμs_{\mu+\lambda}=s_\muを満たすものがただ一つ存在する。この組はμ+λ≤∣S∣\mu+\lambda\le|S|を満たす。
  2. (1)の組(μ,λ)(\mu,\lambda)と整数j,k≥0j,k\ge0について、sj=sks_j=s_kであるための必要十分条件は、j=kj=kであるか、j,k≥μj,k\ge\muかつλ∣j−k\lambda\mid j-kであることである。
  3. 状態列は(1)のμ\mu以降で周期λ\lambdaをもつ。整数μ′≥0\mu'\ge0以降で状態列が周期λ′\lambda'をもつならば、μ≤μ′\mu\le\mu'かつλ∣λ′\lambda\mid\lambda'である。特にλ\lambdaは状態列の最小周期である。
  4. TTが全単射ならば、(1)のμ\muは00であり、状態列は純周期的である。
  5. 出力列は(1)のμ\mu以降で周期λ\lambdaをもち、出力列の最小周期はλ\lambdaを割る。
  6. (S′,T′,O′)(S',T',O')をO′ ⁣:S′→ZO'\colon S'\to Zを満たす有限状態生成器とし、写像φ ⁣:S→S′\varphi\colon S\to S'がT′∘φ=φ∘TT'\circ\varphi=\varphi\circ TとO′∘φ=OO'\circ\varphi=Oを満たすとする。(S′,T′,O′)(S',T',O')の初期状態φ(s0)\varphi(s_0)からの出力列を(zj′)(z'_j)とすると、任意のj≥0j\ge0に対してzj′=zjz'_j=z_jが成り立つ。特に、遷移、出力写像、初期状態が等しい二つの有限状態生成器の出力列は一致する。

証明.(1)を示す。s0,…,s∣S∣s_0,\ldots,s_{|S|}はSSの元の∣S∣+1|S|+1項の列であるから、鳩の巣原理により、ある0≤i<n≤∣S∣0\le i<n\le|S|に対してsi=sns_i=s_nが成り立つ。sn∈{s0,…,sn−1}s_n\in\{s_0,\ldots,s_{n-1}\}を満たす最小の正整数nnをn0n_0とすると、n0≤∣S∣n_0\le|S|であり、n0n_0の最小性によりs0,…,sn0−1s_0,\ldots,s_{n_0-1}は相異なる。sn0=sμs_{n_0}=s_\muを満たすμ<n0\mu<n_0はただ一つであり、λ:=n0−μ\lambda:=n_0-\muと置けば、組(μ,λ)(\mu,\lambda)は二条件を満たし、μ+λ=n0≤∣S∣\mu+\lambda=n_0\le|S|である。逆に組(μ′,λ′)(\mu',\lambda')が二条件を満たすとし、n:=μ′+λ′n:=\mu'+\lambda'と置く。sn=sμ′s_n=s_{\mu'}であるからsn∈{s0,…,sn−1}s_n\in\{s_0,\ldots,s_{n-1}\}である。1≤n′<n1\le n'<nならばs0,…,sn′s_0,\ldots,s_{n'}は相異なるから、sn′∉{s0,…,sn′−1}s_{n'}\notin\{s_0,\ldots,s_{n'-1}\}である。したがってn=n0n=n_0であり、sμ′=sn0s_{\mu'}=s_{n_0}かつμ′<n0\mu'<n_0であるからμ′=μ\mu'=\mu、λ′=λ\lambda'=\lambdaである。

(2)を示す。任意の整数i≥μi\ge\muに対してsi+λ=sis_{i+\lambda}=s_iであることを、iiについての帰納法で示す。i=μi=\muのときは(1)の等式である。si+λ=sis_{i+\lambda}=s_iならばsi+1+λ=T(si+λ)=T(si)=si+1s_{i+1+\lambda}=T(s_{i+\lambda})=T(s_i)=s_{i+1}である。整数j≥0j\ge0に対し、j<μj<\muのときρ(j):=j\rho(j):=jと置き、j≥μj\ge\muのときj−μj-\muをλ\lambdaで割った余りをrrとしてρ(j):=μ+r\rho(j):=\mu+rと置く。j≥μj\ge\muならば、ある整数t≥0t\ge0に対してj=ρ(j)+tλj=\rho(j)+t\lambdaであり、上の等式をtt回適用してsj=sρ(j)s_j=s_{\rho(j)}を得る。j<μj<\muならばsj=sρ(j)s_j=s_{\rho(j)}は定義から成り立つ。ρ(j)\rho(j)とρ(k)\rho(k)は{0,…,μ+λ−1}\{0,\ldots,\mu+\lambda-1\}に属し、s0,…,sμ+λ−1s_0,\ldots,s_{\mu+\lambda-1}は相異なるので、sj=sks_j=s_kはρ(j)=ρ(k)\rho(j)=\rho(k)と同値である。j<μj<\muならばρ(j)<μ\rho(j)<\muであり、j≥μj\ge\muならばρ(j)≥μ\rho(j)\ge\muであるから、ρ(j)=ρ(k)\rho(j)=\rho(k)は、j=k<μj=k<\muであるか、j,k≥μj,k\ge\muかつλ∣j−k\lambda\mid j-kであることと同値である。

(3)を示す。(2)の証明で示した等式により、状態列はμ\mu以降で周期λ\lambdaをもつ。状態列がμ′\mu'以降で周期λ′\lambda'をもつとすると、sμ′+λ′=sμ′s_{\mu'+\lambda'}=s_{\mu'}かつμ′+λ′≠μ′\mu'+\lambda'\ne\mu'であるから、(2)によりμ′≥μ\mu'\ge\muかつλ∣λ′\lambda\mid\lambda'である。λ∣λ′\lambda\mid\lambda'からλ≤λ′\lambda\le\lambda'であるので、λ\lambdaは状態列の最小周期である。

(4)を示す。TTが全単射であり、μ≥1\mu\ge1であると仮定する。T(sμ−1+λ)=sμ+λ=sμ=T(sμ−1)T(s_{\mu-1+\lambda})=s_{\mu+\lambda}=s_\mu=T(s_{\mu-1})であり、TTは単射であるから、sμ−1+λ=sμ−1s_{\mu-1+\lambda}=s_{\mu-1}が成り立つ。他方、μ−1<μ−1+λ≤μ+λ−1\mu-1<\mu-1+\lambda\le\mu+\lambda-1であり、s0,…,sμ+λ−1s_0,\ldots,s_{\mu+\lambda-1}は相異なるから、sμ−1+λ≠sμ−1s_{\mu-1+\lambda}\ne s_{\mu-1}である。この二つは両立しないので、μ=0\mu=0である。

(5)を示す。j≥μj\ge\muならばzj+λ=O(sj+λ)=O(sj)=zjz_{j+\lambda}=O(s_{j+\lambda})=O(s_j)=z_jである。出力列の最小周期をλz\lambda_zとし、出力列がμz\mu_z以降で周期λz\lambda_zをもつとする。ν:=max⁡{μ,μz}\nu:=\max\{\mu,\mu_z\}と置き、λ=aλz+r\lambda=a\lambda_z+r、a≥0a\ge0、0≤r<λz0\le r<\lambda_zと表す。j≥νj\ge\nuならば、j+r≥μzj+r\ge\mu_zであるから周期λz\lambda_zをaa回用いてzj+r=zj+r+aλz=zj+λz_{j+r}=z_{j+r+a\lambda_z}=z_{j+\lambda}であり、j≥μj\ge\muであるからzj+λ=zjz_{j+\lambda}=z_jである。r≥1r\ge1ならば出力列はν\nu以降で周期rrをもち、r<λzr<\lambda_zはλz\lambda_zの最小性と両立しない。したがってr=0r=0であり、λz∣λ\lambda_z\mid\lambdaである。

(6)を示す。sj′:=T′j(φ(s0))s'_j:=T'^j(\varphi(s_0))と置く。φ(sj)=sj′\varphi(s_j)=s'_jをjjについての帰納法で示す。j=0j=0のときは定義である。φ(sj)=sj′\varphi(s_j)=s'_jならばφ(sj+1)=φ(T(sj))=T′(φ(sj))=T′(sj′)=sj+1′\varphi(s_{j+1})=\varphi(T(s_j))=T'(\varphi(s_j))=T'(s'_j)=s'_{j+1}である。したがってzj′=O′(sj′)=O′(φ(sj))=O(sj)=zjz'_j=O'(s'_j)=O'(\varphi(s_j))=O(s_j)=z_jである。S′=SS'=S、T′=TT'=T、O′=OO'=O、φ=id⁡S\varphi=\id_Sの場合が最後の主張である。▨

注意 1.3. 計算機上の擬似乱数生成器では、状態集合と遷移、出力写像、利用者が与える種から初期状態を定める写像、整数演算の規約(法2b2^bでの剰余、剰余の代表元、符号付き整数の扱い)、および出力を受け取る呼出しの順序を固定すれば、同じ種から同じ出力が再現される。種から初期状態への写像または遷移が異なる二つの実装は、同じ種から異なる出力列を与えることがある。出力列の項は呼出しの順に一つずつ渡されるので、複数の計算が一つの生成器を呼び出すとき、呼出しの順序を変えると、同じ出力列のもとでも各計算が受け取る値が変わることがある。

命題 1.4.(S,T,O)(S,T,O)を有限状態生成器とし、OOの終域ZZは22個以上の元をもつR\Rの有限部分集合であるとする。(Ω,F,P)(\Omega,\mathcal F,P)を確率空間、σ\sigmaをSSに値をとる確率変数とし、Zj:=O(Tj(σ))Z_j:=O(T^j(\sigma))(j≥0)(j\ge0)と置く。整数k≥1k\ge1が∣Z∣k>∣S∣|Z|^k>|S|を満たすならば、Z0,…,Zk−1Z_0,\ldots,Z_{k-1}が相互独立であり、かつ各ZjZ_jがZZ上の一様分布に従う、ということは成り立たない。

証明. 写像Φ ⁣:S→Zk\Phi\colon S\to Z^kをΦ(s):=(O(s),O(T(s)),…,O(Tk−1(s)))\Phi(s):=(O(s),O(T(s)),\ldots,O(T^{k-1}(s)))で定める。Φ\Phiの像は高々∣S∣|S|個の元からなり、∣Zk∣=∣Z∣k>∣S∣|Z^k|=|Z|^k>|S|であるから、Φ\Phiの像に属さないw=(w0,…,wk−1)∈Zkw=(w_0,\ldots,w_{k-1})\in Z^kが存在する。(Z0,…,Zk−1)=Φ(σ)(Z_0,\ldots,Z_{k-1})=\Phi(\sigma)であるから、P(Z0=w0,…,Zk−1=wk−1)=0P(Z_0=w_0,\ldots,Z_{k-1}=w_{k-1})=0である。他方、Z0,…,Zk−1Z_0,\ldots,Z_{k-1}が相互独立であり、各ZjZ_jがZZ上の一様分布に従うならば、§E11.7 命題 1.2を一点集合{w0},…,{wk−1}\{w_0\},\ldots,\{w_{k-1}\}へ適用して、P(Z0=w0,…,Zk−1=wk−1)=∣Z∣−k>0P(Z_0=w_0,\ldots,Z_{k-1}=w_{k-1})=|Z|^{-k}>0となる。この二つの値は両立しないので、主張が成り立つ。▨

2 合同法

命題 2.1.m≥2m\ge2を整数とし、Um:={s∈Z∣0≤s≤m−1, gcd⁡(s,m)=1}U_m:=\{s\in\Z\mid 0\le s\le m-1,\ \gcd(s,m)=1\}と置く。aaをgcd⁡(a,m)=1\gcd(a,m)=1を満たす整数とし、s∈Ums\in U_mに対してasasをmmで割った余りをTa(s)T_a(s)と置く。

  1. TaT_aはUmU_mをUmU_mへ写す。任意のs0∈Ums_0\in U_mに対し、有限状態生成器(Um,Ta,id⁡Um)(U_m,T_a,\id_{U_m})のs0s_0からの状態列は純周期的であり、その最小周期はord⁡m(a)\operatorname{ord}_m(a)である。
  2. mmが素数ならば、ord⁡m(a′)=m−1\operatorname{ord}_m(a')=m-1を満たす整数a′a'が存在する。mmが素数でありord⁡m(a)=m−1\operatorname{ord}_m(a)=m-1ならば、任意のs0∈{1,…,m−1}s_0\in\{1,\ldots,m-1\}に対し、状態列の最小周期はm−1m-1であり、s0,…,sm−2s_0,\ldots,s_{m-2}は{1,…,m−1}\{1,\ldots,m-1\}の各元をちょうど一度ずつとる。

証明.(1)を示す。s∈Ums\in U_mならばgcd⁡(as,m)=1\gcd(as,m)=1であり、Ta(s)≡as(modm)T_a(s)\equiv as\pmod mであるからgcd⁡(Ta(s),m)=1\gcd(T_a(s),m)=1であり、Ta(s)∈UmT_a(s)\in U_mである。d:=ord⁡m(a)d:=\operatorname{ord}_m(a)と置く。jjについての帰納法によりsj≡ajs0(modm)s_j\equiv a^js_0\pmod mである。0≤j<k0\le j<kとする。gcd⁡(ajs0,m)=1\gcd(a^js_0,m)=1であるからajs0a^js_0は法mmの逆元をもち、ajs0≡aks0(modm)a^js_0\equiv a^ks_0\pmod mはak−j≡1(modm)a^{k-j}\equiv1\pmod mと同値である。§A4.12 命題 1.2により、ak−j≡1(modm)a^{k-j}\equiv1\pmod mはd∣k−jd\mid k-jと同値である。sj,sk∈{0,…,m−1}s_j,s_k\in\{0,\ldots,m-1\}であるから、sj=sks_j=s_kはsj≡sk(modm)s_j\equiv s_k\pmod mと同値であり、したがってd∣k−jd\mid k-jと同値である。よってs0,…,sd−1s_0,\ldots,s_{d-1}は相異なり、sd=s0s_d=s_0であるから、命題 1.2 (1)の組は(0,d)(0,d)である。命題 1.2 (3)により、状態列は00以降で周期ddをもつ純周期的な列であり、その最小周期はddである。

(2)を示す。mmが素数ならばZ/mZ\Z/m\Zは位数mmの体であり、§E8.7 定理 2.1により、その乗法群は位数m−1m-1の巡回群である。その生成元を剰余類とする整数a′a'をとると、gcd⁡(a′,m)=1\gcd(a',m)=1であり、a′k≡1(modm)a'^k\equiv1\pmod mを満たす最小の正整数kkは生成元の位数m−1m-1に等しいので、ord⁡m(a′)=m−1\operatorname{ord}_m(a')=m-1である。mmが素数ならばUm={1,…,m−1}U_m=\{1,\ldots,m-1\}である。ord⁡m(a)=m−1\operatorname{ord}_m(a)=m-1ならば、(1)により状態列は純周期的で最小周期はm−1m-1であるから、命題 1.2 (3)により命題 1.2 (1)の組は(0,m−1)(0,m-1)である。したがってs0,…,sm−2s_0,\ldots,s_{m-2}はUmU_mの相異なるm−1m-1個の元であり、UmU_mの各元をちょうど一度ずつとる。▨

例 2.2.m=7m=7とする。

  1. 31,…,363^1,\ldots,3^6を77で割った余りは3,2,6,4,5,13,2,6,4,5,1であるから、ord⁡7(3)=6\operatorname{ord}_7(3)=6である。乗数33、初期状態11の状態列は1,3,2,6,4,5,1,…1,3,2,6,4,5,1,\ldotsであり、最小周期は66である。
  2. 21,22,232^1,2^2,2^3を77で割った余りは2,4,12,4,1であるから、ord⁡7(2)=3\operatorname{ord}_7(2)=3である。乗数22の状態列は、初期状態11からは1,2,41,2,4を、初期状態33からは3,6,53,6,5を繰り返し、どちらの最小周期も33である。
  3. s∈{0,…,6}s\in\{0,\ldots,6\}に対して3s3sを77で割った余りをT(s)T(s)と置くと、T(0)=0T(0)=0であり、初期状態00の状態列は0,0,…0,0,\ldotsで、最小周期は11である。0∉U70\notin U_7であるから、この状態列は命題 2.1 (1)の対象ではない。

例 2.3.S={0,…,15}S=\{0,\ldots,15\}とし、s∈Ss\in Sに対して5s+15s+1を1616で割った余りをT(s)T(s)と置く。初期状態00の状態列は

0,1,6,15,12,13,2,11,8,9,14,7,4,5,10,3,0,…0,1,6,15,12,13,2,11,8,9,14,7,4,5,10,3,0,\ldots

である。s0,…,s15s_0,\ldots,s_{15}はSSの各元をちょうど一度ずつとり、s16=s0s_{16}=s_0であるから、命題 1.2 (1)の組は(0,16)(0,16)であり、状態列の最小周期は1616である。

  1. 5≡1(mod4)5\equiv1\pmod 4であるから、任意のs∈Ss\in Sに対してT(s)≡s+1(mod4)T(s)\equiv s+1\pmod 4であり、sj≡j(mod4)s_j\equiv j\pmod 4である。出力写像を「ssを22で割った余り」とする出力列は0,1,0,1,…0,1,0,1,\ldotsで最小周期は22であり、「ssを44で割った余り」とする出力列は0,1,2,3,0,…0,1,2,3,0,\ldotsで最小周期は44である。どちらの最小周期も、命題 1.2 (5)のとおり1616を割り、1616より小さい。
  2. JJを{0,…,15}\{0,\ldots,15\}に値をとり、その集合上の一様分布に従う確率変数とし、X:=sJX:=s_J、Y:=sJ+1=T(X)Y:=s_{J+1}=T(X)と置く。j↦sjj\mapsto s_jとj↦sj+1j\mapsto s_{j+1}はどちらも{0,…,15}\{0,\ldots,15\}からSSへの全単射であるから、XXとYYはどちらもSS上の一様分布に従う。他方、T(0)=1T(0)=1であるからP(X=0, Y=0)=0P(X=0,\ Y=0)=0であり、P(X=0)P(Y=0)=1/256P(X=0)P(Y=0)=1/256と異なるので、XXとYYは独立でない。

3 有限ビット整数の変換

命題 3.1.

  1. 整数n≥1n\ge1とL≥0L\ge0をL=qn+rL=qn+r、q≥0q\ge0、0≤r<n0\le r<nと表す。0≤j<n0\le j<nを満たす整数jjに対し、{0,…,L−1}\{0,\ldots,L-1\}のうちnnで割った余りがjjである元の個数は、j<rj<rならばq+1q+1、j≥rj\ge rならばqqである。
  2. 整数b≥1b\ge1に対してM:=2bM:=2^bと置き、整数1≤n≤M1\le n\le MをM=qn+rM=qn+r、0≤r<n0\le r<nと表す。XXを{0,…,M−1}\{0,\ldots,M-1\}に値をとり、その集合上の一様分布に従う確率変数とし、XXをnnで割った余りをRRとする。0≤j<n0\le j<nを満たす整数jjに対して、j<rj<rならばP(R=j)=(q+1)/MP(R=j)=(q+1)/Mであり、j≥rj\ge rならばP(R=j)=q/MP(R=j)=q/Mである。RRが{0,…,n−1}\{0,\ldots,n-1\}上の一様分布に従うための必要十分条件はr=0r=0であり、r=0r=0はnnが22の冪であることと同値である。

証明.(1)を示す。nnで割った余りがjjである非負整数は、整数t≥0t\ge0によるj+tnj+tnであり、j+tn≤L−1j+tn\le L-1はtn≤qn+(r−1−j)tn\le qn+(r-1-j)と同値である。j<rj<rならば0≤r−1−j<n0\le r-1-j<nであるから、この不等式はt≤qt\le qと同値であり、元の個数はq+1q+1である。j≥rj\ge rならばj≤n−1j\le n-1から−n≤r−1−j<0-n\le r-1-j<0であるから、この不等式はt≤q−1t\le q-1と同値であり、元の個数はqqである。

(2)を示す。P(R=j)P(R=j)は{0,…,M−1}\{0,\ldots,M-1\}のうちnnで割った余りがjjである元の個数をMMで割った値であるから、(1)をL=ML=Mに適用して、確率の式を得る。r=0r=0ならば任意のjjに対してP(R=j)=q/M=1/nP(R=j)=q/M=1/nである。r≥1r\ge1ならばn−1≥rn-1\ge rであるから、P(R=0)=(q+1)/M≠q/M=P(R=n−1)P(R=0)=(q+1)/M\ne q/M=P(R=n-1)であり、RRは一様分布に従わない。r=0r=0はn∣2bn\mid 2^bと同値であり、素因数分解の一意性により、n∣2bn\mid2^bはn=2in=2^iを満たす整数0≤i≤b0\le i\le bが存在することと同値である。n≤Mn\le Mであるから、後者はnnが22の冪であることと同値である。▨

例 3.2.M=16M=16とする。n=3n=3ならば16=5⋅3+116=5\cdot3+1であるから、P(R=0)=6/16P(R=0)=6/16、P(R=1)=P(R=2)=5/16P(R=1)=P(R=2)=5/16である。n=6n=6ならば16=2⋅6+416=2\cdot6+4であるから、0≤j≤30\le j\le3ではP(R=j)=3/16P(R=j)=3/16、j=4,5j=4,5ではP(R=j)=2/16P(R=j)=2/16である。M=264M=2^{64}、n=3n=3ならば264≡1(mod3)2^{64}\equiv1\pmod3であるからr=1r=1であり、P(R=0)P(R=0)はP(R=1)=P(R=2)P(R=1)=P(R=2)より2−642^{-64}だけ大きい。

定理 3.3. 整数b≥1b\ge1に対してM:=2bM:=2^bと置き、整数1≤n≤M1\le n\le Mに対してq:=⌊M/n⌋q:=\lfloor M/n\rfloor、L:=qnL:=qn、θ:=L/M\theta:=L/Mと置く。確率空間(Ω,F,P)(\Omega,\mathcal F,P)上で{0,…,M−1}\{0,\ldots,M-1\}に値をとる確率変数の列(Xk)k≥1(X_k)_{k\ge1}は相互独立であり、各XkX_kは{0,…,M−1}\{0,\ldots,M-1\}上の一様分布に従うとする。Ω∞:=⋂k≥1{Xk≥L}\Omega_\infty:=\bigcap_{k\ge1}\{X_k\ge L\}と置く。ω∈Ω∖Ω∞\omega\in\Omega\setminus\Omega_\inftyに対してN(ω):=min⁡{k≥1∣Xk(ω)<L}N(\omega):=\min\{k\ge1\mid X_k(\omega)<L\}と置き、XN(ω)(ω)X_{N(\omega)}(\omega)をnnで割った余りをY(ω)Y(\omega)と置く。ω∈Ω∞\omega\in\Omega_\inftyに対してはN(ω):=1N(\omega):=1、Y(ω):=0Y(\omega):=0と置く。00=10^0=1と約束する。

  1. 整数k≥1k\ge1と0≤j<n0\le j<nに対して、P(N=k, Y=j)=(1−θ)k−1q/MP(N=k,\ Y=j)=(1-\theta)^{k-1}q/Mが成り立つ。
  2. P(Ω∞)=0P(\Omega_\infty)=0である。
  3. YYは{0,…,n−1}\{0,\ldots,n-1\}上の一様分布に従い、NNとYYは独立である。
  4. 1/2<θ≤11/2<\theta\le1であり、N−1N-1は幾何分布Geom⁡(θ)\operatorname{Geom}(\theta)に従う。NNは可積分であり、E[N]=M/L<2E[N]=M/L<2が成り立つ。

証明.B:={L,…,M−1}B:=\{L,\ldots,M-1\}と置き、0≤j<n0\le j<nに対して、{0,…,L−1}\{0,\ldots,L-1\}のうちnnで割った余りがjjである元の集合をCjC_jと置く。L=qnL=qnであるから、命題 3.1 (1)により∣Cj∣=q|C_j|=qである。整数k≥1k\ge1と0≤j<n0\le j<nに対して

Ak,j:={X1∈B,…,Xk−1∈B, Xk∈Cj}A_{k,j}:=\{X_1\in B,\ldots,X_{k-1}\in B,\ X_k\in C_j\}

と置く。相互独立性の定義によりX1,…,XkX_1,\ldots,X_kは相互独立であり、BBとCjC_jは有限集合であるから Borel 集合である。§E11.7 命題 1.2により

P(Ak,j)=P(X1∈B)⋯P(Xk−1∈B) P(Xk∈Cj)=(M−LM)k−1qM=(1−θ)k−1qMP(A_{k,j})=P(X_1\in B)\cdots P(X_{k-1}\in B)\,P(X_k\in C_j)=\left(\frac{M-L}{M}\right)^{k-1}\frac{q}{M}=(1-\theta)^{k-1}\frac qM

である。

n≤Mn\le Mであるからq≥1q\ge1であり、L≥n≥1L\ge n\ge1であるからθ>0\theta>0である。0≤1−θ<10\le1-\theta<1であるから、§E11.5 補題 1.1により

∑k≥1∑j=0n−1P(Ak,j)=nqM∑k≥1(1−θ)k−1=θ⋅1θ=1\sum_{k\ge1}\sum_{j=0}^{n-1}P(A_{k,j})=n\frac qM\sum_{k\ge1}(1-\theta)^{k-1}=\theta\cdot\frac1\theta=1

である。事象Ak,jA_{k,j}(k≥1, 0≤j<n)(k\ge1,\ 0\le j<n)は互いに素であり、その和集合はΩ∖Ω∞\Omega\setminus\Omega_\inftyである。実際、各XkX_kは{0,…,L−1}∪B\{0,\ldots,L-1\}\cup Bに値をとり、{0,…,L−1}\{0,\ldots,L-1\}はC0,…,Cn−1C_0,\ldots,C_{n-1}の互いに素な和集合であるから、ω∈Ω∖Ω∞\omega\in\Omega\setminus\Omega_\inftyはAN(ω),Y(ω)A_{N(\omega),Y(\omega)}にだけ属し、Ω∞\Omega_\inftyの元はどのAk,jA_{k,j}にも属さない。したがってP(Ω∖Ω∞)=1P(\Omega\setminus\Omega_\infty)=1であり、(2)が成り立つ。{N=k, Y=j}\{N=k,\ Y=j\}は、(k,j)≠(1,0)(k,j)\ne(1,0)ならばAk,jA_{k,j}に等しく、(k,j)=(1,0)(k,j)=(1,0)ならばA1,0∪Ω∞A_{1,0}\cup\Omega_\inftyに等しい。よってNNとYYは確率変数であり、P(Ω∞)=0P(\Omega_\infty)=0からP(N=k, Y=j)=P(Ak,j)P(N=k,\ Y=j)=P(A_{k,j})であるので、(1)が成り立つ。

(3)を示す。(1)と§E11.5 補題 1.1により、0≤j<n0\le j<nに対して

P(Y=j)=qM∑k≥1(1−θ)k−1=qMθ=qL=1nP(Y=j)=\frac qM\sum_{k\ge1}(1-\theta)^{k-1}=\frac q{M\theta}=\frac qL=\frac1n

であり、k≥1k\ge1に対してP(N=k)=n(1−θ)k−1q/M=θ(1−θ)k−1P(N=k)=n(1-\theta)^{k-1}q/M=\theta(1-\theta)^{k-1}である。したがってP(N=k, Y=j)=P(N=k)P(Y=j)P(N=k,\ Y=j)=P(N=k)P(Y=j)である。 Borel 集合D1,D2⊂RD_1,D_2\subset\Rに対し、NNとYYは整数値であるから、確率の可算加法性により

P(N∈D1, Y∈D2)=∑k∈D1∩N≥1 ∑j∈D2∩{0,…,n−1}P(N=k)P(Y=j)=P(N∈D1)P(Y∈D2)P(N\in D_1,\ Y\in D_2)=\sum_{k\in D_1\cap\NN}\ \sum_{j\in D_2\cap\{0,\ldots,n-1\}}P(N=k)P(Y=j)=P(N\in D_1)P(Y\in D_2)

であり、§E11.7 命題 1.2によりNNとYYは独立である。

(4)を示す。M−LM-LはMMをnnで割った余りであるからM−L<n≤qn=LM-L<n\le qn=Lであり、M<2LM<2Lからθ>1/2\theta>1/2である。L≤ML\le Mからθ≤1\theta\le1である。整数i≥0i\ge0に対してP(N−1=i)=θ(1−θ)iP(N-1=i)=\theta(1-\theta)^iであるから、N−1N-1はGeom⁡(θ)\operatorname{Geom}(\theta)に従う。§E11.5 命題 1.5によりE[N−1]=(1−θ)/θ<∞E[N-1]=(1-\theta)/\theta<\inftyであるから、非負確率変数N−1N-1は可積分である。§E11.4 命題 1.2によりN=(N−1)+1N=(N-1)+1は可積分であり、

E[N]=1−θθ+1=1θ=MLE[N]=\frac{1-\theta}\theta+1=\frac1\theta=\frac ML

である。θ>1/2\theta>1/2からE[N]<2E[N]<2である。▨

例 3.4.M=16M=16とする。n=3n=3ならばq=5q=5、L=15L=15、θ=15/16\theta=15/16、E[N]=16/15E[N]=16/15である。n=6n=6ならばq=2q=2、L=12L=12、θ=3/4\theta=3/4、E[N]=4/3E[N]=4/3である。n=9n=9ならばq=1q=1、L=9L=9、θ=9/16\theta=9/16、E[N]=16/9E[N]=16/9である。n=1n=1とn=16n=16ではどちらもL=16L=16、θ=1\theta=1であり、N=1N=1がほとんど確実に成り立つ。一般にb≥1b\ge1とn=2b−1+1n=2^{b-1}+1に対してはq=1q=1、θ=1/2+2−b\theta=1/2+2^{-b}であるから、任意のc>1/2c>1/2に対して1/2+2−b<c1/2+2^{-b}<cとなるbbがあり、定理 3.3 (4)の下界1/21/2を1/21/2より大きい定数で置き換えることはできない。

注意 3.5.定理 3.3の仮定は、入力の列(Xk)(X_k)の相互独立性と一様性である。出力写像の終域が{0,…,M−1}\{0,\ldots,M-1\}である有限状態生成器の出力を入力に用いる場合、初期状態を確率変数としても、命題 1.4により、Mk>∣S∣M^k>|S|を満たすkk個の連続する出力は相互独立で一様とはならず、この仮定は満たされない。他方、一周期の中で{0,…,M−1}\{0,\ldots,M-1\}の各元がちょうど一度ずつ現れる出力列では、命題 3.1 (1)により、一周期のうちLL未満の出力をnnで割った余りは、各jjについてちょうどqq回現れる。例 2.3の状態列を出力列としn=3n=3とすると、一周期の1616個の出力のうち1515を除く1515個の余りは、0,1,20,1,2をそれぞれ55回ずつとる。

4 整数から一様格子への変換

命題 4.1. 整数b≥1b\ge1に対してM:=2bM:=2^bと置き、

Gb:={kM | k∈Z, 0≤k≤M−1},Hb:={2k+12M | k∈Z, 0≤k≤M−1}G_b:=\left\{\frac kM\ \middle|\ k\in\Z,\ 0\le k\le M-1\right\},\qquad H_b:=\left\{\frac{2k+1}{2M}\ \middle|\ k\in\Z,\ 0\le k\le M-1\right\}

と置く。F=F(2,p,emin⁡,emax⁡)F=F(2,p,e_{\min},e_{\max})をemax⁡≥0e_{\max}\ge0を満たす浮動小数点数系とする。

  1. 写像x↦x/Mx\mapsto x/Mは{0,…,M−1}\{0,\ldots,M-1\}からGbG_bへの全単射であり、写像x↦(x+1/2)/Mx\mapsto(x+1/2)/Mは{0,…,M−1}\{0,\ldots,M-1\}からHbH_bへの全単射である。Gb⊂[0,1)G_b\subset[0,1)の最小元は00、最大元は1−1/M1-1/Mであり、Hb⊂(0,1)H_b\subset(0,1)の最小元は1/(2M)1/(2M)、最大元は1−1/(2M)1-1/(2M)である。XXが{0,…,M−1}\{0,\ldots,M-1\}上の一様分布に従うならば、X/MX/Mと(X+1/2)/M(X+1/2)/MはそれぞれGbG_b上とHbH_b上の一様分布に従う。
  2. 整数c≥1c\ge1とkkが1≤k<2c1\le k<2^c、k<2pk<2^p、emin⁡≤−ce_{\min}\le-cを満たすならば、k/2c∈Fk/2^c\in Fである。
  3. p≥bp\ge bかつemin⁡≤−be_{\min}\le-bならばGb⊂FG_b\subset Fであり、p≥b+1p\ge b+1かつemin⁡≤−(b+1)e_{\min}\le-(b+1)ならばHb⊂FH_b\subset Fである。

証明.(1)を示す。二つの写像は狭義単調増加であり、GbG_bとHbH_bの定義によりそれぞれの像はGbG_bとHbH_bである。最小元と最大元はx=0x=0とx=M−1x=M-1の値である。0≤k≤M−10\le k\le M-1に対してP(X/M=k/M)=P((X+1/2)/M=(2k+1)/(2M))=P(X=k)=1/MP(X/M=k/M)=P((X+1/2)/M=(2k+1)/(2M))=P(X=k)=1/Mである。

(2)を示す。2t≤k<2t+12^t\le k<2^{t+1}を満たす整数ttをとり、e:=t−ce:=t-cと置く。1≤k<2c1\le k<2^cから0≤t≤c−10\le t\le c-1であり、−c≤e≤−1-c\le e\le-1であるからemin⁡≤e≤emax⁡e_{\min}\le e\le e_{\max}である。k<2pk<2^pからt≤p−1t\le p-1であるので、m:=k2p−1−tm:=k2^{p-1-t}は整数である。2t≤k<2t+12^t\le k<2^{t+1}の各辺に2p−1−t2^{p-1-t}を掛けて2p−1≤m≤2p−12^{p-1}\le m\le2^p-1を得る。k/2c=m2e+1−pk/2^c=m2^{e+1-p}であるから、k/2ck/2^cは指数eeの正規化数であり、FFに属する。

(3)を示す。0∈F0\in Fである。1≤x≤M−11\le x\le M-1ならば、p≥bp\ge bのときx<2b≤2px<2^b\le2^pであるから、(2)をc=bc=b、k=xk=xに適用してx/M∈Fx/M\in Fを得る。0≤x≤M−10\le x\le M-1ならば(x+1/2)/M=(2x+1)/2b+1(x+1/2)/M=(2x+1)/2^{b+1}であり、1≤2x+1<2b+11\le2x+1<2^{b+1}である。p≥b+1p\ge b+1のとき2x+1<2b+1≤2p2x+1<2^{b+1}\le2^pであるから、(2)をc=b+1c=b+1、k=2x+1k=2x+1に適用して(x+1/2)/M∈F(x+1/2)/M\in Fを得る。▨

例 4.2. binary64 形式の有限数の全体であるF=F(2,53,−1022,1023)F=F(2,53,-1022,1023)を考え、fl⁡\operatorname{fl}をFFの最近接丸めとする。M:=264M:=2^{64}と置き、XXを{0,…,M−1}\{0,\ldots,M-1\}に値をとり、その集合上の一様分布に従う確率変数とする。

  1. G53⊂FG_{53}\subset FかつH52⊂FH_{52}\subset Fである。
  2. 整数xxが1≤x<2531\le x<2^{53}を満たすならば、x/M∈Fx/M\in Fであり、fl⁡(x′/M)=x/M\operatorname{fl}(x'/M)=x/Mを満たすx′∈{0,…,M−1}x'\in\{0,\ldots,M-1\}はxxだけである。特にP(fl⁡(X/M)=x/M)=1/MP(\operatorname{fl}(X/M)=x/M)=1/Mである。
  3. 整数kkが252<k<2532^{52}<k<2^{53}を満たすとし、yk:=k2−53y_k:=k2^{-53}と置く。yk∈Fy_k\in Fであり、∣x−211k∣≤210−1|x-2^{11}k|\le2^{10}-1を満たす任意の整数xxに対してfl⁡(x/M)=yk\operatorname{fl}(x/M)=y_kである。特にP(fl⁡(X/M)=yk)≥(211−1)/MP(\operatorname{fl}(X/M)=y_k)\ge(2^{11}-1)/Mであり、fl⁡(X/M)\operatorname{fl}(X/M)はその値の集合の上の一様分布に従わない。
  4. fl⁡((M−1)/M)=1\operatorname{fl}((M-1)/M)=1である。特に、X/M<1X/M<1がつねに成り立つにもかかわらず、P(fl⁡(X/M)=1)≥1/MP(\operatorname{fl}(X/M)=1)\ge1/Mである。

証明.(1)を示す。FFはp=53p=53、emin⁡=−1022e_{\min}=-1022、emax⁡=1023≥0e_{\max}=1023\ge0の浮動小数点数系である。b=53b=53はp≥bp\ge bとemin⁡≤−be_{\min}\le-bを満たし、b=52b=52はp≥b+1p\ge b+1とemin⁡≤−(b+1)e_{\min}\le-(b+1)を満たすので、命題 4.1 (3)によりG53⊂FG_{53}\subset FかつH52⊂FH_{52}\subset Fである。

(2)を示す。1≤x<2531\le x<2^{53}であるから、命題 4.1 (2)をc=64c=64、k=xk=xに適用してx/M∈Fx/M\in Fを得る。0∈F0\in Fであるから、0≤x′<2530\le x'<2^{53}を満たす任意の整数x′x'に対してx′/M∈Fx'/M\in Fであり、§E20.1 定理 2.2 (1)によりfl⁡(x′/M)=x′/M\operatorname{fl}(x'/M)=x'/Mである。したがって、0≤x′<2530\le x'<2^{53}を満たす整数x′x'のうちfl⁡(x′/M)=x/M\operatorname{fl}(x'/M)=x/Mを満たすものはxxだけである。整数x′x'が253≤x′≤M−12^{53}\le x'\le M-1を満たすとする。命題 4.1 (2)をc=11c=11、k=1k=1に適用して2−11∈F2^{-11}\in Fであり、x′/M≥253/M=2−11x'/M\ge2^{53}/M=2^{-11}である。y∈Fy\in Fがy<2−11y<2^{-11}を満たすならば∣y−x′/M∣>∣2−11−x′/M∣|y-x'/M|>|2^{-11}-x'/M|であるから、fl⁡(x′/M)≥2−11>x/M\operatorname{fl}(x'/M)\ge2^{-11}>x/Mである。したがってfl⁡(x′/M)=x/M\operatorname{fl}(x'/M)=x/Mを満たすx′∈{0,…,M−1}x'\in\{0,\ldots,M-1\}はxxだけであり、P(fl⁡(X/M)=x/M)=P(X=x)=1/MP(\operatorname{fl}(X/M)=x/M)=P(X=x)=1/Mである。

(3)を示す。§E20.1 補題 1.2 (1)をe=−1e=-1に適用すると

F∩[1/2,1)={m2−53∣m∈Z, 252≤m≤253−1}F\cap[1/2,1)=\{m2^{-53}\mid m\in\Z,\ 2^{52}\le m\le2^{53}-1\}

であるから、yk∈Fy_k\in Fである。y∈F∖{yk}y\in F\setminus\{y_k\}とする。y∈[1/2,1)y\in[1/2,1)ならばy−yky-y_kは2−532^{-53}の零でない整数倍であり、y<1/2y<1/2ならば∣y−yk∣>yk−1/2≥2−53|y-y_k|>y_k-1/2\ge2^{-53}、y≥1y\ge1ならば∣y−yk∣≥1−yk≥2−53|y-y_k|\ge1-y_k\ge2^{-53}であるから、いずれの場合も∣y−yk∣≥2−53|y-y_k|\ge2^{-53}である。整数xxが∣x−211k∣≤210−1|x-2^{11}k|\le2^{10}-1を満たすとし、z:=x/Mz:=x/Mと置く。ykM=211ky_kM=2^{11}kであるから∣z−yk∣=∣x−211k∣/M≤(210−1)2−64<2−54|z-y_k|=|x-2^{11}k|/M\le(2^{10}-1)2^{-64}<2^{-54}であり、y∈F∖{yk}y\in F\setminus\{y_k\}に対して

∣y−z∣≥∣y−yk∣−∣z−yk∣>2−53−2−54=2−54>∣z−yk∣|y-z|\ge|y-y_k|-|z-y_k|>2^{-53}-2^{-54}=2^{-54}>|z-y_k|

であるから、fl⁡(z)=yk\operatorname{fl}(z)=y_kである。このような整数xxは211−12^{11}-1個あり、252<k<2532^{52}<k<2^{53}から0<x<M0<x<Mを満たすので、P(fl⁡(X/M)=yk)≥(211−1)/MP(\operatorname{fl}(X/M)=y_k)\ge(2^{11}-1)/Mである。他方、(2)をx=1x=1に適用すると、fl⁡(X/M)\operatorname{fl}(X/M)は値1/M1/Mを確率1/M1/Mでとる。(211−1)/M>1/M(2^{11}-1)/M>1/Mであるから、fl⁡(X/M)\operatorname{fl}(X/M)はその値の集合の上の一様分布に従わない。

(4)を示す。z:=(M−1)/M=1−2−64z:=(M-1)/M=1-2^{-64}と置く。§E20.1 補題 1.2 (1)をe=0e=0とe=−1e=-1に適用すると、1=252⋅2−521=2^{52}\cdot2^{-52}はFFに属し、11未満のFFの最大元は1−2−531-2^{-53}である。y∈Fy\in Fがy<1y<1ならば∣y−z∣≥2−53−2−64>2−64=∣1−z∣|y-z|\ge2^{-53}-2^{-64}>2^{-64}=|1-z|であり、y>1y>1ならば∣y−z∣>∣1−z∣|y-z|>|1-z|であるから、fl⁡(z)=1\operatorname{fl}(z)=1である。X≤M−1X\le M-1であるからX/M<1X/M<1がつねに成り立ち、X=M−1X=M-1ならばfl⁡(X/M)=fl⁡(z)=1\operatorname{fl}(X/M)=\operatorname{fl}(z)=1であるから、P(fl⁡(X/M)=1)≥P(X=M−1)=1/MP(\operatorname{fl}(X/M)=1)\ge P(X=M-1)=1/Mである。▨

注意 4.3.命題 4.1 (3)は、実数x/Mx/Mと(x+1/2)/M(x+1/2)/MがFFに属することを述べる。計算機上の変換は、整数xxから浮動小数点数への変換、1/21/2の加算、MMによる除算のような演算の列である。各演算の結果がその厳密な結果の最近接丸めであるとき、厳密な中間結果がすべてFFに属するならば、§E20.1 定理 2.2 (1)により計算値は格子値に等しい。たとえば(x+1/2)/M(x+1/2)/Mをx+1/2x+1/2を先に計算して求める場合、中間結果x+1/2=(2x+1)/2x+1/2=(2x+1)/2はx=M−1x=M-1で[2b−1,2b)[2^{b-1},2^b)に属するので、それがFFに属するには命題 4.1 (3)の条件に加えてemax⁡≥b−1e_{\max}\ge b-1が要る。

前提記事