Giuga numbers: Difference between revisions
m (typo) |
m (→{{header|Free Pascal}}: cheating like julia version) |
||
Line 205: | Line 205: | ||
=={{header|Pascal}}== |
=={{header|Pascal}}== |
||
==={{header|Free Pascal}}=== |
==={{header|Free Pascal}}=== |
||
OK.Cheating to find square free numbers like julia in distance 6<BR> |
|||
changing main routine at the end of [[Factors_of_an_integer#using_Prime_decomposition]] <BR> |
|||
That means always factors 2,3 and minimum one of 5,7,11.<br> |
|||
Mostly, about 70%, the highest primefactor remains in pfRemain.<br> |
|||
<br> |
|||
<lang pascal> |
|||
<lang pascal>program Giuga; |
|||
function ChkGuiga(n:Uint64;pPrimeDecomp :tpPrimeFac):boolean;inline; |
|||
{$IFDEF FPC} |
|||
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON} |
|||
{$ELSE} |
|||
{$APPTYPE CONSOLE} |
|||
{$ENDIF} |
|||
uses |
|||
sysutils |
|||
{$IFDEF WINDOWS},Windows{$ENDIF} |
|||
; |
|||
//###################################################################### |
|||
//prime decomposition |
|||
type |
|||
tprimeFac = packed record |
|||
pfRemain : Uint64; |
|||
pfMaxIdx : Uint32; |
|||
pfpotPrimIdx : array[0..9] of 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; |
|||
p := pfRemain; |
|||
If p >1 then |
|||
Begin |
|||
str(p,s); |
|||
chk *= p; |
|||
result += '*'+s; |
|||
end; |
|||
str(chk,s); |
|||
result += '_chk_'+s+'<'; |
|||
end; |
|||
end; |
|||
function IsSquarefreeDecomp(var res:tPrimeFac;n:Uint64):boolean;inline; |
|||
//factorize only not prime/semiprime and squarefree |
|||
var |
|||
pr,i,q :NativeUInt; |
|||
Begin |
|||
with res do |
|||
Begin |
|||
pfMaxIdx := 0; |
|||
q := n DIV 2; |
|||
if n = 2*q then |
|||
Begin |
|||
pfpotPrimIdx[pfMaxIdx] := 2; |
|||
n := q; |
|||
if (q AND 1) = 0 then |
|||
EXIT(false); |
|||
inc(pfMaxIdx); |
|||
end; |
|||
q := n DIV 3; |
|||
if n = 3*q then |
|||
Begin |
|||
pfpotPrimIdx[pfMaxIdx] := 3; |
|||
n := q; |
|||
q := q div 3; |
|||
if n = 3*q then |
|||
EXIT(false); |
|||
inc(pfMaxIdx); |
|||
end; |
|||
q := n DIV 5; |
|||
if n = 5*q then |
|||
Begin |
|||
pfpotPrimIdx[pfMaxIdx] := 5; |
|||
n := q; |
|||
q := q div 5; |
|||
if q*5=n then |
|||
EXIT(false); |
|||
inc(pfMaxIdx); |
|||
end; |
|||
q := n DIV 7; |
|||
if n = 7*q then |
|||
Begin |
|||
pfpotPrimIdx[pfMaxIdx] := 7; |
|||
n := q; |
|||
q := q div 7; |
|||
if q*7=n then |
|||
EXIT(false); |
|||
inc(pfMaxIdx); |
|||
end; |
|||
q := n DIV 11; |
|||
if n = 11*q then |
|||
Begin |
|||
pfpotPrimIdx[pfMaxIdx] := 11; |
|||
n := q; |
|||
q := q div 11; |
|||
if q*11=n then |
|||
EXIT(false); |
|||
inc(pfMaxIdx); |
|||
end; |
|||
if pfMaxIdx < 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[pfMaxIdx] := pr; |
|||
n := q; |
|||
q := n div pr; |
|||
if pr*q = n then |
|||
EXIT(false); |
|||
inc(pfMaxIdx); |
|||
end; |
|||
inc(i); |
|||
end; |
|||
pfRemain := n; |
|||
end; |
|||
exit(true); |
|||
end; |
|||
function ChkGiuga(n:Uint64;pPrimeDecomp :tpPrimeFac):boolean;inline; |
|||
var |
var |
||
pFr : Uint64; |
pFr : Uint64; |
||
Line 217: | Line 410: | ||
Begin |
Begin |
||
pfR := pfRemain; |
pfR := pfRemain; |
||
idx := pfMaxIdx; |
idx := pfMaxIdx; |
||
//check squarefree. i primefactors -> count of diviors = 2^i |
|||
p := 1 shl idx; |
|||
if pfR <> 1 then |
if pfR <> 1 then |
||
begin |
begin |
||
if (p+p <> pfdivcnt) then |
|||
EXIT(false); |
|||
if sqr(pfR)>=n then |
if sqr(pfR)>=n then |
||
EXIT(false); |
EXIT(false); |
||
Line 232: | Line 421: | ||
idx -= 1; |
idx -= 1; |
||
repeat |
repeat |
||
p := |
p := pfpotPrimIdx[idx]; |
||
result := (((n DIV p)-1)MOD p) = 0; |
result := (((n DIV p)-1)MOD p) = 0; |
||
if Not(result) then |
if Not(result) then |
||
Line 242: | Line 431: | ||
else |
else |
||
begin |
begin |
||
if (p <> pfdivcnt) then |
|||
EXIT(false); |
|||
result := true; |
result := true; |
||
repeat |
repeat |
||
dec(idx); |
dec(idx); |
||
p := |
p := pfpotPrimIdx[idx]; |
||
result := (((n DIV p)-1)MOD p) = 0; |
result := (((n DIV p)-1)MOD p) = 0; |
||
if not(result) then |
if not(result) then |
||
Line 255: | Line 442: | ||
end; |
end; |
||
end; |
end; |
||
const |
const |
||
LMT = |
LMT = 24423128562;//2214408306;// |
||
var |
var |
||
PrimeDecomp :tPrimeFac; |
|||
pPrimeDecomp :tpPrimeFac; |
|||
T0:Int64; |
T0:Int64; |
||
n : NativeUInt; |
n : NativeUInt; |
||
Line 265: | Line 452: | ||
Begin |
Begin |
||
InitSmallPrimes; |
InitSmallPrimes; |
||
T0 := GetTickCount64; |
T0 := GetTickCount64; |
||
Init_Sieve(1); |
|||
n := |
n := 6; |
||
cnt := 0; |
cnt := 0; |
||
repeat |
repeat |
||
//only |
//only multibles of 6 |
||
inc(n, |
inc(n,6); |
||
GetNextPrimeDecomp; |
|||
if IsSquarefreeDecomp(PrimeDecomp,n)then |
|||
pPrimeDecomp:= GetNextPrimeDecomp; |
|||
if ChkGuiga(n,@PrimeDecomp) then |
|||
//if not prime/semiprime |
|||
if pPrimeDecomp^.pfDivCnt >4 then |
|||
if ChkGuiga(n,pPrimeDecomp) then |
|||
begin |
begin |
||
inc(cnt); |
inc(cnt); |
||
writeln(cnt:3,'..',OutPots( |
writeln(cnt:3,'..',OutPots(@PrimeDecomp,n),' ',(GettickCount64-T0)/1000:6:3,' s'); |
||
end; |
end; |
||
until n >= LMT; |
until n >= LMT; |
||
Line 287: | Line 472: | ||
writeln('Tested til ',n,' runtime ',T0/1000:0:3,' s'); |
writeln('Tested til ',n,' runtime ',T0/1000:0:3,' s'); |
||
writeln; |
writeln; |
||
writeln(OutPots( |
writeln(OutPots(@PrimeDecomp,n)); |
||
end.</lang> |
end.</lang> |
||
{{out|@home AMD 5600G fpc3.2.2 -O4 -Xs |
{{out|@home AMD 5600G fpc3.2.2 -O4 -Xs} |
||
<pre> |
<pre> |
||
1..30 : 8 : 2*3*5_chk_30 |
1..30 : 8 : 2*3*5_chk_30< 0.000 s |
||
2..858 : |
2..858 : 1 : 2*3*11*13_chk_858< 0.000 s |
||
3..1722 : |
3..1722 : 1 : 2*3*7*41_chk_1722< 0.000 s |
||
4..66198 : |
4..66198 : 1 : 2*3*11*17*59_chk_66198< 0.000 s |
||
5..2214408306 : |
5..2214408306 : 1 : 2*3*11*23*31*47057_chk_2214408306< 17.723 s |
||
6..24423128562 : |
6..24423128562 : 1 : 2*3*7*43*3041*4447_chk_24423128562< 457.019 s |
||
Found 6 |
Found 6 |
||
Tested til 24423128562 runtime |
Tested til 24423128562 runtime 457.019 s |
||
</pre> |
|||
=={{header|Perl}}== |
=={{header|Perl}}== |
Revision as of 19:16, 5 July 2022
- 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
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.
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
type
tprimeFac = packed record pfRemain : Uint64; pfMaxIdx : Uint32; pfpotPrimIdx : array[0..9] of 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; p := pfRemain; If p >1 then Begin str(p,s); chk *= p; result += '*'+s; end; str(chk,s); result += '_chk_'+s+'<'; end;
end;
function IsSquarefreeDecomp(var res:tPrimeFac;n:Uint64):boolean;inline; //factorize only not prime/semiprime and squarefree var
pr,i,q :NativeUInt;
Begin
with res do Begin pfMaxIdx := 0;
q := n DIV 2; if n = 2*q then Begin pfpotPrimIdx[pfMaxIdx] := 2; n := q; if (q AND 1) = 0 then EXIT(false); inc(pfMaxIdx); end;
q := n DIV 3; if n = 3*q then Begin pfpotPrimIdx[pfMaxIdx] := 3; n := q; q := q div 3; if n = 3*q then EXIT(false); inc(pfMaxIdx); end;
q := n DIV 5; if n = 5*q then Begin pfpotPrimIdx[pfMaxIdx] := 5; n := q; q := q div 5; if q*5=n then EXIT(false); inc(pfMaxIdx); end;
q := n DIV 7; if n = 7*q then Begin pfpotPrimIdx[pfMaxIdx] := 7; n := q; q := q div 7; if q*7=n then EXIT(false); inc(pfMaxIdx); end;
q := n DIV 11; if n = 11*q then Begin pfpotPrimIdx[pfMaxIdx] := 11; n := q; q := q div 11; if q*11=n then EXIT(false); inc(pfMaxIdx); end; if pfMaxIdx < 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[pfMaxIdx] := pr; n := q; q := n div pr; if pr*q = n then EXIT(false); inc(pfMaxIdx); end; inc(i); end; pfRemain := n; end; exit(true);
end;
function ChkGiuga(n:Uint64;pPrimeDecomp :tpPrimeFac):boolean;inline; var
pFr : Uint64; idx: NativeInt; p : Uint32;
begin
with pPrimeDecomp^ do Begin pfR := pfRemain; idx := pfMaxIdx; if pfR <> 1 then begin if sqr(pfR)>=n then EXIT(false); //squarefree result := (((n DIV pfR)-1)MOD pfR) = 0; if result then begin idx -= 1; repeat p := pfpotPrimIdx[idx]; result := (((n DIV p)-1)MOD p) = 0; if Not(result) then EXIT(result); dec(idx); until idx <0; end; end else begin result := true; repeat dec(idx); p := pfpotPrimIdx[idx]; result := (((n DIV p)-1)MOD p) = 0; if not(result) then EXIT(false); until idx<=0; end; end;
end;
const
LMT = 24423128562;//2214408306;//
var
PrimeDecomp :tPrimeFac; T0:Int64; n : NativeUInt; cnt:Uint32;
Begin
InitSmallPrimes; T0 := GetTickCount64;
n := 6; cnt := 0; repeat //only multibles of 6 inc(n,6); if IsSquarefreeDecomp(PrimeDecomp,n)then if ChkGuiga(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> {{out|@home AMD 5600G fpc3.2.2 -O4 -Xs}
1..30 : 8 : 2*3*5_chk_30< 0.000 s 2..858 : 1 : 2*3*11*13_chk_858< 0.000 s 3..1722 : 1 : 2*3*7*41_chk_1722< 0.000 s 4..66198 : 1 : 2*3*11*17*59_chk_66198< 0.000 s 5..2214408306 : 1 : 2*3*11*23*31*47057_chk_2214408306< 17.723 s 6..24423128562 : 1 : 2*3*7*43*3041*4447_chk_24423128562< 457.019 s Found 6 Tested til 24423128562 runtime 457.019 s
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}
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
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]
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