円周率の計算 (12) -- BBP 公式による計算 (1)

フラッシュ ROM 2 K バイト、RAM 128 バイトの ATtiny2313 で円周率Πの値を BBP 公式 (後述) を使って UART 経由で 16 進数で 4000 桁出力するという記事 (ドイツ語のサイト)

4000 Stellen von Pi mit ATtiny2313 – Mikrocontroller.net

があることを知り、興味を持ちました。
フラッシュ ROM 1 K バイト、RAM 32 バイト、内部レジスタ 16 本の ATtiny10 で実現してみたいと思ったのですが、この制限はかなりキツいものです。
そこで、まずは、フラッシュ ROM 1 K バイト、RAM 128 バイト、内部レジスタ 32 本の ATtiny13A 向けにプログラムを作成しました。 (ソフトウェア UART など一部アセンブラ、大部分は C で記述)
その結果、プログラム・サイズ 1020 バイトで小数点以下 16 進数 8191 桁まで正しい値を出力するプログラムができました。
実行時間は内部 9.6 MHz クロックで約 3 時間 40 分です。
1000 桁では約 2 分 26 秒です。
デバッグは Arduino および PC 上で行い、コア部分のプログラムを共用して、PC / Arduino / ATtiny13A それぞれで動作可能です。
今回の記事では、単独で Arduino 上で動作可能なコア部分のソースプログラムと、ATtiny13A 用のオブジェクトの HEX ファイルのみを掲載します。
BBP (Bailey-Borwein-Plouffe) 公式とは、1995 年に Simon Plouffe が発見し、現在では 3 人のイニシャルを冠して呼ばれる、16 進数での表現での任意の桁の位置から円周率Πの値を求められる公式のことです。

arctan 系の公式など、通常の級数の和による公式では、求めたい桁数分のデータを全て保持できる長大な作業領域を用意し、級数の計算を繰り返して結果をそこに蓄積していきます。
それに対して BBP 公式では、目的の桁以前のデータの寄与は非常にコンパクトなサイズの変数に収まり、目的の桁以降で必要な精度分のデータを保持できる作業領域さえあれば計算できます。
そのため、RAM 容量の少ないマイコンでの計算に適しています。
しかし、結果は 16 進表現、あるいは同じことですが 2 進表現でしか得られないので、通常の 10 進表現での結果を得るには、16 進数出力をすべてキャプチャして、外部のプログラムで 10 進表現に変換しなければなりません。
BBP 公式のアルゴリズムの説明の PDF およびプログラムは、BBP の最初の「B」である David H. Bailey の web サイト

BBP Code Directory

に置かれています。
C で記述されたプログラム "piqpr8.c" では、目的の桁以降のデータを保持するのに double 型の変数を使用しており、その精度は仮数部のビット数である 53 ビットになります。
このビット数で、16 進で約 1180 万桁までは正しく求まるようです。
前述の ATtiny2313 用のプログラムも Bailey のプログラムを元にしており、float 型の変数を使っています。
float では仮数部 23 ビットですが、求める桁数が 4000 桁で、16 進 1 桁ずつ BBP アルゴリズムで求めているので問題ないようです。
ATtiny13A 用のプログラムでは、同様に Bailey のプログラムを元にしていますが、プログラムサイズを減らすために浮動小数点演算は使わず、すべて整数演算で行っています。
小数表現を格納するのに 8 バイト (64 ビット) 使用して、8000 桁程度に対しては 16 進 8 桁 (32 ビット) 程度の精度が得られました。
そこで、BBP アルゴリズム 1 回の計算で 16 進 8 桁 を表示し、計算ステップを 8 桁きざみとしています。
また、前述の ATtiny2313 用のプログラムでは、

の計算を C プログラムで符号なし 16 ビット整数で行っているので、m が 16 ビットの数だと計算途中でオーバーフローして正しく結果が求まりません。
BBP 公式の第 4 項の分母が 15 ビット以下の整数となる条件

から、正しく計算可能な桁数 k を求めると、

