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

Task
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


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

ABC

Based on the Algol 68 sample,

HOW TO REPORT is.giuga n:
    \ each prime factor must appear only once, e.g.: for 2:
    \ [ ( n / 2 ) - 1 ] mod 2 = 0  => n / 2 is odd => n isn't divisible by 4
    \ similarly for other primes
    PUT  1, 3, 1, floor( n / 2 ) IN f.count, f, giuga, v
    WHILE f <= v AND giuga = 1:
        IF v mod f = 0:
            PUT f.count + 1 IN f.count
            IF ( ( floor( n / f ) ) - 1 ) mod f <> 0: PUT 0 IN giuga
            PUT floor( v / f ) IN v
        PUT f + 2 IN f
    IF giuga = 1:             \ n is still a candidate, check it is not prime
        IF f.count = 1: FAIL       \ only 1 factor - it is prime so not giuga
    REPORT giuga = 1

PUT 0, -2 IN g.count, n
WHILE g.count < 4:
    PUT n + 4 IN n                          \ assume the numbers are all even
    IF is.giuga n:
        PUT g.count + 1 IN g.count
        WRITE n
Output:
30 858 1722 66198

Ada

-- Rosetta Code Task written in Ada
-- Giuga numbers
-- https://rosettacode.org/wiki/Giuga_numbers
-- Translated from the Nim solution
-- July 2024, R. B. E.

with Ada.Text_IO; use Ada.Text_IO;
with Ada.Integer_Text_IO; use Ada.Integer_Text_IO;
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;

procedure Giuga_Numbers is

  function Is_Giuga (M : Natural) return Boolean is
    N : Natural := M;
    F : Natural := 2;
    L : Integer := Integer (sqrt (Float (N)));
  begin
    loop
      if ((N mod F) = 0) then
        if ((((M / F) - 1) mod F) /= 0) then
          return False;
        end if;
        N := N / F;
        if (F > N) then
          return True;
        end if;
        else
          F := F + 1;
          if (F > L) then
            return False;
          end if;
      end if;
    end loop;
  end Is_Giuga;

  C : Natural := 0;
  N : Natural := 3;
  Limit : constant Integer := 4;

begin
  Put ("The first ");
  Put (Limit, 1);
  Put (" Giuga numbers are: ");
  loop
    if (Is_Giuga (N)) then
      C := C + 1;
      Put (N, 1);
      Put (" ");
    end if;
    N := N + 1;
    exit when C = Limit;
  end loop;
  New_Line;
end Giuga_Numbers;
Output:
The first 4 Giuga numbers are: 30 858 1722 66198

ALGOL 68

As with the Wren and other samples, assumes the first four Giuga numbers are even and also uses the fact that no prime squared can be a divisor of a Giuga number.

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                                      #
    # each prime factor must appear only once, e.g.: for 2:                  #
    # [ ( n / 2 ) - 1 ] mod 2 = 0  => n / 2 is odd => n isn't divisible by 4 #
    # similarly for other primes                                             #
    INT g count := 0;
    FOR n FROM 2 BY 4 WHILE g count < 4 DO # assume the numbers are all even #
        INT  v        := n OVER 2;
        BOOL is giuga := TRUE;
        INT  f count  := 1;
        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;
                v OVERAB f
            FI
        OD;
        IF is giuga THEN
            # n is still a candidate, check it is not prime                  #
            IF f count > 1 THEN
                g count +:= 1;
                print( ( " ", whole( n, 0 ) ) )
            FI
        FI
    OD
END
Output:
 30 858 1722 66198

ALGOL W

Translation of: 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                                      %
    % each prime factor must appear only once, e.g.: for 2:                  %
    % [ ( n / 2 ) - 1 ] mod 2 = 0  => n / 2 is odd => n isn't divisible by 4 %
    % similarly for other primes                                             %

    integer gCount, n;
    gCount :=  0;
    n      := -2;
    while begin n := n + 4;
                gCount < 4
          end
    do begin    % assume the numbers are all even                            %
        integer v, f, fCount;
        logical isGiuga;
        v       := n div 2;
        isGiuga := TRUE;
        fCount  := 1;
        f       := 1;
        while begin f := f + 2;
                    f <= v and isGiuga
              end
        do begin
            if v rem f = 0 then begin
                % have a prime factor                                        %
                fCount  := fCount + 1;
                isGiuga := ( ( n div f ) - 1 ) rem f = 0;
                v := v div f
            end if_v_rem_f_eq_0
        end while_f_le_v_and_isGiuga ;
        if isGiuga then begin
            % n is still a candidate, check it is not prime                  %
            if fCount > 1 then begin
                gCount := gCount + 1;
                writeon( i_w := 1, s_w := 0, " ", n )
            end if_fCount_gt_1
        end if_isGiuga
    end while_gCount_lt_4
