Fortran Advent Calendar 2018 がとても勉強になっているので、途切れるのが惜しいと思い、即席のネタを投下してみました。お許しを。
#浮動小数点数の精度
Fortran の実数型（REAL 型）は実際は浮動小数点数で、FORTRAN II まではマニュアルでは浮動小数点数と呼んでいましたが、FORTRAN IV で型宣言が導入されたときに、今のように REAL 型と呼ばれるようになったようです。
浮動小数点数は実数を有限桁の二進数で表わしているので実際は有理数です。実数は適当に丸められて、近隣の浮動小数点数によって代表されます。平成末年現在、内部的には単精度は 32bit のビット列、倍精度は 64bit のビット列で表現されています。
実数を write 文などで書かせると、特に指定しない限り十進数に変換されて表示されます。この時、形式的に二進表現を正確に十進表現に変換できますが、元々浮動小数点数は適当に丸めてあるので、十進表現では有効桁は６～７桁、倍精度では１４～１５桁程度から先は実数の表現としてはあまり意味がありません。したがって普通有効桁以上は書き出しても意味がないとされます。
しかしながら、書き出した十進表現を読み込んで内部二進表現を回復するという目的なら（丸めの規則が書き出しと読み込みで同一のとき）、有効桁以上に十進表現の数を書き出す意味があります。
では、どの位の桁だけ書き出せばいいのでしょうか？必要最小限の桁数はいくらなのでしょうか？それは D. Goldberg の "What every computer scientist should know about floating-point arithmetic" の付録に証明付きで示してあって、単精度では９桁とのことです。以下で Fortran で単精度で確かめてみたあと、証明をサボって倍精度でも計算してみます。
https://dl.acm.org/citation.cfm?id=103163
（ちなみにネットには機械翻訳っぽい和訳が落ちていますが日本語的に難があります。ビット別冊 (1993) により日本語な和訳があります。）
http://iss.ndl.go.jp/books/R100000039-I001404293-00
「浮動小数点演算について計算機科学者は何を知っておくべきか / David Goldberg ; 西村恕彦」
##単精度の場合
プログラムでは、乱数で生成した実数をはじめは十進８桁で、次に９桁で内部ファイル（文字列）に書き出し、それを再び同一フォーマットで別の実数変数に読み込み、ビット表現を比較します。ビット表現に違いがある時にのみ、変数の値を十進と二進の両方の表現で書き出します。これをループで複数回繰り返します。
program single
use, intrinsic :: iso_fortran_env
implicit none
real(real32) :: x, y
character(len = 40) :: fmt, buff, bit_x, bit_y
integer :: i
write(fmt, *) '(e40.8)'
call wrfmt()
write(fmt, *) '(e40.9)'
call wrfmt()
contains
subroutine wrfmt()
do i = 1, 1000
call random_number(x)
write(bit_x, '(b32.32)') transfer(x, 0_int32)
write(buff, fmt) x
read (buff, fmt) y
write(bit_y, '(b32.32)') transfer(y, 0_int32)
if (bit_x /= bit_y) then
print '(e20.8, e20.9)', x, x
print '(e20.8, e20.9)', y, y
print *, bit_x
print *, bit_y
print *
end if
end do
print *, 'end of wrfmt'
end subroutine wrfmt
end program single
実行結果
確かに８桁出力の場合、読み込んだ時の丸めのために最後の１ビットが合わないことがしばしばありますが、９桁出力の場合は元のビット表現が完全に再現できています。
0.12251496E+00 0.122514956E+00
0.12251496E+00 0.122514963E+00
00111101111110101110100100011111
00111101111110101110100100100000
0.10337310E+00 0.103373095E+00
0.10337310E+00 0.103373103E+00
00111101110100111011010101000110
00111101110100111011010101000111
0.12379099E+00 0.123790994E+00
0.12379099E+00 0.123790987E+00
00111101111111011000011000100010
00111101111111011000011000100001
0.10202292E+00 0.102022916E+00
0.10202292E+00 0.102022924E+00
00111101110100001111000101100100
00111101110100001111000101100101
0.15125697E-01 0.151256975E-01
0.15125697E-01 0.151256965E-01
00111100011101111101000111000110
00111100011101111101000111000101
0.10813946E+00 0.108139455E+00
0.10813946E+00 0.108139463E+00
00111101110111010111100000111000
00111101110111010111100000111001
0.11314926E+00 0.113149256E+00
0.11314926E+00 0.113149263E+00
00111101111001111011101011001100
00111101111001111011101011001101
end of wrfmt
end of wrfmt
続行するには何かキーを押してください . . .
##倍精度の場合
倍精度の場合１７桁書き出せば良いようです。
program double
use, intrinsic :: iso_fortran_env
implicit none
real(real64) :: x, y
character(len = 78) :: fmt, buff, bit_x, bit_y
integer :: i
write(fmt, *) '(e60.16)'
call wrfmt()
write(fmt, *) '(e60.17)'
call wrfmt()
contains
subroutine wrfmt()
do i = 1, 50
call random_number(x)
write(bit_x, '(b64.64)') transfer(x, 0_int64)
write(buff, fmt) x
read (buff, fmt) y
write(bit_y, '(b64.64)') transfer(y, 0_int64)
if (bit_x /= bit_y) then
print '(e39.16, e39.17)', x, x
print '(e39.16, e39.17)', y, y
print *, bit_x
print *, bit_y
print *
end if
end do
print *, 'end of wrfmt'
end subroutine wrfmt
end program double
出力結果
0.2548044275764261E-01 0.25480442757642607E-01
0.2548044275764261E-01 0.25480442757642611E-01
0011111110011010000101111000101110010001010100111010001010100000
0011111110011010000101111000101110010001010100111010001010100001
0.3001758179231283E+00 0.30017581792312825E+00
0.3001758179231283E+00 0.30017581792312831E+00
0011111111010011001101100001010010100010010000011110100110101001
0011111111010011001101100001010010100010010000011110100110101010
0.4031338096905369E-01 0.40313380969053694E-01
0.4031338096905369E-01 0.40313380969053687E-01
0011111110100100101000111111010010011001101101001110000001100010
0011111110100100101000111111010010011001101101001110000001100001
0.1210777067074573E+00 0.12107770670745734E+00
0.1210777067074573E+00 0.12107770670745729E+00
0011111110111110111111101111001011010110100101010100110100110001
0011111110111110111111101111001011010110100101010100110100101110
0.1966502222769279E-01 0.19665022227692792E-01
0.1966502222769279E-01 0.19665022227692789E-01
0011111110010100001000110001000101001101010111110100100101110010
0011111110010100001000110001000101001101010111110100100101110001
0.1698636521764145E+00 0.16986365217641447E+00
0.1698636521764145E+00 0.16986365217641450E+00
0011111111000101101111100001011110010111011100000011101110011110
0011111111000101101111100001011110010111011100000011101110011111
0.4800633698726932E+00 0.48006336987269316E+00
0.4800633698726932E+00 0.48006336987269321E+00
0011111111011110101110010101101110110110011001110001011011010110
0011111111011110101110010101101110110110011001110001011011010111
0.2542020527679349E-01 0.25420205276793492E-01
0.2542020527679349E-01 0.25420205276793489E-01
0011111110011010000001111100000100011001010010010010011000101100
0011111110011010000001111100000100011001010010010010011000101011
0.4448168500370487E+00 0.44481685003704874E+00
0.4448168500370487E+00 0.44481685003704868E+00
0011111111011100011101111110000100010111111001111001101101101010
0011111111011100011101111110000100010111111001111001101101101001
0.4387013564303580E+00 0.43870135643035796E+00
0.4387013564303580E+00 0.43870135643035801E+00
0011111111011100000100111010111011011010101001010001001000001101
0011111111011100000100111010111011011010101001010001001000001110
end of wrfmt
end of wrfmt
続行するには何かキーを押してください . . .
#結論
結局、十進表現でも単精度なら９桁、倍精度なら１７桁で読み書きすれば、テキストファイルでもバイナリで内部表現のまま読み書きしたのと同じ精度が保たれるようです。
IEEE754 に対応した I/O ならば NaN、±Inf も読み書きできるので大丈夫だと思います。Fortran では丸めの種類をコントロールできるので、これは合わせる必要があると思いますが･･･
ニシン
https://www.zukan-bouz.com/syu/%E3%83%8B%E3%82%B7%E3%83%B3