Giuga numbers: Difference between revisions

From Rosetta Code
Content added Content deleted
(→‎{{header|Wren}}: Added a second, very fast, version.)
Line 886: Line 886:


=={{header|Wren}}==
=={{header|Wren}}==
===Version 1 (Brute force)===
Simple brute force but assumes all Giuga numbers will be even, must be square-free and can't be semi-prime.
Simple brute force but assumes all Giuga numbers will be even, must be square-free and can't be semi-prime.


Line 954: Line 955:
The first 4 Giuga numbers are:
The first 4 Giuga numbers are:
[30, 858, 1722, 66198]
[30, 858, 1722, 66198]
</pre>
<br>
===Version 2 (Pari-GP translation)===
{{libheader|Wren-math}}
{{libheader|Wren-rat}}
This is a translation of the very fast Pari-GP code in the talk page. Only takes 0.015 seconds to find the first six Giuga numbers.
<lang ecmascript>import "./math" for Math, Int
import "./rat" for Rat

var giuga = Fn.new { |n|
System.print("n = %(n):")
var p = List.filled(n, 0)
var s = List.filled(n, null)
for (i in 0..n-2) s[i] = Rat.zero
p[2] = 2
p[1] = 2
var t = 2
s[1] = Rat.half
while (t > 1) {
p[t] = Int.isPrime(p[t] + 1) ? p[t] + 1 : Int.nextPrime(p[t] + 1)
s[t] = s[t-1] + Rat.new(1, p[t])
if (s[t] == Rat.one || s[t] + Rat.new(n - t, p[t]) <= Rat.one) {
t = t - 1
} else if (t < n - 2) {
t = t + 1
p[t] = Math.max(p[t-1], (s[t-1] / (Rat.one - s[t-1])).toFloat).floor
} else {
var c = s[n-2].num
var d = s[n-2].den
var k = d * d + c - d
var f = Int.divisors(k)
for (i in 0...((f.count + 1)/2).floor) {
var h = f[i]
if ((h + d) % (d-c) == 0 && (k/h + d) % (d - c) == 0) {
var r1 = (h + d) / (d - c)
var r2 = (k/h + d) / (d - c)
if (r1 > p[n-2] && r2 > p[n-2] && r1 != r2 && Int.isPrime(r1) && Int.isPrime(r2)) {
var w = d * r1 * r2
System.print(w)
}
}
}
}
}
}

for (n in 3..6) {
giuga.call(n)
System.print()
}</lang>

{{out}}
<pre>
n = 3:
30

n = 4:
1722
858

n = 5:
66198

n = 6:
24423128562
2214408306
</pre>
</pre>


Revision as of 17:37, 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

Version 1 (Brute force)

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]


Version 2 (Pari-GP translation)

Library: Wren-math
Library: Wren-rat

This is a translation of the very fast Pari-GP code in the talk page. Only takes 0.015 seconds to find the first six Giuga numbers. <lang ecmascript>import "./math" for Math, Int import "./rat" for Rat

var giuga = Fn.new { |n| 

   System.print("n = %(n):")
   var p = List.filled(n, 0)
   var s = List.filled(n, null)
   for (i in 0..n-2) s[i] = Rat.zero
   p[2] = 2
   p[1] = 2
   var t = 2
   s[1] = Rat.half
   while (t > 1) {
       p[t] = Int.isPrime(p[t] + 1) ? p[t] + 1 : Int.nextPrime(p[t] + 1)
       s[t] = s[t-1] + Rat.new(1, p[t])
       if (s[t] == Rat.one || s[t] + Rat.new(n - t, p[t]) <= Rat.one) {
           t = t - 1
       } else if (t < n - 2) {
           t = t + 1
           p[t] = Math.max(p[t-1], (s[t-1] / (Rat.one - s[t-1])).toFloat).floor
       } else {
           var c = s[n-2].num
           var d = s[n-2].den
           var k = d * d + c - d
           var f = Int.divisors(k)
           for (i in 0...((f.count + 1)/2).floor) {
               var h = f[i]
               if ((h + d) % (d-c) == 0 && (k/h + d) % (d - c) == 0) {
                   var r1 = (h + d) / (d - c)
                   var r2 = (k/h + d) / (d - c)
                   if (r1 > p[n-2] && r2 > p[n-2] && r1 != r2 && Int.isPrime(r1) && Int.isPrime(r2)) {
                       var w = d * r1 * r2
                       System.print(w)
                   }
               }
           }
       }
   }

}

for (n in 3..6) { 

   giuga.call(n)
   System.print()

}</lang>

Output:
n = 3:
30

n = 4:
1722
858

n = 5:
66198

n = 6:
24423128562
2214408306

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=321430"