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:

Zitat von Aphton (Beitrag 1086923)

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:

Zitat von Aphton (Beitrag 1086920)

... 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 :-D 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:

Zitat von Aphton (Beitrag 1086925)

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.

[gammatester]

Zitat:

Zitat von Aphton (Beitrag 1086929)

... die Tau Funktion geht Blockweise bis zu N/2 durch ...

Statt N/2 ist nimmt man doch normalerweise die Wurzel aus N (oder übersehe ich hier etwas Offensichtliches?). Dies spielt bei zusammengesetzten Zahlen oft keine Rolle (wenn dynamisch abgebrochen wird), aber sehr wohl, wenn N einen sehr großen Primfaktor hat.

[Aphton]

Ne, denn die kleineste Primzahl ist ja bekanntlicherweise 2!

Edit: Mir ist übrigens eingefallen, wie man das noch verbessern könnte.
Ich erhöhe beim Primfraktorzerlegen den Int64 i immer um 1. Ich werd es so umprogrammieren, dass beim ersten Durchgang, also wenn i = 2, es wie gewohnt abläuft, aber dann, immer in zweier Schritten inkrementiere, da ja alle anderen geraden Zahlen ein vielfaches von 2 sind.
Ich werd einen weiteren Inkrementor implementieren, der, sobald i = 2, von 1 auf 2 inkrementiert wird und letztendlich am Schluss, wo inc( i ) steht, der Inkrementor drin stehen wird.

Die Werte die I bei dieser Schleife annimmt, sehen so aus (bei Primes.Min = 2)
{2, 3, 5, 7, 9, 11, 13, 15, 17, ...}

Update (Edit):
Habs nun ausgebessert. Hab die Variable "PrimIncIncrementor" eingeführt. Nun dürfte es doppelt so schnell die Primfaktorzerlegung durchführen!
Siehe Beitrag #1

[gammatester]

Zitat:

Zitat von Aphton (Beitrag 1087007)

Ne, denn die kleineste Primzahl ist ja bekanntlicherweise 2!

Versteh ich trotzdem nicht ganz. Für die Primfaktorzerlegung mußt Du doch nur alle Primzahlen bis sqrt(N) prüfen. Warum erstellst Du also (zumindest theoretisch) ein Siebfenster in der Nähe von N/2.

Mal anders: Was liefert Dein Programm für Tau(8937393460516237311) und wie lange braucht es?

Mein mit einem Primzahlgenerator kurz zusammengehacktes liefert das (von Wolfram Alpha bestätigte) Ergebnis in 1.1 s.

[Aphton]

Hmm. Es gibt, so wie es aussieht, ein paar große Probleme. Ich arbeite gerade dran.
Der Code ist so nicht funktionstüchtig! (Also nur teilweise).
Edit: So, nun dürfte er so richtig funktionieren.

Edit: @gammatester
Wow, das wird sowas von in die Knie gezwungen.

Ok, das Problem ist, dass er ewig lang iterieren muss, bis er beim entsprechenden Werten für den Sieb in CreatePrimes angelangt ist =\ Irgendwelche Ideen?

Edit: Nun, das habe ich schon vor Ewigkeiten ausgebessert. Es sind nun unzählige Fehler weggekommen; trotzdem ists momentan noch so, dass es bei sehr großen Zahlen etwas länger dauert, WENN bei der Primfaktorzerlegung der Quotient dann eine weitere Primzahl ist!
Oder einfacher ausgedrückt:
Je größer die Zahl und je kleiner der Tau() Wert dieser Zahl, desto länger braucht es.

90 ms für divisorTau( 26113434792554522 ) = 128
Wie schon gesagt, bei großen Werten, deren Tau klein ist, dauerts sehr lange.

Hilfe

[gammatester]

Zitat:

Zitat von Aphton (Beitrag 1087036)

Je größer die Zahl und je kleiner der Tau() Wert dieser Zahl, desto länger braucht es.

90 ms für divisorTau( 26113434792554522 ) = 128
Wie schon gesagt, bei großen Werten, deren Tau klein ist, dauerts sehr lange.

Hilfe

