I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

# Giuga numbers

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

Determine and show here the first four Giuga numbers.

Stretch

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

References

## AWK

` # syntax: GAWK -f GIUGA_NUMBER.AWKBEGIN {    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    LoopEnd Function Dim As Uinteger n = 3, c = 0, limit = 4Print "The first "; limit; " Giuga numbers are: ";Do    If isGiuga(n) Then c += 1: Print n; "  ";    n += 1Loop 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 []intvar 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]
```

## 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.1e530 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=: [: /:[email protected], */&>@{@((^ [email protected]>:)&.>/)@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. ({.~ <[email protected]:@>:@#) 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 330   giuga 41722 858   giuga 566198   giuga 624423128562 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    endend 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 foundend @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 numbersconst  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 6var  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;//1e8var  SmallPrimes: tPrimes;  T0 : Int64;  cnt:Uint32; procedure InitSmallPrimes;//get primes. #0..65535.Sieving only odd numbersconst  //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. `
@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_numbersuse 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`

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

## 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    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 = 4var giuga = []var n = 6 // can't be 2 or 4while (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)`
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, Intimport "./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 numberint  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
```