１６番会議室「玉石混淆みんなで作るSample蔵」に寄せられたサンプル
"階乗関数・組み合わせ関数"
数値計算を扱ったサンプルは少ないようなので...。
表引きによる高速な階乗関数とそれを利用した組み合わせ関数です。
ユニットにしてあります。動作確認は Delphi3 Pro です。
Mathユニットを使っているのは、二項分布 function BinomialDistribution
の中でPower関数を呼んでいるからで、この関数を削除すればMathユニット
は必要ありません。
(*-----------------------------------------------------------------------
Combi.pas
組合せの計算は、
nＣm = n!/(m!(n-m)!)
で計算できるが、階乗の計算はオーバーフローしやすいので、
nＣm = n!/(m!(n-m)!)
= n*(n-1)....*(n - m + 1)/m!
= n/i * (n - 1)/2 * ... * (n - k + 1)/m
というように変形してから計算したほうが、計算途中でのオーバーフローが起
こりにくく、n = 16000 くらいまではextended型の範囲内で計算可能である。
それ以上は対数を取るしかない。
組み合わせの計算は繰り返しの多い計算なので頻繁に使う場合には計算が重
くなる原因になる。階乗の計算を高速化することでこの欠点を補うことを試み
た。
階乗の計算自体繰り返しの多い計算ではあるが、どの道extended型では1754!
より大きい値は扱えない。つまり階乗関数の返す値は 0! - 1754! の1755通り
しかないのであるから、この値をあらかじめ数値テーブルに入れてしまって表
引きにすることで劇的に高速化できる。
n < 1755 のときには階乗を用いた計算でもオーバーフローしないので、
nＣm = n!/m!/(n-m)!
で高速に計算し、
n > 1754 のときには、
nＣm = n/i * (n - 1)/2 * ... * (n - k + 1)/m
で計算する。
組み合わせを用いた関数のサンプルとして、二項分布関数と超幾何分布関数
を挙げておいた。
------------------------------------------------------------------------*)
(*----------------------------------------------------------------------*)
Unit Combi;
{$N+}{$R-}
(*------------------------*)
(* INTERFACE *)
(*------------------------*)
interface
(* 階乗 *)
function Fact(k : integer) : extended;
(* 階乗の自然対数 *)
function LogFact(k : integer) : extended;
(* 組み合せ *)
function Combination(n, m : integer) : extended;
(* 組み合せの自然対数*)
function LogCombination(n, m : integer) : extended;
(* 超幾何分布 *)
function HyperGeomDist(x, m1, n1, N : integer) : extended;
(* 二項分布 *)
function BinomialDistribution(n, x : integer; p : extended) : extended;
(*------------------------*)
(* IMPLEMENTATION *)
(*------------------------*)
implementation
uses Math,SysUtils;
const
(* 階乗が求められる数値の限界 1754!まではextended型に収まる *)
MaxFact = 1754;
Type
TVector = array[0..0] of extended;
PVector = ^TVector;
var
(* 1754!までの階乗の値を保持しておく数値テーブル *)
FactTable : PVector;
(*-----------------------------------------------------------------------
階乗の計算
2000.04.07
Extended型では 1754! までしか入らないので、階乗の計算値はグローバルな
配列(FactTable)で持ち、表参照をすることにした。
-----------------------------------------------------------------------*)
function Fact(k : integer) : extended;
begin
if k > MaxFact then
Raise Exception.Create('階乗は 1754! が限度です')
else
Result := FactTable^[k];
end;
(*-----------------------------------------------------------------------
階乗の計算(対数)
-----------------------------------------------------------------------*)
function LogFact(k : integer) : extended;
var
i : integer;
begin
if k <= MaxFact then
Result := Ln(Fact(k))
else begin
Result := Ln(Fact(MaxFact));
for i := MaxFact + 1 to k do Result := Result + Ln(i);
end;
end;
(*-----------------------------------------------------------------------
組合せ （n = 17000, m = 8500 くらいでオーバーフローする）
nＣm = n!/(m!(n-m)!)
= n*(n-1)....*(n - m + 1)/m!
= n/i * (n - 1)/2 * ... * (n - k + 1)/m
1997.12.23 改訂
1997.12.26 改訂
2000.04.07 改訂
Fact関数を表参照で高速化したのに伴い、
n<=1754のときは nＣm = n!/(m!(n-m)!) で計算するようにした。
------------------------------------------------------------------------*)
function Combination(n, m : integer) : extended;
var
i : integer;
begin
if n <= MaxFact then
Result := Fact(n)/Fact(m)/Fact(n-m)
else begin
if n - m < m then m := n - m;
Result := 1;
for i := 1 to m do Result := Result * (n - i + 1) / i;
end;
end;
(*-----------------------------------------------------------------------
組合せの自然対数（普通に求めるとオーバーフローするときはこちら）
ln(nＣm) = ln(n!/(m!(n-m)!))
= ln( n*(n-1)....*(n - m + 1)/m!)
= ln(n/i * (n - 1)/2 * ... * (n - k + 1)/m)
= ln(n) - ln(i) + ln(n - 1) - ln(2)..+ ln(n - k + 1) - ln(m)
1997.12.26 追加
1997.12.30 k = 0 のとき、 1 を返していたバグを 0 を返すように訂正。
2000.04.07 Fact関数の改定に伴い変更。n<16000のときはConbination()を呼ぶ。
------------------------------------------------------------------------*)
function LogCombination(n, m : integer) : extended;
var
i : longint;
begin
if n <= 16000 then
Result := Ln(Combination(n,m))
else begin
if n - m < m then m := n - m;
Result := Ln(1);
for i := 1 to m do Result := Result + Ln(n - i + 1) - Ln(i);
end;
end;
(*---------------------------------------------------------
超幾何分布
白黒合わせてＮ個の石に、黒石がn1個含まれている。
そこからm1個の石をランダムに抜き出したとき、
そこに黒石がx個含まれる確率。
次のような分割表で周辺和を固定したときに与えられて表の
ようなデータが得られる確率の求め方と一致する。
x m1 - x | m1
n1-x m2 - n1 + x | m2
----------------------+------
n1 n2 | N
H(x,m1,n1,N) = m1Ｃx * m2Ｃ(n1-x) / NＣn1
で求められる。
-----------------------------------------------------------*)
function HyperGeomDist(x, m1, n1, N : integer): extended;
begin
Result := Combination(m1, x) * Combination(N-m1, n1 - x) / Combination(N, n1);end;
(*----------------------------------------------------------------------
二項分布
確率ｐの事象がｎ回の試行中にｘ回生じる確率
-----------------------------------------------------------------------*)
function BinomialDistribution(n, x : integer; p : extended) : extended;
begin
if x > n then
Result := 0
else
Result := Combination(n, x) * Power(p, x) * Power((1 - p),(n - x));
end;
(*-----------------
ユニットの初期化
------------------*)
var
i : integer;
begin
(* FactTableのメモリ確保と値の初期化 *)
GetMem(FactTable,SizeOf(Extended)*(MaxFact+1));
FactTable^[0] := 1;
for i := 1 to MaxFact do FactTable^[i] := i * FactTable^[i-1];
end. (* of Unit Combi *)
2000/06/05(Mon) 00:47am GFD03044 ひの
Original document by ひの 氏 ID:(GFD03044)
ここにあるドキュメントは NIFTY SERVEの Delphi Users' Forum の16番会議室「玉石混淆みんなで作るSample蔵」に投稿されたサンプルです。これらのサンプルはボーランド株式会社がサポートする公式のものではありません。また、必ずしも動作が検証されているものではありません。これらのサンプルを使用したことに起因するいかなる損害も投稿者、およびフォーラムスタッフはその責めを負いません。使用者のリスクの範疇でご使用下さい。