calculation using _

1. i386 の 4 倍精度の概説

gcc 4.6 から，浮動小数点数として 4 倍精度 (__float128) を利用することができるようになりました．[1]

とはいえ，拡張倍精度 (long double) がそうであるように，単精度や倍精度に比べて実装の仕方が不自然です．次節以降では，プログラムを書くために必要な 4 倍精度の「癖」を概説します．

当然ながら，拡張倍精度との得失を比較する必要があるでしょう．長所としては，以下のことが挙げられます．

表現精度が著しく向上する．(19 桁 → 34 桁)
標準 library の算術関数等が充実している．

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

FPU (i387) で処理できないので，計算に時間がかかる．
数値の入出力に printf(), scanf() 等が使えず，不自由である．

この 4 倍精度は，現時点では「実装の試用段階」だと思います． "__float128" という名前からして「仮設」の感がありますし，printf(), scanf() の対象外という点からもその印象が強まります．

それと，拡張倍精度との共存がややこしいです．私的な願望ですが，将来的には拡張倍精度としての long double は廃止して，__float128 の方を long double として実装して欲しいです． (要は，SPARK のようにして欲しい．[2])

4 倍精度の算術関数を FPU で処理してくれる CPU が比較的安価で出てくれればうれしいのですが，需要は限定されるのでしょう．私自身，それを常用するわけではありませんし． (算術関数がどうなっていたのかは分かりませんが，昔の SPARC の FPU は 4 倍精度が扱えたらしいです． alpha でも扱えたらしいです．)

[1] gcc の資料には次のように書かれています．(抜粋)

supports additional floating types, __float80 and __float128 to support 80bit (XFmode) and 128 bit (TFmode) floating types.

XFmode : "Extended Floating" mode represents an IEEE extended floating point number. This mode only has 80 meaningful bits (ten bytes). Some processors require such numbers to be padded to twelve bytes, others to sixteen; this mode is used for either.

TFmode : "Tetra Floating" mode represents a sixteen byte floating point number all 128 of whose bits are meaningful. One common use is the IEEE quad-precision form at.

[2] SPARK は使ったことがないので，あくまで資料からの想像です．
また，FreeBSD/arm64 でもそうなっているらしいです．

2. 4 倍精度の利用

以下に，これを書いている時点で最新の gcc 4.7 (ports/lang/gcc47) を使って例示します．

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

quadmath.h を include する．
float, double を使わないならば，math.h を include する必要はない．
__complex128 に関連した関数の宣言は quadmath.h の中にあるので，complex.h は不要である．
変数を __float128 で型宣言する．
定数を __float128 で扱う場合，数値の後に Q を付ける．
出力は，先ず quadmath_snprintf() を用いて数値を文字列に変換し，それを出力する．
出力の書式は %Qf, %Qe, %Qg, %Qa のように指定し，必要なら桁を付記する．
入力は，先ず文字列として読み込み，strtoflt128() を用いてそれを数値に変換する．

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

snprintf() 等とは異なり，" " の中には書式以外の要素は書けません．また，複数の変数を変換することもできません．
以下の使用例は誤りです．

quadmath_snprintf(c,80,"sin(pi/4) = %.40Qf\n",y);
quadmath_snprintf(c,80,"%.40Qf %.40Qf",x,y);

数値の入力は次のように行ないます．この例では strtoflt128() の第 1 引数は文字列の定数ですが，変数でも可です．

   __float128 a;

   a = strtoflt128("3.1415926535897932384626433832795029",NULL);

私自身は，第 2 引数は NULL しか使わないでしょう．文字列として分断しているのなら strcat() すれば良いのですし．

利用できる算術関数や定数は quadmath.h を見れば解ります．

/usr/local/lib/gcc47/gcc/i386-portbld-freebsd8.2/4.7.0/include/quadmath.h

/usr/local/lib/gcc12/gcc/x86_64-portbld-freebsd12.3/12.2.0/include/quadmath.h

上記の例ならば，円周率として M_PIq を用いることができます．

