公開日：2021/09/04

17th Season

#17-8　円周率とデュアルコア

　Moscoviumデュアルコア版MCOC115DC0416の実力を試すため選択したのが円周率の計算プログラムだったのですが、実はあまりデュアルコアの真価を発揮できませんでした(↓)。

　　「#17-3　双頭の狼(sun)」：　Moscovium (Dual Core)，calc_pi.asm

　今回このテーマをもう少し掘り下げて再挑戦してみます。

円周率の計算(再考)

　これまで何度か扱ってきたので、過去記事も当たってみてください。
　本HPで円周率を求めるために使って来た公式は、以下となります。



　マチンの公式と、逆正接atan()関数のグレゴリー級数展開です。

　元々、円周率πとatan()関数の関係は、以下式で表されます。

　　atan (1)＝π/4

　　∵tan (π/4)＝sin (π/4)/cos (π/4)＝(1/√2)/(1/√2)＝1

　atan (1)がπ/4となるので、そのまま、グレゴリー級数に代入して×4倍すると、πが計算できます。グレゴリー(テイラー)級数については、「#16-9　数学関数を整備しよう」で解説しました。
　グレゴリー級数展開は、項番kの値が∞まで続きます。コンピュータなどで近似値を求める時は、必要な精度が得られるkまでを計算して求めます。ところが、一般に1/xの値が大きくなってくると、必要な精度を確保するための項番kの値がどんどん大きくなっていきます。
　つまり、1/xの値を小さくするのは、高速に計算するための必須条件なのです。

　そこで登場するのが、atan()関数の加法定理です(↓)。



　a×b＝n²＋1となる整数を選びます。スタートのatan (1)＝π/4はn＝1と見ることが出来るので、例えばa＝1、b＝2が合致します。すると、式が以下の様に変形できます。

　　π/4＝atan (1)＝atan (1/2)＋atan (1/3)

　これも円周率πの計算公式として知られており、恐らく、最初に提唱した数学者の名前を取って、オイラーの公式と呼ばれています。1/x＝1だった最初の式atan (1)から、1/x＝1/2と1/x＝1/3に絶対値が小さくなって、収束が速くなりました。その代わり、atan()関数のグレゴリー級数を2回計算しなければなりません。それでも、1/xの絶対値が小さくなって高速に収束する効果の方が、計算時間短縮に寄与するという訳です。
　後は、この調子で加法定理を次々と適用して行くと、マチンの公式も導出できます。この時、以下の関係式も使用されます。

　atan()関数は、原点に対して点対象なので、

　　atan (x)＝−atan (−x)

　従って、

　　atan (−x)＝−atan (x)

　a＜0，b＜0となるようなa，bを選択して、atan()関数の引数が負数になった場合、引数の絶対値を取って係数の符号を反転させることが出来る、と言う関係式です。


　以下表に、円周率πの計算公式をいくつか挙げておきます。

　以降、上表の「名前」欄に記した名前で表記します。

　本HPで使って来たのは上表の中でマチン(Machin)公式と呼ばれるものです。
　グレゴリー級数の特徴として、1/xの値が小さい方が収束が速く、結果として高速に(短時間で)円周率を求めることが出来ます。従って、各公式の中で、atan()関数の引数がなるべく小さくなるものが有利です。とは言え、各公式は2〜4個のatan()関数を加減算する形になっていて、1/xに相当する数値も2〜4個存在するので、トータルでどれが一番速いの？ってことです。

　手元には、円周率を計算するプログラム(calc_pi.ino)が既にあるので、実際に計算させてみましょう。

実測！

　ここからは、Arduino DUEを使います。プログラム(calc_pi.ino)を改造しておきましたので、これを使います。

　　円周率計算プログラム：calc_pi.ino

　まず、先頭付近のUSE_MACHINマクロを有効化してください(デフォルトではコメントしてあります)。

