HSPCLの開発やパズドラルート解析GPGPU版のOpneCL移植版や数値流体力学で大気シミュレーションなどやりたいことはまだいろいろあるのだが、ここで少し寄り道しようと思う。
「今度こそ円周率小数点以下第1000億～1兆桁くらいまで計算したい！」
昔からキチガイで円周率が無駄に大好きで小学生の頃よく覚えて今も3.14・・・・・・・・・・7982くらいまでならおもいだせるくらい好きなので、何年前かHSPコンテストに高速電子計算機と言う名で円周率を計算するソフトを出したくらいだ
これは高速化にも力を入れていて、HSP標準命令を使わず多倍長計算用プラグインをこしらえて重い処理はプラグインのほうに投げていたし、 円周率を求める公式には世界最速級のソフトで多用されているチュドノフスキーの公式を使った。当然ソースは計算の無駄のないようにした。
だが、これだけでは案の定計算速度が遅かった。
理由は３つある。一つは巨大桁の乗算に筆算法を用いていたこと、もうひとつはGPUではなくCPU1スレッドだけで計算していたこと、もうひとつはSSE4などで高速化を行なわなかったこと。
後者２つは合わせてもせいぜい10～100倍の速度差にしかならないが、乗算アルゴリズムをどうするかというのは、この円周率を計算する分野において非常に重要になってくる。1億桁の乗算でも15万倍近く計算時間に差が出、1兆桁なら多分1億倍くらい（？）処理量が違ってくるのである！
作品を出してからも乗算のアルゴリズムの勉強は続け、どうしてもFFTが理解できずkaratsuba法に逃げた事もあったが、最近FFTが実は流体シミュのCIP法なんかよりもよっぽど簡単だったんじゃないかということに気づきOpenCLに移植するところまではできた。
昔企画していた1億桁×1億桁の乗算はあまりに大変だとのことで一度凍結させたのが、このFFTの乗算アルゴリズムを使うことでいとも簡単に3000万桁×3000万桁の計算がHSPのしかも標準命令でできることがわかってしまったので、もっと上を目指してもいいんじゃないかということでこのプロジェクトを掲げてみた。
でもやはりいきなりは無理だろうから、まずは100万桁～1億桁の円周率を求めるプログラムを作ることにする。
早速だが企画内容はこうだ
使用言語はC++とOpenCL
数値計算部分は主にOpenCLつまりGPUで行なう
円周率を求めるのに使う公式はラマヌジャン系公式のチュドノフスキーの公式
チュドノフスキーの公式のΣの内部の計算手順はBinary Splitting 法を使用
√10005の計算にはニュートン法で
最後の除算もニュートン法で
ニュートン法やバイナリースプリッティングでところどころにでてくる多倍長の乗算はFFTを用いる
FFT自体の計算精度はdouble-double精度でやる
double-double精度の四則演算は当然GPU上でサポートされていないのでソフトウェア的に行なう
FFTの１つ１つの要素は基数を102502425600として計算
これはhttp://www.nextftp.com/swlabo/m0_pctech/hp_ultraprecision/up_815_3.htmを見てdouble-double精度なら100億配列要素までなら基数を10^11にしてギリギリいけると思ったからである。あと基数を102502425600にすればチュドノフスキーの公式の性質上、乗算を配列シフトで済むところがあってここでも利点があるためである。
目指す目標は2つのGPU(Radeon HD7970)を使って１分以内に１億桁を計算しファイルに出力すること。
http://ysrken.blog.fc2.com/blog-entry-38.html
このブログによればCPUだけでも1億桁1分を切るものが出ているようなのでとりあえずはGPGPUでアドバンテージがある分1分という目標は妥当かなと・・・
そして次のHSPコンテストではHSPCLの公開と同時に、高速電子式計算機,改を出そうと思う！
盛大な寄り道だが、これもまた今まで作ってきたプログラムの集大成みたいなものなので個人的にはぜひチャレンジしたい！
よし、まずは周波数間間引きFFTと時間間引きFFTの勉強から始めるとするか。
------------------------------------------------------------------------------------------------------------------------------------------
2019/10追記
FMTをPythonとCUDAで実装しました。ガウスルジャンドルで3億桁まで求めた記事をQiitaに公開しました。1億桁1分は切れなかったよ・・・
https://qiita.com/Red_Black_GPGPU/items/e933e0d846b874b86c32
ソースはこちら
https://github.com/toropippi/FMT_QiitaSample
「今度こそ円周率小数点以下第1000億～1兆桁くらいまで計算したい！」
昔から
これは高速化にも力を入れていて、HSP標準命令を使わず多倍長計算用プラグインをこしらえて重い処理はプラグインのほうに投げていたし、 円周率を求める公式には世界最速級のソフトで多用されているチュドノフスキーの公式を使った。当然ソースは計算の無駄のないようにした。
だが、これだけでは案の定計算速度が遅かった。
理由は３つある。一つは巨大桁の乗算に筆算法を用いていたこと、もうひとつはGPUではなくCPU1スレッドだけで計算していたこと、もうひとつはSSE4などで高速化を行なわなかったこと。
後者２つは合わせてもせいぜい10～100倍の速度差にしかならないが、乗算アルゴリズムをどうするかというのは、この円周率を計算する分野において非常に重要になってくる。1億桁の乗算でも15万倍近く計算時間に差が出、1兆桁なら多分1億倍くらい（？）処理量が違ってくるのである！
作品を出してからも乗算のアルゴリズムの勉強は続け、どうしてもFFTが理解できずkaratsuba法に逃げた事もあったが、最近FFTが実は流体シミュのCIP法なんかよりもよっぽど簡単だったんじゃないかということに気づきOpenCLに移植するところまではできた。
昔企画していた1億桁×1億桁の乗算はあまりに大変だとのことで一度凍結させたのが、このFFTの乗算アルゴリズムを使うことでいとも簡単に3000万桁×3000万桁の計算がHSPのしかも標準命令でできることがわかってしまったので、もっと上を目指してもいいんじゃないかということでこのプロジェクトを掲げてみた。
でもやはりいきなりは無理だろうから、まずは100万桁～1億桁の円周率を求めるプログラムを作ることにする。
早速だが企画内容はこうだ
使用言語はC++とOpenCL
数値計算部分は主にOpenCLつまりGPUで行なう
円周率を求めるのに使う公式はラマヌジャン系公式のチュドノフスキーの公式
チュドノフスキーの公式のΣの内部の計算手順はBinary Splitting 法を使用
√10005の計算にはニュートン法で
最後の除算もニュートン法で
ニュートン法やバイナリースプリッティングでところどころにでてくる多倍長の乗算はFFTを用いる
FFT自体の計算精度はdouble-double精度でやる
double-double精度の四則演算は当然GPU上でサポートされていないのでソフトウェア的に行なう
FFTの１つ１つの要素は基数を102502425600として計算
これはhttp://www.nextftp.com/swlabo/m0_pctech/hp_ultraprecision/up_815_3.htmを見てdouble-double精度なら100億配列要素までなら基数を10^11にしてギリギリいけると思ったからである。あと基数を102502425600にすればチュドノフスキーの公式の性質上、乗算を配列シフトで済むところがあってここでも利点があるためである。
目指す目標は2つのGPU(Radeon HD7970)を使って１分以内に１億桁を計算しファイルに出力すること。
http://ysrken.blog.fc2.com/blog-entry-38.html
このブログによればCPUだけでも1億桁1分を切るものが出ているようなのでとりあえずはGPGPUでアドバンテージがある分1分という目標は妥当かなと・・・
そして次のHSPコンテストではHSPCLの公開と同時に、高速電子式計算機,改を出そうと思う！
盛大な寄り道だが、これもまた今まで作ってきたプログラムの集大成みたいなものなので個人的にはぜひチャレンジしたい！
よし、まずは周波数間間引きFFTと時間間引きFFTの勉強から始めるとするか。
------------------------------------------------------------------------------------------------------------------------------------------
2019/10追記
FMTをPythonとCUDAで実装しました。ガウスルジャンドルで3億桁まで求めた記事をQiitaに公開しました。1億桁1分は切れなかったよ・・・
https://qiita.com/Red_Black_GPGPU/items/e933e0d846b874b86c32
ソースはこちら
https://github.com/toropippi/FMT_QiitaSample