また，c++ であっても，正弦関数は sinq() であって sin() ではないので要注意です．その他の算術関数も同様です．
これを sin() にしたければ，次のようにします． (/usr/include/c++/4.2/cmath が参考になります．)

inline __float128 sin( __float128 x )  {  return sinq(x);  }

この配慮はそのうち不要になるかもしれません． gcc を更新するごとに quadmath.h を確認すると良いでしょう．

自作のプログラム s1.c は以下のようにして compile します． (諸々の option は省略．)

% gcc47 s1.c -o s1 -lquadmath

c++ の場合には，#include <quadmath.h> の部分を次のようにします．これは bug に対する暫定対策であり，そのうち c と同じ記述で通用するようになるはずです．
→ gcc 4.9 ではこの bug は修正されており，普通に記述できます．

extern "C" {
   #include <quadmath.h>
}

compile は以下のようにします． (ソースは s2.C)

% g++47 s2.C -o s2 -Wl,-rpath=/usr/local/lib/gcc47 -lquadmath

option の -Wl,-rpath=/usr/local/lib/gcc47 は，libquadmath.so だけでなく，c++ に必要な libgcc_s.so や libstdc++.so の所在と優先順位を指定しています．この代わりに，以下のように環境変数を設定しても良いです．

% setenv LD_LIBRARY_PATH /usr/local/lib/gcc47

これらは状況に応じて選択して下さい．

3. 参考資料

先ずは gcc の公式ページです． (ここからの link は設定しておりません．)

http://gcc.gnu.org/onlinedocs/gcc/Floating-Types.html
http://gcc.gnu.org/onlinedocs/gccint/Machine-Modes.html
http://gcc.gnu.org/onlinedocs/libquadmath/quadmath_005fsnprintf.html  [1]

ソースを読むのは基本です．

/usr/ports/lang/gcc47/work/gcc-4.7-20120128/libquadmath/
/usr/ports/lang/gcc47/work/gcc-4.7-20120128/libgcc/fp-bit.c   ← 四則演算はこちらです．

参考までに，libm のソースはこちらにあります．

/usr/src/lib/msun/

assembler 出力． (以下の例では s1.s が得られます．)

% gcc47 -S s1.c
(以下，s1.s を抜粋)
        call    __multf3     ← 乗算
        call    sinq         ← sinq()

nm の出力．

% nm s1
(以下，出力を抜粋)
08048970 T __multf3             ← この情報からは分かりませんが，乗算は libgcc から
         U sinq@@QUADMATH_1.0   ← sinq() は libquadmath から

[1] これには「環境により printf("%Qf",1.2Q) のような使い方も可能」とあります． OS の更新の際に試してみると良いでしょう．
あるいは，libc の printf() 関連のソースや man page を見た方が速いかもしれません．

/usr/src/lib/libc/stdio/vfprintf.c
/usr/src/lib/libc/stdio/vfscanf.c

% man vfprintf
% man vfscanf

4. 補足

いつの間にか -mlong-double-128 という option が加えられました． (ややこしいですが，-m128bit-long-double と混同しないで下さい．) これを設定すると，long double で宣言した変数が __float128 と同じになります．

数値の後に L を付けた定数は 128bit の long double として扱われます．
long double で宣言した変数や定数の四則演算も 128bit 精度で行なわれます．

しかし，少なくとも FreeBSD においては以下の点で対応ができておりません．

math.h には，精度に関係なく long double 用の定数は用意されておりません．
(一方，quadmath.h には M_PIq 等の定数が用意されております．)
sinl() 等の算術関数は 128bit の long double に対応しておりません．
printf(), scanf() は 128bit の long double に対応しておりません．

なお，__float128 のための関数を 128bit の long double から使うことは可能です．定数も使えます．しかし，それらが必要なのであれば，素直に __float128 を使った方がソースが読み易いでしょう．

__float128 の利用法

1. i386 の 4 倍精度の概説

2. 4 倍精度の利用

3. 参考資料

4. 補足