Hundredsoft

x64では最低でもSSE2まで使えるので、今更FPUを使うことはあまりありません。

もし、倍精度(仮数部53bit)でも精度が足りない場合、
SSE2を使って４倍精度浮動小数点ライブラリを組むか、
４倍精度に対応した x64向けのコンパイラを用意するのが普通でしょう。

しかし、計算速度を要求されるようなケースで、もし64bit精度で足りるなら、
FPUによる拡張倍精度(仮数部64bit)という選択肢もあります。
（例えば10進で、あと1桁だけ精度が欲しい場合にコンパイラまで変えないですよね）

以前にx86用の拡張倍精度クラスを書きましたが、
今度はx64用のサンプルコード。

x64では、インラインアセンブラが使えないので、組み込み関数(Compiler Intrinsics)を使うか、
MASM ( ml64.exe ) を使ってアセンブリをコーディングするかの選択になりますが、
FPUを使うので、ml64用の.asmファイルを作ることにします。
ターゲットは、x64 のWindowsで Visual Studio で開発しているものとします。
他のコンパイラだと呼び出し規約( MSのFastCall )が異なることがあるので注意。

double をオーバーラップするような Ldouble クラスを作る事にします。
メンバとして、拡張倍精度のデータを１つ持ちます。
必要なのは 10バイトですが、切りの良いところで 16バイト取っています。

まず、ヘッダーファイル「Ldouble.h」です。
基本的な演算をクラスに入れて、三角関数,指数関数等を外部関数として呼び出せるようにしています。



アセンブリで作成するのは、extern "C"の関数群です。
これが「Ldouble.asm」です。



最後に、Ldoubleクラスを使ったサンプルコード「main.cpp」です。




VS2010でBuildする際は、新規プロジェクト(VisualC++,Win32,Win32コンソールアプリケーション,空のプロジェクト)で、
ソース( Ldouble.asm, Ldouble.h, main.cpp )を追加して、
構成マネージャーでプラットフォームをx64にします。
新規でx64が選べない場合は、VC++2010の64bit対応のインストールが足りないかもしれません。
その場合は、ml64.exeも入ってないかもしれないので確認して下さい。

続いて、ソリューションエクスプローラの Ldouble.asm を右クリックして、プロパティ。
構成プロパティ・全般で項目の種類→カスタムビルドツールに変更して下さい。

全般の下に出てくるカスタムビルドツールのコマンドライン欄に、
ml64.exe /DWIN_X64 /Zi /c /Cp /Fl /Fo$(IntDir)\$(InputName).obj Ldouble.asm
出力ファイル欄に
$(IntDir)\$(InputName).obj
を指定して下さい。

これで Build が通るようになります。

例外は考えていません。ご利用の際にはFPU制御ワードを適宜設定して下さい。
(Ldouble.asmのdeffcwを変えて下さい。)

(参考) FPU制御ワードについて

 MSB
  F  E  D  C  B  A  9  8  7  6  5  4  3  2  1  0
 +-----------+-----+-----+-----+----------------+
      NC      丸め  精度   NC      例外マスク
 +-----------+-----+-----+-----+----------------+

  例外マスク bit5〜0
   0:無効操作例外
   1:デノーマル例外
   2:零除算例外
   3:オーバーフロー例外
   4:アンダーフロー例外
   5:不正確結果例外

  精度 bit9,8
   00: 24bit(単精度)
   01: reserved
   10: 53bit(倍精度)
   11: 64bit(拡張倍精度)

  丸め bit11,10
   00: 最近値（四捨五入に近い）
   01: 切り捨て(floor)
   10: 切り上げ(ceil)
   11: 零方向切り捨て

---

Ldoubleへの入力は、double か整数です。
出力は double型 だけです。これは、double へのキャストとして実現されます。
もちろん、LdoubleからLdoubleへのコピーで精度が落ちることはありません。

Ldoubleは、64bit精度を保持していますので、
Ldoubleを使った演算中は常に64bit精度を保ちます。
出力(キャスト) は double(53bit) 精度に丸めた値を返しますが、
Ldouble が持っている値が 53bit に丸められるわけではありません。