が導かれます。
そのため、ATtiny2313 用のプログラムでは 4000 桁までとしているようです。
C プログラムでの計算を 32 ビット演算で行えば、m の値は 16 ビット最大値の 65535 まで計算可能ですが、プログラムサイズが大きくなります。
PC 用のプログラムではプログラムサイズは気にする必要はないので 32 ビット演算で行い、8191 桁まで求めることが可能です。
ATtiny13A 用のプログラムでは、剰余を求める関数をアセンブラで記述することにより、C では不可能なキャリー発生を検出して、16 ビット演算で m の最大値 65535、つまり 8191 桁の計算まで対応しています。
Arduino ではアセンブラで記述したソースファイル (「.S」 サフィックス) は処理の対象にならないので、インライン・アセンブラの形式で記述しなければなりませんが、面倒なので C での記述のままにしているため、桁数は 4095 桁までとなっています。
ATtiny13A のプログラムの実行結果をキャプチャしたものを下に示します。

3.
243F 6A88 85A3 08D3 1319 
8A2E 0370 7344 A409 3822 
299F 31D0 082E FA98 EC4E 
6C89 4528 21E6 38D0 1377 
BE54 66CF 34E9 0C6C C0AC 

. . . . ＜中略＞ . . . .

BFE2 35BD D2F6 7112 6905 
B204 0222 B6CB CF7C CD76 
9C2B 5311 3EC0 1640 E3D3 
38AB BD60 2547 ADF0 BA38 
209C F746 CE75 

最後の「CE75」が小数点以下 8189 〜 8192 桁目で、正しい値は「CE76」です。
プログラム容量の制限のため、16 進数 4 桁を 5 組の計 20 桁を 1 行に出力しているだけで、その他の区切りはありません。
PC 版および Arduino 版では 16 進 100 桁、つまり 5 行ごとに空行 1 行を出力して見やすくしてあります。
ATtiny13A 用のオブジェクトの HEX ファイルを下に示します。

:1000000009C016C015C014C013C012C011C010C062
:100010000FC00EC011241FBECFE9CDBF10E0A0E677
:10002000B0E001C01D92A537B107E1F77ED1E4C170
:10003000E7CF40916000A9E6B0E0E1E6F0E0342FC0
:100040002D9141112095808190E0830F911D820FA9
:10005000911D819320E0A137B20711F0392FF0CF25
:100060000895CF92DF92EF92FF920F931F93CF9359
:10007000DF936C018B0161307105F1F0C0E1D0E0DC
:1000800002C0C01BD10BC017D107D8F791E0E92EF1
:10009000F12C0EC0C0FE05C0C701BE01A80168D189
:1000A0007C01CE01BE01A80163D1EC01D694C794B6
:1000B000C114D10479F702C0EE24FF24C701DF91F7
:1000C000CF911F910F91FF90EF90DF90CF90089507
:1000D0008F929F92AF92BF92CF92DF92EF92FF9258
:1000E0000F931F93CF93DF93CC24DD24C0E0D0E0A7
:1000F0002EC0809174007601E80EF11CC8018C1BA3
:100100009D0BB701AEDF6091730070E0A70130D1A5
:10011000BC0180E090E0CB017727662701E710E083
:100120004701AA24BB24A501940144D1F8013293CC
:1001300022938F01CB0177276627F0E009361F074E
:1001400091F777DF219688E090E0C80ED91E0091E4
:10015000710010917200C017D10758F2780163E066
:10016000EE0CFF1C6A95E1F754E0C52ED12C4AC075
:10017000809174004701880E911C6091730070E0BB
:1001800080E090E0CB0177276627302F321B232FAA
:100190002370220F220F04C0969587957795679557
:1001A0002A95D2F710926A001092690010926C00A2
:1001B00010926B0010926E0010926D001092700001
:1001C00010926F0036953695E601C31BD109CC0F0E
:1001D000DD1FC759DF4FAA24BB2408C0A501940125
:1001E000E9D03A932A93CB0177276627E0E0C93616
:1001F000DE0709F098F71DDF0F5F1F4F88E090E0E2
:10020000E80EF91E2091710030917200C9014096EC
:100210008017910708F0ACCFDF91CF911F910F911C
:10022000FF90EF90DF90CF90BF90AF909F908F9016
:1002300008951F93182F82958F7085D0812F83D0BA
:100240001F9108951F93182F892FF3DF812FF1DF5E
:1002500080E27ED01F910895CF92DF92EF92FF92BD
:100260000F931F93109272001092710014E001E03E
:1002700094E0F92E82E0C82EB5E0DB2EA6E0EA2E4F
:100280001AC05DD010E08091670090916800DADFBD
:100290001F5F153011F453D010E0809165009091EC
:1002A0006600D0DF80917100909172000896909363
:1002B000720080937100809171009091720020E231
:1002C0008030920758F510926200109261001092EF
:1002D0006400109263001092660010926500109204
:1002E00068001092670000937400F09273001092FF
:1002F0006000EEDEF0927400C09273000093600024
:10030000E7DED092740000937300E2DEE0927400A6
:10031000DFDE1F5F153008F4B6CFB3CF1F910F910A
:10032000FF90EF90DF90CF900895C19AB99A07D0CF
:1003300083E30ED08EE20CD08FDF01D008958DE0E4
:1003400007D08AE005C08F708A3008F0895F805D31
:100350000FB680959AE0F89471E57A95F1F708F078
:10036000C19A08F4C19886959A95B1F70FBE089581
:1003700002C0841B950B84179507D8F720E030E066
:1003800014C060FF08C0280F391F18F02417350764
:1003900010F0241B350B880F991F18F08417950750
:1003A00010F0841B950B769567956115710549F7DB
:1003B000C9010895A1E21A2EAA1BBB1BFD010DC0A5
:1003C000AA1FBB1FEE1FFF1FA217B307E407F50705
:1003D00020F0A21BB30BE40BF50B661F771F881FE1
:1003E000991F1A9469F760957095809590959B0177
:0C03F000AC01BD01CF010895F894FFCFCF
:00000001FF

