# Giuga numbers

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

Giuga numbers
You are encouraged to solve this task according to the task description, using any language you may know.
Definition

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

Determine and show here the first four Giuga numbers.

Stretch

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

References

## 11l

Translation of: Python
```F isGiuga(m)
V n = m
V f = 2
V l = sqrt(n)
L
I n % f == 0
I ((m I/ f) - 1) % f != 0
R 0B
n I/= f
I f > n
R 1B
E
f++
I f > l
R 0B

V n = 3
V c = 0
print(‘The first 4 Giuga numbers are:’)
L c < 4
I isGiuga(n)
c++
print(n)
n++```
Output:
```The first 4 Giuga numbers are:
30
858
1722
66198
```

## Arturo

```giuga?: function [n]->
and? -> not? prime? n
-> every? factors.prime n 'f
-> zero? (dec n/f) % f

print.lines select.first:4 1..∞ => giuga?```
Output:
```30
858
1722
66198```

## AWK

```# syntax: GAWK -f GIUGA_NUMBER.AWK
BEGIN {
n = 3
stop = 4
printf("Giuga numbers 1-%d:",stop)
while (count < stop) {
if (is_giuga(n)) {
count++
printf(" %d",n)
}
n++
}
printf("\n")
exit(0)
}
function is_giuga(m,  f,l,n) {
n = m
f = 2
l = sqrt(n)
while (1) {
if (n % f == 0) {
if (((m / f) - 1) % f != 0) { return(0) }
n /= f
if (f > n) { return(1) }
}
else {
if (++f > l) { return(0) }
}
}
}
```
Output:
```Giuga numbers 1-4: 30 858 1722 66198
```

## C++

### Brute force

Based on the Go solution. Takes 26 minutes on my system (Intel Core i5 3.2GHz).

```#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;
}
}
}
```
Output:
```First 5 Giuga numbers:
30
858
1722
66198
2214408306
```

### Faster version

Library: Boost
Translation of: Wren
```#include <boost/rational.hpp>

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <vector>

using rational = boost::rational<uint64_t>;

bool is_prime(uint64_t n) {
if (n < 2)
return false;
if (n % 2 == 0)
return n == 2;
if (n % 3 == 0)
return n == 3;
for (uint64_t p = 5; p * p <= n; p += 4) {
if (n % p == 0)
return false;
p += 2;
if (n % p == 0)
return false;
}
return true;
}

uint64_t next_prime(uint64_t n) {
while (!is_prime(n))
++n;
return n;
}

std::vector<uint64_t> divisors(uint64_t n) {
std::vector<uint64_t> div{1};
for (uint64_t i = 2; i * i <= n; ++i) {
if (n % i == 0) {
div.push_back(i);
if (i * i != n)
div.push_back(n / i);
}
}
div.push_back(n);
sort(div.begin(), div.end());
return div;
}

void giuga_numbers(uint64_t n) {
std::cout << "n = " << n << ":";
std::vector<uint64_t> p(n, 0);
std::vector<rational> s(n, 0);
p[2] = 2;
p[1] = 2;
s[1] = rational(1, 2);
for (uint64_t t = 2; t > 1;) {
p[t] = next_prime(p[t] + 1);
s[t] = s[t - 1] + rational(1, p[t]);
if (s[t] == 1 || s[t] + rational(n - t, p[t]) <= rational(1)) {
--t;
} else if (t < n - 2) {
++t;
uint64_t c = s[t - 1].numerator();
uint64_t d = s[t - 1].denominator();
p[t] = std::max(p[t - 1], c / (d - c));
} else {
uint64_t c = s[n - 2].numerator();
uint64_t d = s[n - 2].denominator();
uint64_t k = d * d + c - d;
auto div = divisors(k);
uint64_t count = (div.size() + 1) / 2;
for (uint64_t i = 0; i < count; ++i) {
uint64_t h = div[i];
if ((h + d) % (d - c) == 0 && (k / h + d) % (d - c) == 0) {
uint64_t r1 = (h + d) / (d - c);
uint64_t r2 = (k / h + d) / (d - c);
if (r1 > p[n - 2] && r2 > p[n - 2] && r1 != r2 &&
is_prime(r1) && is_prime(r2)) {
std::cout << ' ' << d * r1 * r2;
}
}
}
}
}
std::cout << '\n';
}

int main() {
for (uint64_t n = 3; n < 7; ++n)
giuga_numbers(n);
}
```
Output:
```n = 3: 30
n = 4: 1722 858
n = 5: 66198
n = 6: 24423128562 2214408306
```

## 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
```
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.

```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)
}
```
Output:
```The first 5 Giuga numbers are:
[30 858 1722 66198 2214408306]
```

```--for obvious theoretical reasons the smallest divisor of a number bare 1
--must be prime
primeFactors :: Int -> [Int]
primeFactors n = snd \$ until ( (== 1) . fst ) step (n , [] )
where
step :: (Int , [Int] ) -> (Int , [Int] )
step (n , li) = ( div n h , li ++ [h] )
where
h :: Int
h = head \$ tail \$ divisors n --leave out 1

divisors :: Int -> [Int]
divisors n = [d | d <- [1 .. n] , mod n d == 0]

isGiuga :: Int -> Bool
isGiuga n = (divisors n /= [1,n]) && all (\i -> mod ( div n i - 1 ) i == 0 )
(primeFactors n)

