Gute Nacht, liebe Programmierfreunde
Hiermit möchte die Tau(N) Funktion näher bringen.
Tau(N) berechnet die Anzahl der ganzzahligen Divisoren von N - bsp. Tau(20) = 6, denn 20 ist teilbar durch {1, 2, 4, 5, 10, 20}.
Diese Funktion dürfte den Derive Benutzern auch unter dem Namen "divisorTau()" bekannt sein; so habe ich sie auch genannt.
Mathemathischer Hintergrund:
- wenn p = Primzahl ^ k = Positiver Integer, dann gilt:
tau(p^k) = k+1
- wenn ggT(a, b) = 1 dh. wenn a und b =
Primzahlen teilfremd (danke BUG), dann gilt:
tau(ab) = tau(a) * tau(b)
Quelle
Da bei der Tau Funktion Primzahlen bzw. Primfaktorzerlegung benötigt wird, habe ich weiters noch den Sieb des Eratosthenes implementiert und das so gemacht, dass er im Wertebereich von {_From .. _Till} die Primzahlen bestimmen kann!
Delphi-Quellcode:
type
TBoolArray = Array of Boolean;
TPoint64 = record
X, Y: Int64;
end;
TNumberState = (nsIsPrime, nsIsNotPrime, nsOutOfRange);
TPrimeTable = record
private
FTable: TBoolArray;
FMin, FMax: Int64;
function GetNumberState(const Number: Int64): TNumberState;
public
procedure CreateTable(const AMin, AMax: Int64);
property Min: Int64 read FMin write FMin;
property Max: Int64 read FMax write FMax;
property NumberState[const Prime: Int64]: TNumberState read GetNumberState; default;
end;
{TPrimeTable}
function CreatePrimeTable(const Min, Max: Int64): TBoolArray;
var
Size : Int64;
i, j : Int64;
begin
Size := Max - Min + 1;
SetLength( Result, Size );
FillChar( Result[0], Length(Result), True );
i := 1;
while i <= Max div 2 do
begin
inc( i );
if (i >= Min) and (not Result[i-Min]) then
continue;
if (Min = Max) and (not Result[0]) then
Exit;
if i < Min then
j := i * Trunc(Min / i)
else
j := i * 2;
while j < Min do
inc( j, i );
while j <= Max do
begin
Result[j-Min] := False;
inc( j, i );
end;
end;
end;
procedure TPrimeTable.CreateTable(const AMin, AMax: Int64);
begin
FMin := Math.Min( AMin, AMax );
FMax := Math.Max( AMin, AMax );
FTable := CreatePrimeTable( FMin, FMax );
end;
function TPrimeTable.GetNumberState(const Number: Int64): TNumberState;
var
i: Int64;
begin
Result := nsOutOfRange;
i := Number - Min;
if (i >= 0) and (i <= Max-Min) then
if FTable[i] then
Result := nsIsPrime
else
Result := nsIsNotPrime;
end;
function divisorTau(N: Int64): Int64;
var
Primes : TPrimeTable;
exp : Array of TPoint64;
i, Index : Int64;
BlockStart: Int64;
HalfN : Int64;
NoPrimeDividerFound: Boolean;
primIndexIncrementor: Integer;
const
BlockSize = 1000; // for the table
function _GetExponent(const E: Int64): Integer;
var
i: Integer;
begin
Result := -1;
for i := 0 to High(exp) do
if exp[i].X = E then
begin
Result := i;
Exit;
end;
end;
procedure NewExp(const E: Int64; const F: Int64 = 0);
begin
SetLength( exp, Length(exp) + 1 );
Index := High(exp);
exp[Index].X := E;
exp[Index].Y := F;
end;
function IsPrime(const P: Int64): Boolean;
var
SmallTtbl: TPrimeTable;
begin
if (P >= Primes.Min) and (P <= Primes.Max) then
Result := Primes.NumberState[P] = nsIsPrime
else
begin
SmallTtbl.CreateTable( P, P );
Result := SmallTtbl.NumberState[P] = nsIsPrime;
end;
end;
procedure PF();
begin
Primes.Min := 0;
halfN := Trunc( SQRT( N ) );
SetLength( exp, 0 );
repeat
BlockStart := 2;
NoPrimeDividerFound := True;
repeat
if Primes.Min <> BlockStart then
Primes.CreateTable( BlockStart, BlockStart + Min( BlockSize, HalfN ) );
i := Primes.Min;
if i = 2 then
primIndexIncrementor := 1
else
begin
if i and 1 = 0 then
inc( i );
primIndexIncrementor := 2;
end;
while i <= Primes.Max do
begin
if Primes.NumberState[i] = nsIsPrime then
if N mod i = 0 then
begin
NoPrimeDividerFound := False;
Index := _GetExponent( i );
if Index = -1 then
NewExp( i );
repeat
inc( exp[Index].Y );
N := N div i;
until N mod i > 0;
// die folgende notwendige Zeile bringt alles zum Stocken
if (N > 1) and IsPrime( N ) then
begin
NewExp( N, 1 );
N := 1;
Exit;
end;
if i > 2 then
break;
end;
inc( i, primIndexIncrementor );
primIndexIncrementor := 2;
end;
if not NoPrimeDividerFound then
break;
inc( BlockStart, BlockSize );
until (N <= 1) or (BlockStart + BlockSize > halfN);
until NoPrimeDividerFound or (N <= 1);
end;
begin
PF;
if N > 1 then
Result := 0
else
begin
Result := 1;
i := 0;
while i < Length(exp) do
begin
Result := Result * ( exp[i].Y + 1 );
Inc( i );
end;
end;
end;
Ich habe diesen Code für die Projekt Euler Aufgabe #12 programmiert und dachte mir, dass jemand anderer auch noch Nutzen daraus ziehen könnte.
Ich garantiere keine 100%'ige Richtigkeit; eines kann ich aber sagen - die Aufgabe habe ich innerhalb von 5 Sekunden gelöst =)
Bin für Verbesserungsvorschläge und konstruktive Kritik offen.
Edit:
Mir fällt gerade ein, das man da eventuelle etwas verbessern könnte.
Die erste (äußerste) (repeat) Schleife könnte man eventuell wegbringen, wenn mir einer bestätigen könnte, dass die Primzahlen, die beim Primfaktorzerlegen zum Dividieren verwendet werden, in aufsteigender Reihenfolge drankommen und keine davon sich wiederholt.
Konkret: wenn soetwas z.B.: nicht geschehen kann:
2*2*2 * 3*3*3 * 5*5*5 *
2*2*2
Also hier wiederholt sich die Primzahlenpotenz (2^3). Wenn das nie eintritt, kann man die repeat Schleife sparen, denn sie sorgt dafür, dass die Anfangsprimzahlen wieder zum Dividieren verwendet werden. Und weiter müsste ich dann auch kein Array of TPoint64 verwenden und könnte mir auch noch die _GetExponent Funktion (ermittelt den Index) sparen.
Gute Nacht