Tau Funktion (+Ressourcensparend; +Erweiterter Sieb von Eratosthenes)

[Aphton]

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

[sx2008]

Sollte man folgenden Codeabschnitt nicht mit einer For-Schleife anstelle des While-Ersatzkonstrukts programmieren?

Delphi-Quellcode:

i := Primes.Min;
while i <= Primes.Max do
begin
  if Primes.PrimeState[i] = psIsPrime then
    ....

  inc( i );
end;

// also vereinfacht so
for i:=Primes.Min to Primes.Max do
begin
  if Primes.PrimeState[i] = psIsPrime then
    ....
end;

Ich sehe da noch weitere While-Schleifen, die sich mit einer For-Schleife ersetzen lassen.

[Aphton]

Das sind Int64 Werte (Min & Max)!
Edit: Int64 Variablen können in Turbo Delphi Explorer (2006) bei der For Schleife nicht verwendet werden.
Danke für die Aufmerksamkeit!

[sx2008]

Zitat von Aphton:

Das sind Int64 Werte (Min & Max)!
Edit: Int64 Variablen können in Turbo Delphi Explorer (2006) bei der For Schleife nicht verwendet werden.
Danke für die Aufmerksamkeit!

Schade. Der Code wird so schwerer zu lesen als nötig.
Wenn der Bereich von Min & Max grössere Dimensionen annimmt, dann müssste man
vielleicht noch TBoolArray als packed array deklarieren.

[Aphton]

Hmmm, wenn du mir bestätigen könntest, dass packed Array die Performance nicht verschlechtert?
Imho gibts bei packed Arrays kein Aligning, was zu nem kleineren Speicherverbrauch führt, jedoch leidet das Zugreifen auf die Werte darunter!

[BUG]

Zitat von Aphton:

... wenn ggT(a, b) = 1 dh. wenn a und b = Primzahlen, ...

4 ist keine Primzahl (2*2)
13 ist eine Primzahl (13)
ggT(4, 13) = 1

ggT(a, b) = 1 <=> a und b teilerfremd

[sx2008]

Natürlich wird der Zugriff etwas langsamer.
Ich kann jetzt nur grob schätzen, dass der Zugriff vielleicht 30% länger dauert und der gesamte Algorithmus ~10% langsamer. wird
Das müsste man testen und ich hab keinen Compiler auf'm Rechner, der Funktionen in Records mag.
Aber wenn die Werte Min & Max schon auf Werte > 2^31 ausgelegt sind, dann muss man Kompromisse eingehen.

[BUG]

Andererseits muss diese Menge an Speicher ja auch verwaltet werden.
Bei einer Array-Größe von 2^32 sind das selbst bei einem packed Array schon 4GB, da stößt du mit Delphi im Moment auf die Grenzen.
Wenn du auf Bitebene arbeiten würdest, könntest du da mit 512MB auskommen.

[Aphton]

@BUG
Danke, du hast recht. a und b sind teilfremd und müssen nicht zwingend Primzahlen sein.
Weiters - die Indices des Arrays werden getestet, ob es sich um eine Primzahl handelt. Deswegen auch Int64. Dabei muss das Array nicht > 2^32 sein.
Man kann ja Min und Max so wählen, dass die Differenz < 2^32 bzw Max Arbeitsspeicher ist.

Deshalb auch Ressourcensparend: die Tau Funktion geht Blockweise bis zu N/2 durch und belegt nie einen Speicher > 1000 Bytes (PrimeTable - Array of Boolean)

[himitsu]

Zitat von Aphton:

Hmmm, wenn du mir bestätigen könntest, dass packed Array die Performance nicht verschlechtert?

Zitat:

Delphi-Quellcode:

TBoolArray = Array of Boolean;
TPoint64 = record
  X, Y: Int64;
end;

exp : Array of TPoint64;

Bei diesen beiden Records macht packed absolut keinen Unterschied.
Also, da sich absolut nix an der Ausrichtung ändert, kann sich auch nichts verbesser/verschlechtern.

Ein Array ist immer packed (es sei denn ein enthaltener Record ist packed und eine Ausrichtung könnte das Array optimieren).
Und bei einem Record greift eine Ausrichtung nur, wenn unterschiedlich große Basistypen im Record verbaut wurden ... also bei gleich großen Basistypen kann packed nichts verändern.

PS: Array of Boolean/Byte/.../ShortInt ist ein Sonderfall, denn hier ist ein Arrayfeld jeweils 1 Byte und es ist somit immer packed, da packed die Ausrichtung temporär auf 1 heruntersetzt.