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
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
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
#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
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
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
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
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
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
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
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
-- -- 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
... 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
#!/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
≪ 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
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)
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
Categories:
- Programming Tasks
- Solutions by Programming Task
- 11l
- ABC
- Ada
- ALGOL 68
- ALGOL W
- Arturo
- AWK
- BASIC
- BASIC256
- FreeBASIC
- Gambas
- PureBasic
- C++
- Boost
- Delphi
- SysUtils,StdCtrls
- EasyLang
- Euler
- Go
- Haskell
- J
- Java
- Jq
- Julia
- Lua
- Mathematica
- Wolfram Language
- Maxima
- Nim
- PARI/GP
- Pascal
- Free Pascal
- PascalABC.NET
- Perl
- Phix
- PL/M
- Python
- Quackery
- Raku
- Ring
- RPL
- Ruby
- Rust
- Scala
- Sidef
- Wren
- Wren-math
- Wren-rat
- XPL0