ATtiny13A の PB1/MISO (6 番ピン) がロジック・レベルのシリアル出力となっています。　通信速度は 38.4 kbps の設定です。
RS232C インターフェース IC などを介して PC 側と接続します。
デフォルトのフューズ設定と一致する 9.6 MHz の内部クロックを使いますが CLKDIV8 ビットはデフォルトの「0」から「1」に書き換えてクロックを 8 分周しない設定にする必要があります。
具体的にはフューズ Low バイト = 0x7A、High バイト = 0x1F です。
ソフトウェア UART の実現には ChaN さんの「suart.S」を利用しています。
ボーレート・タイミング用の定数は、公称 9.6 MHz のクロック周波数に合わせてありますから、実際の発振周波数によってはボーレートに誤差が生じ、文字化けする場合があります。
HEX ファイル上で対処するには、アドレス 0x358 から 2 バイト 0x71, 0xE5 (ldi r23,0x51) のコードを (チェックサムも含めて) 書き換えます。
Arduino 用のスケッチを下に示します。
Arduino IDE 0022 で作成しています。　他のバージョンでは試していません。
「pi_BBP.pde」

/**************************************************/
/* pi_BBP.pde: Pi calculation by BBP algorithm    */
/*             up to 4096 hexadecimal digits      */
/*             using integer arithmetics only     */
/*                                                */
/*             based on "piqpr8.c" (2006-09-08)   */
/*             wrote by David H. Bailey           */
/*                                                */
/* 2012/10/18 Created by pcm1723                  */
/**************************************************/

#if (!defined(USE_AVR) & !defined(USE_X86))
  #define USE_ARDUINO (1)
#endif

#if (USE_ARDUINO)
  #define USE_SERIAL (1)
  #define N_digits (1000)
#else  
  #define N_digits (8192)
//  #define N_digits (1000)
#endif  

#if (USE_SERIAL)
  #define PRINTLN Serial.println
#endif  

#include <stdint.h>