end.
Output:

Same as Algol 68

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

Delphi

Works with: Delphi version 6.0


function IsGiugaNumber(N: integer): boolean;
var IA: TIntegerDynArray;
var I,V: integer;
begin
Result:=False;
if IsPrime(N) then exit;
GetPrimeFactors(N,IA);
for I:=0 to High(IA) do
	begin
	V:=N div IA[I] - 1;
	if (V mod IA[I])<>0 then exit;
	end;
Result:=True;
end;

procedure ShowGiugaNumbers(Memo: TMemo);
var I,Cnt: integer;
begin
Cnt:=0;
for I:=4 to High(integer) do
if IsGiugaNumber(I) then
	begin
	Inc(Cnt);
	Memo.Lines.Add(IntToStr(I));
	if Cnt>=4 then break;
	end;
end;
Output:
30
858
1722
66198

Elapsed Time: 4.991 Sec.


EasyLang

Translation of: Python
func giuga m .
   n = m
   for f = 2 to sqrt n
      while n mod f = 0
         if (m div f - 1) mod f <> 0
            return 0
         .
         n = n div f
         if f > n
            return 1
         .
      .
   .
.
n = 3
while cnt < 4
   if giuga n = 1
      cnt += 1
      print n
   .
   n += 1
.

Euler

Uses the C style for loop procedure from the Sieve of Eratosthenes task

