パソコンによる 円 周 率 小 数 点 以 下 5 兆 桁 の 計 算 文 部 科 学 省 科 学 技 術 政 策 研 究 所 シンポジウム 近 未 来 への 招 待 状 ナイスステップな 研 究 者 2010からのメッセージ 文 部 科 学 省 旧 文 部 省 庁 舎 6 階 第 2 講 堂 平 成 23 年 6 月 30 日 旭 松 食 品 株 式 会 社 情 報 システム 課 近 藤 茂 1
円 周 率 とは 円 周 率 = 円 周 直 径 = 3.14 円 周 率 πの 語 源 はギリシア 語 の 周 りを 意 味 する περιφερειαの 頭 文 字 をとったと 思 われる 円 周 率 は 無 理 数 かつ 超 越 数 円 周 率 は 小 学 校 5 年 生 の 算 数 の 教 科 書 に 登 場 する 円 周 率 は 幾 何 学 解 析 学 力 学 統 計 学 などに 使 用 される 公 式 に 多 く 含 まれている 2
円 周 率 を 計 算 から 求 める 初 めて 円 周 率 を 計 算 から 求 めたのは アルキメデスと 言 われる 直 径 1の 円 からその 円 周 はその 円 の 内 接 外 接 する 正 多 角 形 の 周 囲 の 長 さの 間 にあると 考 え 正 96 角 形 まで 計 算 して 10 3 < π < 71 10 3 70 (3.14084 ) (3.14285 ) を 見 出 し 円 周 率 を3.14まで 計 算 した ( 紀 元 前 3 世 紀 ) ルドルフにおいては 正 2^62 多 角 形 まで 計 算 して35 桁 まで 求 めた (1621 年 ) π=3.1415926535 8979323846 2643383279 69399(35 桁 ) 2^62=4,611,686,018,427,387,904 3
手 計 算 による 時 代 17 世 紀 になって ニュートン ライプニッツによる 微 分 積 分 学 が 確 立 されると 無 限 級 数 を 使 用 したπを 求 める 式 が 発 見 される ようになった マチンが 以 下 の 公 式 を 発 見 し 100 桁 まで 求 めた いままでの 公 式 より 収 束 が 早 い -> 早 く 計 算 ができる π = 16 ta 1 4 ta 5 1 1 1 239 シャンクスはこの 式 を 使 用 して707 桁 (1874 年 )まで 求 めた その 後 1945 年 にファーガソンによる 計 算 から 計 算 間 違 いが 指 摘 され527 桁 までしか 正 しくない 事 がわかったが 70 年 以 上 この 値 が 信 用 されていた その 後 コンピュータの 登 場 で 計 算 桁 数 は 飛 躍 的 に 伸 びていった 4
計 算 機 による 時 代 年 代 計 算 者 コンピュータ 計 算 機 を 使 用 桁 数 公 式 計 算 時 間 1949(S24) ライトウィズナー 他 ENIAC( 真 空 管 ) 2,037 マチン 70h 1973(S48) ジャン ギュー 他 CDC7600 1,001,250 ガウス 23h 1985(S60) ウィリアム ゴスパー Symbolics3670 17,526,200 ラマヌジャン - 1989(H01) チュドノフスキー 兄 弟 IBM3090 1,011,196,691 ラマヌジャン 他 - 1997(H09) 金 田 高 橋 他 HITACHI SR2201 51,539,607,552 ルジャンドル 29h 1999(H11) 金 田 高 橋 他 HITACHI SR8000 206,158,430,208 ルジャンドル 37h 2002(H14) 金 田 後 他 HITACHI SR8000/MPP 1,241,177,304,180 高 野 他 432h 2009(H21) 高 橋 T2K( 筑 波 大 学 ) 2,576,980,377,524 ルジャンドル 29h 2009(H21) ファブリスベラール i7-940( 自 作 パソコン) 2,699,999,990,000 チュドノフスキー 2760h 5
円 周 率 の 計 算 のきっかけ 学 生 時 代 (1974 年 )にFACOM230-25 用 の FORTRANを 学 習 する 際 たまたま 目 に 止 まった 雑 誌 それが 学 習 コンピュータ だった この 中 に 円 周 率 を1000 桁 求 めるプログラムが 掲 載 されて おり これを 基 に230で 円 周 率 を 計 算 させた 事 が 円 周 率 の 計 算 のきっかけとなった 以 来 いままで 37 年 間 この 数 字 の 並 びを 追 い 続 けてきた 多 分 これからも 追 い 続 けていくと 思 う 6
なぜ5 兆 桁 になったか 2009 年 4 月 筑 波 大 学 のスーパコンピュータで 記 録 される 計 算 桁 数 は 約 2 兆 5 千 億 桁 2009 年 12 月 Fabriceのパソコンで 記 録 される 計 算 桁 数 は 約 2 兆 7 千 億 桁 2010 年 4 月 に 自 身 が1 兆 桁 の 計 算 に 成 功 最 初 10 兆 桁 を 考 えたが ハードウェアが 足 りない 計 算 時 間 が 半 年 以 上 等 により 断 念 結 果 当 時 の 記 録 の 約 2 倍 の5 兆 桁 を 挑 戦 する 桁 数 と 決 めた 7
5 兆 桁 の 計 算 にこだわった 事 個 人 の 所 有 する1 台 のパソコンを 使 用 して その 持 っている 計 算 能 力 を 最 大 限 利 用 して 計 算 を 行 う 8
5 兆 桁 の 計 算 履 歴 計 算 は2010 年 5 月 4 日 に 開 始 し 91 日 後 の8 月 3 日 に 終 了 し ました その 間 に 計 算 が 停 止 したのは1 回 です 原 因 について は 以 下 の 理 由 により 詳 しくわかりませんでした 別 のマシンで 再 現 できなかった 事 チェックポイントからの 再 計 算 で 同 じ 場 所 を 問 題 なく 通 過 してしまった 事 昨 年 の 猛 暑 にも 耐 えて 全 ての 部 品 が 長 期 間 連 続 稼 働 で 故 障 せずに 動 作 してくれた 事 は 幸 運 だったと 思 っています 9
主 計 算 に 使 用 した 公 式 チュドノフスキーの 公 式 ( 1) (6)!(13591409 + 545140134) π = 426880 10005 3 0 (3!)(!) 640320 3 = =0 : 3.1415926535 8973420766 8453591578 2983407622 3326091570 6590894145 =1 : 3.1415926535 8979323846 2643383587 3506884758 6634599637 4315654905 =2 : 3.1415926535 8979323846 2643383279 5028841971 6767885484 6287912727 =3 : 3.1415926535 8979323846 2643383279 5028841971 6939937510 5820984947 1 を1つ 進 めるごとに 約 14 桁 精 度 が 上 がる を352,568,346,744まで 進 めると 精 度 が5 兆 桁 に 達 する 10
計 算 の 高 速 化 1 大 きな 桁 の 計 算 に 関 しての 高 速 化 について 説 明 します 円 周 率 の 計 算 において 多 桁 x 多 桁 の 乗 算 がキーになります 筆 算 のように 乗 算 を 行 うと 桁 数 の2 乗 に 比 例 した 計 算 時 間 が かかります 計 算 時 間 の 短 縮 にはいろいろな 手 法 がありますが FFTを 利 用 して 計 算 する 事 で 高 速 化 を 図 ります C=AxBの 計 算 の 概 要 を 示 します 1. A,Bの 各 桁 を 数 列 と 見 なしFFT 変 換 を 行 う 変 換 後 の 値 を Ak,Bkとする 2. 項 毎 にAkとBkを 乗 算 し Ckを 得 る 3. Ckを 逆 FFT 変 換 し 桁 上 げ 処 理 をすればCが 得 られる FFT(Fast Fourier Trasform)は 離 散 Fourier 変 換 を 高 速 に 行 う 方 法 です 計 算 工 数 がO(^2)からO(log)に 減 少 する 11
計 算 の 高 速 化 2 チュドノフスキーの 公 式 のΣの 部 分 はbiary splittigの 手 法 を 使 用 し て 計 算 します 最 後 の 一 項 になるまでは 除 算 抜 きで 計 算 できます SL = L / 2 1 = 0 1 l= 0 l= 0 B B 2l 2l + 1 C 2lC 2l + 1 ( A2 C 2 + 1 + B2 A2 + 1) A = 13591409 + 545140134 B = - (2 + 1) (6 + 1) (6 + 5) C = 10939058860032000^3 =2 : 3.141592653589793238462643383587350688475866345996374315654905 =4 : 3.141592653589793238462643383279502884197169399375105820984947 =7 : 3.141592653589793238462643383279502884197169399375105820974944 5923078164062862089986280348253421170657 残 り1 項 まで 計 算 を 進 め 円 周 率 は 最 終 的 に C 0 π = 10005 A0 で 求 めます 12
計 算 の 高 速 化 3 C 0 10005及 び の 計 算 はニュートン 法 を 使 用 して 計 算 します A0 平 方 根 平 方 根 はその 逆 数 を 計 算 し 元 の 数 を 乗 算 して 平 方 根 を 計 算 する C = 1 / SQR(A) SQR(A) = A * (1 / SQR(A)) 漸 化 式 a +1 = a + 0.5 * a *(1-A * a 2 ) A = 2 A = 0.7として ( 計 算 を 繰 り 返 す 毎 に 精 度 が2 倍 になっていく -> 収 束 が 早 い) 1 : 0.707 2 : 0.707106757 3 : 0.7071067811 8654628345 1619907 4 : 0.7071067811 8654752440 0844362101 5823014625 2648639901 0579459709 5 : 0.7071067811 8654752440 0844362104 8490392848 3593768847 4036588317 除 算 除 算 は 除 数 の 逆 数 を 求 めてその 値 と 被 除 数 を 乗 算 して 求 める C = A B = A * (1 / B) 漸 化 式 a +1 = a *(2-A * a ) 13
各 計 算 の 検 証 いくつかある 中 の 代 表 例 として 1. 計 算 する 上 位 桁 に0を 付 けて 計 算 する 例 えば12x34なら 0012x0034として 計 算 し 計 算 結 果 の 上 位 桁 に 予 想 される 数 の 0があるか 確 認 する 12x34=0012x0034=00000408 2. 計 算 結 果 の 余 りと 計 算 する 値 の 余 りを 計 算 結 果 の 演 算 と 同 じ 計 算 をしてその 余 りと 等 しい 事 を 確 認 する 余 りの 計 算 に 使 用 する 除 数 は 素 数 を 使 用 する (A + B) mod p=((a mod p) + (B mod p)) mod p (A - B) mod p=((a mod p) - (B mod p)) mod p (A x B) mod p=((a mod p) x (B mod p)) mod p pは(2^61)-1=2,305,843,009,213,693,951( 素 数 ) 14
15 検 証 計 算 に 使 用 した 公 式 ) 6 8 1 5 8 1 4 8 2 1 8 4 ( 16 1 0 + + + + = = π BBP(Bailey-Borwei-Plouffe)の 公 式 + + + + + + + + + = = 9 10 1 7 10 2 5 10 2 3 10 2 1 10 2 3 4 1 1 4 2 1024 1) ( 2 1 2 2 6 8 5 0 6 π Fabrice Bellardの 公 式 上 記 の 公 式 は 収 束 は 早 くないが 2 進 展 開 における それまでの 桁 を 計 算 しなくても 目 的 の 桁 を 計 算 でき る
計 算 結 果 の 検 証 1.16 進 の 計 算 結 果 の 検 証 10 進 5 兆 桁 に 相 当 する16 進 4,152,410,118,610 桁 の 手 前 32 桁 値 と ピンポイントで 任 意 の 桁 が 計 算 ができるBBPとBellardの 両 式 で 計 算 した32 桁 の 値 が 一 致 しているかを 確 認 した E4 7C587FB338 E55505A1930 5CCB829FC8 : 4,152,410,118,610 2.16 進 -10 進 変 換 の 検 証 以 下 の 式 を 計 算 し 一 致 する 事 を 確 認 した A = Floor(π10 * 10 ^ N) mod p B = Floor(π16 * 10 ^ N) mod p A = Bになる 事 を 確 認 する π10は10 進 の 円 周 率 π16は16 進 の 円 周 率 N = 5,000,000,000,000(10 進 の 計 算 桁 数 ) P = 64ビットの 素 数 16
パソコンの 諸 元 CPU : X5680(6core,3.33GHz) x 2 Memory : 96GB(8GBx12) Motherboard : Z8PE-D12 OS : Widows Server 2008 R2 VGA:マザーボード 内 蔵 HDD: 1TB( 起 動 用 ) 1TB x 1 6TB( 結 果 保 存 ) 2TB x 3 24TB( 計 算 用 ) 2TB x 16 17
ハードウェアの 選 択 1 限 られた 予 算 入 手 先 などから 部 品 を 集 めテストをして 現 在 の 仕 様 になりました 1.CPU 2 個 のCPUが 連 携 できる 事 からXeoを 選 択 2.メモリ メモリ 容 量 をできる 限 り 増 やす 為 1 枚 当 たり8GBの 容 量 を 採 用 3.マザーボード 速 度 と 安 定 性 の 両 方 を 鑑 みZ8PE-D12とした 4.OS ディスクアクセスの 早 い Server 版 のOSとした 18
ハードウェアの 選 択 2 5 兆 桁 の 計 算 に 必 要 なメモリは24TBです PCではこの 容 量 のメモリは 搭 載 できません ハードディスクとRaidカードの 選 択 は 重 要 です C P U チ ッ プ セ ッ ト RAIDカ ー ド 6.4GT/S 8GB/S 600MB/S x 8 ハ ー ド デ ィ ス ク 100MB/S 19
ハードウェアの 選 択 3 5.ハードディスク メーカや 同 じ 型 式 でもそのロットを 変 えて 安 定 度 と 速 度 の 両 面 から 総 合 評 価 をしてST32000641ASを 選 んだ 接 続 台 数 についてもテストし 16 台 が 最 適 と 判 断 した 6.Raidカード ハードディスクを8 台 接 続 した 場 合 でも 台 数 に 比 例 した 転 送 速 度 を 引 き 出 せた9260-8iを 選 んだ 20
ハードディスクのアクセス 速 度 ハードディスクのアクセス 速 度 は メモリに 比 較 すると 相 当 遅 い また 機 械 的 動 作 を 伴 うので ランダムアクセス は 特 に 遅 い と 言 うことから この 速 度 が 計 算 時 間 に 大 きな 影 響 を 与 える アクセス 速 度 を 上 げる 為 に... 1.ドライブは 複 数 同 時 にアクセスする 方 法 2.OSやRaidに 頼 らないアクセス 方 法 21
5 兆 桁 の 計 算 結 果 5 兆 桁 の 手 前 100 桁 は 以 下 の 通 り 2597691971 6538537682 7963082950 0909387733 3987211875 6399906735 0873400641 7497120374 4023826421 9484283852 0から9までの 出 現 頻 度 は 0 : 499,998,976,328 1 : 499,999,966,055 2 : 500,000,705,108 3 : 500,000,151,332 4 : 500,000,268,680 5 : 499,999,494,448 6 : 499,998,936,471 7 : 500,000,004,756 8 : 500,001,218,003 9 : 500,000,278,819 円 周 率 の 値 の 並 びは 乱 数 のように 見 える 22
ギネスブックへの 申 請 2010 年 12 月 1 日 に 申 請 を 行 い 2011 年 1 月 13 日 付 で 認 定 書 をもらいました 23
ご 静 聴 ありがとうございました 24