　　`define　USE_MACHIN

　後はcalc_pi()関数の中にコメントしてある公式のどれか1つだけを有効にして、コンパイル＆書き込んでください。

　　　　// Machin
　　　　m_atan_1_x (A, 16,5);
　　　　m_atan_1_x (A, -4,239);

　　//　// Euler
　　//　m_atan_1_x (A, 1*4,2);
　　//　m_atan_1_x (A, 1*4,3);

　　//　// Euler (2)
　　//　m_atan_1_x (A, 4*4,5);
　　//　m_atan_1_x (A, -1*4,70);
　　//　m_atan_1_x (A, 1*4,99);

　　//　// Gauss
　　//　m_atan_1_x (A, 12*4,18);
　　//　m_atan_1_x (A, 8*4,57);
　　//　m_atan_1_x (A, -5*4,239);

　　//　// Stormer
　　//　m_atan_1_x (A, 6*4,8);
　　//　m_atan_1_x (A, 2*4,57);
　　//　m_atan_1_x (A, 1*4,239);

　公式毎に、小数点以下1,000桁の円周率を計算して、計算に要した時間をまとめたのが以下表となります。

※計算時間が短かった順に並べてあります。また、計算した1,000桁の円周率が、下位数桁を除いて期待値と一致することを確認して備考欄に「ok」と記しました。

　マチン(Machin)公式も良い線行っていました(約282ms)が、この中で一番速かったのはKlingenstierna公式で約273msでした。備考欄に変なコメントが付いているKlingenstierna公式は、実は以前作ったプログラム(calc_pi.ino)のままだと、正しい円周率が計算できませんでした。

　プログラムのm_atan_1_x()関数内に、以下の様にxx＝x²で割り算する箇所(↓)があったのですが、この処理で計算を間違っていました。

　　// B=quo...rem=B/sq(x)
　　quo=int16 ((B[idx] + rem)/xx);
　　rem=(B[idx] + rem - quo*xx)*BASE;
　　B[idx]=quo;

　問題だったのは、Klingenstierna公式のatan (1/515)の計算で、1/x＝1/515より、xx＝x²＝515²＝265,225となります。多倍長固定小数を構成する16ビット符号付整数配列B[ ]の上位桁から順にxxで割り算して行きます。この時、上位桁で生じた剰余を下位桁に伝搬させながら割り算を実行しています(↓)。



　上位桁をxxで割り算した剰余remは、0≦rem＜xxですが、下位桁に伝搬するとき、rem×BASEとして桁合わせした後で足し込みます。BASE＝10,000ですので、rem×BASEの最大値は以下となります。

　　MAX (rem×BASE)＝(265,225−1)×10,000＝2,652,240,000＝0x9e15_f880

　莫大な数値になります。ビット数は約31.3 → 32ビット！にもかかわらず、「int32　rem；」で宣言していたのが仇となりました。はい。符号無で32ビットになる数値を符号付32ビット(int32)では計算できません。
　と、言う事で、変数remの宣言をuint32に変更したところ、正しく円周率を計算できるようになりました。

　実は、更に大きなxを使う、Stormer (2)公式があったのですが、こちらはx＝682および12,943ということで、全くダメ。32ビット演算では足りません。Stormer (2)公式に関しては、今回のスコープから除外しました。

デュアルコアに最適な公式は？

　前節で上位3つにランクインした公式について、更に、各atan()項に関してバラバラに計算を実行して、計算に要した時間を計測しました。単にatan()関数値を計算しただけなので、計算した数値自体は無視して、計算に要した時間だけに着目します(下表)。

　上表の結果を一言で言い表すなら、「1/xの値が小さいほど、収束が速い」という基本のまんま、です。問題は、これらのatan()項を如何にしてデュアルコアCPUに振り分けるか？と言う事です。

• Machin公式：
• 最初に使っていた公式です。
• 従って、既にデュアルコア版で計算済みです。
• CPU1：　atan (1/239)項で、約65ms。
• CPU2：　atan (1/5)項で、約217ms。
• この場合、遅い方のCPU2で演算時間が決まるので、217msに、計算結果を統合する処理分を追加した時間で終了することになります。
• デュアルコア化による速度比を簡易的に求めると、(217＋65)/217≒1.30倍速程度と見積もることが出来ます。
• 実測でも似たような結果でした。ガックシ……。
• atan (1/5)項とatan (1/239)項の計算時間が上表の様に圧倒的に不均衡なため、デュアルコアに分散した効果がほとんど見られなかったのが敗因です(↓)。
• 
  
  
  
• Gauss公式：
• 以下の様に分散するのがベストでしょう。
• CPU1：　atan (1/57)項とatan (1/239)項で、約154ms。
• CPU2：　atan (1/18)項で、約123ms。
• この場合、遅い方のCPU1で演算時間が決まるので、154msに、計算結果を統合する処理分を追加した時間で終了することになります。
• デュアルコア化による速度比を簡易的に求めると、(154＋123)/154≒1.80倍速程度と見積もることが出来ます。
• Klingenstierna公式：
• 以下の様に分散するのがベストでしょう。
• CPU1：　atan (1/239)項とatan (1/515)項で、約123ms。
• CPU2：　atan (1/10)項で、約151ms。
• この場合、遅い方のCPU2で演算時間が決まるので、151msに、計算結果を統合する処理分を追加した時間で終了することになります。
• デュアルコア化による速度比を簡易的に求めると、(123＋151)/151≒1.81倍速程度と見積もることが出来ます。

　と、言う事で、デュアルコアに分散した場合でも、Klingenstierna公式の方が、少しだけ速度比が高い(約×1.81倍速)様ですし、そもそも、全体の計算時間も一番速かったので、Klingenstierna公式を使うことに決定します。
　Machin公式の様な、約×1.30倍速なんて寂しい結果には、ならないはずです。


　以下に、Klingenstierna公式を使ってデュアルコアに対応したバージョンのプログラムをリンクしておきます。

　　完成版プログラム：calc_pi.asm

※必要なインクルードファイルなどはこちらからダウンロードできます。


　以下リンクより、必要なデザインデータをダウンロードして、MCOC115DC0416のデザインをコンパイル(論理合成・配置配線)してください。上のアセンブラをアセンブルしたMIFファイルを所定のディレクトリに格納しておくことをお忘れ無き様。手順はこちらです。



　FPGAボード(Sipeed Lichee Tang Primer)に書き込んだら、準備完了です。

実験だ

　「#17-3　双頭の狼(sun)」で実施したのと同じ実験なので、さらりと結果だけをまとめておきます(以下表)。

※計算時間の単位は秒。

　今回、旧Machin公式版も含めて、最新版のプログラムで計測し直しています。プログラム(calc_pi.ino)中の「`define　USE_MACHIN」を無効化してあると、Klingenstierna公式を使って計算します。アセンブラ(calc_pi.asm)の方は、条件付きアセンブルには対応していないので、デフォルトでKlingenstierna公式を使うようになっています。

