exflib Getting Started with FORTRAN90/95

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

• FORTRAN90/95 コンパイラ
• ソース・プログラム
• libexfloat.a
• exflib.F90

libexfloat.a と exflib.F90 の入手方法を説明します．

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/exflib.F90 .
% 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/exflib.F90 .
% cp exflib-x64-bin-20180620/exflib-x86_64-intelmac/libexfloat.a .

これで準備は完了です．

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

! declare.f90
PROGRAM declare

  USE exflib           ! モジュールの利用を宣言
  TYPE(exfloat) :: x   ! 多倍長数型

END PROGRAM declare

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

• Intel Compiler ifort の場合

  
  % ifort -fast -free -o declare.out -DEXFLIB_EXFLOAT_PRECISION10=100 exflib.F90 declare.f90 libexfloat.a
  % ./declare.out    ←何も起きず，すぐにプロンプトが表示されます
  %
  
  

• GCC gfortran (Linux, MacOSX) の場合

  
  % gfortran -O2 -o declare.out -DEXFLIB_EXFLOAT_PRECISION10=100 exflib.F90 declare.f90 libexfloat.a
  % ./declare.out  ←何も起きず，すぐにプロンプトが表示されます
  % 
  
  

• GCC gfortran (on MinGW)の場合

  
  % gfortran -O2 -o declare.out -DEXFLIB_EXFLOAT_PRECISION10=100 exflib.F90 declare.f90 exfloat.dll
  % ./declare.out  ←何も起きず，すぐにプロンプトが表示されます
  % 
  
  

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

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

! show_digits.f90
PROGRAM show_digits

  USE exflib
  TYPE(exfloat) :: x

  ! コンパイル時に指定した計算精度を表示
  WRITE(*,*) exflib_exfloat_precision10

  ! 実際の計算精度を表示
  WRITE(*,*) exflib_exfloat_precision

END PROGRAM show_digits

これをコンパイルして実行すると，次のようになります．

(Intel Compiler ifort の場合)
% ifort -fast -free -o digits.out -DEXFLIB_EXFLOAT_PRECISION10=100 exflib.F90 digits.f90 libexfloat.a

(GCC gfortran の場合)
% gfortran -O2 -o digits.out -DEXFLIB_EXFLOAT_PRECISION10=100 exflib.F90 digits.f90 libexfloat.a

% ./digits.out
100         ← 指定精度 exflib_exfloat_precision10 の値
115.595520  ← 計算精度 exflib_exfloat_precision の値
% 

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

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

• exfloat 型の変数
• 整数型の変数，リテラル整数値 (即値)
• 文字列によるリテラルの実数値

を代入することができます．最後の「リテラルの実数値」は， 実数値をシングル・クォーテーションで囲んだものです．
単・倍精度の変数およびリテラル値を代入することはできません．

! init.f90
PROGRAM init

  USE exflib
  TYPE(exfloat) :: x, y
  INTEGER :: 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'   ! これも有効

END PROGRAM init

シングル・クォーテーションで囲まれた文字列中の式は， 実行時に適切な精度で評価されて，結果が左辺の変数に代入されます． 文字列中に空白を含めることはできません． この文字列中では次の3つの定数が利用可能です．

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

また，配列，2次元配列(行列) の宣言と初期化も同じようにできます．

! array.f90
PROGRAM array

  USE exflib
  TYPE(exfloat) :: a(10,20), b(10)
  INTEGER :: i, j

  DO i = 1, 10
    DO j = 1, 20
      a(i,j) = i
    END DO
  END DO

  DO i = 1, 10
    b(i) = i
  END DO

END PROGRAM array

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

// dynamic_allocation.f90
PROGRAM dynamic_allocation

  USEexflib
  INTEGER, PARAMETER :: ROW = 10
  INTEGER, PARAMETER :: COL = 20
  TYPE(exfloat), ALLOCATABLE :: a(:), b(:)
  INTEGER :: i, j

  ALLOCATE( a( ROW*COL ) )

  DO i = 1, ROW
    DO j = 1, COL
      a( (i-1)*COL + (j-1) + 1 ) = i;
    END DO
  END DO
  

  ALLOCATE( b( 10000 ) )

  DO i = 1, 10000
    b(i) = i
  END DO

  ! after use
  DEALLOCATE( a )
  DEALLOCATE( b )

END PROGRAM dynamic_allocation

出力の形式は，次の3通りがあります．

