乱数の原器としての円周率

乱数の原器としての円周率

三好和憲

11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111

1

.

はじめに

真の乱数，すなわちあらゆるパターンが等しい確率で かつ独立に出現し得るような無限数列を得ることが可能 であろうか.人間がこれを自ら生成すること，たとえば サイコロを完全に均等に作りまた均等な配位で同等の条 件で(すなわち重力ポテンシャルなどを全く同じにして) 振ることなどは原理的に不可能である. 自分で作り出さなくても世の中に真に確率的な現象が 存在するならば，これを目に見える形に取り出すことを 考えればよい.放射性物質の崩療やツェナーダイオード の熱雑音を利用する物理乱数はこの思想にしたがったも のである.しかし仮に現象そのものが確率的であったと してもこれを取り出す電気回路への環境の影響を完全に 排除することは難しく，得られた出力が本当に確率的で あるかどうかはわからない. 物理法則とならんで宇宙の創世者から与えられたもの として数がある.無理数の無限小数による表現では数字 が循環しな L 、から O. 10011000111000011110000011111000000111111... _..."--t--''''--r-'~'''--v--''一一一一

¥ / ¥ / ¥ / ¥ / ¥ /

2 コ 3 コ 4 :z 5:z 6 コ のような意図的に構成したものでな L 、かぎり数字の列の 中にはあらゆるパターンが等しい確率で出現するものと 期待できる. 数を(たとえば、12 と)特定し基数を(たとえば通常 用いている 10 と)指定すれば数字の列は(現実に計算で 表 1 は.f2の値の小数部分における各数字の出現頻度 を 10進 3000万桁までの範囲でまとめたものである .χ2 値 は 3.19 (1 000万桁まで)を除けばすべて 5%から 95% の 範囲に分布しており乱数的であるといってよい.しかし 基数を 10以外の任意の数に選んでも同じことが L 、える保 証はない. 、12を正則違分数で表現すると，

v

'

2=[1 : 2

,

2

,

2

,…]

とすべての部分商が 2 である. 、/すにかぎらず 2 次の無 理数はすべて正則違分数の部分商が循環することが知ら れており，はっきりした規則性をもっている.連分数表 現における規則性と数値表現における数字分布との関係 は不明であるが，明白な規則性をもっ数の方が数値表現 の『かき混ぜ方』が不足していると言えそうである. 、/互のような単純な無理数よりも高尚な数として思い 浮かぶのは超越数である.世界初の電子計算機 ENIAC によって 1949年円周率 π と自然対数の底 e の値がそれぞ れ 2， 000桁以上計算されたのも {Reitwiesner'50} 超越 数の乱数性に着目した J.von Neumann の示唆による ものであった.ところで同じ超越数でも自然対数の底は 正則連分数で表わすと e =[2 : 1 ， 2， 1 ， 1 ， 4， 1 ， 1 ，.・ .2n ， I ， I ， 2{n+1} ， I ， 1 …] とごく単純な規則性を示している.これに対して円周率 の方は きるかは別として)直ちに無限の先まで確定するわけで， と全く不規則である. 司 EE 」 q 4

,

I

,

l η4

,

a 守 ー

,

ー

,

ヨ J

I

,

の4

,

1 1

,

l

,

n 4 n y η ， h

I

,

F吋，， l

,

マ 4 司、， F i L 一一 π その中の数字が本当の意味で独立であろうはずはない. 数値表現の聞に相関がない場合に部分商 n が出現する しかしまた一方で数字の列の中の特定のパターンがある 理論的な確率は log{(n+1}2/n{n+2}}/log 2 に等し 特定のパターンをどこか別の所定の箇所に生成する確率 いが円周率の部分商の分布は実際にこの値に近いことが が高いと L 寸理論的な根拠もない.したがって数字の列 確かめられている.このことからも円周率は真の乱数の は分布としては独立とみなせよう. 候補または乱数の原器として有力であるといえよう.後 に述べるように円周率の値は古典的な統計検定や筆者の 提唱するポテンシャル検定の結果においても模範的な乱 みよしかずのり 工学院大学電子工学科 〒 160 新宿区西新宿 1-24-2 1991 年 12 月号 数の性質を示している. {5) 5

7

3

© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

円唱.帥-r唱伺@rpoon NJ・dFOG的FO納付 前向.AWFnJWNPOON ト唱.aF唱ONFO凶F AFF.円AWAwmooo-F F附.的喝oe。。@ AW0.暗唱ト"。。的 暗唱.0F岨円@曲HOOA 噌 nJw.qFA司@かかかF 噛噛.唱-rFNPoop hFN.国400om ""・唱FFOON 岡田.寸前噌'aF F咽.噛唱。。帥 ap嗣.喧A可FON F4.A 噌 AVFOF 円ト.AF@Fm 寸0.0F‘FAFF 国内.曲N' 岨 0. 寸 M問的 OF.噌ON O唱.AF@ N N a F

円周率の計算

円周率の値を高精度に計算することは乱数とは離れ でも計算機の性能の検証，誇示の目的で ENIAC 以 来行なわれてきた.手計算の時代から記録レースが続 いていたことからも明らかなように円周率には数の世 界の中で独特の地位が与えられている. 100-200桁程 度なら円周率の値を暗唱している人は多いが e や 、12:を 100 桁といわず50桁で、も暗唱している人が何人

2

.

唱NoaFAF'N A W N H A F h F 4 N 暗唱唱曲AFAFF OFト固AF4F N帥ぜ島AFAV nN帥AFAFト 。。唱品かA司 @。噛aF品副円 筒 hF創刊00創刊 F寸400-F 帥トかか4 N寸。ON ド000F FドhF唱 節的 ON トNOV Fd 的 内 NN 怜 FF 。" L F V N F

•

N円N唱AFAFN Nト4hFAF白可制 的。帆‘Fhr‘FF @唱曲円崎酔AWF AF帥maFかか 044AFAFト 崎寸。。。怖 か制寸ooaq eト晶FAFAFF 目。。。。F 晶 NFO帥 NAF 串 aFF ω000F O A F A F A W AFJ 喝_AFF J A w h v 寸国寸 OON 40F 唱曲刷 @“岡 崎 F ド いるであろうか. 手計算『こせよ計算機によるにせよ最近に至るまで円 周率を計算する手法は逆正接の加法公式を利用するも のであった. すなわち 1671 年に James Gregory の 求めた級数展開 NaFF向。。岨円ド0トトavaFNG唱曲HM問。。門前トFFOO円FAFトドaFAFN 国制品創刊O嗣制的pap噛aFJWNJ円F向。副削刊トF嗣00凶N帥NN固AFJWN ドO曲HNOONN帥FトかかFAW噛制NOON唱hFA 噌OOONAF円陣AF--F NhFFFO的F@AFO@島4FmNトFO凶F門的FOO的FNAV 闘争 AV4F 唱寸NF。。F由納nAFAFAF唱NF000F寸ト噛AVAFAF由O制00OF O帥唱F。@帥附句hFAFト唱@AFAFAFトmn4aFAF作曲叫円00@ 門跡的F。凶@FMHAFAF4"。F。。闘的AFAFAFAFJW噛寸前。。的 トドOF。喧島崎凶nhFAFnmA可AFAF島問的帥白色伊跡的曲Nト。。寸 附4@。。倒的ト。。。N@的問。。創刊Oト@酢aFFFトNOON 唱舗@AFか凶帥-F。。FA可制。。OFFAFF。。F唱円寸。OF 的NAFAFe唱N凶円。的帥的。。的@凶円F。mF帥00 自 前噌酔AFFN曲。。倒的mNOω帥鍋maFF唱POON AF向島aFNO。。FOOFOF@mOOF唱ド崎酔 nhphF喧OドaFAW"噌F帥刷ドドhF寸的‘FAFAW F。。NNOON拘凶。Nm唱AFF帥00制 ""。-rFOOF唱F。FoahF400F NN凶作om。かぜoomNhF寸 kFO則的FNOかF4 め F唱FN OAF40F。。FNahFOF 帥噌ト凶暗唱血FA可血FA司 -FNNaFFhvF肉付 。，トaFFF@

i

i

m 唱曲Hω

‘'

ム Y­ n 一 。‘一 -E 且 Z 一 一+ n 一 、 lLn

斗一

2 a u 句ヤハ】一← 一一

z

n

a

c

r

A

を 1706年と公表された Machin の公式 AF創刊ト島酔AVN 円ト国OO的N nAW円000N 帥AFNOO凶P Aw-「 F000F 曲目aFAVAFト ト門NAFAF寸 FOト@酔拘 前回噛@AFF 寸制aF聞か 00帥AF4 qト色FhFF 唱OFOF A可申 O 的 NOON 凶00F 由崎寸 o a F F e a F A 司A 司 』 FF h p

.

.

倒閣端隔週QMW緩Q

同、

"

母島由000門 N ト COO 帥 N noooooN 円唱唱。。凶F ト跡。かかか AV由凶naVAFap aFト噌血VAFAW N e a a F h F " 噛内唱‘FAFF AVF固かか かか伺hF4 o n a h F F AFm ・­ F n h F q ?ド血 FF N 凶酔 @凶喧 @‘ FF @。 F A司的 FM 問 。 F 。 。。000000向。」F 0000000帥NOLF a0000000制OLF 。。。0000帥FOLF 00000000FOト 白000000由O」F 。。。0000凶OLF 0000000寸OLF aooooooNO 』 F aoooo。OFOLF a00000的O 』 F O0000ONOLF 00000OFOト 白0000的O」F 。。。。ONO 』F aooooFO 』「 。。。。的OLF 。。OONOLF a。。。FO 』F aoo帥O 』 F aoO 剣 O 』F OO。FOLF

i

t

-」F 岡田 Ha のような π/4 を与える逆正接加法公式と組み合せて 計算するものであった.

Abraham

Sharp は Gregory

級数に x=

1 /

JTを代入して 1699年に72桁を計算し

たが，彼に Gregory の級数展開を利用して円周率を 計算するよう示唆したのはハレ一等星で有名な英国王 立天文台長 Edmund Halley である.また上記の巧 妙な逆正援の加法公式を発見し自身でも 100 桁を計算 した John Machin は London 大学の天文学の教授 であった.別段天文学で特に高精度の円周率の値を必 要としているわけではない.円周率のロマンと天文学 者のロマンとに何らかの接点があったのであろうか. 逆正接加法公式を利用する場合，計算する桁数を 2 倍にすると必要な演算行程はその自乗の 4 倍になり計 算時間も 4 倍かかる.これは手計算でも機械計算でも 事情は同じであるが，計算機では桁数を増すと多倍長 計算の基数をより大きくとる必要が生じるので計算時 間は自乗よりもう少し余分にかかる. ところが近年のスーパーコンピュータではベクトル の高速演算が可能になると同時に実装される記憶装置 の容量も飛躍的に増大し，これにより高速フーリエ変 換を用いて多倍長の乗算が高速に行なえるようになっ た. Gauss-Legendre の公式によれば円周率は次のよ うに得られる. l!'

=

4

Arctan..L -Arctan~

4 - -

5

-

-

-

-

2

3

9

t寸，

xo=

1

an=l. bn

= '

-

.

ー、(2 '

5

7

4

(6)

表 2 円周率計算の歴史(最近のものは代表的なもののみ掲載)

c

a

l

c

u

l

a

t

e

d

b

y

m

a

c

h

i

n

e

u

s

e

d

d

a

t

e

p

r

e

c

i

s

i

o

n

ti回e

f

O

r

l

D

u

l

a

R

e

i

t

w

i

e

s

n

e

r

e

t

a

1

E

N

I

A

C

1

9

4

9

2

0

3

7

7

0

h

M

a

c

h

i

n

Nicho1son 量 Jeene1

N

O

R

C

1

9

5

4

3

0

9

2

1

3

m

M

a

c

h

i

n

F

e

1

t

o

n

P

e

g

a

s

u

s

1

9

5

7

7

4

8

0

3

3

h

K

l

i

n

g

e

n

s

t

i

e

r

n

a

G

e

n

u

y

s

1

8

1

0

1

7

0

4

1

9

5

8

1

0

0

0

0

lh4 0圃 Machin

S

h

a

n

k

s

L 胃rench

1

8

1

0

1

7

0

9

0

1

9

6

1

1

0

0

2

6

5

8

h4

3

m

Stoer田er

G

u

i

l

l

o

u

d

L

F

i

l

i

a

t

r

e

1

8

1

0

1

7

0

3

0

1

9

6

6

2

5

0

0

0

0

4

1

h

5

5

.

G

a

u

s

s

Gui l1 0ud 量 Oicha闘pt

C

O

C

6

6

0

0

1

9

6

7

5

0

0

0

0

0

2

8

h

l

0

m

G

a

u

s

s

Gui l1 0ud 量 80uyer

C

O

C

7

6

0

0

1

9

7

3

1

0

0

0

0

0

0

2

3

h

1

8

m

G

a

u

s

s

M

i

y

o

s

h

i

L

K

a

n

a

d

a

F

A

C

O

M

1

0

1

2

0

0

1

9

8

1

2

0

0

0

0

3

6

1

3

7

h

K

l

i

n

g

e

n

s

t

i

e

r

n

a

Ta岡ura

L

K

a

n

a

d

a

I

I

l

T

A

C

1

0

1

2

8

0

1

1

1

9

8

2

4

1

9

4

2

8

8

2

h

2

1

m

G

a

u

s

s

-

L

e

g

e

n

d

r

e

U

s

h

i

r

o

L

K

a

n

a

d

a

I

¥

I

T

A

C

S

8

1

0

/

2

0

1

9

8

3

1

0

0

1

3

3

9

5

2

4

h

G

a

u

s

s

-

L

e

g

e

n

d

r

e

G

o

s

p

e

r

S

y

m

b

o

l

i

c

s

3

6

7

0

1

9

8

5

1

7

0

0

0

0

0

0

6

m

o

n

t

h

R

a

m

u

n

u

j

a

n

8

a

i

l

e

y

C

R

A

Y

-

2

1

9

8

6

2

9

3

6

0

1

1

1

2

8

h

8

o

r

w

e

i

n

-

8

o

r

w

e

i

n

K

a

n

a

d

a

e

t

a

1

N

E

C

S

X

-

2

1

9

8

7

1

3

4

2

1

7

7

2

8

3

7

h

G

a

u

s

s

-

L

e

g

e

n

d

r

e

K

a

n

a

d

a

e

t

a

l

H

I

T

A

C

S

8

2

0

/

8

0

E

1

9

8

9

1

0

7

3

7

4

1

7

9

9

9

9

h

5

5

m

G

a

u

s

s

-

L

e

g

e

n

d

r

e

!γ

内旬

ω

…+1♂

1戸=空主守乎乎和

f子子b/c主

tlc+叫l=t九k 一 X/c (何ak 一 ak +1)戸2 Xk+l=2xk (

k

;

;

;

;

0

)

このとき出竺ぜ→ π

守'k この手続きは Newton 法のように 2 次の収束をする ので求める桁数を 2 倍にするには途中の計算もまた倍の 精度で行なう必要はあるがループをたった 1 回余分に回 るだけでよい.したがって多倍長の乗算が主記憶上で実 行できるだけの記憶装置の容量さえあれば計算する桁数 を増やすことに障害はない. 1982年以来次々に記録が更 新され{金田他 '83， '84等)現在では 10億桁以上が求め られている. 桁数の多い計算には不向きであるが円周率を与える異 色な公式として Ramanujan により発見された次の関係 がある.

よ=五互ら(物)

!

(I

10H26390n)

π 兜Ol ，j叫 (n!

)

'

39がn 数学的な輿味からこの公式自体を検証する目的で 1985 年には 1700万桁の計算も行なわれている.以上をまとめ たものが衰 2 である. 10億桁の数字列を数値計算での適当な精度，たとえば 8 桁ずつに分割すれば 1 億2500万個の乱数標本が得られ ることになる.これだけではスーパーコンピュータを用 いた大規模なシミュレーションには乱数の個数としてま だ不足である.しかし乗算合同法による擬似乱数などと 1991 年 12 月号 異なり円周率の各桁の聞には軽重の差はなく相関も認め られないことから分割する箇所をずらしたりは倍)全 体を素数倍して増殖することができる.これにより円周 率の数字列から実用的な個数の乱数を供給することも十 分可能である.

3

.

古典的統計検定

円周率の値がどの程度乱数らしいのか，特に擬似乱数 と比べてどのくらい質的に異なるのかは興味深い.超越 数の数字列は周期を持たないので擬似乱数発生器の検定 に用いられるスベクトル検定 (Coveyou他 '67) を適用 することはできない.ここでは古典的な検定として 1 万 桁の π の値の検定 (Pathria'62) と同様に数字の頻度， 系列相関， ラ個ずつ区切ったポーカ一手の検定を行なっ た. 表 l と何様にして円周率の小数部分における各数字の 出現頻度を 10進 10億桁までの範囲でまとめたものが表 3 である. 32 レンジのうち f 値がその期待値である 9を越 えるものはわずか 2 個しかなく全体に小さいほうに偏っ ており、/互に比べてややお行儀が良すぎるきらいはあ る.近年の円周率記録レースの発端となった 200 万桁ま ででは偶然 f 値が期待値に一致している. 表 41土 5 億桁， 10億桁の範囲それぞれを 10 のプロック に分け全体および各プロック内でそれぞれ数字の頻度， 系列相関， 5 個ずつ区切ったポーカ一手検定を行ない f 値をまとめたものである.頻度検定(自由度 9 )では全

(7)

5

1

5

© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

回.『圃{伺} 円周率の数字の出現頻度 表 3 gnHvnHvn 尚 νan 岬aaaznd'ntuntuphvnwu'EL-n 兎 U''EbnHU 内_J “自民 vnMUnAMm ，・ n ， zo 内uwaA3n ，， vhum ，・ AHu--RV"hdv 内 HMm 尻“ ntu 内 JUe ワ台国内 'unKMenv 円， E冒命‘ u'n ， s 内ベ und-仇 KMAHuntuwn ，.， hdwAHM の河口 WOAun''" ， E 内 Huv'EAn 〆 un ぺ H''Ba-A 叫_unλM 合唱 H・ 4EAphd' 十ゐ aaEnhun 司 uw-、ぺ四・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・ ・ 8az 伶白内 Daazaaznuqu"t 戸 Baa-nf 巧 tRunB 守 'nlnonLS 唖 aaτEuaaznbRuwnfaazaa-qu'nt7 ・内 ba4 ・

-E

-nB-aAτna 。 -ahun41&a 唖内 unuq ，M'i1 ‘ FO--aaτauznonua-nundqu 肉 aa 唖 'in4 内 4n'U 噌 iaonU 向。­ 帽‘ aaa 向 ZUF 内 UAnu--EE'avu ・ s ・'の yu ，‘， .-AHνaa--AUd-AHWA 叫 U 内 wuwpnu 噌 EAam 唱の令"の 4HV の yu の udv の_5u 曹の屯 HV の xu'Ea' 戸 hvORun ，・内 hv 戸川 un''­ -・・ 2_・ azup 同 uv 肉 Hu'hHvnHvn 判 u ・ mmu"'' 唱， a・ "hvnnu 内 Hun 点 H-nHunJ “ nhu'hdw 唱 a--"''nHuv 陶 hunHVP 同 dvnruphunAun ，・向 Ju----aa ，“ F同“曹の WW 《 wuw 《叫 unHW 内 HVAHV 内HU 内 HURHWA 吋uwnwu の wu'AHun 叫 M'SAn 吋 U の_FunHV 向 HMOAUORuw 町九日 'nHv-''A-an 宅内 HVRHunHvnHMAHvnHU 曹の wun 叫 U 内 WMAHV 内 wunHvnwu 白 HvnHVAHvn 叫 un 叫_unHvnHu-­ 剛唱 EA-n'uaa--Fhuaku のHVaa--n 叫 ua “‘ nHvn 吋 uvnHvnuunHvnHvnHV 内 WH 曹の叫 ωnHunHV­ -E. ‘唱 a 且.噌 EA 倫，“ a5uaaEaAvnwuwphHV 白HunHun"" の曽 WAHV 内 Hu-­ 'E& 。，“ nAun5uaaτhnunHu-­

.,

A-8-25312282086157308450924806735239"

司 -EE ・内 J“ F_同 unHunHUAwuaaz 戸 hvn ‘un ， zphu ・ 24nkua 仏 Ea 民 unRU 円九 uw 噌 EAh5u 曹の Fuamuauun 叫_uwhHU 自民 uw'az-内 HuwaRU 内 5uwnAun 九日曹の丸 uv­ 帽唱'且 'aJuaa 宮内 uvauunHM 内 uunwunMM'n 明日 'aa 官宅 EAOAuwβhH-a 長 u'aa 宮崎， z-EE--hHwaa--の 4U4EA-aa ‘‘‘ uvnhu'"' ・ βhu 内対 wamu­ --aA'F 内 U 《wunwun" “内 WU 《吋 ωnHU 曹の WUAHM 向叫 eaHMm 叫uwn 叫 uaud ・ -AUAH" の_FU'Ea. 。，“ 45wp 同 d'Fhv 内 ZM'aa-­ .， 4aa ・ n 国 nBAHvnEnuquaHV 内百 nEnEnBnunununU 内 Ununununu-1 ・ aa 骨 an--且 MqψEun 司 4aqunB 内 UOW 内 U 《 U 《 uunwnU 内 Unu­ -a 咽_Aqda4a 峰。 onuEU 内 unununuaU 《 u­ '41 ・・。 49Ava 宅 PD 。 onu­ -EA 田 ' n ，‘， nAun ，“内 hvph “， ng-nkunHW 内 Avn'''hdvnAvnHuv 自 H 》内 HVAHuaHwaHU"5 ・ aa--aaτnHuphup 内 un 叫_unHuw の λMn ，， .EA-aa--内 HU 《 wuv'Ea-­ 自 -EA 偽点 uw 内 WHwn 河 u-aκun ， zphun ，目。，“拘 hun ぺ u ・向 Hu--aaap 同U 自民 u'mhuRHUmnu 内 Huw 内 WU 円， enwuwaλv ・ a ・ -nhvnnvn ， E'Ea-a ， UAHV 戸 hu­ -'EA-a&Z 内対日 wn 刊 HwnwunHM ・ E ・-内 JU 白民 un ‘u 仰険 uwa 翼 UW 崎 U 向 2“ F 同 uaaτ 内 J “唱 E--anunAuvpnvaaτn'unHunAV 内魚 HW''A-nHu-­ .‘，‘ aaznHunHuwAHuw 内 wunHuwnHunH 切除 HV 内 Huh-unRUnAuwm 民 uwnkun 叫 M'n 叫 un 輔 wnHvn 喝 UAHHvnHM' 同 unhu-唱_B-bnr “伊川 ω 内 wuvAHV 内 HUAHUAHuvnHvm 叫 umwuvm 叫 unwun 河u'n 叫 dv 内 wuvn 吋 UAHV 内叫 HV 内 HvnHVA 叫 HW 内 wdve --の，“ aax ・ "hda 民 MAHuvaa 宮内切 uaa 否仇切 ωAWU 《叫 dv 内切 dw 内切 dw 内 HWAMd 向 HWAMu-《 wu'Aud ・ --'E 且・・ Ea. ， EA-。ノ M 向_F “ aaz 内， e 内司 uaaa-nHvnwunHU 曹の Hunwun" “自 'lA' の，“。，“ aaa 開 phu"'' 内

MU

-ｭ

tu-naEu 。 Daa 官 nuqdwt ム nI 。 equooa 峰。 oqu 斗 λ て lQ 凶て UA ・唱 4 ・ 14

体に χq砲が小さく， 2 項系列相関(自由度如)や ポーカー手検定(自由度 6 }では x2_{値のばらつき} が大きい.この傾向は2∞万桁， 8∞万桁の検定で も見られたものである.

4

.

ポテンシャル検定

表 4 円周率の古典統計の f 値

(

)内の数値は上側確率 ふ v

-,

nu -n 刊、，， p a 6 -L H = -F a g -a u --L 医局 Ju -nu

,

zt -n v ETIlt -B I ' F n u -n w u F'BA2 F a ｭ -A P T --FA ・ ' e J u -eu ，，、 i ' E 1 3 ' a ' -vJ 、，， -c n B

-nE

-aLV ・

-u

e

A

-nu 、・

-e

d

-Ta ， s 、 -S A ETZEE--e

,.

FtE ・

-c

b

-o m --aAHU

-b

n

+,

E'14 t- ーーーー+ーーー句ーー ---t ーーーーーーーーー『ーーー+ーーー---ーーー-ーーー+

Iwho1el 7.42

1

1

0

1

.

8

1

8.98

G

.

1

7

1

86.42

2.04(>.90) 1

2 1 1

1

.

3

0

1

1

0

1

.

3

6

1 1

1

.

05(<.10) I

3

2.08 (

>

.

9

9

)

1

7

1

.

6

8

(

>

.

9

0

)

1

1

.

6

7

(

>

.

9

0

)

1

4

G.19

1

94.78

118.91(<.005)1

5

5.11

1113.54(<.05} 1 0.42(>.995) 1

6

8.50

1

86.51

3.98

7

4.10

1

92.95

1

.

1

2

(

>

.

9

0

)

1

8

5.63

1102.50

6.13

9 1 13.02

1

1

3

1

.

4

1

(

<

.

0

0

5

)

1 5.86

1 1

0

1 12.68

1

1

1

1

.

24(<.10) 1 5.93

シミュレーションを行なううえで常に真の乱数 が求められていると L 、うわけではなく，それぞれ の目的に応じた検定を行なって棄却されなければ 合格とするのでかまわない. 擬似乱数の格子構造 (Marsaglia '68) は多く のシミュレーションで有害と考えられる.これを 検出する理論的な検定法としてはスベクトル検定 法があるが，前述のように超越数の数字の列のよ うな周期を持たない標本には適用できないため阿 ←一一←一一一一一ート一一一一一ート一一一一一一+

周率を疑似乱数と直接比較することができない (1)

First

500 冊 il1ion

decima1 digits

筆者の提唱したポテンシャル検定法(三三好 '83， '87)は乱数列を用いて同等な粒子を 3 次元空間内 に配置するとき，粒子系の 2 体相関によるポテン シャルエネルギーが格子構造を鋭敏に検出するこ とならびに近傍および途方の寄与が均衡している ことを用いた統計的検定法であり，乱数発生器の 周期と無関係に行なえる利点がある. N個の同等な粒子を 1 辺2Lの単位立方体の内部 に配置しこの立方体が無限に周期構造をなすもの とするとき，粒子系の 2 体相関によるポテンシャ ルエネルギーは単位立方体当たり，無次元化して 次のように書かれる. ここで Xi は i 番目の粒子 の 3 次元座標を意味する. Lr ・ 4 ・ m m 司 EFIz' F-nHM 向HU -J u -A n 3 -nH 、.，，--- -aE-く〉 側冒 HH---ft ，・ 1 ---nqd 勾，・ nHU'EAnHU' ・ An 拠υnqJV 勾 'h 勾，， mud -rrA-FUrBn303pbqa 内 dnEaqa4qa -alv ・---'ua ・ du-qtuvphunHV-aA 勾'uFnuvan 岨zan 噌 phdq'han 唱 -nυ ，，‘、 -'aA -n y --T ・ 'lEEl-?BEa--EISE--laaBEt ・ 'EE--' ・・ l'lIBEt--'l 問、，』'­ m n H u w ｭ 『 nwd ・-m-EAm 一四 an 守内‘ uwn4u'EA 内JLHawd 勾， en ，・ nHvnudwO 拠υ 』笥 u ・ -FhAw--AOAunhuFhd 勾，.句、 u'gA 勾''nHu-EA ---AF'A 一---P ・A ・ -nqunHuan 崎zphd ・， A 司 ''q ‘ uwphd 向 wdnwd 。， h -eJu-勾 tod 買 U7sRuqdaonuqJvaonu -B I t -4IBI--+l'l

,

l'IB--El'IB--lil--E --‘.，，、，.，、，，， FFnHunHUF 内 JU --nwun 吋dnwJW -vJ 、

.

.

,---.

-cnB-、ノ、ノ、ノ -nH=-''t't 、，，、 司自 FU--内，ゐ勾'・ nnu 内 'uan ーの JuphJW 内 JU 内tuFhunA リ -uuFT ・ --nwd 内 JugF 内 dvan 噌白瓜 u"hvphJvan 噌唱 EAaAunwJW 』伺 H2 ・---­ -aUAM-aaa7 ・ qaqJnnaqnORunE? “司 t -FA ， E 、 -F r a -4 ・・ Bt'ill ・ E21'lss'IB--zBI'l ，， a11SIZE--1 ・・ 'EE

-K

-e

-c t o -i -nvm-o'I 勾_ιqJaqRUFO マ toon3nu -'BAHU--nM ・ E ・--・ nunH-司" lt'll-TIll--Illit--IIllit---tl

U

N

=

ー+ーー---ーーーー町ーー+ーーーーーーーーーーーーー+ーーーーーーーーーーーーー+ 3

r

r{P(x-xtl ー 1H ，E ð(x'-xJ} ー 1}

¥

¥

-

-

.

i

j _-_ _

-dxdx'

(

2

)

First 1 bi11ion decimal digits

4 1t"~~

Ix-x'l

(v=(2L)8=N) で与えられるので，厳密にはN(Nー 1)/2個の I/r は 一 3 ~明 I 6log(2+ イす)ー π九 一五七れ雨戸可T

-

4L } 距離の逆数の単位立方体内での期待値は (61og

(2+

、IT- 1t" )/4L に等しく，分散は

2(1)Ji三

_7-Fjω+

=Arctan----l

_{1---.; z2+ 1}

一戸

log(2+.;T-1t"

l

y

キ O.

5

0

2

3

3

8

9

9

x

ら

1991 年12月号 互いに独立ではないが粒子数Nを大きくしたときUNは 漸近的に正規分布N (0， σ刊にしたがうとしてよい. ここで、_σ02_，主

σ0

2 =

(J-

_Y

担

N-I)

_σ2(

_斗与

_(1.

_9143xL叩と

_l

、_lπf 辺、rI 1V である. U

_N

を計算するとき，粒子簡の距離 !x;-xJ! は 着目した粒子を中心とする単位立方体内部で考える. 粒子数Nを 108

=

₁₀₀₀

_,

₁

₆

8

=

₄₀₉₆

_,

₂

₅

8

=

_{15625の 3 通} りに選び円周率 5 億桁の小数部の数字を上位から8桁ず つ区切って粒子の座標とし，それぞれ 17500， 41∞，

5

0

0

(9)

₅

₇

₇

© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

個の標本を得た .(24Nx 標本数が 5 億 桁に達していないのはディスク上に保 存できるデータセット容量の上限値と CPU時間の関係である.

)

標本全体とこれを 10個に分割したプ ロッグの各々について，標本平均

U

=μ (U

N

) および標本標準偏差 s.d.= σ (U

_N

) を検定した結果を表 5 に掲げ る.表で z は標本平均÷標本標準偏差 ×イ房王子才ヌ， prob. は標準正規分 布の確率密度関数を Ixl から∞まで積 分したものの 2 倍， chi.sq. は標本分 散 ÷σ。2X 標本サイズである. 乱数が少数の超平面上に落ちる格子 構造があると U

N

は負の値をとり，ま た粒子がすべて相異なる格子点上に配 置されると U_Nは正の値をとる.後者 はU

_N

の患大値を与えNが奇数の 3乗 のときはほぼ0.53N また偶数の 3 乗の ときはほぼ0.90N となる， 粒子数N を大きくするほど格子構造 の検出能力は高くなるが， 32 ビット機 で悪名高い乗数 65539 による乗算合同 法擬似乱数から粒子系を構成すると N =1000でも μ (UN) がー 190にもなり正 規分布N (0， σ ぷ)からほど遠いこと がはっきりする. 5. まとめ 円周率の数字列は種々の検定で良好 な結果を示している.また厳密な検定 ではないが円周率の値を整数部の 3 も 含めて頭から I1贋に数字をとって整数を {乍るときこれが素数になるのは 160 桁 までの範囲では 1 桁( 3)

,

2 桁 (31

),

6 桁 (314159) ， 38桁(略) だけであ る.これも素数定理から導かれる素数 の分布とよく合致している. 現時点では，乱数の原器として 1 つ だけ選ぶとすれば類似品がないこと， 再現性があること，などから円周率で あろう.凝似乱数や物理乱数との比較 表 S 円周率のポテンシャル検定 +1+ItEI--l

,

EIF+ --曲中- -・『*-FqF5574q46??ZTｭ ---・・・・・・ E ---・・・・ヲ 'RV-J-avhu 『， h4 ・ av 噌・" -・・-反 uvau--。， av' 句句'』 'e-・ τdAv- -h-56778777777-

eE71111111111

-E

t+ ・・ +It--lit--' ・ 11 ・ T -? ? ? ? 3 t 3 5 1 7 0 - 03-・ -88346844 ・ 360-向VRV 圃 d"---0 ・ 4---7478 守 778778-tt 17 ・ -s-44444644444 目・・ 4 ・申 -FVu' -+Ea'It---'ttItl+ ・・ ----xχxzχ%xvhχχx-tt = = e -b -5 8 5 5 8 5 4 7 7 7 2 -。-。 -34713216754-tt 'r-r-・・・・・・・・・・・ -nn -n p -p -4 ? 4 8 0 0 7 5 0 6 6 -ｭ t o -7 4 8 8 e V 6 5 4 2 -e c e--+I+ 』 111111IBI--?" ・ r ltr--35107184455 甲・・ t ・ oa--7?8?1054296-pp vEO-目 27Z4 守 534621-筒 VE--x-36818204671 一 05 p ・ 7t-E ・・・・・・・・・・・ -t dt1--00102000oo--ｭ t g -t t ndt+1+Itliti--111+aaｭ -r E d -5 2 6 3 3 3 6 4 2 2 9 -d -F 3 -7 8 8 4 6 O z -Z 7 4 -t t nd't-u-45?76?65230-nn -n-a-t-81020501?13-aa paLm-n-16773841438-ee -tP1-a-170T42057Ez--t dgme 目・-・・・・・・・・・・・ -ff n a -m -0 0 2 0 3 0 0 0 o o -1 1 1L ・ 4-2FFFFFFnn 絹 -Y--gg fcfO+1+11111111tlt+11 。 l 。 o-F 』・・ t o e k -r ・ ro-e-t-gg -r ・ o-。-。 12345678?OF* b 。 bo-L-ht-** m -m O M b -H ｭ u h U 2 -4 NTN4+I+11111111111+ Numb ・ r ()f 1 nd ・ p ・ nd ・ nt P..rtlel ・..

=

40守晶

Th oratleal .tand ・ rd d ・ vlatlon

=

122.5016

Numb ・ 2・。 f "ampl ・.

"

'

4100

1.03046400 d ・ el 同 aL dl 冒 It.. or ・ þroe ・ g ・・ d

+日ーーーーーー+ーーーーーーーーー -iトーーーーーー -iトーーーーーー+ー---ーー+-ーーー『ーー『同+

1 bloek Iyl 岡・ an(U) x 1 prob.1 ...d. 1 ehl- ・ q. I

+日ーーーーー -iトー--ー』ーーーーー+ー-ー』ーーー+ー---ーー+ーー-ーーー+ーーーーーーーーー+ 1 whol ・ 0.6797931 0.3493172.69Y.lf24.621 4243.1 1 -4.0832751-0.6590150.99 Y. 1125.461 与 30. 唱 E 2 四 1.4099431-0.2202182.57 Y. lf29.661 459.3 牢 1 3 2.7264651 O.4365166.24Y.1126.461 437.0 4 1 -6.79709fl-l.0915127.51 Y. lf2ι.091 434.4 5 0.9143281 0.1494188.13Y.1123.951 41 守 .8 6 4.0328581 0.6951148.70XI117.481 377.0 7 112.1.7564912.13051 3.31Y.llf8.571 384.1 8 1 -4.0058141-0.6211153.45Y.1130.591 465. 宇牢 I 9 1 -0.2035591-0.0336197.32Y.1122.551 410

,

3 10 3.1483081 0.5098161.02Y.1125.041 427.2 +ーーー四ーーー+ー--ーー『ーーーー+-司ーーーー田+ーー四ー田町+ーーー--ー+ーー皿ーーーーーー+

*

I ・ Ignlfleant at 10 p ・ re ・ nt l ・ v ・ t Numb..r of I nd ・ þ"nd ・ nt partleL ・..

=

15625

Th..or ・ tle..l standard d ・ vlatlon

=

299.1031

Numb ・ r 01 S8刷 pl ・.

=

500

187500000 d ・ e I 同 al dlglt. ar ・ proe ・..・ d

+ー一ーーーーー+ーーーーーーー ---iトーーーーー』ー+由ーーー-ー+ーーーー』ー+ーーー--ーーーー+

I bloek Iy' 冊。 an(U) x 1 prob.1 ...d. 1 ehl-sq.

+--ーーーー-+ー---ー+四ー四『ーーー+目『ーーーー+ーーーーーー+ーーーーー-ー--+ 1 whol ・ 1-24.6270521-1.82641 6.78Y.1301.511 508.1 1 -8.8799021-0.2022183.98Y.1310.551 53.9 2 1 39.6500991 0.8727138.28XI321.261 57.7 3 1 11.5831971 0.2646179.13χ1309.571 53.6 4 I-S9.3749251-1.4235115.46Y.I294.931 48.6 5 1 11.4333071 0.2577179.66Y.1313.681 55.0 6 1-32.7605091-0.8154/41.48Y.1284.081 45.1 7 1-16.02626 守 1-0.4501/65.26 Y. 1263.201 44.8 8 1-94.0778071-2.01811 4.36 1. 132 守 .631 60.7 9 1-35.6865421-0.9458134.42Y./266.811 39.8 10 1-60.1291651-1.4221115.50Y.1296.971 50.0 +ーー『ーーーー+ーーーーーーーーーー+司ーーー『ー-+-ーー四ーー+ーー一ーーー+ー-ーーー白骨--+

に利用するのに最適であり，実際のシミュレーションを 多数個の擬似乱数を用いて行なう場合でも円周率の数字 列を用いたテスト結果と照合することで結果の信頼性を 高めることができょう.

参考文献

[

1

] Coveyou

,

R. R. and MacPherson

,

R. D.

Fourier a

n

a

l

y

s

i

s

o

f

uniform random number

generators

,

J

.

A

s

s

.

Comput. Mach.

,

14

,

1

0

0

-

1

1

9

(

1

9

6

7

)

.

[2] Kanada

,

Y.

,

Tamura

,

Y.

,

Yoshino

,

S

.

and

Usihro

,

Y.: C

a

l

c

u

l

a

t

i

o

n

o

f

πto

10

,

013

,

395

decimals based on the Gauss.Legendre a

l

g

o

r

i

ｭ

thms and Gauss arctangent relatian

,

TR84-01

,

Computer Centre

,

the University of Tokyo

(

1

9

8

3

)

[3] Marsaglia

,

G.

Random numbers f

a

l

I

mainly on t

h

e

planes

,

P

r

o

c

.

N

a

t

.

Acad. Sci.

,

81

,

2

5

-

2

8

(1

9

6

8

)

.

[4] Miyoshi

,

K. and Nakayama

,

K.: Calculaｭ

t

i

o

n

ofπto

2

,

000

,

000 decimals

,

P

r

o

c

.

Annual

Conference o

f

the I

n

f

o

.

P

r

o

c

.

S

o

c

.

Japan.

,

23

,

9

4

1

-

9

4

2

(

1

9

8

1

)

[ 5 ] 三好和憲:ポテンシャル法による古L数の検定.数 理解析研究所講究録，

498

,

1

9

1

-

1

9

8

(

1

9

8

3

)

[

6

] Miyoshi

,

K.: On a S

t

a

t

i

s

t

i

c

a

l

Test Based

on t

h

e

P

o

t

e

n

t

i

a

l

Energy f

o

r

t

h

e

System o

f

Non.Periodic Random Number Cenerators

,

Tensor

,

N. S

.

48

,

2

3

6

-

2

3

9

(1

9

87

)

[

7] Pathria

,

R.

K

.

:

A s

t

a

t

i

s

t

i

c

a

l

study o

f

randomness among the f

i

r

s

t

10

,

000 d

i

g

i

t

s

o

f

11:,

Math. Comp.

,

16

,

1

8

8

-

1

9

7

(1

9

6

2

)

[8] Reitwiesner

,

G

.

:

An ENIAC determination

o

f

πand

e

t

o

more than 2

0

0

0

decimal p

l

a

c

e

s

.

MTAC

,

4

,

1

1

-

1

5 (

1

9

5

0

)

[9] Tamura

,

Y. and Kanada

,

Y

.

:

C

a

l

c

u

l

a

t

i

o

n

。fπto

4

,

194

,

293 d

ecimals based on t

h

e

Gaｭ

uss-Legendre algorithms

,

TR83-01

,

Computer

Centre

,

the University o

f

Tokyo

(1

9

8

3

)

1991 年 12 月号 11 月号/発売中/定価 930 円

最新ネットワーク

コミュニケーション

一一パソコン LAN 構築の実際一一 日本国 LAN 事情 蓑固定彦 パソコン LAN 入門 吉野益弘・鈴木洋 応用分散システムとは何か 伊吹公夫 インテリジェントネットワークとは 野口正ー 簡単なネットワークでも実用性は大である 青山光伸 箱の中のイーサネット 芝野耕司 〈トレンド〉 オフ.ジェクト指向ソフトウェア技術の l1i新動向 原因賢一・字郵宮公富Q ・深海良彰

対称性とその破れ

物理学における対称・非対称と対称性の破れ 伸滋文 素粒子におけるカ学的対称伎の破れ菅本舗夫 時間のおりかえし点 小野健ー 超伝導，趨流動にお t坊対称性の破れ 大見暫巨 「結晶J ・料学・文化 小川泰 宇宙における対称性の破れ大規模構造 池内 了 時間の矢 ニ間瀬敏史 生物の発生における対羽市Eの破れ 大隅典子・土居洋文 新しい対称性量ヲ常手への入門 上野喜三雄他 圃最新事i 好僻発売中

ザ・ UNIX

戸J" 隼人箸/A5/定価 1751 円 ー価格表示は，税込み価格となっています.

サイエンス社

東京都千代田区神田須田町2-4 安部徳ピル

電話 (ω)3256・1091 (代〉振替東京 7

-2387

(11)

5

7

9

© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.