solution :: [Int]
solution = take 4 \$ filter isGiuga [2..]
```
Output:
```[30,858,1722,66198]
```

## 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):

```giguaP=: {{ (1<y)*(-.1 p:y)**/(=<.) y ((_1+%)%]) q: y }}"0
```
```   1+I.giguaP 1+i.1e5
30 858 1722 66198
```

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:

```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
}}
```
Example use:
```   giuga 1

giuga 2

giuga 3
30
giuga 4
1722 858
giuga 5
66198
giuga 6
24423128562 2214408306
```

## Julia

```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)
```
Output:
```30
858
1722
66198
```

```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)
```
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.

```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.
```
@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.

```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');
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');
with Facs[2] do
begin
fMul := 2*3*7;fPrime := 7;fPrimIdx:= 3;
end;
InsertNextPrimeFac(Facs,3);
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.
```
@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

```#!/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;
```
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}
```

### Faster version

Translation of: Wren
```--
-- demo\rosetta\Giuga_number.exw
-- =============================
--
with javascript_semantics
requires("1.0.2") -- (is_prime2 tweak)

procedure giuga(integer n)
printf(1,"n = %d:",n)
sequence p = repeat(0,n),
s = repeat(0,n)
p[1..2] = 1
s[1] = {1,2}
integer t = 2
while t>1 do
integer pt = p[t]+1
p[t] = pt
pt = get_prime(pt)
integer {c,d} = s[t-1]
{c,d} = {c*pt+d,d*pt}
s[t] = {c,d}
if c=d
or c*pt+(n-t)*d<=d*pt then
t -= 1
elsif t < (n - 2) then
t += 1
{c,d} = s[t-1]
p[t] = max(p[t-1], is_prime2(floor(c/(d-c)),true))
else
{c,d} = s[n-2]
integer dmc = d-c,
k = d*d-dmc
sequence f = factors(k,1)
for i=1 to floor(length(f)/2) do
integer h = f[i]
if remainder(h+d,dmc) == 0
and remainder(k/h+d,dmc) == 0 then
integer r1 = (h + d) / dmc,
r2 = (k/h + d) / dmc,
pn2 = get_prime(p[n-2])
if  r1 > pn2
and r2 > pn2
and r1 != r2
and is_prime(r1)
and is_prime(r2) then
printf(1," %d",d * r1 * r2)
end if
end if
end for
end if
end while
printf(1,"\n")
end procedure

papply({3,4,5,6},giuga)
```
Output:
```n = 3: 30
n = 4: 1722 858
n = 5: 66198
n = 6: 24423128562 2214408306
```

You can (almost) push things a little further on 64-bit:

```if machine_bits()=64 then giuga(7) end if
```

and get

```n = 7: 432749205173838 550843391309130318 14737133470010574
```

It took about 3 minutes for those to show, but then carried on in a doomed quest for an hour before I killed it.

## Python

Translation of: FreeBASIC
```#!/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
```

## Quackery

`dpfs` is Distinct Prime Factors.

``` [ [] swap
dup times
[ dup i^ 2 + /mod
0 = if
[ nip dip
[ i^ 2 + join ]
[ dup i^ 2 + /mod
0 = iff
nip again ] ]
drop
dup 1 = if conclude ]
drop ]                    is dpfs  ( n --> [ )

[ dup dpfs
dup size 2 < iff
[ 2drop false ]
done
true unrot
witheach
[ 2dup / 1 -
swap mod 0 != if
[ dip not
conclude ] ]
drop ]                    is giuga ( n --> b )

[] 0
[ 1+ dup giuga if
[ tuck join swap ]
over size 4 = until ]
drop
echo```
Output:
`[ 30 858 1722 66198 ]`

## Raku

```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) %% \$_ };
}
}
```
Output:
`First 4 Giuga numbers: 30 858 1722 66198`

## Rust

```use prime_tools ;

fn prime_decomposition( mut number : u32) -> Vec<u32> {
let mut divisors : Vec<u32> = Vec::new( ) ;
let mut divisor : u32 = 2 ;
while number != 1 {
if number % divisor == 0 {
divisors.push( divisor ) ;
number /= divisor ;
}
else {
divisor += 1 ;
}
}
divisors
}

fn is_giuga( num : u32 ) -> bool {
let prime_factors : Vec<u32> = prime_decomposition( num ) ;
! prime_tools::is_u32_prime( num ) &&
prime_factors.into_iter( ).all( |n : u32| (num/n -1) % n == 0 )
}

fn main() {
let mut giuga_numbers : Vec<u32> = Vec::new( ) ;
let mut num : u32 = 2 ;
while giuga_numbers.len( ) != 4 {
if is_giuga( num ) {
giuga_numbers.push( num ) ;
}
num += 1 ;
}
println!("{:?}" , giuga_numbers ) ;
}
```
Output:
```[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.

```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
n = (n/2).truncate
while (n%3 == 0) {
if (last == 3) {
factors.clear()
return
}
last = 3
n = (n/3).truncate
}
while (n%5 == 0) {
if (last == 5) {
factors.clear()
return
}
last = 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
n = (n/k).truncate
} else {
k = k + inc[i]
i = (i + 1) % 8
}
}
}

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 }) {
}
n = n + 4 // can't be divisible by 4
}
System.print("The first %(limit) Giuga numbers are:")
System.print(giuga)
```
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.

```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()
}
```
Output:
```n = 3:
30

n = 4:
1722
858

n = 5:
66198

n = 6:
24423128562
2214408306
```

## 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;
];
]```
Output:
```30 858 1722 66198
```