公開日:2021/09/04
title.png

17th Season


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

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

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

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

円周率の計算(再考)

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

a106_calcpi1.png

 マチンの公式と、逆正接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()関数の加法定理です(↓)。

a178_atanadd.png

 a×b=n2+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()関数の引数が負数になった場合、引数の絶対値を取って係数の符号を反転させることが出来る、と言う関係式です。


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

名前 公式
Machin 4×atan (1/5)−atan (1/239)=π/4
Klingenstierna 8×atan (1/10)−atan (1/239)−4×atan (1/515)=π/4
Gauss 12×atan (1/18)+8×atan (1/57)−5×atan (1/239)=π/4
Stormer 6×atan (1/8)+2×atan (1/57)+atan (1/239)=π/4
Stormer (2) 44×atan (1/57)+7×atan (1/239)−12×atan (1/682)+24×atan (1/12,943)=π/4
Euler atan (1/2)+atan (1/3)=π/4
Euler (2) 4×atan (1/5)−atan (1/70)+atan (1/99)=π/4

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

 本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桁の円周率を計算して、計算に要した時間をまとめたのが以下表となります。

名前 公式 DUE[μs] 備考
Klingenstierna 8×atan (1/10)−atan (1/239)−4×atan (1/515)=π/4 272,679 ok, uint32 rem
Gauss 12×atan (1/18)+8×atan (1/57)−5×atan (1/239)=π/4 276,668 ok
Machin 4×atan (1/5)−atan (1/239)=π/4 282,167 ok
Stormer 6×atan (1/8)+2×atan (1/57)+atan (1/239)=π/4 323,079 ok
Euler (2) 4×atan (1/5)−atan (1/70)+atan (1/99)=π/4 378,204 ok
Euler atan (1/2)+atan (1/3)=π/4 831,876 ok
Stormer (2) 44×atan (1/57)+7×atan (1/239)−12×atan (1/682)+24×atan (1/12,943)=π/4 NG

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

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

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

  // 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=x25152=265,225となります。多倍長固定小数点数を構成する16ビット符号付整数配列B[ ]の上位桁から順にxxで割り算して行きます。この時、上位桁で生じた剰余を下位桁に伝搬させながら割り算を実行しています(↓)。

a106_calcpi3.png

 上位桁を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()関数値を計算しただけなので、計算した数値自体は無視して、計算に要した時間だけに着目します(下表)。

名前 数式 DUE[μs] 備考
Klingenstierna 4×(-4)×atan (1/515) 57,209 Σ=122,512
Klingenstierna 4×(-1)×atan (1/239) 65,303
Klingenstierna 4×8×atan (1/10) 150,503  
Gauss 4×(-5)×atan (1/239) 65,302 Σ=153,512
Gauss 4×8×atan (1/57) 88,210
Gauss 4×12×atan (1/18) 123,466  
Machin 4×(-1)×atan (1/239) 65,303  
Machin 4×4×atan (1/5) 217,027  

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

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


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

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

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


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

btn_moscovium.png

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

実験だ

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

円周率計算時間(calc_pi) 5,116桁 速度比 備考
ボード マイコン [MHz] 実測 @1MHz @1MHz
Machin公式版
ESP32 Xtensa LX6 240 1.579 379 1.75 Machin,32bit,MUL,DIV
Arduino DUE SAM3X8E 84 7.355 618 1.07 Machin,32bit,MUL,DIV
STM32 STM32F103C 72 7.757 558 1.19 Machin,32bit,MUL,DIV
FPGAボード MCOC115DC0416 24 39.270 942 −1.42 Machin,16bit,MUL,DIV,DUAL
MCOC115CP0408 24 50.703 1,217 −1.84 Machin,16bit,MUL,DIV
Klingenstierna公式版
ESP32 Xtensa LX6 240 1.525 366 1.81 Klingenstierna,32bit,MUL,DIV
Arduino DUE SAM3X8E 84 7.105 597 1.11 Klingenstierna,32bit,MUL,DIV
STM32 STM32F103C 72 7.493 540 1.23 Klingenstierna,32bit,MUL,DIV
FPGAボード MCOC115DC0416 24 27.577 662 1.00 Klingenstierna,16bit,MUL,DIV,DUAL
MCOC115CP0408 24 49.017 1,176 −1.78 Klingenstierna,16bit,MUL,DIV

※計算時間の単位は秒。

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

 以下雑感です。


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


btn_arduino.png

btn_prev.png btn_next.png

home