円周率の計算 (14) -- BBP 公式による計算 (3)

今回は PC (intel x86 プロセッサ) 用の C ソース・プログラムを示します。
前々回の「pi_BBP.pde」ファイルと、下の「pi_BBP.c」ファイルを同じフォルダに置き、「pi_BBP.c」をコンパイルするだけです。
整数演算だけなので、math ライブラリ (libm) をリンクする必要はありません。
cygwin 上の x86 用 gcc を使用しています。　他の環境では試していません。
「pi_BBP.c」

/************************************************/
/* pi_BBP.c : Pi calculation by BBP algorithm   */
/*                                              */
/* 2012/10/18 Created by pcm1723                */
/************************************************/
 
#define USE_X86    (1)
#define USE_PRINTF (1)
#define USE_ASM    (1)
#define PRINTLN()  printf("\r\xa")

#include "pi_BBP.pde"

#include <stdint.h>
#include <stdlib.h>
 
uint16_t mulmod(uint16_t aa, uint16_t b, uint16_t ak)
{
  uint32_t r; // result
  uint32_t a;
  
  a = aa;
  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()

int main()
{
  printf("pi = 3.\r\xa");
  setup();
  return(0);
}

やっていることは、単に mulmod() 関数をこのファイル中で定義されているものに置き換えているだけです。
mulmod() 関数の定義でも、単にワーク用の変数を 32 ビットとして 8191 項までの計算に対応させているだけです。
浮動小数点演算ハードウェアを持たない ATtiny13A のために、すべての演算を整数型で行っているのですが、x86 で実行させると、かえって遅くなり、浮動小数点演算を使用する Bailey の元のコード「piqpr8.c」を利用したプログラムの方が実行時間が短くなります。
例をあげると、Celeron M430 (1.73 GHz, 65 nm Yonah-1M コア, FSB 533 MHz) の PC 上で小数点以下 8192 桁分を求める実行時間が

• 浮動小数点演算: 約 26 秒
• 整数演算 : 約 49 秒

と、ほぼ倍の時間となっています。