exflib Getting Started with C++

多倍長計算ライブラリ exflib の使い方を説明します.
ご覧になりたい項目をクリックしてください. 再度クリックすると,その項目を閉じます.

ソース・プログラム中の文字の色は

を表します.

 1. 必要なファイルの準備

プログラムのコンパイルには,次のものが必要です。 libexfloat.a と exfloat.h の入手方法を説明します.

Linux, MacOSX の場合

ターミナルを開いて次を実行します. PROJ はソース・プログラムを作成するディレクトリです.

% mkdir PROJ

ここからファイルをダウンロードし, PROJ に保存してください.保存した後,次を実行します.

(Linux の場合)
% cd PROJ
% tar jxf exflib-x64-bin-20180620.tar.bz2
% cp exflib-x64-bin-20180620/exflib-x86_64-linux/exfloat.h .
% cp exflib-x64-bin-20180620/exflib-x86_64-linux/libexfloat.a .

(MacOSX の場合)
% cd PROJ
% tar jxf exflib-x64-bin-20180620.tar.bz2
% cp exflib-x64-bin-20180620/exflib-x86_64-intelmac/exfloat.h .
% cp exflib-x64-bin-20180620/exflib-x86_64-intelmac/libexfloat.a .

これで準備は完了です.

 2. 多倍長型の変数の型宣言とコンパイル

exflib では,多倍長精度の実数型 exfloat を利用できます. 次のプログラムは,変数 x を exfloat 型と宣言するだけで, 他に何もしないプログラムです.
// declare.cxx
#include "exfloat.h"    // インターフェースの利用を宣言

int main()
{
  exfloat x;   // 多倍長数型
  return 0;
}
上の内容を declare.cxx という名前のファイルで保存します. 保存先のディレクトリは PROJ にします.
コンパイルは次のようにします. 生成された実行ファイル(プログラム) declare.out を実行しても, 何も表示されません.

コンパイル時のオプション -DPRECISION は 多倍長数の精度を10進法で表します. 上の例では (10進) 100桁を指定しています. 50桁以上,2万桁以下で指定できます.

計算精度は,およそ 19.6 桁 (正確には 64 log 2 10 桁) 刻みで保持しています. 例えば,100桁を指定した場合でも,105桁を指定した場合でも,同じ精度で実行されます. 計算精度を表示するには,次のようにします.

// show_digits.cxx
#include "exfloat.h"
#include <iostream>
using namespace std;

int main()
{
  exfloat x;

  // コンパイル時に指定した計算精度を表示
  cout << exfloat::precision10 << endl;

  // 実際の計算精度を表示
  cout << exfloat::precision << endl;

  return 0;
}
これをコンパイルして実行すると,次のようになります.

(Intel Compiler icc の場合)
% icc -O2 -o digits.out -DPRECISION=100 digits.cxx -lexfloat

(GCC g++ の場合)
% g++ -O2 -o digits.out -DPRECISION=100 digits.cxx -lexfloat

% ./digits.out
100         ← 指定精度 exfloat::precision10 の値
115.595520  ← 計算精度 exfloat::precision の値
% 

10進100桁を指定すると, 実際には 115.6 桁程度の桁数が確保されていることがわかります.

 3. 初期化と代入

exfloat 型の変数に代入できるのは,次のとおりです. を代入することができます.最後の「リテラルの実数値」は, 実数値をダブル・クォーテーションで囲んだものです.
単・倍精度の変数およびリテラル値を代入することはできません.

// init.cxx
#include "exfloat.h"

