円周率の計算 (9)

8 月 20 日の記事で検討した spigot アルゴリズムの変更点を、8 月 15 日の記事の Arduino 用のスケッチに適用し、乗除算部分をインライン・アセンブラによる 8 ビット演算に置き換えたプログラムを下に示します。 (Spigot 関数の部分のみ)
全体のスケッチのコンパイル後のバイナリ・スケッチ・サイズは 5574 バイトです。

void Spigot( mp8_t num ) // MP work 
{
  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
#if defined(ARDUINO)
  register uint16_t temp2;
#endif
  
  out = 2;
  for( i = 0 ; i < nterms ; i++ ) {
    num[i] = out; // initalize array
  }
  PRINT_STR  ("Spigot (nterms = ");
  PRINT_DEC  (nterms);
  PRINTLN_STR(")");
  PRINT_STR  ("pi = ");
  for( n = nterms ; n >= (nterms-2*Ndigits_Spigot); n -= n_step ) {
    temp = 0;
    for( i = n - 1 ; i >= 2 ; i-- ) {
      denom = 2 * i - 1;
#if defined(ARDUINO)
//             product,  multiplicand, lier, carry 
      mul_add_u8(temp2, (uint8_t)temp,    i, 0);
      mul_add_u8(temp,         num[i], base, 0);
      temp2 += temp;
      temp = div_u16_u8(temp2, denom); // temp = rem8:quo8
      num[i] = (temp >> 8); // remainder
      temp  &= 0xff;        // quotient
#else
      temp = temp * i + num[i] * base;
      num[i] = temp % denom;
      temp /= denom;
#endif      
    } // for (i ...
#if defined(ARDUINO)
    temp = div_u16_u8(temp, base); // temp = rem8:quo8
    i = out + (0xff & temp); // quitient
    out = (temp >> 8); // remainder
    temp = div_u16_u8(i, 10); // temp = rem8:quo8
    if (nterms != n) {
      PRINT_DEC(0xff & temp); // i / 10
    } // if (nterms != n) ...
    PRINT_DEC(temp >> 8); // i % 10
    if (0 == (div_u16_u8((nterms-n),(n_step * 5)))>>8) {
      PRINT_STR(((nterms) == n) ? ".\r\xa" : " ");
    } // if (0 == ...
#else
    i = (out + temp / base);
    if (nterms != n) {
      PRINT_DEC(i / 10);
    } // if (nterms != n) ...
    PRINT_DEC(i % 10);
    if (0 == ((nterms-n) % (n_step * 5))) {
      PRINT_STR(((nterms) == n) ? ".\r\xa" : " ");
    } // if (0 == ...
    out = temp % base;
#endif    
  } // for (n ...
  PRINTLN_STR("");
  PRINTLN_STR("");
} // void Spigot()

Spigot アルゴリズムの分の出力結果を下に示します。
整数部が 100 進法の「31」から「3」になったため、小数点以下の有効桁数が 38 から 37 に減りました。

Spigot (nterms = 128)
pi = 3.
1415926535 8979323846 2643383279 5028841624