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[edit]
F isGiuga(m)
V n = m
V f = 2
V l = sqrt(n)
L
I n % f == 0
I ((m I/ f) - 1) % f != 0
R 0B
n I/= f
I f > n
R 1B
E
f++
I f > l
R 0B
V n = 3
V c = 0
print(‘The first 4 Giuga numbers are:’)
L c < 4
I isGiuga(n)
c++
print(n)
n++- Output:
The first 4 Giuga numbers are: 30 858 1722 66198
ALGOL 68[edit]
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[edit]
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[edit]
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[edit]
# 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[edit]
BASIC256[edit]
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[edit]
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[edit]
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[edit]
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++[edit]
Brute force[edit]
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[edit]
#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
Euler[edit]
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[edit]
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[edit]
--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[edit]
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
}}
giuga 1
giuga 2
giuga 3
30
giuga 4
1722 858
giuga 5
66198
giuga 6
24423128562 2214408306
Java[edit]
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]
Julia[edit]
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[edit]
using Primes
function getgiugas(numberwanted, verbose = true)
n, found, nfound = 6, Int[], 0
starttime = time()
while nfound < numberwanted
if n % 5 == 0 || n % 7 == 0 || n % 11 == 0
for (p, e) in eachfactor(n)
(e != 1 || rem(n ÷ p - 1, p) != 0) && @goto nextnumber
end
verbose && println(n, " (elapsed: ", time() - starttime, ")")
push!(found, n)
nfound += 1
end
@label nextnumber
n += 6 # all mult of 6
end
return found
end
@time getgiugas(2, false)
@time getgiugas(6)
- Output:
30 (elapsed: 0.0) 858 (elapsed: 0.0) 1722 (elapsed: 0.0) 66198 (elapsed: 0.0009999275207519531) 2214408306 (elapsed: 18.97099995613098) 24423128562 (elapsed: 432.06500005722046) 432.066249 seconds (235 allocations: 12.523 KiB)
Nim[edit]
import std/math
func isGiuga(m: Natural): bool =
var n = m
var f = 2
var l = int(sqrt(n.toFloat))
while true:
if n mod f == 0:
if (m div f - 1) mod f != 0:
return false
n = n div f
if f > n:
return true
else:
inc f
if f > l:
return false
var n = 3
var c = 0
const Limit = 4
stdout.write "The first ", Limit, " Giuga numbers are: "
while true:
if n.isGiuga:
inc c
stdout.write n, " "
if c == Limit: break
inc n
echo()
- Output:
The first 4 Giuga numbers are: 30 858 1722 66198
Pascal[edit]
Free Pascal[edit]
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[edit]
Generating recursive squarefree numbers of ascending primes and check those numbers.
2*3 are set fixed. 2*3*5 followed 2*3*7 than 2*3*11. Therefor the results are unsorted.
program Giuga;
{
30 = 2 * 3 * 5.
858 = 2 * 3 * 11 * 13.
1722 = 2 * 3 * 7 * 41.
66198 = 2 * 3 * 11 * 17 * 59.
2214408306 = 2 * 3 * 11 * 23 * 31 * 47057.
24423128562 = 2 * 3 * 7 * 43 * 3041 * 4447.
432749205173838 = 2 * 3 * 7 * 59 * 163 * 1381 * 775807.
14737133470010574 = 2 * 3 * 7 * 71 * 103 * 67213 * 713863.
550843391309130318 = 2 * 3 * 7 * 71 * 103 * 61559 * 29133437.
244197000982499715087866346 = 2 * 3 * 11 * 23 * 31 * 47137 * 28282147 * 3892535183.
554079914617070801288578559178 = 2 * 3 * 11 * 23 * 31 * 47059 * 2259696349 * 110725121051.
1910667181420507984555759916338506 = 2 * 3 * 7 * 43 * 1831 * 138683 * 2861051 * 1456230512169437. }
{$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
{$ELSE}
{$APPTYPE CONSOLE}
{$ENDIF}
uses
sysutils
{$IFDEF WINDOWS},Windows{$ENDIF}
;
const
LMT =14737133470010574;// 432749205173838;//24423128562;//2214408306;
type
tFac = packed record
fMul :Uint64;
fPrime,
fPrimIdx,
fprimMaxIdx,dummy :Uint32;
dummy2: Uint64;
end;
tFacs = array[0..15] of tFac;
tPrimes = array[0..62157] of Uint32;//775807 see factor of 432749205173838
// tPrimes = array[0..4875{14379}] of Uint32;//sqrt 24423128562
// tPrimes = array[0..1807414] of Uint32;//29133437
// tPrimes = array[0..50847533] of Uint32;// 1e9
// tPrimes = array[0..5761454] of Uint32;//1e8
var
SmallPrimes: tPrimes;
T0 : Int64;
cnt:Uint32;
procedure InitSmallPrimes;
//get primes. #0..65535.Sieving only odd numbers
const
//MAXLIMIT = (trunc(sqrt(LMT)+1)-1) shr 1+4;
MAXLIMIT = 775807 DIV 2+1;//(trunc(sqrt(LMT)+1)-1) shr 1+4;
var
pr : array of byte;
pPr :pByte;
p,j,d,flipflop :NativeUInt;
Begin
SmallPrimes[0] := 2;
setlength(pr,MAXLIMIT);
pPr := @pr[0];
p := 0;
repeat
repeat
p +=1
until pPr[p]= 0;
j := (p+1)*p*2;
if j>MAXLIMIT then
BREAK;
d := 2*p+1;
repeat
pPr[j] := 1;
j += d;
until j>MAXLIMIT;
until false;
SmallPrimes[1] := 3;
SmallPrimes[2] := 5;
j := 3;
d := 7;
flipflop := (2+1)-1;//7+2*2,11+2*1,13,17,19,23
p := 3;
repeat
if pPr[p] = 0 then
begin
SmallPrimes[j] := d;
inc(j);
end;
d += 2*flipflop;
p+=flipflop;
flipflop := 3-flipflop;
until (p > MAXLIMIT) OR (j>High(SmallPrimes));
setlength(pr,0);
end;
procedure OutFac(var F:tFacs;maxIdx:Uint32);
var
i : integer;
begin
write(cnt:3,' ');
For i := 0 to maxIdx do
write(F[i].fPrime,'*');
write(#8,' = ',F[maxIdx].fMul);
writeln(' ',(GetTickCount64-T0)/1000:10:3,' s');
//readln;
end;
function ChkGiuga(var F:tFacs;MaxIdx:Uint32):boolean;inline;
var
n : Uint64;
idx: NativeInt;
p : Uint32;
begin
n := F[MaxIdx].fMul;
idx := MaxIdx;
repeat
p := F[idx].fPrime;
result := (((n DIV p)-1)MOD p) = 0;
if not(result) then
EXIT;
dec(idx);
until idx<0;
inc(cnt);
end;
procedure InsertNextPrimeFac(var F:tFacs;idx:Uint32);
var
Mul : Uint64;
i,p : uint32;
begin
with F[idx-1] do
begin
Mul := fMul;
i := fPrimIdx;
end;
while i<High(SmallPrimes) do
begin
inc(i);
with F[idx] do
begin
if i >fprimMaxIdx then
BREAK;
p := SmallPrimes[i];
if p*Mul>LMT then
BREAK;
fMul := p*Mul;
fPrime := p;
fPrimIdx := i;
IF (Mul-1) MOD p = 0 then
IF ChkGiuga(F,idx) then
OutFac(F,idx);
end;
// max 6 factors //for lmt 24e9 need 7 factors
if idx <5 then
InsertNextPrimeFac(F,idx+1);
end;
end;
var
{$ALIGN 32}
Facs : tFacs;
i : integer;
Begin
InitSmallPrimes;
T0 := GetTickCount64;
with Facs[0] do
begin
fMul := 2;fPrime := 2;fPrimIdx:= 0;
end;
with Facs[1] do
begin
fMul := 2*3;fPrime := 3;fPrimIdx:= 1;
end;
i := 2;
//search the indices of mx found factor
while SmallPrimes[i] < 11 do inc(i); Facs[2].fprimMaxIdx := i;
while SmallPrimes[i] < 71 do inc(i); Facs[3].fprimMaxIdx := i;
while SmallPrimes[i] < 3041 do inc(i); Facs[4].fprimMaxIdx := i;
while SmallPrimes[i] < 67213 do inc(i); Facs[5].fprimMaxIdx := i;
while SmallPrimes[i] < 775807 do inc(i); Facs[6].fprimMaxIdx := i;
{
writeln('Found ',cnt,' in ',(GetTickCount64-T0)/1000:10:3,' s');
//start with 2*3*7
with Facs[2] do
begin
fMul := 2*3*7;fPrime := 7;fPrimIdx:= 3;
end;
InsertNextPrimeFac(Facs,3);
//start with 2*3*11
writeln('Found ',cnt,' in ',(GetTickCount64-T0)/1000:10:3,' s');
with Facs[2] do
begin
fMul := 2*3*11;fPrime := 11;fPrimIdx:= 4;
end;
InsertNextPrimeFac(Facs,3);
}
InsertNextPrimeFac(Facs,2);
writeln('Found ',cnt,' in ',(GetTickCount64-T0)/1000:10:3,' s');
writeln;
end.
- @TIO.RUN:
1 2*3*5 = 30 0.000 s 2 2*3*7*41 = 1722 0.810 s 3 2*3*7*43*3041*4447 = 24423128562 0.871 s 4 2*3*11*13 = 858 1.057 s 5 2*3*11*17*59 = 66198 1.089 s 6 2*3*11*23*31*47057 = 2214408306 1.152 s Found 6 in 1.526 s Real time: 1.682 s CPU share: 99.42 % @home: Limit:432749205173838 start 2*3*7 Found 0 in 0.000 s 1 2*3*7*41 = 1722 147.206 s 2 2*3*7*43*3041*4447 = 24423128562 163.765 s 3 2*3*7*59*163*1381*775807 = 432749205173838 179.124 s Found 3 in 197.002 s start 2*3*11 4 2*3*11*13 = 858 197.002 s 5 2*3*11*17*59 = 66198 219.166 s 6 2*3*11*23*31*47057 = 2214408306 244.468 s real 5m10,271s Limit :14737133470010574 start 2*3*7 Found 0 in 0.000 s 1 2*3*7*41 = 1722 1330.819 s 2 2*3*7*43*3041*4447 = 24423128562 1567.028 s 3 2*3*7*59*163*1381*775807 = 432749205173838 1788.203 s 4 2*3*7*71*103*67213*713863 = 14737133470010574 2051.552 s Found 4 in 2129.801 s start 2*3*11 5 2*3*11*13 = 858 2129.801 s 6 2*3*11*17*59 = 66198 2305.752 s 7 2*3*11*23*31*47057 = 2214408306 2591.984 s Found 7 in 3654.610 s real 60m54,612s
Perl[edit]
#!/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[edit]
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[edit]
-- -- 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[edit]
... 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[edit]
#!/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[edit]
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[edit]
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[edit]
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...
Ruby[edit]
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[edit]
use prime_tools ;
fn prime_decomposition( mut number : u32) -> Vec<u32> {
let mut divisors : Vec<u32> = Vec::new( ) ;
let mut divisor : u32 = 2 ;
while number != 1 {
if number % divisor == 0 {
divisors.push( divisor ) ;
number /= divisor ;
}
else {
divisor += 1 ;
}
}
divisors
}
fn is_giuga( num : u32 ) -> bool {
let prime_factors : Vec<u32> = prime_decomposition( num ) ;
! prime_tools::is_u32_prime( num ) &&
prime_factors.into_iter( ).all( |n : u32| (num/n -1) % n == 0 )
}
fn main() {
let mut giuga_numbers : Vec<u32> = Vec::new( ) ;
let mut num : u32 = 2 ;
while giuga_numbers.len( ) != 4 {
if is_giuga( num ) {
giuga_numbers.push( num ) ;
}
num += 1 ;
}
println!("{:?}" , giuga_numbers ) ;
}
- Output:
[30, 858, 1722, 66198]
Wren[edit]
Version 1 (Brute force)[edit]
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)[edit]
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[edit]
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