結論

gamma関数を使いましょう．

integer(int32) :: n = 10
print *, gamma(dble(n + 1)) ! 3628800.0000000000

概要

Fortranで階乗を計算する方法をまとめました．

「Fortranは数値計算用のプログラミング言語だと言われているのに階乗を計算する関数がないじゃないか！」としばしば指摘されます．この記事を読めば，理由がわかります．

ついでに，他の方法で階乗を計算する方法をまとめました．

`gamma`関数を使う

Fortranには階乗を計算する関数，いうなればfactrial関数がありません．理由は簡単で，階乗の概念を複素数全体に拡張したgamma関数が実装されているからです．

ガンマ関数と階乗には，次の関係があります．

n! = \Gamma(n+1)

そのため，Fortranで整数nの階乗を計算したい場合，gamma関数の引数にn+1を渡せばよいということになります．ただし，gamma関数は引数の型が実数に制限され，戻り値も実数なので，整数を渡して整数を受け取るには型変換が必要です．

int(gamma(dble(n + 1)))

factorial関数として定義するなら，下記のようになるでしょう．

integer(int64) function factorial(n)
    implicit none
    integer(int32), intent(in) :: n

    factorial = int(gamma(dble(n + 1)), int64)
end function factorial

整数の数列を作ってproductを計算する

Fortranには，配列の全要素の積を計算するproduct関数があります．product関数と配列構成子を利用することで，比較的短い行数で階乗を計算できます．

integer(int32) :: n = 10
integer(int32) :: i
integer(int64) :: factorial

factorial = product([(i, i=1, n)])

この計算方法は，Fortranの内部手続しか用いないので，parameterとしてコンパイル時に評価できます．コンパイル時に評価され，実行時には値を参照するだけなので，Fortranで階乗の計算結果を得る最速の方法です．他言語ユーザがFortranとの実行速度の比較を行う際は，この方法を利用しないと公平ではありません．

integer(int32), parameter :: n = 10
integer(int32) :: i
integer(int64), parameter :: factorial = product([(i, i=1, n)])

カウンタ用変数iを別途宣言していますが，これはFortran 2023で[(i, integer(int32) :: i=1, n)]と書けるようになります．

再帰関数を用いる方法

よくある再帰関数を用いる方法です．負の数が渡された場合の対処は省略しています．

Fortranの関数は，標準で再帰呼び出しが許可されておらず，recursive属性を付与することで再帰呼び出しを可能にしていました．しかし，Fortran 2018では標準で再帰呼び出しが可能になります．-stand f18や-std=f0218，-f2018を付けてコンパイルすれば，recursive属性は不要です．

recursive function factorial(n) result(factorial_n)
    implicit none
    integer(int32), intent(in) :: n
    integer(int64) :: factorial_n

    if (n > 0) then
        factorial_n = n*factorial(n - 1)
        return
   end if

   factorial_n = 1
end function

まとめ

Fortranに階乗を計算する関数がないのは，より一般的なgamma関数があるため．
階乗はコンパイル時に評価できる．
Fortran 2018では，関数は標準で再帰呼び出しできる．

Fortranによる階乗の計算

結論

概要

gamma関数を使う

整数の数列を作ってproductを計算する

再帰関数を用いる方法

まとめ

`gamma`関数を使う