// use 4 elements for fraction data (4 * 16 = 64 bits)
#define N_frac (4)

// pi calculation digit increment (8 hex digits per 1 BPP calc)
#define N_step (8)

// hexadecimal digit position
uint16_t id;

// parameters for series() function
uint8_t  coef, sub_flag;
uint8_t  m;

// fraction part by integer array
uint16_t frac[N_frac];
uint16_t sum[N_frac];

//
// prototype for mulmod() function
//
uint16_t mulmod(uint16_t a, uint16_t b, uint16_t ak);

//
// add/sub frac[] data to/from sum[] data
//
void add_sub_frac( void )
{
  uint8_t  carry;
  uint16_t temp;
  uint8_t  ff;
  register uint8_t *s_p, *f_p;
  uint8_t  i;
  
  carry = sub_flag;
  f_p = (uint8_t *)(&(frac[0]));
  s_p = (uint8_t *)(&(sum [0]));
  i = (2 * N_frac);
  do {
// make 1's complement if subtract  
    ff = *(f_p++);
    if (sub_flag) { ff ^= 0xff; }
    temp = ((uint16_t)carry + ff + (*(s_p)));
    *(s_p++) = (uint8_t)temp;
    carry  = (temp >> 8);
  } while (--i);
} // void add_sub_frac()

//
// (a * b) mod ak
//
#if !defined(USE_ASM)
uint16_t mulmod(uint16_t a, uint16_t b, uint16_t ak)
{
  uint16_t r; // result
  
  while (a >= ak) { // modulo ak
    a -= ak;
  } // while (a >= ak) ...
  r = 0;
  while (b) { // 
    if (0x01 & b) { // test b0
      r += a;
      if (r >= ak) { // modulo ak
        r -= ak;
      } // if (r >= ak) ...
    } // if (0x01 & b) ...
    a <<= 1; // mult by 2
    if (a >= ak) { // modulo ak
      a -= ak;
    } // if (a >= ak) ...
    b >>= 1; // divide by 2
  } // while (b >= 1) ...
  return(r);
} // uint16_t mulmod()
#endif // #if !defined(USE_ASM)

//
//  expm = (16^p) mod ak. 
//
uint16_t expm (uint16_t p, uint16_t ak)
{
  uint16_t r; // result
  uint16_t b; 
  
  if (1 == ak) { // early return
    return(0);
  }
  r = 1; // intialize result
  b = 16;
  while (b >= ak) { // b = (b mod ak)
    b -= ak;
  } // while (b >= ak) ... 
  while (p) {
    if (0x01 & p) { // (1 == b0)
      r = mulmod(r, b, ak);
    } // if (0x01 & p) ...
    b = mulmod(b, b, ak);
    p >>= 1; // divide by 2
  } // while (p >= 1) ...
  return(r);
} // uint16_t expm()


//  This routine evaluates the series  sum_k 16^(id-k)/(8*k+m) 
//  using the modular exponentiation technique.

volatile void series( void )
{
  uint8_t  j, jj;
//  uint16_t p, tt, temp, ak;
  uint16_t k;
  uint16_t ak;
  uint16_t t;
  uint32_t prem; // partial remainder
  uint16_t *f_p; // fraction data array pointer
  
  prem = 0;
/*  Sum the series up to id. */
// first term (has integer part) of the series
  for (k = 0; k < id; k++){
    ak = 8 * k + m;      // denominator
    t = expm (id - k, ak);   // calc ((16^(id-k)) mod ak)
    t = mulmod(t, coef, ak); // calc ((t * coef)  mod ak)
    prem   = t;
    prem <<= 16;
    f_p = (&(frac[N_frac])); // calc fraction
    do {
      *(--f_p)   = (prem / ak);
      prem  %= ak;
      prem <<= 16; // update partial remainder
    } while (f_p != (&(frac[0]))); 
    add_sub_frac();
    /*  Compute a few terms where k >= id. */
  } // for (k = 0;
// second term (fraction part only) of the series
  for (k = id; k <= id + (N_frac * 4); k++){
    ak = 8 * k + m; // denominator
    prem  = ((uint32_t)coef << 16); // initially 0x40000 or 0x2000 or 0x10000
    { register uint8_t sft;
      sft    = (0x03 &(k-id)); //
      prem >>= (sft << 2); // shift right count = (16^(0x03 & (k-id))
    }
    for (j = 0; j < N_frac; j++) { frac[j] = 0; } // clear result area
    jj = ((uint8_t)(k - id) >> 2); // offset to skip leading zero word(s)
    for (f_p = (&(frac[N_frac-jj])); f_p > (&(frac[0])) ; ) {
      *(--f_p) = (prem / ak);
      prem   %= ak;
      prem  <<= 16; // update partial remainder
    }
    add_sub_frac();
  } // for (k = id; ...
} // void series()