begin
   new for; new n; new gCount;
   for   <- ` formal init; formal test; formal incr; formal body;
              begin
                label again;
                init;
again:          if test then begin body; incr; goto again end else 0
              end
            '
          ;
   gCount <- 0;
   for( ` n <- 2 ', ` gCount < 4 ', ` n <- n + 4 '
      , ` begin
             new v; new f; new isGiuga; new fCount;
             v       <- n % 2;
             isGiuga <- true;
             fCount  <- 1;
             for( ` f <- 3 ', ` f <= v and isGiuga ', ` f <- f + 2 '
                , ` if v mod f = 0 then begin
                       fCount  <- fCount + 1;
                       isGiuga <- [ [ n % f ] - 1 ] mod f = 0;
                       v <- v % f
                    end else 0
                  '
                );
             if isGiuga then begin
                if fCount > 1 then begin
                   gCount <- gCount + 1;
                   out n
                end else 0
             end else 0
          end
        '
      )
end $

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

Java

Algorithm uses the mathematical facts that a Giuga number must be square-free and cannot be a semi-prime.

It does not assume that a Giuga number is even, and exhaustively tests all possible candidates up to approximately 147,000 in around 20 milliseconds.

import java.util.ArrayList;
import java.util.Collections;
import java.util.List;

public final class GiugaNumbers {
	
	public static void main(String[] aArgs) {
        primes = List.of( 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59 );
        
        List<Integer> primeCounts = List.of( 3, 4, 5 );
        for ( int primeCount : primeCounts ) {        
        	primeFactors = new ArrayList<Integer>(Collections.nCopies(primeCount, 0));
        	combinations(primeCount, 0, 0);
        }
        
        Collections.sort(results);
        System.out.println("Found Giuga numbers: " + results);
    }
	
	private static void checkIfGiugaNumber(List<Integer> aPrimeFactors) {
		final int product = aPrimeFactors.stream().reduce(1, Math::multiplyExact);
		
		for ( int factor : aPrimeFactors ) {
			final int divisor = factor * factor;
			if ( ( product - factor ) % divisor != 0 ) {
				return;
			}
		}
		
		results.add(product);		
	}

    private static void combinations(int aPrimeCount, int aIndex, int aLevel) {
        if ( aLevel == aPrimeCount ) {
        	checkIfGiugaNumber(primeFactors);
            return;
        }
        
        for ( int i = aIndex; i < primes.size(); i++ ) {
            primeFactors.set(aLevel, primes.get(i));
            combinations(aPrimeCount, i + 1, aLevel + 1);
        }
    }
 
    private static List<Integer> primes;
    private static List<Integer> primeFactors;
    private static List<Integer> results = new ArrayList<Integer>();

}
Output:
Found Giuga numbers: [30, 858, 1722, 66198]

jq

Works with jq, the C implementation of jq

Works with gojq, the Go implementation of jq

Works with jaq, the Rust implementation of jq

The algorithm used here is naive but suffices to find the first four Giuga numbers in a reasonable amount of time whether using the C, Go, or Rust implementations. On my machine, gojq is fastest at about 12s.

# `div/1` is defined firstly to take advantage of gojq's infinite-precision
# integer arithmetic, and secondly to ensure jaq returns an integer.
def div($j):
  (. - (. % $j)) / $j | round;  # round is for the sake of jaq

# For convenience
def div($i; $j): $i|div($j);

def is_giuga:
  . as $m
  | sqrt as $limit
  | {n: $m, f: 2}
  | until(.ans;
      if (.n % .f) == 0
      then if ((div($m; .f) - 1) % .f) != 0 then .ans = 0
           else .n = div(.n; .f)
           | if .f > .n then .ans = 1 end
           end
      else .f += 1
      | if .f > $limit then .ans = 0 end
      end)
  | .ans == 1 ;

limit(4; range(4; infinite) | select(is_giuga))
Output:
30
858
1722
66198

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)

Lua

Translation of: ALGOL 68
do --[[ find the first 4 Giuga numbers, composites n such that all their
        distinct prime factors f exactly divide ( n / f ) - 1

        each prime factor must appear only once, e.g.: for 2:
        [ ( n / 2 ) - 1 ] mod 2 = 0  => n / 2 is odd => n isn't divisible by 4
        similarly for other primes
     ]]

    local gCount, n = 0, 2
    while gCount < 4 do
        n = n + 4                           -- assume the numbers are all even
        local fCount, f, isGiuga, v = 1, 1, true, math.floor( n / 2 )
        while f <= ( v - 2 ) and isGiuga do
            f = f + 2
            if v % f == 0 then                          -- have a prime factor
                fCount  = fCount + 1
                isGiuga = ( math.floor( n / f ) - 1 ) % f == 0
                v       = math.floor( v / f )
            end
        end
        if isGiuga then       -- n is still a candidate, check it is not prime
            if fCount > 1 then
                gCount = gCount + 1
                io.write( " ", n )
            end
        end
    end
end
Output:
 30 858 1722 66198

Mathematica /Wolfram Language

Translation of: Julia
IsGiuga[n_Integer] := Module[{factors},
  factors = FactorInteger[n];
  AllTrue[factors, Function[{f},
    Mod[Quotient[n, f[[1]]] - 1, f[[1]]] == 0 && f[[1]] != n]]
]

GetGiuga[N_Integer] := Module[{giugaNumbers = {}, i = 4},
  While[Length[giugaNumbers] < N,
    If[IsGiuga[i], AppendTo[giugaNumbers, i]];
    i++
  ];
  giugaNumbers
]

Print[GetGiuga[4]]
Output:
{30, 858, 1722, 66198}


Maxima

/* Predicate function that checks wether an integer is a Giuga number or not */
giugap(n):=if not primep(n) then block(ifactors(n),map(lambda([x],mod((n/x)-1,x)=0),map(first,%%)),
    if length(unique(%%))=1 and apply(lhs,unique(%%))=0 then true)$

/* Function that returns a list of the first len Giuga integers */
giuga_count(len):=block(
    [i:1,count:0,result:[]],
    while count<len do (if giugap(i) then (result:endcons(i,result),count:count+1),i:i+1),
    result)$

/* Test case */
giuga_count(4);
Output:
[30,858,1722,66198]

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 

PARI/GP

Translation of: Julia
is_giuga(n) = {
  my(factors = factor(n));
  for (i = 1, #factors[,1],
    if (factors[i,1] == n, return(0));
    if (Mod(n \ factors[i,1] - 1, factors[i,1]) != 0, return(0));
  );
  return(1);
}

get_giuga(N) = {
  my(giugaNumbers = [], i = 4);
  while(#giugaNumbers < N,
    if (is_giuga(i), giugaNumbers = concat(giugaNumbers, i));
    i++;
  );
  giugaNumbers
}

print(get_giuga(4))
Output:
[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

PascalABC.NET

uses School;

begin
  var n := 1;
  var cnt := 0;
  while cnt < 4 do
  begin
    if n.Factorize.All(f -> (f <> n) and (n div f - 1).Divs(f)) then
    begin
      cnt += 1;
      Print(n);
    end;
    n += 1;
  end;
end.
Output:
30 858 1722 66198 

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.

PL/M

Works with: 8080 PL/M Compiler

... under CP/M (or an emulator)

Only finds 3 Giuga numbers as the fourth is larger than 65535, the largest integer supported by the 8080 PL/M compiler.

100H: /* FIND SOME GIUGA NUMBERS, COMPOSITES N SUCH THAT ALL THEIR DISTINCT */
      /*                       PRIME FACTORS F EXACTLY DIVIDE ( N / F ) - 1 */

   /* CP/M BDOS SYSTEM CALL AND I/O ROUTINES                                */
   BDOS: PROCEDURE( FN, ARG ); DECLARE FN BYTE, ARG ADDRESS; GOTO 5; END;
   PR$CHAR:   PROCEDURE( C ); DECLARE C BYTE;    CALL BDOS( 2, C );  END;
   PR$STRING: PROCEDURE( S ); DECLARE S ADDRESS; CALL BDOS( 9, S );  END;
   PR$NL:     PROCEDURE;   CALL PR$CHAR( 0DH ); CALL PR$CHAR( 0AH ); END;
   PR$NUMBER: PROCEDURE( N ); /* PRINTS A NUMBER IN THE MINIMUN FIELD WIDTH  */
      DECLARE N ADDRESS;
      DECLARE V ADDRESS, N$STR ( 6 )BYTE, W BYTE;
      V = N;
      W = LAST( N$STR );
      N$STR( W ) = '$';
      N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
      DO WHILE( ( V := V / 10 ) > 0 );
         N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
      END;
      CALL PR$STRING( .N$STR( W ) );
   END PR$NUMBER;

   /* TASK                                                                  */

   /* FIND THE FIRST THREE GIUGA NUMBERS (THE FOURTH IS > 65535)            */
   /* EACH PRIME FACTOR CAN ONLY APPEAR ONCE, E.G.,: FOR 2:                 */
   /* (( N / 2 ) - 1) MOD 2 = 0  => N / 2 IS ODD => N NOT DIVISIBLE BY 4    */
   /* SIMILARLY FOR OTHER PRIMES                                            */
   DECLARE ( N, V, FCOUNT, F ) ADDRESS;
   DECLARE IS$GIUGA            BYTE;
   N     = 2;
   DO WHILE N < 65000;                   /* ASSUME THE NUMEBRS ARE ALL EVEN */
      V        = N / 2;
      IS$GIUGA = 1;
      FCOUNT   = 1;
      F        = 1;
      DO WHILE ( F := F + 2 ) <= V AND IS$GIUGA;
         IF V MOD F = 0 THEN DO;
            /* HAVE A PRIME FACTOR                                          */
            FCOUNT   = FCOUNT + 1;
            IS$GIUGA = ( ( N / F ) - 1 ) MOD F = 0;
            V = V / F;
         END;
      END;
      IF IS$GIUGA THEN DO;
         IF FCOUNT > 1 THEN DO;
            /* N IS NOT PRIME, SO IS GIUGA                                  */
            CALL PR$CHAR( ' ' );CALL PR$NUMBER( N );
         END;
      END;
      N = N + 4;
   END;

EOF
Output:
 30 858 1722

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 1722 66198
done...

RPL

Works with: HP version 49
≪ DUP → n
  ≪ FACTORS 1 SF
     1 OVER SIZE FOR j
        DUP j GET
        n OVER / 1 -
        IF SWAP MOD THEN 1 CF END
     2 STEP DROP 
     1 FS?
≫ ≫ 'GIUGA?' STO

≪ { } 4
   WHILE OVER SIZE 4 < REPEAT
      IF DUP GIUGA? THEN LOOK SWAP OVER + SWAP END
      2 +
   END DROP
≫ 'TASK' STO
Output:
1: {30 858 1722 66198}

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]

Scala

Translation of: Java
object GiugaNumbers {

  private var results: List[Int] = List()
  def main(args: Array[String]): Unit = {
    val primes: List[Int] = List(2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59)
    
    val primeCounts: List[Int] = List(3, 4, 5)
    for (primeCount <- primeCounts) {
      var primeFactors: List[Int] = List.fill(primeCount)(0)
      combinations(primes, primeCount, 0, 0, primeFactors)
    }

    val sortedResults = results.sorted
    println(s"Found Giuga numbers: $sortedResults")
  }

  private def checkIfGiugaNumber(primeFactors: List[Int]): Unit = {
    val product: Int = primeFactors.reduce(Math.multiplyExact)
    for (factor <- primeFactors) {
      val divisor: Int = factor * factor
      if ((product - factor) % divisor != 0) {
        return
      }
    }
    results :+= product
  }

  private def combinations(primes: List[Int], primeCount: Int, index: Int, level: Int, primeFactors: List[Int]): Unit = {
    if (level == primeCount) {
      checkIfGiugaNumber(primeFactors)
      return
    }

    for (i <- index until primes.length) {
      val updatedPrimeFactors = primeFactors.updated(level, primes(i))
      combinations(primes, primeCount, i + 1, level + 1, updatedPrimeFactors)
    }
  }  
}
Output:
Found Giuga numbers: List(30, 858, 1722, 66198)

Sidef

4..Inf `by` 2 -> lazy.grep {|n|
    n.factor.all {|p| p `divides` (n/p - 1) }
}.first(4).say
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