int main()
{
  exfloat x, y;
  int i;

  x = y;      // 多倍長数の代入
  x = i;      // 整数型変数の代入
  x = 1;      // リテラルの整数値の代入

  // x = 0.1;     // リテラルの単・倍精度数の代入は禁止

  // 文字列での代入が可能
  // ダブル・クォーテーションで囲む
  x = "0.1";
  x = "1/10";     // 文字列中に四則を含めることも可能
  x = "1e-10";    // 指数表現も可能
  // x = "1.0D-1"; // 倍精度表現には対応しない
  x = "#PI/2";    // 数学定数の利用例
  x = "12.34*(5/6)*(7e-8)*#E";   // これも有効

  return 0;
}
ダブル・クォーテーションで囲まれた文字列中の式は, 実行時に適切な精度で評価されて,結果が左辺の変数に代入されます. 文字列中に空白を含めることはできません. この文字列中では次の3つの定数が利用可能です. また,配列,2次元配列(行列) の宣言と初期化も同じようにできます.
// array.cxx
#include "exfloat.h"

int main()
{
  exfloat a[10][20], b[10];

  for(int i = 0; i < 10; i++){
    for(int j = 0; j < 20; j++){
      a[i][j] = i;
    }
  }

  for(int i = 1; i < 10; i++)
    b[i] = i;

  return 0;
}
配列の要素数が大きすぎると、エラー (Segmentation Fault) がでる場合が あるようです。 その場合は「動的メモリ確保」を利用してください。 使用後は必ずメモリ領域を解放してください。 行列などの2次元配列を自然につかう方法は、ご相談ください。
// dynamic_allocation.cxx
#include "exfloat.h"

int main()
{
  const int row = 10;
  const int col = 20;

  // exfloat a[200];
  exfloat *a = new exfloat [ row*col ];

  for(int i = 0; i < row; i++)
    for(int j = 0; j < col; j++)
      a[ i*col + j ] = i;

  exfloat *b = new exfloat [ 10000 ];      

  for(int i = 0; i < 10000; i++)
    b[i] = i;

  // after use
  delete [] a;
  delete [] b;

  return 0;
}

 4. 出力

出力の形式は,次の3通りがあります. 出力の精度を要求しないのであれば,倍精度経由の方法が高速です.
// output.cxx
#include "exfloat.h"
#include <iostream>
#include <iomanip>
using namespace std;

int main()
{
  exfloat x = "#PI";

  // 倍精度数に変換しての出力,高速
  cout << static_cast<double>( x ) << endl;

  // 固定小数点表記,10進法で小数点以下50桁
  cout << setiosflags(ios::fixed) << setprecision(50) << x << endl;

  // 浮動小数点表記,10進法で整数部(1桁)も含めて50桁
  cout << setiosflags(ios::scientific) << setprecision(50) << x << endl;

  return 0;
}

 5. 単・倍精度数の入力

単・倍精度の変数を代入するには,明示的な型変換 exflib_cast が必要です. exflib_cast を使う場合は,多倍長型の数の精度の指定に関わらず, 代入される値は単・倍精度数と同じ精度になります.

// cast.cxx
#include "exfloat.h"
#include <iostream>
#include <iomanip>
using namespace std;

int main()
{
  exfloat x, y, z;

  float f;
  double d;

  y = "0.1";         // 多倍長数として 0.1 を代入

  // 単・倍精度を代入するには,明示的に型変換する
  x = static_cast<exfloat>( 0.1f );

  // 精度の確認  
  z = x - y;
  cout << setiosflags(ios::scientific) << setprecision(50) << z << endl;

  x = static_cast<exfloat>( 0.1 );

  // 精度の確認  
  z = x - y;
  cout << setiosflags(ios::scientific) << setprecision(50) << z << endl;

  // 次も可能だが,x の精度は上と同じになる
  f  = 0.1;
  d = 0.1;
  x = static_cast<exfloat>( f );   // x は 0.1に対し,10進約7桁の精度
  x = static_cast<exfloat>( d );   // x は 0.1に対し,10進約16桁の精度

  return 0;
}
上のプログラムのコンパイル・実行結果は次のとおりです.

(Intel Compiler icc の場合)
% icc -O2 -o cast.out -DPRECISION=100 cast.cxx -lexfloat

