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

[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.