2項係数プログラム [fortran]

こんにちは、さくらえびです。
本稿では、2項係数のプログラムをご紹介します。
2項係数では、以前ものの記事でも取り上げた階乗プログラムを使用します。
式の展開部分に出てくる係数です。
下記の展開などで出てくる係数部分(x,yの前についてくる数字)です。
\[
(x-y)^2 = x^2 – 2xy + y^2
\]
\[
(x-y)^3 = x^3 -3x^2y + 3xy^2 -y^3
\]
さて、この部分ですが、下記の数式で書けてしまいます。
\[
nC_k = \frac{n(n-1)…(n-k)}{k(k-1)…2} = \frac{n!}{k!(n-k)!}
\]
これをプログラムすれば良いので、下記のようになります。

function bin_coef(n,k)
implicit none
integer*4 n,k
double precision bin_coef,Fact
bin_coef = Fact(n)/(Fact(k)*Fact(n-k))
end function
function Fact(n)
implicit none
integer*4 n,i
double precision Fact,tmp
if (n == 0 .or. n == 1) then
Fact = 1.d0
endif
tmp = 1.d0
do i=n,1,-1
tmp = tmp*dble(i)
enddo
Fact = tmp
end function

このように、関数の中で関数を使うことで、
簡単な関数の作り込みを十分にしておけば、使いまわすことができます。
CAD方面で使われるBern-Stein基底(ベジェとか？)なんかもこれで書けたりします。
追々、Bern-Stein基底についても触れようと思います。
関数の動作テストの時間が取れず、
なかなか記事が書けてませんが、随時更新します。

by さくらえび