Fortran 95 - ネイビア数の計算（多倍長演算ライブラリ FMLIB 使用）！

Fortran 95 で多倍長演算ライブラリ FMLIB を使用してネイピア数を計算してみました。

0. 前提条件

• LMDE 3 (Linuz Mint Debian Edition 3; 64bit) での作業を想定。
• GCC 6.3.0 でのビルド（コンパイル＆リンク）を想定。
• Fortran 用多倍長演算ライブラリ FMLIB がインストール済みであること。
  （参照：「Fortran - 多倍長演算ライブラリ FMLIB のインストール！」）

1. ネイピア数について

当ブログ過去記事を参照のこと。（但し、多倍長演算用のライブラリは使用せず、自前で実装している）

• C++ - ネイピア数（自然対数の底）e 計算！
• Ruby - ネイピア数（自然対数の底）e 計算！
• Python - ネイピア数（自然対数の底）e 計算！

2. ソースコードの作成

File: napier.f95

!****************************************************
! ネイピア数（自然対数の底 e）の計算
!
! date          name            version
! 2018.11.17    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
!****************************************************
!
program napier
  use FMZM
  implicit none
  integer, parameter  :: KETA = 1000  ! 小数点以下の桁数
  integer, parameter  :: N    =  451  ! 計算項数
  type(IM)            :: k            ! 階乗部分
  type(FM)            :: np           ! ネイビア数
  character(16)       :: fmt_s        ! 結果出力フォーマット
  character(KETA + 2) :: s            ! 結果文字列
  integer             :: i            ! ループインデックス

  call FM_SET(KETA + 2)
  k  = TO_IM('1')
  np = TO_FM('1')
  do i = 1, N
    k = k * i
    np = np + 1.0d0 / k
  end do
  !call FMPRINT(np)  ! <= FMLIB 標準の書式で出力
  write (fmt_s, '("F", I0, ".", I0)') KETA + 2, KETA
  call FM_FORM(fmt_s, np, s)
  print '(A)', s

  stop
end program napier

• Gist - Fortran 95 source code to compute Napier’s constant with FMLIB.

3. FMLIB の準備

コンパイルに使用するので、同じディレクトリ内に FMLIB のオプジェクト／モジュールファイルのリンクを貼っておく。（煩わしくないのなら、コンパイル時にフルパスで指定してもよい。また、 Makefile を作成してビルドする方法でも良いだろう）

4. コンパイル＆リンク

$ gfortran napier.f95 -c -O3
$ gfortran fmsave.o fm.o fmzm90.o napier.o -o napier

5. 実行

$ ./napier
2.71828182845904523536028747135266249775724709369995957496696762772407663035354759
  45713821785251664274274663919320030599218174135966290435729003342952605956307381
  32328627943490763233829880753195251019011573834187930702154089149934884167509244
  76146066808226480016847741185374234544243710753907774499206955170276183860626133
  13845830007520449338265602976067371132007093287091274437470472306969772093101416
  92836819025515108657463772111252389784425056953696770785449969967946864454905987
  93163688923009879312773617821542499922957635148220826989519366803318252886939849
  64651058209392398294887933203625094431173012381970684161403970198376793206832823
  76464804295311802328782509819455815301756717361332069811250996181881593041690351
  59888851934580727386673858942287922849989208680582574927961048419844436346324496
  84875602336248270419786232090021609902353043699418491463140934317381436405462531
  52096183690888707016768396424378140592714563549061303107208510383750510115747704
  1718986106873969655212671546889570350354

（実際は1行で出力される）

---

次回は、 FMLIB を使用して円周率を計算してみます。（Chudnovsky の公式を使用）

以上。