Mega Code Archive
Compute binomial coefficient [fast]
<![CDATA[
{+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
(Schnelle) Berechnung des Binomialkoeffizienten
(Fast) binomial coefficient computing
------------------------------------------------
Der folgende Algorithmus berechnet den Binomialkoeffizienten "n über k".
Ziel war es einen schnellen Algorithmus, sowie eine Datenstruktur,
die mit riesigen Zahlen umgehen kann, zu entwerfen.
Diese Unit nutzt zwei weitere Units: Funct, MyBigInt. In Funct ist die
Funktion zur Berechnung des ggTs drin. MyBigInt ist eine Klasse, die die
Darstellung von großen Zahlen ermöglicht.
Alle drei Units können auf der folgenden Seite heruntergeladen werden:
http://www.hitstec.de/archiv.php?delphi=3
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}
//*** Code: CalcBinomKoeff *******************************************
unit CalcBinomKoeff;
interface
uses Sysutils, Funct, MyBigInt;
function CalculateBinomKoeff(n, k: Integer): string;
implementation
function CalculateBinomKoeff(n, k: Integer): string;
const
zaehler = 0;
nenner = 1;
Pos = 0;
var
m, i, j, iGgT: Integer;
arBruch: array [0..1] of array of Integer;
mbResult: TMyBigInt;
begin
// "n-k" soll kleiner als "k" sein, gilt wegen Symmetrie des BinomKoeff. im Nenner
if (n - k) > k then k := n - k;
m := n - k;
// Initialisierung der Variablen
for i := Low(arBruch) to High(arBruch) do SetLength(arBruch[i], m + 1);
// Setzen der gebliebenen Faktoren, "nur" n-k Faktoren sind wichtig
for i := 1 to m do
begin
arBruch[zaehler][i] := i + k;
arBruch[nenner][i] := i;
end;
arBruch[zaehler][Pos] := 1; // max{m+1: arBruch[zaehler][i]=1 für alle i=1...m}
arBruch[nenner][Pos] := 2; // max{m+1: arBruch[nenner][i]=1 für alle i=1...m}
while arBruch[nenner][Pos] <= m do
begin
for i := m downto arBruch[nenner][Pos] do
begin
for j := m downto arBruch[zaehler][Pos] do
begin
// Bestimmung des ggTs
iGgT := ggT(arBruch[nenner][i], arBruch[zaehler][j]);
if iGgT > 1 then
begin
// Kürzen
arBruch[zaehler][j] := Round(arBruch[zaehler][j] / iGgT);
arBruch[nenner][i] := Round(arBruch[nenner][i] / iGgT);
if arBruch[nenner][i] = 1 then break;
end;
end;
// Verschieben der Position im Zaehler
for j := arBruch[zaehler][Pos] to m do
if arBruch[zaehler][j] = 1 then arBruch[zaehler][Pos] := j + 1
else
break;
end;
// Verschieben der Position im Nenner
for i := arBruch[nenner][Pos] to m do
if arBruch[nenner][i] = 1 then arBruch[nenner][Pos] := i + 1
else
break;
end;
mbResult := TMyBigInt.Create(1);
try
// Faktoren im Zaehler aufmultiplizieren
for i := arBruch[zaehler][Pos] to m do mbResult.Multiply(mbResult, arBruch[zaehler][i]);
Result := mbResult.ToString;
finally
FreeAndNil(mbResult);
end;
end;
end.
//*** Code: Funct ****************************************************
unit Funct;
interface
function ggT(x, y: Integer): Integer;
implementation
function ggT(x, y: Integer): Integer;
begin
if y > x then Result := ggT(y, x)
else if y = 0 then Result := x
else
Result := ggT(y, x mod y);
end;
end.
//*** Code: MyBigInt *************************************************
unit MyBigInt;
interface
uses Sysutils, Math;
const
Base = 10;
type
TMyBigInt = class
private
Len: Integer;
Value: AnsiString;
procedure Trim;
procedure Shift(k: Integer);
procedure MultiplyAtom(Multiplier1: TMyBigInt; Multiplier2: Integer);
public
constructor Create(iValue: Integer = 0);
procedure Add(Addend1, Addend2: TMyBigInt);
procedure Multiply(Multiplier1, Multiplier2: TMyBigInt); overload;
procedure Multiply(Multiplier1: TMyBigInt; Multiplier2: Integer); overload;
function ToString: string;
procedure CopyFrom(mbCopy: TMyBigInt);
end;
implementation
constructor TMyBigInt.Create(iValue: Integer = 0);
var
sTmp: ShortString;
i: Integer;
begin
inherited Create;
sTmp := IntToStr(abs(iValue));
Len := Length(sTmp);
SetLength(Value, Len);
for i := 1 to Len do Value[i] := Chr(StrToInt(sTmp[Len - i + 1]));
end;
procedure TMyBigInt.Add(Addend1, Addend2: TMyBigInt);
{ zwei TMyBigInt miteinander addieren }
var
i, iCarry, iTemp: Integer;
begin
// Länge der Wert-Strings angleichen
iTemp := max(Addend1.Len, Addend2.Len);
SetLength(Value, iTemp);
for i := Len + 1 to iTemp do Value[i] := #0; // Für den Fall Addend1/Addend2=Self
Len := iTemp;
// Berechnung von Übertrag und Summe
iCarry := 0;
for i := 1 to Len do
begin
iTemp := iCarry;
if i <= Addend1.Len then iTemp := iTemp + Ord(Addend1.Value[i]);
if i <= Addend2.Len then iTemp := iTemp + Ord(Addend2.Value[i]);
Value[i] := Char(iTemp mod Base);
iCarry := iTemp div Base;
end;
if iCarry > 0 then
begin
Len := Len + 1;
SetLength(Value, Len);
Value[Len] := Char(iCarry);
end;
end;
procedure TMyBigInt.Multiply(Multiplier1, Multiplier2: TMyBigInt);
{ zwei TMyBigInt miteinander multipliziren }
var
mbResult, mbTemp: TMyBigInt;
i: Integer;
begin
mbResult := TMyBigInt.Create;
try
mbTemp := TMyBigInt.Create;
try
for i := 1 to Multiplier2.Len do
begin
// Multiplizieren nach der "Schulmethode"
mbTemp.MultiplyAtom(Multiplier1, Ord(Multiplier2.Value[i]));
mbTemp.Shift(i - 1);
mbResult.Add(mbResult, mbTemp);
end;
finally
FreeAndNil(mbTemp);
end;
CopyFrom(mbResult);
finally
FreeAndNil(mbResult);
end;
end;
procedure TMyBigInt.Multiply(Multiplier1: TMyBigInt; Multiplier2: Integer);
{ TMyBigInt und einen Integer multiplizieren }
var
mbTemp: TMyBigInt;
begin
mbTemp := TMyBigInt.Create(Multiplier2);
try
Multiply(Multiplier1, mbTemp);
finally
end;
end;
function TMyBigInt.ToString: string;
{ Zahl in einen String umwandeln }
var
i: Integer;
begin
Trim;
Result := '';
for i := Len downto 1 do Result := Result + IntToStr(Ord(Value[i]));
end;
procedure TMyBigInt.CopyFrom(mbCopy: TMyBigInt);
{ von mbCopy kopieren }
begin
Value := mbCopy.Value;
Len := mbCopy.Len;
end;
procedure TMyBigInt.Trim;
{ führende Nullen entfernen }
var
i, p: Integer;
begin
p := Len;
for i := Len downto 1 do
begin
if not (Value[i] in ['0']) then break;
p := i - 1;
end;
if p < Len then
begin
SetLength(Value, p);
Len := p;
end;
end;
procedure TMyBigInt.Shift(k: Integer);
{ von hinten mit k Nullen auffüllen, also mit Base^k multiplizieren }
var
i: Integer;
begin
if k = 0 then Exit;
SetLength(Value, Len + k);
for i := Len downto 1 do Value[i + k] := Value[i];
for i := 1 to k do Value[i] := #0;
Len := Len + k;
end;
procedure TMyBigInt.MultiplyAtom(Multiplier1: TMyBigInt; Multiplier2: Integer);
{ Multiplikation mit einer Ziffer }
var
i, iCarry, iTemp: Integer;
begin
// Multiplikation mit 1
if Multiplier2 = 1 then
begin
CopyFrom(Multiplier1);
Exit;
end;
SetLength(Value, Multiplier1.Len);
Len := Multiplier1.Len;
iCarry := 0;
for i := 1 to Len do
begin
iTemp := Ord(Multiplier1.Value[i]) * Multiplier2 + iCarry;
Value[i] := Char(iTemp mod Base);
iCarry := iTemp div Base;
end;
if iCarry > 0 then
begin
Len := Len + 1;
SetLength(Value, Len);
Value[Len] := Char(iCarry);
end;
end;
end.
]]>