総論

概要

gcc を用いて，より高精度の数値計算を行なう方法を紹介します．

ここで「高精度」というのは，科学技術計算における物理現象の予測精度とか大次元の連立一次方程式を高精度で解く方法とか精度保証などではなく，浮動小数点数の表現精度と演算精度のことです．

また，「総論」というのはここの一連のページの中の総論という意味です． 題材はあくまで私が実際に使っているものに限定しており，世の中にある諸々の alternative を統べるという意味ではありません．

対象は次のようなものです．

1. OS は FreeBSD の 5 系以降 (/i386 および /pc98)
2. CPU は i386
3. gcc は version 3 以降
4. 題材は gcc 標準の拡張倍精度 (long double)，4 倍精度 (__float128)
   および非標準の library である qlib

御存知のとおり，c などで使われている浮動小数点数の表現精度は有限です． i386 + gcc の環境では，浮動小数点数には単精度，倍精度，拡張倍精度，4 倍精度 [1] (float, double, long double [2], __float128) の 4 種類があり，各々の表現精度は 10 進数で約 7, 16, 19, 34 桁です． (5 節 を参照．)

ほとんどの場合には倍精度で事足りると思われますが，そうでないこともあり得ます． (そもそも，事足りている人はここを読んでいないでしょうし．) 倍精度以上の精度が必要になる事態としては，次のような場合が想定されます．

1. 原理的に精度が必要な計算をしたい．
   例えば，値が 20 桁以上も違う物理量 (例えば化学種のモル分率) の差が必要になるとか，ランク落ちに近い状態の連立一次方程式を解くような場合．[3]
2. 解の誤差が計算結果を定性的に変えていないか確認したい．
   例えば，対称な問題を FEM で解く場合，その連立一次方程式の解は完全に対称ではなく誤差の分だけ非対称になります． その誤差の源や影響を評価したい場合など． (対称性のある問題は，その対称性を考慮して計算領域を限定した方が精度良く解ける． しかし，何らかの物理的な要因でその対称性が破れるような現象に注目する場合，物理的な非対称性と誤差によるそれが混在しては困る．)
3. 収束の悪い級数の値を求めたい．
   例えば，近似解法の精度を検討するために，解析解があればそれと比較するのが普通だと思います．[4] その解析解が精度良く計算できれば良いのですが，それが無限級数で与えられており，かつ計算結果が各々の項よりも 10 桁，20 桁小さいとか悩ましい場合もあります． (20 桁小さいならば，その計算結果は正しいはずがない！)[5]

これらを解決するためには，定式化や計算法の改良，あるいは精度保証付き数値計算法を追求するのも良いでしょう． ただ，浮動小数点数の表現精度，演算精度を高めれば解決する場合もありますし，それで済めば楽なので，選択肢の一つとしてお勧めいたします．[6] 以下ではそのやり方を紹介します．

[1] 4 倍精度 (__float128) を利用できるのは gcc 4.6 以降からです． (対象となる CPU は i386, x86_64, ia64, etc. です．)
  厳密に言えば __float80 も使えるようになりましたが，情報が少なく long double との差異も分からないので，言及しません．

[2] 環境を i386 + gcc に限定しない場合，long double は必ずしも同じものを指していないので要注意です．
   例えば SPARK + gcc では long double はサイズが 128 bit ，表現精度は約 34 桁です． (これは「拡張倍精度」でなく「4 倍精度」と呼ばれています．)

[3] これらの場合には定式化を再考するべきでしょう． そうするとしても，新旧の定式化の結果が比較できれば検証が楽です．

[4]「近似」とは，「粗雑」なことではなく精度を任意に制御できる状態である． (by Van Dyke)
   適切に精度を評価していないような「近似」計算は時として危険です．

[5] その無限級数の収束性を高めるために Euler 変換や Shunks 変換などを使うにしても，その変換自体の正当性を確認するために，結局は精度の良い計算結果が必要です．

[6] これを書いた動機は，FreeBSD/i386 + gcc という環境下にて表現精度に関係した良い資料が見つけられなかったことです．

高精度化のために拡張倍精度 (long double) を使う長所としては，以下のことが挙げられます．

