Giuga numbers

From Rosetta Code
Task
Giuga numbers
You are encouraged to solve this task according to the task description, using any language you may know.
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




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

ALGOL 68

BEGIN # find some Giuga numbers, composites n such that all their distinct   #
      #                        prime factors f exactly divide ( n / f ) - 1  #

    # find the first four Giuga numbers                                      #
    INT g count := 0;
    FOR n FROM 2 WHILE g count < 4 DO
        INT  v        := n;
        BOOL is giuga := TRUE;
        INT  f count  := 0;
        IF NOT ODD v THEN
            # 2 is a prime factor of n                                       #
            is giuga := ( ( n OVER 2 ) - 1 ) MOD 2 = 0;
            IF is giuga THEN
                f count +:= 1;
                WHILE v MOD 2 = 0 DO v OVERAB 2 OD
            FI
        FI;
        IF is giuga THEN
            FOR f FROM 3 BY 2 WHILE f <= v AND is giuga DO
                IF v MOD f = 0 THEN
                    # have a prime factor                                    #
                    f count +:= 1;
                    is giuga := ( ( n OVER f ) - 1 ) MOD f = 0;
                    IF is giuga THEN
                        # n is still a candidate                             #
                        v OVERAB f;
                        WHILE v > 1 AND v MOD f = 0 DO v OVERAB f OD
                    FI
                FI
            OD;
            IF is giuga THEN
                # n is still a candidate, check it is not prime              #
                is giuga := f count > 1
            FI
        FI;
        IF is giuga THEN
            g count +:= 1;
            print( ( " ", whole( n, 0 ) ) )
        FI
    OD
END
Output:
 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

BASIC

BASIC256

n = 3
c = 0
limit = 4
print "The first"; limit; " Giuga numbers are: ";
do
	if isGiuga(N) then c += 1: print n; "  ";
	n += 1
until c = limit
end

function isGiuga(m)
	n = m
	f = 2
	l = sqr(n)
	while True
		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
	end while
end function
Output:
Same as FreeBASIC entry.

FreeBASIC

Function isGiuga(m As Uinteger) As Boolean
    Dim As Uinteger n = m, 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

Gambas

Public Sub Main() 
  
  Dim n As Integer = 3, c As Integer = 0, limit As Integer = 4 
  
  Print "The first "; limit; " Giuga numbers are: "; 
  Do 
    If isGiuga(n) Then 
      c += 1
      Print n; "  "; 
    Endif
    n += 1 
  Loop Until c = limit
  
End

Function isGiuga(m As Integer) As Boolean
  
  Dim n As Integer = m, f As Integer = 2, l As Integer = Sqr(n) 

  Do 
    If n Mod f = 0 Then
      Dim q As Integer = (m / f)
      If (q - 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
Output:
Same as FreeBASIC entry.

PureBasic

Procedure.b isGiuga(m.i)
  Define.i n = m,	f = 2, l = Sqr(n)
  While #True
    If Mod(n, f) = 0:
      If Mod(((m / f) - 1), f) <> 0: ProcedureReturn #False: EndIf
      n = n / f
      If f > n: ProcedureReturn #True: EndIf
    Else
      f + 1
      If f > l: ProcedureReturn #False: EndIf
    EndIf
  Wend
EndProcedure

OpenConsole()
Define.i n = 3, c = 0, limit = 4
Print("The first " + Str(limit) + " Giuga numbers are: ")
Repeat
  If isGiuga(N):
    c + 1
    Print(Str(n) + "  ")
  EndIf
  n + 1
Until c = limit

Input()
CloseConsole()
Output:
Same as FreeBASIC entry.

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

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]

Haskell

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

Ad hoc faster version

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)

Nim

Translation of: FreeBASIC
import std/math

func isGiuga(m: Natural): bool =
  var n = m
  var f = 2
  var l = int(sqrt(n.toFloat))
  while true:
    if n mod f == 0:
      if (m div f - 1) mod f != 0:
        return false
      n = n div f
      if f > n:
        return true
    else:
      inc f
      if f > l:
        return false

var n = 3
var c = 0
const Limit = 4
stdout.write "The first ", Limit, " Giuga numbers are: "
while true:
  if n.isGiuga:
    inc c
    stdout.write n, " "
    if c == Limit: break
  inc n
echo()
Output:
The first 4 Giuga numbers are: 30 858 1722 66198

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');
  //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_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

Ring

see "working..." + nl
see "The first 4 Giuga numbers are:" + nl
load "stdlibcore.ring"

Comp = []
num = 0
n = 1
while true
      n++
      if not isPrime(n) 
         Comp = []
         for p = 1 to n
             if isPrime(p) AND (n % p = 0)
                 add(Comp,p)
             ok
         next
         flag = 1
         for ind = 1 to len(Comp)
             f = Comp[ind]
             res = (n/f)- 1
             if res % f != 0
                flag = 0
                exit
             ok
         next
         if flag = 1
            see "" + n + " "
            num++
         ok
         if num = 4
            exit
         ok
      ok
end
see nl + "done..." + nl
Output:
working...
The first 4 Giuga numbers are:
30 858 
done...
working...
The first 4 Giuga numbers are:
30 858 1722 66198
done...

Ruby

require 'prime'

giuga = (1..).lazy.select do |n| 
  pd = n.prime_division
  pd.sum{|_, d| d} > 1 &&  #composite
  pd.all?{|f, _| (n/f - 1) % f == 0} 
end

p giuga.take(4).to_a
Output:
[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
    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)
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
Retrieved from "https://rosettacode.org/wiki/Giuga_numbers?oldid=343066"