#if defined(USE_SERIAL)
void setup_serial( void )
{
  Serial.begin(38400);
} // void setup_serial()
#endif

void  setup( void )
{
  int8_t i;
  uint8_t h_cnt1;
#if !defined(USE_AVR)  
  uint8_t h_cnt2;
#endif  
  
#if defined(USE_SERIAL)  
  setup_serial();
  Serial.print  ("pi by BBP algorithm (");
  Serial.print  (N_digits, DEC);
  Serial.println(" hexadecimal digits)");
  Serial.println("");
  Serial.print("pi = 3.");
#endif 
  h_cnt1 = 5-1; // item counter per line
#if !defined(USE_AVR)  
  h_cnt2 = 5-1; // block counter
#endif  
  for (id = 0; id < N_digits; id += N_step) {
    for (i = 0; i < N_frac; i++) { // clear result
      sum[i] = 0;
    } // for (i = 0; ...
  
    m = 1; coef = 4; sub_flag = 0; 
    series(); // +4 * S1
    
    m = 4; coef = 2; sub_flag = 1; 
    series(); // -2 * S4
    
    m = 5; coef = 1; 
    series(); // -1 * S5
    
    m = 6;
    series(); // -1 * S6

// output 8 hex digits (two elements of uint16_t array)    
    for (i = N_frac-1; i >= (N_frac-2); i--) {
      if (5 <= (++h_cnt1)) {
        h_cnt1 = 0;
        PRINTLN();
#if !defined(USE_AVR)        
        if (5 <= (++h_cnt2)) {
          h_cnt2 = 0;
          PRINTLN();
        } // if (5 <= (++h_cnt1)) ...
#endif        
      } // if (5 <= (++h_cnt2)) ...
#if defined(USE_SERIAL)
      if (0x1000 > sum[i]) { Serial.print("0"); }
      if (0x0100 > sum[i]) { Serial.print("0"); }
      if (0x0010 > sum[i]) { Serial.print("0"); }
      Serial.print(sum[i], HEX);
      Serial.print(" ");
#elif defined(USE_PRINTF)
      printf("%.4X ", sum[i]);
#elif defined(USE_SUART)
      xmit_hex4s(sum[i]);
#endif // #if defined(USE_SERIAL) ... #elif ...
    } // for (i = 0; ...
  } // for (id = 0; ....
#if defined(USE_SERIAL)
  Serial.println();
  Serial.println();
  Serial.print  ("(Elapsed time = ");
  Serial.print  (millis());
  Serial.println(" ms)");
#endif
} // void setup()

#if defined(USE_ARDUINO)
void loop( void )
{
} // void loop()
#endif

シリアルの通信速度は 38.4 kbps です。
「Serial Monitor」ボタンを押してシリアル・コンソールを起動して結果を見ることができます。
前述のように、桁数は 4095 桁までに制限されています。
桁数を決めているのは N_digits の定義で、デフォルトでは 1000 になっているので必要に応じて書き換えます。
16 MHz クロックの Arduino で 1000 桁の計算に約 2 分かかります。
PC 用や ATtiny13A 用のプログラムでは、上のプログラムをコア部分としてインクルードしています。
PC 用や ATtiny13A 用のソースプログラムは次回以降に掲載します。