1. gcc に標準で備わっているので，比較的容易に利用できる．
2. 計算時間は double を使う場合と大差ない．

また短所としては，以下のことが挙げられます．

1. 倍精度に比べて 10 進で 3 桁程しか向上しない．
2. 算術関数は，拡張倍精度に対応していないものもある．
3. i386 に関しては，現時点では入念に整備されているとは言い難い状態である．

拡張倍精度は，次のようにして利用することができます．

1. 変数を long double で型宣言する．
2. 定数を long double で扱う場合，数値の後に L を付ける． (実行文だけでなく，数値を読み込む場合も．)
3. 入出力の書式を %Lf または %Le で指定する．
4. 拡張倍精度 (変数，定数，関数) の四則演算を行なうならば，先ず fpsetprec(FP_PE) を実行する．

#include <stdio.h>
#ifdef __FreeBSD__
#include <floatingpoint.h>                    /*  必ず include する．  */

long double sinl(long double);     /*  sinl() を prototype 宣言する．  */
#endif

int main(void)
{
   long double x,y;                       /*  型は long double である．*/

#ifdef __FreeBSD__
   fpsetprec(FP_PE);                /*  これを演算の実行文の前におく． */
#endif

   x = 0.25L*3.14159265358979323846L;       /*  定数はこのようにする． */
   y = sinl(x);                  /*  sinl() は long double 用の sin()  */

   printf("sin(pi/4) = %Le\n",y);    /* 入出力の書式は L を付加する．  */

   return 0;
}

一見簡単そうですが，細かい注意事項がいくつかあります． 詳細は long double をご覧下さい．

高精度化のために 4 倍精度 (__float128) を使う長所としては，以下のことが挙げられます．

1. 表現精度が著しく向上する．(34 桁)
2. gcc 4.6 以降には標準で備わっているので，比較的容易に利用できる．
3. 標準 library の算術関数等が充実している．

また短所としては，以下のことが挙げられます．

1. 倍精度，拡張倍精度に比べて演算に時間がかかる．
2. 数値の入出力が他の浮動小数点数とは別扱いであり，不自由である．

4 倍精度は，次のようにして利用することができます．

1. 変数を __float128 で型宣言する．
2. 定数を __float128 で扱う場合，数値の後に Q を付ける．

#include <stdio.h>
#include <quadmath.h>                               /*  必ず include する．  */

int main(void)
{
   char c[80];
   __float128 x,y;                               /*  型は __float128 である．*/

   x = 0.25Q*3.1415926535897932384626433832795029Q;    /*  定数はこのように  */
   y = sinq(x);                         /*  sinq() は __float128 用の sin()  */

   quadmath_snprintf(c,80,"%.40Qf",y);     /* 入出力の書式は Q を付加する．  */
   printf("sin(pi/4) = %s\n",c);

   return 0;
}

一見簡単そうですが，細かい注意事項がいくつかあります． 詳細は __float128 をご覧下さい．

i386 にて long double 以上の精度で計算をしたい場合，より高精度な浮動小数点数を定義し，その演算をソフトウェアにてエミュレートするしか方策はありません． 幸いなことに，Cephes からそれを実現するための library が配布されています．

高精度化のために qlib を使う長所としては，以下のことが挙げられます．

1. 新たに定義される浮動小数点数 qfloat の精度は 10 進で 100 桁以上ある． (5 節 を参照．)
2. 算術関数が比較的充実している．

また短所としては，以下のことが挙げられます．

1. 計算時間は double を使う場合に比して激増する．
2. FreeBSD の ports になっておらず，自力で環境を整備しなければならない．

この library は c および c++ から利用することができます．

c++ 言語から使う場合，演算子オーバーローディングの機能が使えるので，四則演算も含めて普通の c 言語のように書けます．

#include <stdio.h>
#include "qfloat.h"        // 必ず include する．

extern qfloat qpi;         // M_PI のような定数は変数として用意されている．

int main(void)
{
   qfloat x,y;             // 型は qfloat である．

   x = 0.25*qpi;
   y = sin(x);                          // sin() は qlib の qsin() が使われる．

   qprint(30,"sin(pi/4) =",y,"\n");     // 数値を出力する．

   return 0;
}

この library の作成や利用法の詳細は qlib をご覧下さい．