• 一旦，倍精度型に変換
• 固定小数点表記
• 浮動小数点表記

出力の精度を要求しないのであれば，倍精度経由の方法が高速です．

! output.f90
PROGRAM output

  USE exflib
  TYPE(exfloat) :: x
  CHARACTER(60) :: str

  x = '#PI'

  ! 倍精度数に変換しての出力
  WRITE(*,*) DBLE( x )

  ! 固定小数点表記，10進法で小数点以下50桁
  str = exflib_format('F.50', x)
  WRITE(*,*) TRIM( str )

  ! 浮動小数点表記，10進法で整数部(1桁)も含めて50桁
  str = exflib_format('E.50', x)
  WRITE(*,*) TRIM( str )

END PROGRAM output

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

! cast.f90
PROGRAM cast

  USE exflib
  TYPE(exfloat) :: x, y, z
  CHARACTER(60) :: str

  REAL*4 :: float
  REAL*8 :: double

  y = '0.1'         ! 多倍長数として 0.1 を代入

  ! 単・倍精度を代入するには，明示的に型変換する
  x = exflib_cast( 0.1 )

  ! 精度の確認  
  z = x - y
  str = exflib_format('E.50', z)
  WRITE(*,*) TRIM( str )

  x = exflib_cast( 0.1D0 )

  ! 精度の確認  
  z = x - y
  str = exflib_format('E.50', z)
  WRITE(*,*) TRIM( str )

  ! 次も可能だが，x の精度は上と同じになる
  float  = 0.1E0
  double = 0.1D0
  x = exflib_cast( float  )   ! x は 0.1に対し，10進約7桁の精度
  x = exflib_cast( double )   ! x は 0.1に対し，10進約16桁の精度

END PROGRAM cast

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

(Intel Compiler ifort の場合)
% ifort -fast -free -o cast.out -DEXFLIB_EXFLOAT_PRECISION10=100 exflib.F90 cast.f90 libexfloat.a

(GCC gfortran の場合)
% gfortran -O2 -o cast.out -DEXFLIB_EXFLOAT_PRECISION10=100 exflib.F90 cast.f90 libexfloat.a

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

これより，単精度数 0.1E0 は，精確には

0.100000001490...

です ($0.1+1\times {10}^{-9}+5\times {10}^{-10}+\cdots$ )．一方， 倍精度数 0.1D0 は

0.100000000000000005551115123126...

で ($0.1+5\times {10}^{-18}+5\times {10}^{-19}+\cdots$ )， いずれも 0.1 より真に大きい値であることがわかります

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

! file_write.f90
PROGRAM file_write

  USE exflib
  TYPE(exfloat) :: x, y
  CHARACTER(60) :: str1, str2

  x = '0.1'
  y = '#PI'
  str1 = exflib_format('F.50', x)
  str2 = exflib_format('E.50', y)

  OPEN(8, FILE='file.txt')
  WRITE(8,*) TRIM( str1 ), ' ', TRIM( str2 )
  CLOSE(8)

END PROGRAM file_write

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

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

! file_read.f90
PROGRAM file_read

  USE exflib
  TYPE(exfloat) :: x, y
  CHARACTER(80) :: str1, str2, str

  OPEN(9, FILE='file.txt')

  ! 文字列として読み込み，その文字列を exfloat 型変数に代入
  READ(9,*) str1, str2
  x = str1
  y = str2
  
  CLOSE(9)

  ! 読み込んだ値を画面に表示して確認
  str = exflib_format('F.60', x)
  WRITE(*,*) 'x=', str
  str = exflib_format('F.60', y)
  WRITE(*,*) 'y=', str

END PROGRAM file_read

