exflib Getting Started with C++

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

プログラムのコンパイルには，次のものが必要です。

C++ コンパイラ
ソース・プログラム
libexfloat.a
exfloat.h

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 .

これで準備は完了です．

exflib では，多倍長精度の実数型 exfloat を利用できます．次のプログラムは，変数 x を exfloat 型と宣言するだけで，他に何もしないプログラムです．

// declare.cxx
#include "exfloat.h"    // インターフェースの利用を宣言

int main()
{
  exfloat x;   // 多倍長数型
  return 0;
}

上の内容を declare.cxx という名前のファイルで保存します．保存先のディレクトリは PROJ にします．
コンパイルは次のようにします．生成された実行ファイル(プログラム) declare.out を実行しても，何も表示されません．

Intel Compiler icc の場合


% icc -O2 -o declare.out -DPRECISION=100 declare.cxx -lexfloat
% ./declare.out    ←何も起きず，すぐにプロンプトが表示されます
%

GCC g++ (Linux, MacOSX) の場合


% g++ -O2 -o declare.out -DPRECISION=100 declare.cxx -lexfloat
% ./declare.out  ←何も起きず，すぐにプロンプトが表示されます
%

MinGW (Windows) の場合
MinGW では実行時に libgcc_s_dw2-1.dll, libstdc++-6.dll に関するエラーが出る場合があるようです．この場合は次を試してください．


% g++ -O2 -o declare.out -DPRECISION=100 declare.cxx exfloat.dll -static-libgcc -static-libstdc++
% ./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 桁程度の桁数が確保されていることがわかります．

exfloat 型の変数に代入できるのは，次のとおりです．

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つの定数が利用可能です．

#PI : 円周率
#E : 自然対数の底
#G : Euelr の定数 (0.57722...)

また，配列，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;
}

出力の形式は，次の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;
}

単・倍精度の変数を代入するには，明示的な型変換 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 \times 10^{-9} + 5 \times 10^{-10} + \dots

)．一方，倍精度数 0.1 は

0.100000000000000005551115123126...

で (

0.1 + 5 \times 10^{-18} + 5 \times 10^{-19} + \dots

)，いずれも 0.1 より真に大きい値であることがわかります

ファイルへの出力も，画面への出力と同様に倍精度型・文字列を利用します．

// 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;
}

// 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;
}

// 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. サンプル・プログラム

// 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 \times 10^{374}

であることがわかります．これは倍精度型では表現可能な大きさを越えており， static_cast<double> に変換すると inf と表示されます．

// 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!

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

桁数の精度は，ファイル内で exfloat.h をインクルードする前に #define でも指定できます．

// exfloat.h をインクルードする前に桁数を 10 進法で宣言する
#define PRECISION 100
#include "exfloat.h"

int main()
{
  cout << exfloat::precision << endl;
  return 0;
}

プログラムはひとつのファイルで構成してください．分割コンパイルも可能ですが，計算精度の管理が繁雑になります．
多倍長型の変数毎に，異なる精度を設定することはできません．全ての変数で同じ精度になります．
多倍長数の精度をプログラムの実行中に(動的に)変更することはできません．コンパイル時に指定します．
複素数型も利用できますが，ここでは説明しません．附属のドキュメントをご覧ください．

exflib Getting Started with C++

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

Linux, MacOSX の場合

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

3. 初期化と代入

4. 出力

5. 単・倍精度数の入力

6. ファイルの入出力

7. 算術，組込み函数

8. 比較

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

200! (階乗) を計算する

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

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