この方法はとてもうまく問題を解決してくれるかもしれません。
計算精度ボトムの一部だけをLdoubleによる高精度計算に置き換えて、
その後キャストするだけで、今までの計算ルーチンに戻せます。
例えば、１つの double型のスタック変数を Ldouble にしただけ(1文字修正)で問題が解決するかもしれません。
とても自然に組み込めるのは良いのですが、コードを一見しただけだと変化（影響）がわかりにくいので、
コメントは丁寧に書く必要があります。

---

さて、Build後にできたexeを実行すると次の結果が得られました。

as:= 1.0000000000000000e+001
bs=: 1.2300000000000000e+002

+  : 1.3300000000000000e+002: 1.3300000000000000e+002
-  : -1.1300000000000000e+002: -1.1300000000000000e+002
abs: 1.1300000000000000e+002: 1.1300000000000000e+002
if : -1.1300000000000000e+002: -1.1300000000000000e+002
*  : 1.2300000000000000e+003: 1.2300000000000000e+003
/  : 1.2300000000000001e+001: 1.2300000000000001e+001
%  : 3.0000000000000000e+000: 3.0000000000000000e+000
pi : 3.1415926535897931e+000: 3.1415926535897931e+000
sin: 4.9999999999999994e-001: 5.0000000000000000e-001
cos: 8.6602540378443871e-001: 8.6602540378443860e-001
sin: 4.9999999999999994e-001: 5.0000000000000000e-001
cos: 8.6602540378443871e-001: 8.6602540378443860e-001
tan: 5.7735026918962573e-001: 5.7735026918962573e-001
ata: 4.8234790710102493e-001: 4.8234790710102499e-001
at2: 9.8279372324732905e-001: 9.8279372324732905e-001
sqr: 1.4142135623730951e+000: 1.4142135623730951e+000
exp: 1.6439050426651380e+005: 1.6439050426651380e+005
ln : 1.2010000000000000e+001: 1.2010000000000000e+001
pow: 3.3234387002712369e+000: 3.3234387002712369e+000
int: 3.0000000000000000e+000: 3.0000000000000000e+000
dec: 3.2343870027123689e-001: 3.2343870027123672e-001

Machin's pi: 3.1415926535897936e+000: 3.1415926535897931e+000
Leibniz's pi(5000000項での真値): 3.1415924535897932384646433832795
Leibniz's pi(5000000)  double(3.1415924535897797e+000) time = 0.036[sec]
Leibniz's pi(5000000) Ldouble(3.1415924535897930e+000) time = 0.073[sec]

最後のLeibniz's piを除いては、基本演算の確認です。
左側がdouble演算による結果で、右側がLdouble(FPU)による結果です。
Leibniz's piはライプニッツの公式によるπの計算ですが、
わざと演算量を増やして、処理時間を比較しています。
5,000,000項で打ち切って、演算結果と演算に要した時間を表示します。

演算時間に関しては、double演算のおよそ倍の時間が掛っていますが、
double演算では累積誤差により 14 桁の精度しかありませんが、FPUでは 16 桁の精度があります。
三角関数などを使うと、FPUの方が10倍くらい処理時間が掛ることもあります。
doubleに比べて、11 bit多くの精度を出すため処理時間が掛かるようです。
しかし除算については、double演算とFPUで処理時間は同等でした。

逆アセンブラを見てみると、double演算ではSSE2が使われていますが並列化までは行われていません。
うまく並列化（最適化）できれば、double演算はもう少し早くなります。
(Release-Buildではdouble演算側に多少のベクトル化が行われていました。)

正確さと処理時間のトレードオフですので、どちらが良いというものではありません。
目的に即した解決策が見つけられなければ、このような方法もありかと思います。

上記、掲載ソースを圧縮したものをここに置いておきます。

---

Tags: プログラムメモ

現在では10兆桁を超えるπが求められているが、πの特定の桁だけを算出する方法が知られていて、
πの部分算出の記録では500兆桁を超えている。

この部分算出の元になったのが、次の公式です。