! arithmetic.f90
PROGRAM arithmetic

  USE exflib
  TYPE(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 = x ** y

END PROGRAM arithmetic

exfloat 型に対しては，以下の組込み函数を実装しています．

ABS, SQRT, SIN, COS, TAN,
EXP, POW (**), LOG, LOG10, ASIN, ACOS, ATAN,
SINH, COSH, TANH, MOD, FRACTION, INT, FLOOR, CEILING

多倍長数と倍精度数の演算は実装されていません． 多倍長数と「リテラル実数」の演算をおこなうには， 一旦，多倍長数に代入する必要があります．

! arithmetic2.f90
PROGRAM arithmetic2

  USE exflib
  TYPE(exfloat) :: x, y, z

  x = '#PI'

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

  ! 一旦，多倍長型の変数に代入する必要がある
  y = '0.1'
  z = x + y

END PROGRAM arithmetic2

! compare.f90
PROGRAM compare

  USE exflib
  TYPE(exfloat) :: x, y, r

  x = '0.1'
  r = '1e-20'
  y = x + r

  ! 多倍長数同士の比較
  IF ( x < y ) THEN
    WRITE(*,*) 'x is less than y'
  ELSE IF ( x .EQ. y ) THEN
    WRITE(*,*) 'x is equal to y'
  ELSE
    WRITE(*,*) 'x is greater than y'
  ENDIF  

  ! 単・倍精度数との比較も可能だが，注意が必要
  IF ( y < 0.1D0 ) THEN
    WRITE(*,*) 'y is less than 0.1D0'
  ELSE
    WRITE(*,*) 'y is greater than or equal to 0.1D0'
  ENDIF  

END PROGRAM compare

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

 x is less than y
 y is less than 0.1D0

なお，利用可能な比較演算子は以下のとおりです．

== (.EQ.)   /= (.NE.)   < (.LT.)   <= (.LE.)   > (.GT.)   >= (.GE.)

! fact.f90
PROGRAM fact

  USE exflib
  
  TYPE(exfloat) :: f
  INTEGER*8 :: n
  CHARACTER(500) :: str   ! 出力用

  f = 1

  DO n = 1, 200

    f = n * f

    str = exflib_format('F.3', f);
    WRITE(*,*) n, TRIM(str)

    str = exflib_format('E.20', f);
    WRITE(*,*) n, TRIM(str)

    WRITE(*,*) n, DBLE( f )

  END DO

END PROGRAM fact

200! は 100 桁を超える整数なので，例えば 計算精度を 500 桁にするほうがよいでしょう．

(Intel Compiler ifort の場合)
% ifort -fast -free -o fact.out -DEXFLIB_EXFLOAT_PRECISION10=500 exflib.F90 fact.f90 libexfloat.a

(GCC gfortran の場合)
% gfortran -O2 -o fact.out -DEXFLIB_EXFLOAT_PRECISION10=500 exflib.F90 fact.f90 libexfloat.a

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

200! がおよそ $7.89\times {10}^{374}$ であることがわかります． これは倍精度型では表現可能な大きさを越えており， DBLE で変換すると Infinity と表示されます．

! newton.f90
! ------------------------------
! Find zero of f(x)=0
! ------------------------------
FUNCTION f(x)
  USE exflib
  TYPE(exfloat), INTENT(IN) :: x
  TYPE(exfloat) :: f

  f = x**2 - 2

END FUNCTION f

! ------------------------------
! Derivative of f(x)
! ------------------------------
FUNCTION df(x)
  USE exflib
  TYPE(exfloat), INTENT(IN) :: x
  TYPE(exfloat) :: df

  df = 2*x

END FUNCTION df

! ------------------------------
! Newton iteration
! ------------------------------
PROGRAM newton
  USE exflib
  INTEGER, PARAMETER :: MAX_ITR = 20
  REAL*8, PARAMETER :: INITIAL_VALUE = 1.0D0
  REAL*8, PARAMETER :: TOLERANCE = 1.0D-95
  CHARACTER(120) :: str   ! 出力用

  INTEGER :: i
  TYPE(exfloat) :: f, df, x(MAX_ITR)

  ! set initial value
  x(1) = exflib_cast( INITIAL_VALUE )

  DO i = 2, MAX_ITR

    x(i) = x(i-1) - f( x(i-1) ) / df( x(i-1) )

    str = exflib_format('F.100', x(i) )
    WRITE(*,*) i, TRIM( str )

    ! if x(i) is close to zero, then exiting DO-loop
    IF ( ABS( f(x(i)) ) < TOLERANCE ) THEN
      WRITE(*,*) 'Converge!'
      EXIT
    END IF
   
  END DO

END PROGRAM newton

実行結果

% ./newton.out
2 1.5000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
3 1.4166666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667
4 1.4142156862745098039215686274509803921568627450980392156862745098039215686274509803921568627450980392
5 1.4142135623746899106262955788901349101165596221157440445849050192000543718353892683589900431576443402
6 1.4142135623730950488016896235025302436149819257761974284982894986231958242289236217849418367358303566
7 1.4142135623730950488016887242096980785696718753772340015610131331132652556303399785317871612507104752
8 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276416016
Converge!

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

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