Giuga numbers
- 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
C++
Brute force
Based on the Go solution. Takes 26 minutes on my system (Intel Core i5 3.2GHz). <lang cpp>#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; } }
}</lang>
- Output:
First 5 Giuga numbers: 30 858 1722 66198 2214408306
Faster version
<lang cpp>#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);
}</lang>
- Output:
n = 3: 30 n = 4: 1722 858 n = 5: 66198 n = 6: 24423128562 2214408306
FreeBASIC
<lang freebasic>Function isGiuga(m As Uinteger) As Boolean
Dim As Uinteger n = m Dim As Uinteger f = 2, l = Sqr(n) Do If n Mod f = 0 Then If ((m / f) - 1) Mod f <> 0 Then Return False n /= f If f > n Then Return True Else f += 1 If f > l Then Return False End If 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</lang>
- Output:
The first 4 Giuga numbers are: 30 858 1722 66198
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. <lang go>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)
}</lang>
- Output:
The first 5 Giuga numbers are: [30 858 1722 66198 2214408306]
J
We can brute force this task building a test for giuga numbers and checking the first hundred thousand integers (which takes a small fraction of a second):
<lang J>giguaP=: {{ (1<y)*(-.1 p:y)**/(=<.) y ((_1+%)%]) q: y }}"0</lang>
<lang J> 1+I.giguaP 1+i.1e5 30 858 1722 66198</lang>
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:
<lang J>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
}}</lang>
Example use:<lang J> giuga 1
giuga 2
giuga 3
30
giuga 4
1722 858
giuga 5
66198
giuga 6
24423128562 2214408306</lang>
Julia
<lang ruby>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)
</lang>
- Output:
30 858 1722 66198
Ad hoc faster version
<lang ruby>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)
</lang>
- Output:
30 (elapsed: 0.0) 858 (elapsed: 0.0) 1722 (elapsed: 0.0) 66198 (elapsed: 0.0009999275207519531) 2214408306 (elapsed: 18.97099995613098) 24423128562 (elapsed: 432.06500005722046) 432.066249 seconds (235 allocations: 12.523 KiB)
Pascal
Free Pascal
OK.Cheating to find square free numbers like julia in distance 6
That means always factors 2,3 and minimum one of 5,7,11.
<lang pascal>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.</lang>
- @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.
<lang pascal>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. </lang>
- @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
<lang 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;</lang>
- 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.
Python
<lang 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</lang>
Raku
<lang perl6>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) %% $_ }; }
}</lang>
- Output:
First 4 Giuga numbers: 30 858 1722 66198
Wren
Version 1 (Brute force)
Simple brute force but assumes all Giuga numbers will be even, must be square-free and can't be semi-prime.
Takes only about 0.05 seconds to find the first four Giuga numbers but finding the fifth would take many hours using this approach, so I haven't bothered. <lang ecmascript>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)</lang>
- 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. <lang ecmascript>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()
}</lang>
- Output:
n = 3: 30 n = 4: 1722 858 n = 5: 66198 n = 6: 24423128562 2214408306
XPL0
<lang 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; ];
]</lang>
- Output:
30 858 1722 66198