この式からd桁目を求めるために、16^dを乗じて、小数点以下を取りだすことを考える。


ここでSを次のように定義する。

すると


これが1995年にSimon Plouffe(プラウフ)により発見された、BBP(Bailey-Borwein-Plouffe)公式と呼ばれるもので、
これをプログラムしたものが、
http://www.experimentalmath.info/bbp-codes/ にあります。(piqpr8.c)
ここでは右向きバイナリ法を使って剰余(mod)を求めているが、剰余計算にモンゴメリ法を使うことを考える。

モンゴメリ法については過去にも
http://www.hundredsoft.jp/win7blog/log/eid93.html や
http://www.hundredsoft.jp/win7blog/log/eid98.html でサンプルを載せていますが、
高速なモンゴメリ法では、法は奇数であることが望ましい。

S1,S5の法は、8k+1,8k+5であるので奇数であるが、
S4,S6は偶数となるので、このままでは使えない

整数論の定理の中に次のものがある。
ならば

が成り立つ。

これをS4の第一項に当てはめると、

であるから、


同様にS6の第一項は、


これで法が奇数となり、高速なモンゴメリ法を使うことができる。
プログラムは次のようになる。 (defaultではπの100万桁(Hex)目を算出する。)

#include <stdio.h>
#include <math.h>
#include <time.h>

int expmong (int p, int ak);

int main()
{
  double pid, s1, s2, s3, s4;
  double series (int m, int n);
  void ihex (double x, int m, char c[]);
  int id = 1000000;
#define NHX 16
  char chx[NHX];

  clock_t start,end;

/*  id is the digit position.  Digits generated follow immediately after id. */
  start = clock();

  s1 = series (1, id);
  s2 = series (4, id);
  s3 = series (5, id);
  s4 = series (6, id);
  pid = 4. * s1 - 2. * s2 - s3 - s4;
  pid = pid - (int) pid + 1.;
  ihex (pid, NHX, chx);

  end = clock();
  printf (" position = %i\n fraction = %.15f \n hex digits =  %10.10s\n time = %.3f[sec]",
  id, pid, chx, (double)(end-start)/CLOCKS_PER_SEC);

  return 0;
}

void ihex (double x, int nhx, char chx[])

/*  This returns, in chx, the first nhx hex digits of the fraction of x. */

{
  int i;
  double y;
  char hx[] = "0123456789ABCDEF";

  y = fabs (x);

  for (i = 0; i < nhx; i++){
    y = 16. * (y - floor (y));
    chx[i] = hx[(int) y];
  }
}

#define eps 1e-17

double series (int m, int id)

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

{
  int k, ak, p;
  double s, t;

  s = 0.;

/*  Sum the series up to id. */

  if (m == 4){
    for (k = 0; k < id; k++){
      ak = 2 * k + 1;
      p = expmong (id - k - 1, ak);
      t = (p * (4 % ak)) % ak;
      s = s + t / ak;
      s = s - (int) s;
    }
  }else if (m == 6){
    for (k = 0; k < id; k++){
      ak = 4 * k + 3;
      p = expmong (id - k - 1, ak);
      t = (p * (8 % ak)) % ak;
      s = s + t / ak;
      s = s - (int) s;
    }
  }else{
    for (k = 0; k < id; k++){
      ak = 8 * k + m;
      t = expmong (id - k, ak);
      s = s + t / ak;
      s = s - (int) s;
    }
  }

/*  Compute a few terms where k >= id. */

  for (k = id; k <= id + 100; k++){
    ak = 8 * k + m;
    t = pow (16., (double) (id - k)) / ak;
    if (t < eps) break;
    s = s + t;
    s = s - (int) s;
  }
  return s;
}

