Fortranの逆三角関数を使って円周率を求める

概要

Fortranでは，円周率が定数として定義されていません．そのため，プログラム内で円周率を利用するには，ユーザが定義しなければなりません．
逆三角関数を使って円周率πをプログラム中で定義する方法を調べました．

使用環境

コンパイラ	バージョン
intel Parallel Studio XE Composer Edition for Fortran	19.0.0.117
PGI Visual Fortran for Windows	18.7
gfortran	7.3.0 (for Ubuntu on WSL)

更新履歴

• 2019年2月9日 複数コンパイラでの検証結果を追記

方法

円周率は逆三角関数を利用すると，何通りかの方法で求められます．

\begin{eqnarray}
\pi &=& 4 \tan^{-1}1\\
    &=& \cos^{-1}-1\\
    &=& 2\sin^{-1}1
\end{eqnarray}

このうち，FORTRANでは第1式がよく使われているようです．

著者は第2式を最初に知り，また他の式と比較してもかけ算がない分簡単だと感じています．

数値実験

プログラム

これら3式の精度に違いがあるのか，次のプログラムを実行して確認してみました．ちなみに，PGI Fortranはreal128を持っていないので，real64のみを用いています．

program main
    use,intrinsic :: iso_fortran_env
    implicit none
    
    real(real128),parameter :: PI = 3.14159265358979323846264338327950288_real128

    print '(2F23.19)', 4.0_real64*atan( 1.0_real64), PI - 4.0_real64*atan( 1.0_real64)
    print '(2F23.19)',            acos(-1.0_real64), PI -            acos(-1.0_real64)
    print '(2F23.19)', 2.0_real64*asin( 1.0_real64), PI - 2.0_real64*asin( 1.0_real64)

    print '(2F40.36)', 4.0_real128*atan( 1.0_real128), PI - 4.0_real128*atan( 1.0_real128)
    print '(2F40.36)',             acos(-1.0_real128), PI -             acos(-1.0_real128)
    print '(2F40.36)', 2.0_real128*asin( 1.0_real128), PI - 2.0_real128*asin( 1.0_real128)
end program main

基準となるのは，40桁の円周率を直接real128に代入した定数です．ただし，末尾3桁ほどは打ち切られます．

実行結果

Intel Fortran

   3.14159265358979323846264338327950
  3.1415926535897931160  0.0000000000000001225
  3.1415926535897931160  0.0000000000000001225
  3.1415926535897931160  0.0000000000000001225
  3.141592653589793238462643383279502797  0.000000000000000000000000000000000000
  3.141592653589793238462643383279502797  0.000000000000000000000000000000000000
  3.141592653589793238462643383279502797  0.000000000000000000000000000000000000

PGI Fortran

    3.141592653589793
  3.1415926535897931160  0.0000000000000000000
  3.1415926535897931160  0.0000000000000000000
  3.1415926535897931160  0.0000000000000000000

gfortran

   3.14159265358979323846264338327950280
  3.1415926535897931160  0.0000000000000001225
  3.1415926535897931160  0.0000000000000001225
  3.1415926535897931160  0.0000000000000001225
  3.141592653589793238462643383279502797  0.000000000000000000000000000000000000
  3.141592653589793238462643383279502797  0.000000000000000000000000000000000000
  3.141592653589793238462643383279502797  0.000000000000000000000000000000000000

考察

実行結果を見ると，real64の場合は小数点以下15桁まで正確に得られています．そして，3式の結果に違いはありません．これはいずれのコンパイラでも同じです．

PGI Fortranだけは，誤差が0になっています．これはPGIの計算結果の精度が高いわけではなく，定数がreal64で宣言されているためです．

real128のときは，小数点以下32桁まで一致します．そして，real64と同様，3式の結果に違いはありません．

まとめ

逆三角関数を利用して円周率を求める3通りの式を確認しました．得られる値に違いはありませんでした．