![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
R で階乗関数 n! = 1 x 2 x ... x n を計算する関数として，ガンマ関数を使う方法が示されています。
n! = gamma(n+1) # +1 を忘れないこと！
なるほどです。私は今まで，こんなのを使ってました。
prod(1:n)
n が 0 だと，困るんですけどね。
参考：特定の引数値に対して関数の動作を変更するにはこうします
> temp <- function (n) ifelse(n == 0, 1, prod(1:n)) # n=0 の時は 1 を返す > temp(0) [1] 1 > temp(5) [1] 120
引数のチェックなどを含めて，以下のように発展させましょう。
fact <- function(x) ifelse (x != floor(x) | x < 0, NA, ifelse(x == 0 , 1, prod(1:x)))
> fact(-10) [1] NA > fact(0) [1] 1 > fact(1.2) [1] NA > fact(4) [1] 24 > fact(170) [1] 7.257416e+306 > fact(171) [1] Inf
参考：R で n! 用の特別な関数を用意していない真意は、gamma 関数を使う方が一般性があり、計算も早いからだと思います。prod(1:n) は当然無駄な計算をしていますし、ご提案の関数 fact はベクトル化されていないことも障害になりそうです。R の関数はそれが意味のある場合にはベクトル化するのが原則です。結論は 階乗の計算には gamma 関数を使おう です。
> fact(1:5) [1] 1 1 1 1 1 # ベクトル化されていないので警告が出る Warning message: Numerical expression has 5 elements: only the first used in: 1:x > system.time(for (i in 1:100000) fact(100)) # 計算時間の比較 [1] 50.12 0.94 52.04 0.00 0.00 > system.time(for (i in 1:100000) gamma(101)) [1] 0.42 0.03 0.47 0.00 0.00 # 100倍近く早い
＝＝＝＝＝
fact が遅いのは，いろいろチェックをしているからと思います。
> system.time(for (i in 1:100000) prod(1:100)) [1] 2.0703125 0.3203125 2.6718750 0.0000000 0.0000000
> system.time(for (i in 1:100000) gamma(101)) [1] 0.7734375 0.1093750 0.9375000 0.0000000 0.0000000
あまり差がないと言うべきでしょう
> fact <- function(x) ifelse (x != floor(x) | x < 0, NA, ifelse(x == 0 , 1, prod(1:x))) > system.time(for (i in 1:100000) fact(100)) [1] 95.710938 4.351562 106.625000 0.000000 0.000000
> fact2 <- function(x) ifelse (x != floor(x) | x < 0, NA, ifelse(x == 0 , 1, gamma(x+1))) > system.time(for (i in 1:100000) fact2(101)) [1] 93.117188 4.671875 105.640625 0.000000 0.000000
ほとんどが，引数チェックのオーバーヘッドです。
> system.time(for (i in 1:100000) fact2(11)) [1] 93.328125 4.351562 104.125000 0.000000 0.000000 > system.time(for (i in 1:100000) fact(10)) [1] 93.43750 4.59375 105.32812 0.00000 0.00000
> system.time(for (i in 1:100000) fact2(3)) [1] 93.046875 4.320312 106.046875 0.000000 0.000000 > system.time(for (i in 1:100000) fact(2)) [1] 93.750000 3.960938 104.625000 0.000000 0.000000
さらに，実際はほとんどが，関数起動のオーバーヘッドのように思えます。
注：gamma を使うのに反対しているわけではありません。
ベクトル化は以下のようにすればいいでしょうか
> fact <- function(x)
+ {
+ fact2 <- function(x) {
+ ifelse (x != floor(x) | x < 0, NA, ifelse(x == 0 , 1, prod(1:x)))
+ }
+ sapply(x, fact2)
+ }
> fact(1:10)
[1] 1 2 6 24 120 720 5040 40320 362880 3628800
参考：なるほどあらゆる計算機言語で論理判断は時間を食うのでできるだけ回避すべきだというのを思い出しました(どうしてなのか未だに理解できないでいますが)。ベクトル化の方法は参考になりました。ところで、ほとんど重箱の隅をつつく類のコメントですが、x == floor(x) は必ずしもいつも安全な整数判定条件ではありませんね。とはいえ、gamma 関数自身も最後は計算結果を整数化しているものと想像しています。また gamma 関数の計算は恐らく内部的に FORTRAN コードを使っているので、計算量だけからいえば、prod(1:n) の計算の手間とは比較になりません。
> x= <- 1-1e-16; x == floor(x); floor(x) [1] FALSE [1] 0 > x <- 1-1e-17; x == floor(x); floor(x) [1] TRUE [1] 1
R core チームの苦労が多少分かるような気がします。ところで何だか関数紹介のページらしく無くなってきたような気がしてきましたが、いかにも Wiki らしくて(?)、こんなのも面白いかも知れません。
1.9.0 では 内部で gamma(n+1) を呼び出す factorial(n) が追加されました。ただ、gamma を呼び出すだけでもオーバーヘッドはあるようで... prod(1:n) よりは速いようですが。