Leider habe ich keinen Schimmer, was Dein Code genau macht. Ich habe ihn laufen lassen und selbst mit sqrt(N) statt N/2 endet er nie für große N. Ich weiß nur, daß auf meinem Uralt-1.7GHz-Pentium 4 zu Hause die Zeiten für meinen Code wie folgt sind

Code:

Mit Prime-Generator:
tau(8937393460516237311) = 30
Start: 19:41:37.38     Stop: 19:41:40.07     Diff: 2.69 sec

Mit simple nextprime32
tau(8937393460516237311) = 30
Start: 19:43:07.79     Stop: 19:43:32.34     Diff: 24.55 sec

Wenn Interesse besteht, kann ich ihn auch posten.

[Aphton]

Ja mach das bitte.

[gammatester]

Das Testprogramm benutzt die Unit mp_prime aus meiner MPArith-Sammlung, bzw als Direktlink: http://home.netsurf.de/wolfgang.ehrh...2011-01-04.zip

Delphi-Quellcode:

{Testprogram zu DP-Praxis 2010-03-09: Tau Funktion von Aphton}
program t_tau;

{$i STD.INC}

{$ifdef APPCONS}
  {$apptype console}
{$endif}

uses
  mp_prime;

{$ifndef HAS_INT64}
type
  int64 = longint;
{$endif}


type
  TPrimeFac = record
                p: int64;
                e: integer;
              end;

  TFactList = array[1..64] of TPrimeFac;


{---------------------------------------------------------------------------}
procedure factor(n: int64; var pcn: integer; var FLN: TFactList);
  {-Primfaktorzerlegung von n mit Primgenerator}
var
  sieve: TSieve;
  cp: int64;
begin
  prime_sieve_init(sieve,2);
  pcn := 0;
  repeat
    cp := prime_sieve_next(sieve);
    if cp=1 then break;
    if n mod cp = 0 then begin
      {Potenzen von cp anspalten}
      inc(pcn);
      with FLN[pcn] do begin
        p := cp;
        e := 1;
        n := n div cp;
        while (n<>1) and (n mod cp = 0) do begin
          inc(e);
          n := n div cp;
        end;
      end;
    end;
  until cp*cp > n;
  if cp<=1 then begin
    writeln('Überlauf prime_sieve_next');
    halt;
  end
  else if n<>1 then begin
    {Rest n ist prim}
    inc(pcn);
    with FLN[pcn] do begin
      p := n;
      e := 1;
    end;
  end;
  prime_sieve_clear(sieve);
end;


{---------------------------------------------------------------------------}
procedure factor2(n: int64; var pcn: integer; var FLN: TFactList);
  {-Primfaktorzerlegung von n mit nextprime32}
var
  cp: int64;
begin
  pcn := 0;
  cp  := 1;
  repeat
    cp := nextprime32(cp+1);
    if cp<=1 then break;
    if n mod cp = 0 then begin
      {Potenzen von cp anspalten}
      inc(pcn);
      with FLN[pcn] do begin
        p := cp;
        e := 1;
        n := n div cp;
        while (n<>1) and (n mod cp = 0) do begin
          inc(e);
          n := n div cp;
        end;
      end;
    end;
  until cp*cp > n;
  if cp<=1 then begin
    writeln('Überlauf prime_sieve_next');
    halt;
  end
  else if n<>1 then begin
    {Rest n ist prim}
    inc(pcn);
    with FLN[pcn] do begin
      p := n;
      e := 1;
    end;
  end;
end;


{---------------------------------------------------------------------------}
function tau(n: int64): longint;
  {-Tau-Funktion = sigma0(n)}
var
  i,pcn: integer;
  fln: TFactList;
  t: longint;
begin
  factor2(n,pcn,fln);
  t := 1;
  for i:=1 to pcn do t := t*(1+fln[i].e);
  tau := t;
end;


var
  n: int64;
  t: longint;
begin
{$ifdef HAS_INT64}
  n := 8937393460516237311;
{$else}
  n := 2080899072;
{$endif}
  t := tau(n);
  writeln('tau(',n,') = ',t);
end.