(GCC g++ の場合)
% g++ -O2 -o cast.out -DPRECISION=100 cast.cxx -lexfloat

% ./cast.out
 1.49011611938476562500000000000000000000000000000000e-9
 5.55111512312578270211815834045410156250000000000000e-18
% 

これより,単精度数 0.1f は,精確には
0.100000001490...
です ( 0.1 + 1 × 10 -9 + 5 × 10 -10 + ).一方, 倍精度数 0.1 は
0.100000000000000005551115123126...
で ( 0.1 + 5 × 10 -18 + 5 × 10 -19 + ), いずれも 0.1 より真に大きい値であることがわかります

 6. ファイルの入出力

ファイルへの出力も,画面への出力と同様に倍精度型・文字列を利用します.
// file_write.cxx
#include "exfloat.h"
#include <iostream>
#include <iomanip>
#include <fstream>
using namespace std;

int main()
{
  exfloat x, y;

  x = "0.1";
  y = "#PI";

  ofstream fout("file.txt");
  fout << setiosflags(ios::scientific) << setprecision(50) << x << ' ' << y << endl;
  fout.close();

  return 0;
}
1行に複数の値を出力するときは, 文字列の間の空白 ' ' を忘れないようにしてください.

このプログラムをコンパイル・実行すると, 新しくファイル file.txt が生成されます. 次のプログラムは,この file.txt からデータを読み込む例です.

// file_read.cxx
#include "exfloat.h"
#include <iostream>
#include <fstream>
using namespace std;

int main()
{
  exfloat x, y;

  ifstream fin("file.txt");

  // 文字列として読み込み,その文字列を exfloat 型変数に代入
  char str1[60], str2[60];  
  fin >> str1 >> str2;
  x = str1;
  y = str2;
  
  fin.close();

  // 読み込んだ値を画面に表示して確認
  cout << setprecision(ios::fixed) << setprecision(60) << "x=" << x << endl;
  cout << setprecision(ios::fixed) << setprecision(60) << "y=" << y << endl;

  return 0;
}

 7. 算術,組込み函数

// arithmetic.cxx
#include "exfloat.h"

int main()
{
  exfloat x, y, z, w;

  x = "#PI";
  y = "#E";
  w = y;

  z = x + y * w;       // 多倍長数同士の演算

  z = x + 2*y + w/3;   // 整数との四則演算も可能

  y = exp( x );        // 組込み函数の例
  z = pow( x, y );

  return 0;
}
exfloat 型に対しては,以下の組込み函数を実装しています.
fabs, sqrt, sin, cos, tan, exp, pow, log, log10, asin, acos, atan, atan2, sinh, cosh, tanh, mod, floor, ceiling
多倍長数と倍精度数の演算は実装されていません. 多倍長数と「リテラル実数」の演算をおこなうには, 一旦,多倍長数に代入する必要があります.
// arithmetic2.cxx
#include "exfloat.h"

int main()
{
  exfloat x, y, z;

  x = "#PI";

  // 多倍長数と倍精度数の演算は,不可
  // z = x + 0.1

  // 一旦,多倍長型の変数に代入する必要がある
  y = "0.1";
  z = x + y;

  return 0;
}

 8. 比較

// compare.cxx
#include "exfloat.h"
#include <iostream>
using namespace std;

int main()
{
  exfloat x, y, r;

  x = "0.1";
  r = "1e-20";
  y = x + r;

  // 多倍長数同士の比較
  if ( x < y ) {
    cout << "x is less than y" << endl;
  } else if ( x ==  y ) {
    cout << "x is equal to y" << endl;
  } else {
    cout << "x is greater than y" << endl;
  }

  // 単・倍精度数との比較も可能だが,注意が必要
  if ( y < 0.1 ) {
    cout << "y is less than 0.1" << endl;
  } else {
    cout << "y is greater than or equal to 0.1" << endl;
  }

  return 0;
}
上述のとおり,倍精度数 0.1 は 0.1 より真に大きい値です. したがって「0.1 より大きい」数 (y が実際そうである) であっても 「0.1 より小さい」かもしれません. 実際,上のプログラムの実行結果は次のとおりです.

 x is less than y
 y is less than 0.1

