円周率の計算 (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 サイトを参照してください。

円周率.jp

http://www14.ocn.ne.jp/~kk62526/pi/Spigot.html

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 用のプログラムは次回掲載します。