　以下雑感です。

• 速度比：
• MCOC115DC0416(デュアルコア、Klingenstierna)を基準(1.00)として、正値はMCOC115DC0416より速い、負値はMCOC115DC0416の方が速い。
• Moscovium @24MHzのデュアルコアと言う事で、処理速度的にはMoscovium @48MHzに匹敵しているが、@1MHz換算では24MHzで掛け算していることに留意。
• 同じプラットフォーム(マイコンボード)で比べると、Klingenstierna公式の方がMachin公式よりもわずかに速い。
• MCOC115DC0416(デュアルコア)にフォーカスして、MCOC115CP0408(シングルコア)と比べてみると：
• Machin公式では約×1.29倍(≒50.703/39.270)程度の速度比しか出なかった。
• Klingenstierna公式では約×1.78倍(≒49.017/27.577)と、ようやく、デュアルコアに分散処理させた効果を見ることが出来た。
• 高速に計算できるKlingenstierna公式版にフォーカスすると、
• Moscovium(デュアルコア版)は、Cortex-M3、32ビット(DUE @84MHz)と比べてサイクル性能で約1.11倍(≒662/597)。実測値でも約3.88倍(≒27.577/7.105)程度の差まで詰めている。
• Moscovium @24MHzのデュアルコアと言う事で、実質Moscovium @48MHz相当なのだけれど、@84MHzで32ビットMCUのDUE(Cortex-M3)との速度差がこれなら、大健闘。
• これまで何度か、ESP32(Xtensa LX6)とDUE(Cortex-M3)を比べてきた。
• ずっと、「Xtensa LX6速いな」と、思っていた。
• しかし、こうして16ビットCPUコアMoscoviumとの比較をしてみた結果、Xtensa LX6が速いのではなく、「Cortex-M3が遅い」のだという結論に到達した。
• だって、Moscoviumコアは、そんなべらぼうに速いCPUじゃないよ。
• 教科書レベルの標準的なコアだから。
• 商用の32ビットMCUが演算性能でこんなに肉薄されてはいかんでしょう？
• Thumb自体のアイディアは良かったと思うけど、ネイティブARMを捨てたCortex-M系のMCUは、元々のThumbの良さを忘れて暴走しているのではなかろうか？

　検討段階ではデュアルコアに分散することで、シングルコア時の約×1.81倍速で計算できると見積もったのに対して、実測では約×1.78倍速でした。検討段階では円周率1,000桁の測定を元にしたのですが、実測では5,116桁に増えたことで各atan()項の計算時間差が広がったことと、分散できないオーバヘッド分が考慮されていなかった事が誤差の原因だと思われます。想定の範囲内📏です。