Giuga numbers: Difference between revisions

{{out|@home AMD 5600G ( 4,4 Ghz ) fpc3.2.2 -O4 -Xs}}

</pre>

====alternative version====

Generating recurursive squarefree numbers of ascending primes and check those numbers.<BR>

2*3 are set fixed. 2*3*5 followed 2*3*7 than 2*3*11. Therefor the results are unsorted.

<lang pascal>program Giuga;

{$IFDEF FPC}

{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}

{$ELSE}

{$APPTYPE CONSOLE}

{$ENDIF}

uses

sysutils

{$IFDEF WINDOWS},Windows{$ENDIF}

;

const

LMT = 24423128562;//24423128562;//2214408306;//

type

tFac = packed record

fMul :Uint64;

fPrime,

fPrimIdx :Uint32;

end;

tFacs = array[0..31] of tFac;

tPrimes = array[0..14379] of Uint32;

// tPrimes = array[0..50847533] of Uint32;

// tPrimes = array[0..5761454] of Uint32;

var

SmallPrimes: tPrimes;

T0 : Int64;

cnt:Uint32;

procedure InitSmallPrimes;

//get primes. #0..65535.Sieving only odd numbers

const

// MAXLIMIT = (100*1000*1000-1) shr 1;

MAXLIMIT = (trunc(sqrt(LMT)+1)-1) shr 1;

var

pr : array of byte;

pPr :pByte;

p,j,d,flipflop :NativeUInt;

Begin

SmallPrimes[0] := 2;

setlength(pr,MAXLIMIT);

pPr := @pr[0];

p := 0;

repeat

repeat

p +=1

until pPr[p]= 0;

j := (p+1)*p*2;

if j>MAXLIMIT then

BREAK;

d := 2*p+1;

repeat

pPr[j] := 1;

j += d;

until j>MAXLIMIT;

until false;

SmallPrimes[1] := 3;

SmallPrimes[2] := 5;

j := 3;

d := 7;

flipflop := (2+1)-1;//7+2*2,11+2*1,13,17,19,23

p := 3;

repeat

if pPr[p] = 0 then

begin

SmallPrimes[j] := d;

inc(j);

end;

d += 2*flipflop;

p+=flipflop;

flipflop := 3-flipflop;

until (p > MAXLIMIT) OR (j>High(SmallPrimes));

setlength(pr,0);

writeln(p);

writeln(j);

end;

procedure OutFac(var F:tFacs;maxIdx:Uint32);

var

i : integer;

begin

write(cnt:3,' ');

For i := 0 to maxIdx do

write(F[i].fPrime,'*');

write(#8,' = ',F[maxIdx].fMul);

writeln(' ',(GetTickCount64-T0)/1000:10:3,' s');

//readln;

end;

function ChkGiuga(var F:tFacs;MaxIdx:Uint32):boolean;inline;

var

n : Uint64;

idx: NativeInt;

p : Uint32;

begin

n := F[MaxIdx].fMul;

idx := MaxIdx;

repeat

p := F[idx].fPrime;

result := (((n DIV p)-1)MOD p) = 0;

if not(result) then

EXIT;

dec(idx);

until idx<0;

inc(cnt);

end;

procedure InsertNextPrimeFac(var F:tFacs;idx:Uint32);

var

i,p,Mul : uint32;

begin

i := F[idx-1].fPrimIdx+1;

if i<High(SmallPrimes) then

repeat

with F[idx] do

begin

p := SmallPrimes[i];

fPrimIdx := i;

fPrime := p;

Mul := F[idx-1].fMul;

fMul := Mul*p;

if fMul>LMT then

BREAK;

IF (Mul-1) MOD p = 0 then

IF ChkGiuga(F,idx) then

OutFac(F,idx);

end;

InsertNextPrimeFac(F,idx+1);

inc(i);

until i>High(SmallPrimes);

end;

var

{$ALIGN 32}

Facs : tFacs;

Begin

InitSmallPrimes;

T0 := GetTickCount64;

with Facs[0] do

begin

fMul := 2;fPrime := 2;fPrimIdx:= 0;

end;

with Facs[1] do

begin

fMul := 2*3;fPrime := 3;fPrimIdx:= 1;

end;

InsertNextPrimeFac(Facs,2);

T0 := GetTickCount64-T0;

writeln('Found ',cnt);

writeln;

end.</lang>

{{out|@TIO.RUN}}

<pre>

78135 // limit of search prime

14380 // count of primes

1 2*3*5 = 30 0.000 s

2 2*3*7*41 = 1722 12.321 s

3 2*3*7*43*3041*4447 = 24423128562 12.536 s

4 2*3*11*13 = 858 14.444 s

5 2*3*11*17*59 = 66198 14.937 s

6 2*3*11*23*31*47057 = 2214408306 15.326 s

Found 6

Real time: 41.668 s CPU share: 99.30 %

//@home real 0m6,431s