// tのモンゴメリリダクション
int reduction(__int64 t, int n, int n2, int nb)
{
  static int masktbl[] = {
      0,
      0x00000001,0x00000003,0x00000007,0x0000000f,0x0000001f,0x0000003f,0x0000007f,0x000000ff,
      0x000001ff,0x000003ff,0x000007ff,0x00000fff,0x00001fff,0x00003fff,0x00007fff,0x0000ffff,
      0x0001ffff,0x0003ffff,0x0007ffff,0x000fffff,0x001fffff,0x003fffff,0x007fffff,0x00ffffff,
      0x01ffffff,0x03ffffff,0x07ffffff,0x0fffffff,0x1fffffff,0x3fffffff,0x7fffffff,0xffffffff,
  };
  __int64 c = t * n2;
  c &= masktbl[nb];    // mod Rの意味,(nb)bit以上を0クリア
  c = c * n + t;
  c >>= nb;            // 1/r の意味
  if (c >= n) c -= n;
  return c;
}

// モンゴメリ冪剰余計算
int expmong (int p, int ak)
{
/*  expm = 16^p mod ak. */

  int nb = 0;          // 有効Bit数
  int r;               // r: nより大きな２のべき
  for (nb=0, r=1; r<ak; nb++, r<<=1) ;

  // r2 = r^2 mod ak
  int r2 = ((__int64)r * r) % ak;

  // n*n2 ≡ -1 mod r
  int n2 = 0; /* 求めるN' */
  {
    int t=0, vi=1, ni=nb;
    while (ni--){
      if ((t & 1) == 0){
        t += ak;
        n2 += vi;
      }
      t >>= 1;
      vi <<= 1;
    }
  }

  // 冪剰余 16^p mod ak, バイナリ法の下位桁から計算
  int q = reduction(16 * r2, ak, n2, nb);
  int x = reduction(r2, ak, n2, nb);
  while (p > 0){
    if (p & 1){
      x = reduction((__int64)x * q, ak, n2, nb);
    }
    q = reduction((__int64)q * q, ak, n2, nb);
    p >>= 1;
  }
  return reduction(x, ak, n2, nb);
}

改変は、series() のm=4,m=6の処理とexpm()→expmong()です。

intel系x86コードでの64bit整数演算は早くないので、面白い結果が出ました。

BBP計算速度テスト (x86用,x64用ターゲットを変えてBuildしている)

---

【右向きバイナリ法(WindowsXP 32bit)】
(1回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 2.763[sec]
(2回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 2.804[sec]
(3回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 2.763[sec]

【モンゴメリ法(WindowsXP 32bit)】
(1回目)
position = 1000000
fraction = 0.423429797567801
hex digits = 6C65E52CB4
time = 3.114[sec]
(2回目)
position = 1000000
fraction = 0.423429797567801
hex digits = 6C65E52CB4
time = 3.134[sec]
(3回目)
position = 1000000
fraction = 0.423429797567801
hex digits = 6C65E52CB4
time = 3.204[sec]

【結果】
およそ15%、モンゴメリ法が時間が掛る

---

【右向きバイナリ法(Windows7 64bit)】
(1回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 1.788[sec]
(2回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 1.801[sec]
(3回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 1.786[sec]

【モンゴメリ法(Windows7 64bit)】
(1回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 1.255[sec]
(2回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 1.258[sec]
(3回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 1.223[sec]

【結果】
およそ30%、モンゴメリ法が早い

---

VC2010でBuildした限りでは、
32bitではバイナリ法の方が早く、64bitではモンゴメリ法の方が早くなりました。

バイナリ法のコードも整数演算化してみたのですが、
部分的なコード変更だと、かえって成績が悪くなってしまいました。
バイナリ法のコンパイラオプションでSSE2を有効にしてみてもあまり効果はなかったのですが、
Intelコンパイラを使うと、もう少し早くなるかもしれません。

このコードでも、２〜３千万桁までは計算できますが、
expmong()の変数 r2 と series()の変数 p の型を__int64にすれば、１億桁まで計算できるようになります。
但し、x86(32bit)だと64bit整数演算が増えるとそれだけ処理が遅くなります。
x64で１億桁の計算は、私のPCでは150秒ほどでした。
拡張倍精度(仮数部64bit)と長整数(80bit程度)を組み合わせれば、このアルゴリズムでも100億桁ぐらいは行けると思いますが、これだけの計算量になると普通のPCでは難しいです。（電気代も上がるし...）

---

Tags: プログラムメモ