Giuga numbers: Difference between revisions

From Rosetta Code
Content added Content deleted
m (→‎alternative version: extreme cheating: limit the used factors in the specific column 3041 is highest factor at column 5)
Line 189: Line 189:
30 858 1722 66198</lang>
30 858 1722 66198</lang>


These numbers have some interesting properties but there's an issue with guaranteeing correctness of more sophisticated approaches.
These numbers have some interesting properties but there's an issue with guaranteeing correctness of more sophisticated approaches. That said, here's a translation of the pari/gp implementation on the talk page:

<lang J>divisors=: [: /:~@, */&>@{@((^ i.@>:)&.>/)@q:~&__

giuga=: {{
r=. i.0
p=. (2) 0 1} s=. 1r2,}.(2>.y-1+t=.1)$0
while. t do.
p=. p t}~ 4 p:t{p
s=. s t}~ (s{~t-1)+1%t{p
if. (1=t{s) +. 1 >: (t{s)+(y-t+1)%t{p do.
t=. t-1
elseif. t<y-3 do.
p=. p (t+1)}~ (p{~t) >. (%-.)s{~t
t=. t+1
else.
'c d'=. 2 x: s{~y-3
dc=. d-c
k=. (d^2)-dc
for_h. ({.~ <.@-:@>:@#) f=. divisors k do.
if. 0=dc|h+d do.
if. 0=dc|dkh=. d+k%h do.
py3=. p{~y-3
if. py3 < r1=. (h+d)%dc do.
if. py3 < r2=. dkh%dc do.
if. r1~:r2 do.
if. 1 p: r1 do.
if. 1 p: r2 do.
r=. r, d*r1*r2
end.
end.
end.
end.
end.
end.
end.
end.
end.
end.
r
}}</lang>

Example use:<lang J> giuga 1

giuga 2

giuga 3
30
giuga 4
1722 858
giuga 5
66198
giuga 6
24423128562 2214408306</lang>


=={{header|Julia}}==
=={{header|Julia}}==

Revision as of 16:11, 14 July 2022

Giuga numbers is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
Definition

A Giuga number is a composite number n which is such that each of its distinct prime factors f divide (n/f - 1) exactly.

All known Giuga numbers are even though it is not known for certain that there are no odd examples.

Example

30 is a Giuga number because its distinct prime factors are 2, 3 and 5 and:

  • 30/2 - 1 = 14 is divisible by 2
  • 30/3 - 1 = 9 is divisible by 3
  • 30/5 - 1 = 5 is divisible by 5


Task

Determine and show here the first four Giuga numbers.

Stretch

Determine the fifth Giuga number and any more you have the patience for.

References




C++

Based on the Go solution. Takes 26 minutes on my system (Intel Core i5 3.2GHz). <lang cpp>#include <iostream>

// Assumes n is even with exactly one factor of 2. bool is_giuga(unsigned int n) { 

   unsigned int m = n / 2;
   auto test_factor = [&m, n](unsigned int p) -> bool {
       if (m % p != 0)
           return true;
       m /= p;
       return m % p != 0 && (n / p - 1) % p == 0;
   };
   if (!test_factor(3) || !test_factor(5))
       return false;
   static constexpr unsigned int wheel[] = {4, 2, 4, 2, 4, 6, 2, 6};
   for (unsigned int p = 7, i = 0; p * p <= m; ++i) {
       if (!test_factor(p))
           return false;
       p += wheel[i & 7];
   }
   return m == 1 || (n / m - 1) % m == 0;

}

int main() { 

   std::cout << "First 5 Giuga numbers:\n";
   // n can't be 2 or divisible by 4
   for (unsigned int i = 0, n = 6; i < 5; n += 4) {
       if (is_giuga(n)) {
           std::cout << n << '\n';
           ++i;
       }
   }

}</lang>

Output:
First 5 Giuga numbers:
30
858
1722
66198
2214408306

FreeBASIC

<lang freebasic>Function isGiuga(m As Uinteger) As Boolean 

   Dim As Uinteger n = m
   Dim As Uinteger f = 2, l = Sqr(n)
   Do
       If n Mod f = 0 Then
           If ((m / f) - 1) Mod f <> 0 Then Return False
           n /= f
           If f > n Then Return True
       Else    
           f += 1
           If f > l Then Return False
       End If
   Loop

End Function

Dim As Uinteger n = 3, c = 0, limit = 4 Print "The first "; limit; " Giuga numbers are: "; Do 

   If isGiuga(n) Then c += 1: Print n; "  ";
   n += 1

Loop Until c = limit</lang>

Output:
The first 4 Giuga numbers are: 30  858  1722  66198

Go

Translation of: Wren

I thought I'd see how long it would take to 'brute force' the fifth Giuga number and the answer (without using parallelization, Core i7) is about 1 hour 38 minutes. <lang go>package main

import "fmt"

var factors []int var inc = []int{4, 2, 4, 2, 4, 6, 2, 6}

// Assumes n is even with exactly one factor of 2. // Empties 'factors' if any other prime factor is repeated. func primeFactors(n int) { 

   factors = factors[:0]
   factors = append(factors, 2)
   last := 2
   n /= 2
   for n%3 == 0 {
       if last == 3 {
           factors = factors[:0]
           return
       }
       last = 3
       factors = append(factors, 3)
       n /= 3
   }
   for n%5 == 0 {
       if last == 5 {
           factors = factors[:0]
           return
       }
       last = 5
       factors = append(factors, 5)
       n /= 5
   }
   for k, i := 7, 0; k*k <= n; {
       if n%k == 0 {
           if last == k {
               factors = factors[:0]
               return
           }
           last = k
           factors = append(factors, k)
           n /= k
       } else {
           k += inc[i]
           i = (i + 1) % 8
       }
   }
   if n > 1 {
       factors = append(factors, n)
   }

}

func main() { 

   const limit = 5
   var giuga []int
   // n can't be 2 or divisible by 4
   for n := 6; len(giuga) < limit; n += 4 {
       primeFactors(n)
       // can't be prime or semi-prime
       if len(factors) > 2 {
           isGiuga := true
           for _, f := range factors {
               if (n/f-1)%f != 0 {
                   isGiuga = false
                   break
               }
           }
           if isGiuga {
               giuga = append(giuga, n)
           }
       }
   }
   fmt.Println("The first", limit, "Giuga numbers are:")
   fmt.Println(giuga)

}</lang>

Output:
The first 5 Giuga numbers are:
[30 858 1722 66198 2214408306]

J

We can brute force this task building a test for giuga numbers and checking the first hundred thousand integers (which takes a small fraction of a second):

<lang J>giguaP=: {{ (1<y)*(-.1 p:y)**/(=<.) y ((_1+%)%]) q: y }}"0</lang>

<lang J> 1+I.giguaP 1+i.1e5 30 858 1722 66198</lang>

These numbers have some interesting properties but there's an issue with guaranteeing correctness of more sophisticated approaches. That said, here's a translation of the pari/gp implementation on the talk page:

<lang J>divisors=: [: /:~@, */&>@{@((^ i.@>:)&.>/)@q:~&__

giuga=: {{ 

 r=. i.0
 p=. (2) 0 1} s=. 1r2,}.(2>.y-1+t=.1)$0
 while. t do.
   p=. p t}~ 4 p:t{p
   s=. s t}~ (s{~t-1)+1%t{p
   if. (1=t{s) +. 1 >: (t{s)+(y-t+1)%t{p do.
     t=. t-1
   elseif. t<y-3 do.
     p=. p (t+1)}~ (p{~t) >. (%-.)s{~t
     t=. t+1
   else.
     'c d'=. 2 x: s{~y-3
     dc=. d-c
     k=. (d^2)-dc
     for_h. ({.~ <.@-:@>:@#) f=. divisors k do.
       if. 0=dc|h+d do.
         if. 0=dc|dkh=. d+k%h do.
           py3=. p{~y-3
           if. py3 < r1=. (h+d)%dc do.
             if. py3 < r2=. dkh%dc do.
               if. r1~:r2 do.
                 if. 1 p: r1 do.
                   if. 1 p: r2 do.
                     r=. r, d*r1*r2
                   end.
                 end.
               end.
             end.
           end.
         end.
       end.
     end.
   end.
 end.
 r

}}</lang>

Example use:<lang J> giuga 1 

  giuga 2

  giuga 3

30 

  giuga 4

1722 858 

  giuga 5

66198 

  giuga 6

24423128562 2214408306</lang>

Julia

<lang ruby>using Primes

isGiuga(n) = all(f -> f != n && rem(n ÷ f - 1, f) == 0, factor(Vector, n))

function getGiuga(N) 

   gcount = 0
   for i in 4:typemax(Int)
       if isGiuga(i)
           println(i)
           (gcount += 1) >= N && break
       end
   end

end

getGiuga(4)

</lang>

Output:
30      
858 
1722
66198

Ad hoc faster version

<lang ruby>using Primes

function getgiugas(numberwanted, verbose = true) 

   n, found, nfound = 6, Int[], 0
   starttime = time()
   while nfound < numberwanted
       if n % 5 == 0 || n % 7 == 0 || n % 11 == 0
           for (p, e) in eachfactor(n)
               (e != 1 || rem(n ÷ p - 1, p) != 0) && @goto nextnumber
           end
           verbose && println(n, "  (elapsed: ", time() - starttime, ")")
           push!(found, n)
           nfound += 1
       end
       @label nextnumber
       n += 6 # all mult of 6
   end
   return found

end

@time getgiugas(2, false) @time getgiugas(6)

</lang>

Output:
30  (elapsed: 0.0)
858  (elapsed: 0.0)
1722  (elapsed: 0.0)
66198  (elapsed: 0.0009999275207519531)
2214408306  (elapsed: 18.97099995613098)
24423128562  (elapsed: 432.06500005722046)
432.066249 seconds (235 allocations: 12.523 KiB)

Pascal

Free Pascal

OK.Cheating to find square free numbers like julia in distance 6
That means always factors 2,3 and minimum one of 5,7,11.

<lang pascal>program Giuga;

{$IFDEF FPC} 

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

{$ELSE} 

 {$APPTYPE CONSOLE}

{$ENDIF} uses 

 sysutils

{$IFDEF WINDOWS},Windows{$ENDIF} 

 ;

//###################################################################### //prime decomposition only squarefree and multiple of 6

type 

 tprimeFac = packed record
                pfpotPrimIdx : array[0..9] of Uint64;
                pfMaxIdx  : Uint32;
              end;
 tpPrimeFac = ^tprimeFac;
 tPrimes = array[0..65535] of Uint32;

var 

 {$ALIGN 8}
 SmallPrimes: tPrimes;
 {$ALIGN 32}

procedure InitSmallPrimes; //get primes. #0..65535.Sieving only odd numbers const 

 MAXLIMIT = (821641-1) shr 1;

var 

 pr : array[0..MAXLIMIT] of byte;
 p,j,d,flipflop :NativeUInt;

Begin 

 SmallPrimes[0] := 2;
 fillchar(pr[0],SizeOf(pr),#0);
 p := 0;
 repeat
   repeat
     p +=1
   until pr[p]= 0;
   j := (p+1)*p*2;
   if j>MAXLIMIT then
     BREAK;
   d := 2*p+1;
   repeat
     pr[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 pr[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));

end;

function OutPots(pD:tpPrimeFac;n:NativeInt):Ansistring; var 

 s: String[31];
 chk,p: NativeInt;

Begin 

 str(n,s);
 result := s+' :';
 with pd^ do
 begin
   chk := 1;
   For n := 0 to pfMaxIdx-1 do
   Begin
     if n>0 then
       result += '*';
     p := pfpotPrimIdx[n];
     chk *= p;
     str(p,s);
     result += s;
   end;
   str(chk,s);
   result += '_chk_'+s+'<';
 end;

end;

function IsSquarefreeDecomp6(var res:tPrimeFac;n:Uint64):boolean;inline; //factorize only not prime/semiprime and squarefree n= n div 6 var 

 pr,i,q,idx :NativeUInt;

Begin 

 with res do
 Begin
   Idx := 2;

   q := n DIV 5;
   if n = 5*q then
   Begin
     pfpotPrimIdx[2] := 5;
     n := q;
     q := q div 5;
     if q*5=n then
       EXIT(false);
     inc(Idx);
   end;

   q := n DIV 7;
   if n = 7*q then
   Begin
     pfpotPrimIdx[Idx] := 7;
     n := q;
     q := q div 7;
     if q*7=n then
       EXIT(false);
     inc(Idx);
   end;

   q := n DIV 11;
   if n = 11*q then
   Begin
     pfpotPrimIdx[Idx] := 11;
     n := q;
     q := q div 11;
     if q*11=n then
       EXIT(false);
     inc(Idx);
   end;

   if Idx < 3 then
     Exit(false);

   i := 5;
   while i < High(SmallPrimes) do
   begin
     pr := SmallPrimes[i];
     q := n DIV pr;
     //if n < pr*pr
     if pr > q then
       BREAK;
     if n = pr*q then
     Begin
       pfpotPrimIdx[Idx] := pr;
       n := q;
       q := n div pr;
       if pr*q = n then
         EXIT(false);
       inc(Idx);
     end;
     inc(i);
    end;
    if n <> 1 then
    begin
      pfpotPrimIdx[Idx] := n;
      inc(Idx);
    end;
    pfMaxIdx := idx;
 end;
 exit(true);

end;

function ChkGiuga(n:Uint64;pPrimeDecomp :tpPrimeFac):boolean;inline; var 

 p : Uint64;
 idx: NativeInt;

begin 

 with pPrimeDecomp^ do
 Begin
   idx := pfMaxIdx-1;
   repeat
     p := pfpotPrimIdx[idx];
     result := (((n DIV p)-1)MOD p) = 0;
     if not(result) then
       EXIT;
     dec(idx);
   until idx<0;
 end;

end;

const 

 LMT = 24423128562;//2214408306;//

var 

 PrimeDecomp :tPrimeFac;
 T0:Int64;
 n,n6 : UInt64;
 cnt:Uint32;

Begin 

 InitSmallPrimes;

 T0 := GetTickCount64;
 with PrimeDecomp do
 begin
   pfpotPrimIdx[0]:= 2;
   pfpotPrimIdx[1]:= 3;
 end;
 n := 0;
 n6 := 0;
 cnt := 0;
 repeat
   //only multibles of 6
   inc(n,6);
   inc(n6);
   //no square factor of 2
   if n AND 3 = 0 then
     continue;
   //no square factor of 3
   if n MOD 9 = 0 then
     continue;

   if IsSquarefreeDecomp6(PrimeDecomp,n6)then
     if ChkGiuga(n,@PrimeDecomp) then
     begin
       inc(cnt);
       writeln(cnt:3,'..',OutPots(@PrimeDecomp,n),'  ',(GettickCount64-T0)/1000:6:3,' s');
     end;
 until n >= LMT;
 T0 := GetTickCount64-T0;
 writeln('Found ',cnt);
 writeln('Tested til ',n,' runtime ',T0/1000:0:3,' s');
 writeln;
 writeln(OutPots(@PrimeDecomp,n));

end.</lang>

@home AMD 5600G ( 4,4 Ghz ) fpc3.2.2 -O4 -Xs:
  1..30 :2*3*5_chk_30<   0.000 s
  2..858 :2*3*11*13_chk_858<   0.000 s
  3..1722 :2*3*7*41_chk_1722<   0.000 s
  4..66198 :2*3*11*17*59_chk_66198<   0.000 s
  5..2214408306 :2*3*11*23*31*47057_chk_2214408306<  17.120 s 
  6..24423128562 :2*3*7*43*3041*4447_chk_24423128562<  450.180 s
Found 6
Tested til 24423128562 runtime 450.180 s

24423128562 :2*3*7*43*3041*4447_chk_24423128562
TIO.RUN (~2.3 Ghz )takes ~4x runtime ? ( 2214408306 DIV 2 ) in 36 secs :-(

alternative version

Generating recursive squarefree numbers of ascending primes and check those numbers.
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; { 30 = 2 * 3 * 5. 858 = 2 * 3 * 11 * 13. 1722 = 2 * 3 * 7 * 41. 66198 = 2 * 3 * 11 * 17 * 59. 2214408306 = 2 * 3 * 11 * 23 * 31 * 47057. 24423128562 = 2 * 3 * 7 * 43 * 3041 * 4447. 432749205173838 = 2 * 3 * 7 * 59 * 163 * 1381 * 775807. 14737133470010574 = 2 * 3 * 7 * 71 * 103 * 67213 * 713863. 550843391309130318 = 2 * 3 * 7 * 71 * 103 * 61559 * 29133437. 244197000982499715087866346 = 2 * 3 * 11 * 23 * 31 * 47137 * 28282147 * 3892535183. 554079914617070801288578559178 = 2 * 3 * 11 * 23 * 31 * 47059 * 2259696349 * 110725121051. 1910667181420507984555759916338506 = 2 * 3 * 7 * 43 * 1831 * 138683 * 2861051 * 1456230512169437. } {$IFDEF FPC} 

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

{$ELSE} 

 {$APPTYPE CONSOLE}

{$ENDIF} uses 

 sysutils

{$IFDEF WINDOWS},Windows{$ENDIF} 

 ;

const 

 LMT =14737133470010574;// 432749205173838;//24423128562;//2214408306;

type 

 tFac = packed record
                fMul     :Uint64;
                fPrime,
                fPrimIdx,
                fprimMaxIdx,dummy :Uint32;
                dummy2: Uint64;
              end;
 tFacs  = array[0..15] of tFac;
 tPrimes = array[0..62157] of Uint32;//775807 see factor of 432749205173838

// tPrimes = array[0..4875{14379}] of Uint32;//sqrt 24423128562 // tPrimes = array[0..1807414] of Uint32;//29133437 // tPrimes = array[0..50847533] of Uint32;// 1e9 // tPrimes = array[0..5761454] of Uint32;//1e8 var 

 SmallPrimes: tPrimes;
 T0 : Int64;
 cnt:Uint32;

procedure InitSmallPrimes; //get primes. #0..65535.Sieving only odd numbers const 

 //MAXLIMIT = (trunc(sqrt(LMT)+1)-1) shr 1+4;
 MAXLIMIT = 775807 DIV 2+1;//(trunc(sqrt(LMT)+1)-1) shr 1+4;

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);

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 

 Mul : Uint64;
 i,p : uint32;

begin 

 with F[idx-1] do
 begin
   Mul := fMul;
   i := fPrimIdx;
 end;

 while i<High(SmallPrimes) do
 begin
   inc(i);
   with F[idx] do
   begin
     if i >fprimMaxIdx then
       BREAK;
     p := SmallPrimes[i];
     if p*Mul>LMT then
       BREAK;
     fMul := p*Mul;
     fPrime := p;
     fPrimIdx := i;
     IF (Mul-1) MOD p = 0 then
       IF ChkGiuga(F,idx) then
         OutFac(F,idx);
   end;

// max 6 factors //for lmt 24e9 need 7 factors 

   if idx <5 then
     InsertNextPrimeFac(F,idx+1);
 end;

end;

var 

 {$ALIGN 32}
 Facs : tFacs;
 i : integer;

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;
 i := 2;
 //search the indices of mx found factor
 while SmallPrimes[i] < 11 do inc(i); Facs[2].fprimMaxIdx := i;
 while SmallPrimes[i] < 71 do inc(i); Facs[3].fprimMaxIdx := i;
 while SmallPrimes[i] < 3041 do inc(i); Facs[4].fprimMaxIdx := i;
 while SmallPrimes[i] < 67213 do inc(i); Facs[5].fprimMaxIdx := i;
 while SmallPrimes[i] < 775807 do inc(i); Facs[6].fprimMaxIdx := i;

{ 

 writeln('Found ',cnt,' in ',(GetTickCount64-T0)/1000:10:3,' s');
 //start with 2*3*7
 with Facs[2] do
 begin
   fMul := 2*3*7;fPrime := 7;fPrimIdx:= 3;
 end;
 InsertNextPrimeFac(Facs,3);
 //start with 2*3*11
 writeln('Found ',cnt,' in ',(GetTickCount64-T0)/1000:10:3,' s');
 with Facs[2] do
 begin
   fMul := 2*3*11;fPrime := 11;fPrimIdx:= 4;
 end;
 InsertNextPrimeFac(Facs,3);

} 

 InsertNextPrimeFac(Facs,2);
 writeln('Found ',cnt,' in ',(GetTickCount64-T0)/1000:10:3,' s');
 writeln;

end. </lang>

@TIO.RUN:
  1  2*3*5 = 30           0.000 s
  2  2*3*7*41 = 1722           0.810 s
  3  2*3*7*43*3041*4447 = 24423128562           0.871 s
  4  2*3*11*13 = 858           1.057 s
  5  2*3*11*17*59 = 66198           1.089 s
  6  2*3*11*23*31*47057 = 2214408306           1.152 s
Found 6 in      1.526 s
Real time: 1.682 s CPU share: 99.42 %

 @home: 
Limit:432749205173838
start 2*3*7
Found 0 in      0.000 s
  1  2*3*7*41 = 1722         147.206 s
  2  2*3*7*43*3041*4447 = 24423128562         163.765 s
  3  2*3*7*59*163*1381*775807 = 432749205173838         179.124 s
Found 3 in    197.002 s
start 2*3*11
  4  2*3*11*13 = 858         197.002 s
  5  2*3*11*17*59 = 66198         219.166 s
  6  2*3*11*23*31*47057 = 2214408306         244.468 s

real  5m10,271s

Limit :14737133470010574
start 2*3*7
Found 0 in      0.000 s
  1  2*3*7*41 = 1722        1330.819 s
  2  2*3*7*43*3041*4447 = 24423128562        1567.028 s
  3  2*3*7*59*163*1381*775807 = 432749205173838        1788.203 s
  4  2*3*7*71*103*67213*713863 = 14737133470010574        2051.552 s
Found 4 in   2129.801 s
start 2*3*11
  5  2*3*11*13 = 858        2129.801 s
  6  2*3*11*17*59 = 66198        2305.752 s
  7  2*3*11*23*31*47057 = 2214408306        2591.984 s
Found 7 in   3654.610 s

real  60m54,612s

Perl

<lang perl>#!/usr/bin/perl

use strict; # https://rosettacode.org/wiki/Giuga_numbers use warnings; use ntheory qw( factor forcomposites ); use List::Util qw( all );

forcomposites 

 {
 my $n = $_;
 all { ($n / $_ - 1) % $_ == 0 } factor $n and print "$n\n";
 } 4, 67000;</lang>
Output:
30
858
1722
66198

Phix

with javascript_semantics
constant limit = 4
sequence giuga = {}
integer n = 4
while length(giuga)<limit do
    sequence pf = prime_factors(n)
    for f in pf do
        if remainder(n/f-1,f) then pf={} exit end if
    end for
    if length(pf) then giuga &= n end if
    n += 2
end while
printf(1,"The first %d Giuga numbers are: %v\n",{limit,giuga})
Output:
The first 4 Giuga numbers are: {30,858,1722,66198}

Python

Translation of: FreeBASIC

<lang python>#!/usr/bin/python

from math import sqrt

def isGiuga(m): 

   n = m
   f = 2
   l = sqrt(n)
   while True:
       if n % f == 0:
           if ((m / f) - 1) % f != 0:
               return False
           n /= f
           if f > n:
               return True
       else:
           f += 1
           if f > l:
               return False


if __name__ == '__main__': 

   n = 3
   c = 0
   print("The first 4 Giuga numbers are: ")
   while c < 4:
       if isGiuga(n):
           c += 1
           print(n)
       n += 1</lang>


Raku

<lang perl6>my @primes = (3..60).grep: &is-prime;

print 'First four Giuga numbers: ';

put sort flat (2..4).map: -> $c { 

   @primes.combinations($c).map: {
       my $n = [×] 2,|$_;
       $n if all .map: { ($n / $_ - 1) %% $_ };
   }

}</lang>

Output:
First 4 Giuga numbers: 30 858 1722 66198

Wren

Simple brute force but assumes all Giuga numbers will be even, must be square-free and can't be semi-prime.

Takes only about 0.05 seconds to find the first four Giuga numbers but finding the fifth would take many hours using this approach, so I haven't bothered. <lang ecmascript>var factors = [] var inc = [4, 2, 4, 2, 4, 6, 2, 6]

// Assumes n is even with exactly one factor of 2. // Empties 'factors' if any other prime factor is repeated. var primeFactors = Fn.new { |n| 

   factors.clear()
   var last = 2
   factors.add(2)
   n = (n/2).truncate
   while (n%3 == 0) {
       if (last == 3) {
           factors.clear()
           return
       }
       last = 3
       factors.add(3)
       n = (n/3).truncate     
   }
   while (n%5 == 0) {
       if (last == 5) {
           factors.clear()
           return
       }
       last = 5
       factors.add(5)
       n = (n/5).truncate
   }
   var k = 7
   var i = 0
   while (k * k <= n) {
       if (n%k == 0) {
           if (last == k) {
               factors.clear()
               return
           }
           last = k
           factors.add(k)
           n = (n/k).truncate
       } else {
           k = k + inc[i]
           i = (i + 1) % 8
       }
   }
   if (n > 1) factors.add(n)

}

var limit = 4 var giuga = [] var n = 6 // can't be 2 or 4 while (giuga.count < limit) { 

   primeFactors.call(n)
   // can't be prime or semi-prime
   if (factors.count > 2 && factors.all { |f| (n/f - 1) % f == 0 }) {
       giuga.add(n)
   }
   n = n + 4 // can't be divisible by 4

} System.print("The first %(limit) Giuga numbers are:") System.print(giuga)</lang>

Output:
The first 4 Giuga numbers are:
[30, 858, 1722, 66198]

XPL0

<lang XPL0>func Giuga(N0); \Return 'true' if Giuga number int N0; int N, F, Q1, Q2, L; [N:= N0; F:= 2; L:= sqrt(N); loop [Q1:= N/F; 

       if rem(0) = 0 then      \found a prime factor
               [Q2:= N0/F;
               if rem((Q2-1)/F) # 0 then return false;
               N:= Q1;
               if F>N then quit;
               ]
       else    [F:= F+1;
               if F>L then return false;
               ];
       ];

return true; ];

int N, C; [N:= 3; C:= 0; loop [if Giuga(N) then 

               [IntOut(0, N);  ChOut(0, ^ );
               C:= C+1;
               if C >= 4 then quit;
               ];
       N:= N+1;
       ];

]</lang>

Output:
30 858 1722 66198
Retrieved from "https://rosettacode.org/wiki/Giuga_numbers?oldid=321429"