2001年世界記録への挑戦過程

　　　　　　　　2003/1/10 日立製作所 & 早稲田大学　後保範 ( Ushiro Yasunori )

1. 2001年世界記録兆戦の概要と結果

　主計算はSR8000/MPP 128ノードを使用して、16進1兆桁を計算し、検証計算は計算機利用
上の都合で64ノードを使用して16進5千億桁を計算して、両者の結果を比較した。
私の精度に関するプログラムミスで、3597億余桁以降が不一致となり、2001年の挑戦は
諦めざるを得なかった。

2. 2001年世界記録挑戦の経過

2.1 利用公式

　収束の速い級数公式である Ramanujan型の級数の公式に DRM法を適用して計算した。

2.1.1 主計算に使用した公式

　　円周率計算の級数公式で最も収束が速い(高速)なRamanujan型の級数の公式で
　Chudnovskyが発見した下記公式を使用した。最後に全桁の平方根と除算が必要なため、
　DRM法による級数計算部分と、最終的にπを計算する部分の2段階に分割した。
　(1) DRM法による級数計算
　　　X = 13591409 - 1･1･5/(9･1³･10672303³) ･ { (13591409 + 54514013)
　　　　　　+ 3･7･11/(9･2³･10672303³) ･ { (13591409 + 54514013･2)
　　　　　　- 5･13･17/(9･3³･10672303³) ･ { (13591409 + 54514013･3)
　　　　　　+ 7･19･23/(9･4³･10672303³) ･ { (13591409 + 54514013･4)
　　　　　　- ･･･ }
　(2) πの値の計算
　　　π = 53360･640320^1/2 / X

2.1.2 検証計算に使用した公式

　　Ramanujan型の級数の公式でChudnovsky公式についで高速な下記公式を使用した。
　DRM法による級数計算部分と、最終的にπを計算する部分の2段階に分割した。
　　　X = 1103 + 1･1･3/(2･1³･198⁴) ･ { (1103 + 26390)
　　　　　　+ 3･5･7/(2･2³･198⁴) ･ { (1103 + 26390･2)
　　　　　　+ 5･9･11/(2･3³･198⁴) ･ { (1103 + 26390･3)
　　　　　　+ 7･13･15/(2･4³･198⁴) ･ { (1103 + 26390･4)
　　　　　　+ ･･･ }
　(2) πの値の計算
　　　π = 9801 / (2･2^1/2･ X)

2.2 計算経過

2.2.1 16進1兆桁の主計算

　SR8000/MPP 128ノードを使用して計算した。
　(a)　DRM法による級数計算
　　開始日時： 2001/9/16 04:15
　　終了日時： 2001/9/22 19:02
　　使用時間： 343652 s = 95h 28m
　　中断回数： 4回
　(b)　16進πの値計算
　　開始日時： 2001/10/06 00:39
　　終了日時： 2001/10/06 17:41
　　使用時間： 61296 s = 17h 2m
　　中断回数： 0回
　(c)　合計計算時間
　　343652 + 61296 = 404948 s = 112h 29m

2.2.2 16進5千億桁桁の検証計算

　SR8000/MPP 64ノードを使用して計算した。
　(a) DRM法による級数計算
　　開始日時： 2001/10/9 21:49
　　終了日時： 2001/10/14 19:04
　　使用時間： 354443 s = 98h 27m
　　中断回数： 2回
　(b) 16進πの値計算
　　開始日時： 2001/10/15 01:56
　　終了日時： 2001/10/17 09:43
　　使用時間： 51594 s = 14h 20m
　　中断回数： 1回
　(c) 合計計算時間
　　354443 + 51594 = 406037 s = 112h 47h

3. 2001年の計算結果

3.1 主計算と検証計算の比較結果

　　2001年10月19日に主計算と検証計算の16進の結果比較をおこなったところ、
　359,712,809,450桁以降異なることが判明した。
　この時点では、どちらの計算が誤っているかは不明であった。
　両計算の359,712,809,441桁以降の20桁は下記の通り。
　　
　　主計算(誤）　： 8234700EE3 3B42B4CB6C　
　　検証計算(正）： 8234700EE0 B6CC5EE098

3.2 他の区間公式との比較結果

　　2001年10月31日に別の区間演算公式で16進4千億1桁から20桁を計算し、両方の結果と
　比較したところ、主計算の結果が不正であることが判明した。
　両計算の4千億1桁以降の20桁は下記の通り。
　　
　　主計算(誤）　： A7B86E9407 87D3E3E5BF
　　検証計算(正）： 9B8A5C4793 B384EE4BA0

3.3 1年後判明事項

　　約2ヶ月の机上での不正原因調査の結果2001年12月6日に、主計算、検証計算共に結果
　不正が発生していると予測された。その時、結果不正となる桁数も下記のように20～30桁
　の差で推定できた。
　　主計算の結果不正桁推定値　： 359,712,809,427 + 20～30桁
　　検証計算の結果不正桁推定値： 455,568,772,158 + 20～30桁
　2002/10/4の16進1兆桁の計算により検証計算の結果不正発生桁は
　　455,568,772,187桁と確認(推定値+29桁）できた。
　455,568,772,161以降の40桁は下記の通り。
　　
　　2001年検証計算(誤) ： 5F414578B0 A0685B5232 0EFD98D81D B3BDD2E000
　　2002年世界記録(正) ： 5F414578B0 A0685B5232 0EFD98FFF5 99BD47FF3A