なお,利用可能な比較演算子は以下のとおりです.
==    !=    <    <=    >    >=

 9. サンプル・プログラム

200! (階乗) を計算する

// fact.cxx
#include "exfloat.h"
#include <iostream>
#include <iomanip>
using namespace std;

int main()
{
  exfloat f;

  f = 1;

  for(int n = 2; n <= 200; n++){

    f *= n;

    cout << n << ' '
         << setiosflags(ios::fixed) << setprecision(3)
         << f << resetiosflags(ios::fixed) << endl;    

    cout << n << ' '
         << setiosflags(ios::scientific) << setprecision(20)
         << f << resetiosflags(ios::scientific) << endl;    

    cout << n << ' ' << static_cast<double>( f ) << endl;

  }

  return 0;
}
200! は 100 桁を超える整数なので,例えば 計算精度を 500 桁にするほうがよいでしょう.

(Intel Compiler icc の場合)
% icc -O2 -o fact.out -DPRECISION=500 fact.cxx -lexfloat

(GCC g++ の場合)
% g++ -O2 -o fact.out -DPRECISION=500 fact.cxx -lexfloat

% ./fact.out
.....
200 788657867364790503552363213932185062295135977687173263294742533244359449963403342920304284011984623904177212138919638830257642790242637105061926624952829931113462857270763317237396988943922445621451664240254033291864131227428294853277524242407573903240321257405579568660226031904170324062351700858796178922222789623703897374720000000000000000000000000000000000000000000000000.000
200 7.88657867364790503552e+374
200 inf

200! がおよそ 7.89 × 10 374 であることがわかります. これは倍精度型では表現可能な大きさを越えており, static_cast<double> に変換すると inf と表示されます.

Newton 法で 2 の平方根を求める

// newton.cxx
// ------------------------------
// Find zero of f(x)=0
// ------------------------------
#include "exfloat.h"
#include <iostream>
#include <iomanip>
using namespace std;

exfloat f(const exfloat &x)
{
  return x*x - 2;
}

// ------------------------------
// Derivative of f(x)
// ------------------------------
exfloat df(const exfloat &x)
{
  return 2*x;
}

// ------------------------------
// Newton iteration
// ------------------------------
int main()
{
  const int MAX_ITR = 20;
  const exfloat INITIAL_VALUE = 1;
  const double TOLERANCE = 1e-95;

  exfloat x[MAX_ITR];

  // set initial value
  x[0] = INITIAL_VALUE;

  for(int i = 1; i < MAX_ITR; i++){

    x[i] = x[i-1] - f( x[i-1] ) / df( x[i-1] );

    cout << i << ' ' << setprecision(100) << x[i] << endl;

    // if x[i] is close to zero, then exiting FOR-loop
    if ( fabs( f(x[i]) ) < TOLERANCE ){
      cout << "Converge!" << endl;
      break;
    }
   
  }

  return 0;
}
実行結果

% ./newton.out
1 1.500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
2 1.416666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667
3 1.414215686274509803921568627450980392156862745098039215686274509803921568627450980392156862745098039
4 1.41421356237468991062629557889013491011655962211574404458490501920005437183538926835899004315764434
5 1.414213562373095048801689623502530243614981925776197428498289498623195824228923621784941836735830357
6 1.414213562373095048801688724209698078569671875377234001561013133113265255630339978531787161250710475
7 1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641602
Converge!

次の反復結果と合致する桁を赤字で表示しています. 二乗収束する様子がわかります.

 10. 全般にわたる注意事項


fujiwara [atmark] acs.i.kyoto-u.ac.jp