円周率の計算 (6)
「8 ビット縛り」の 8 ビット・マイコンとして、次は手持ちの STM8S Discovery (STmicroelectronics 製、STM8 コア) を選びました。
STM32F103 を使ったフラッシュ・ライタ/デバッグ・インターフェース部分と、切り離し可能な STM8S ターゲットボード部分とで構成されており、STM32F103 側を切り離して Versaloon 化などをした後に残った STM8S ターゲット・ボード部分は、一部では「燃えないゴミ」とも呼ばれています。
「燃えないゴミ」扱いされていますが、STM8S には乗算命令、除算命令があり、円周率の計算には好都合です。
ただし、8 x 8 = 16 ビットの乗算命令 (MUL) は、4 サイクルと高速ですが、16 / 8 = 商 16 ビット、余り 8 ビットの除算命令 (DIV) 、および 16 / 16 = 商 16 ビット、余り 16 ビットの除算命令 (DIVW) は最大 17 クロックと、「普通」の性能です。
まずは、短いプログラムで実現できる spigot アルゴリズムを、アセンブリ言語でインプリメントしました。
spigot アルゴリズムの説明については、下の web サイトを参照してください。
8 月 13 日版のプログラムでは、この「円周率.jp」で変数名を見やすく書き直した版のプログラムを元にしています。
ただし、円周率.jp のプログラムには、配列 f[] から numerator[] への変換し忘れが1か所あります。
8 月 13 日版のプログラムでは、C で記述してあるので、全ての演算が 16 ビット幅 (以上) に「プロモーション」されているので全く問題は生じないのですが、「8 ビット縛り」のプログラムを作成している途中で、8 ビット演算では不都合が生じる場合があることが分かりました。
それは、内側の i に関するループの、
temp /= denom;
の計算で、n の最初の繰り返しで、i が 7 以下の場合に temp の値が 0x100 を超えてしまうことです。
denom で割った商が 0x100 を超えると言うことは、商 8 ビットを求める除算ルーチンでは結果が正しく得られないことになります。
また、
temp = temp * i + num[i] * base;
の計算で、右辺の temp が 8 ビットにおさまっていない場合には、(temp * i) の計算も 16 ビット x 8 ビットで実行する必要が生じます。
これらの不都合は、n の繰り返しの 2 回目以降、最後まで発生しません。 不都合が生じるのは最初の 1 回だけです。
さらに、base = 100 として、num[] に格納する数値を 100 進数として扱っていますが、その初期値を (base / 5) = 20 から (base / 50) = 2 に変更すると、 temp の値が 0x100 を超えるのは i = 1 に限られることが分かりました。
したがって、i = 1 を特別扱いすれば、8 ビット演算だけで処理できることになります。
ただし、初期値 20 では、結果が「31」、「41」、「59」という系列で得られますが、初期値 2 では、「3」、「14」、「15」という系列になり、それにともなって、有効桁数も 38 から 37 に減少します。
まず、i = 1 を特別扱いするために、i に関するループを i = 2 で抜けるようにします。
i = 1 での処理の様子をたどると、
temp = temp * i + num[i] * base = temp + num[1] * base = temp + (base / 50) * base = temp + 2 * base = temp + 200
となります。
また、
denom = (2 * i) - 1 = (2 * 1) - 1 = 1
ですから、denom による除算は 1 による除算となり、どんな被除数でも割り切れ、
num[i] = temp % denom num[1] = temp % 1 = 0
となり、一方
temp /= denom
temp /= 1
は 1 で割ることになり、temp の数値は変化しません。
これで、i = 1 の処理を終え、i に関するループが終了します。
次に出力に移るわけですが、
i = (out + temp/base)
の式によって、下位桁からの補正を加えた最終的な出力値を得ています。
この式の temp に先ほどの i = 1 の場合の temp の式を代入すると、
i = (out + temp / base) = (out + (temp + (base / 50) * base) / base) = (out + (temp / base + base / 50)) = (out + temp / base + 2) = ((out + 2) + temp / base)
となります。
i = 1 での temp の計算で、200 を加えることにより、8 ビットの限界をはみ出してしまうのですが、その効果は出力値に 2 を足すだけのことと同じであることが上の式から分かります。
そこで、ゼロ・クリアされた状態から出発する「out」の値を、あらかじめ +2 しておいて、temp には 200 を足さないようにしても同じ結果になります。
temp に 200 を足さないとすれば、その後の temp /= denom でも temp は変化しませんから、実質的には i = 1 で処理すべきことは何もなくなり、i に関するループを i = 2 で終了するだけで十分です。
その代わり、「out」の初期値を (base / 50) = 2 に設定しておくことが必要です。
また、num[1] の値は、n の繰り返しの最初の回で、初期値からゼロに書き換えられた後は、毎回ゼロに設定され続けるので、(num[1] * base) の項は、初回以外はゼロであり、i = 1 で処理しないことと同等になります。
この変更を加えた C プログラムを下に示します。
void Spigot( mp8_t num ) { const uint8_t base = 100; // n-ary number const uint8_t nterms = 128; // number of terms const uint8_t n_step = 4; // calculation interval register uint8_t i; // loop index register uint8_t n; // loop index register uint16_t temp; // temporary work / carry register uint8_t out; // output digits / carry adjust register uint8_t denom; // denominator of series printf("Spigot (nterms = %d)\r\xa", nterms); printf("pi = "); for( i = 0 ; i < nterms ; i++ ) { num[i] = base / 50; // initalize array } out = num[1]; for( n = nterms ; n >= (nterms-2*Ndigits_Spigot); n -= n_step ) { temp = 0; for( i = n-1 ; i > 1 ; i-- ) { denom = 2 * i - 1; temp = temp * i + num[i] * base; num[i] = temp % denom; temp /= denom; } i = (out + temp / base); if (nterms != n) { printf("%01d", i / 10); } // if printf("%01d", i % 10); if (0 == ((nterms-n) % (n_step * 5))) { printf("%s", ((nterms == n) ? ".\r\xa" : " ")); } // if out = temp % base; } printf("\r\xa"); printf("\r\xa"); } // void Spigot()
STM8S 用のプログラムは次回掲載します。