Giuga numbers: Difference between revisions

From Rosetta Code
Content added Content deleted
m (→‎{{header|Free Pascal}}: forgotten to insert the lang <lang pascal> append used time)
m (→‎{{header|Free Pascal}}: n is now alsways even.remove this check.Now 6.th guiga number found < 9 min.)
Line 136: Line 136:
=={{header|Pascal}}==
=={{header|Pascal}}==
==={{header|Free Pascal}}===
==={{header|Free Pascal}}===
changing main routine at the end of [[Factors_of_an_integer#using_Prime_decomposition]]
changing main routine at the end of [[Factors_of_an_integer#using_Prime_decomposition]] <BR>
Mostly, about 70%, the highest primefactor remains in pfRemain.<br>
<lang pascal>
<lang pascal>
function ChkGuiga(n:Uint64;pPrimeDecomp :tpPrimeFac):boolean;inline;
var
pFr : Uint64;
idx: NativeInt;
p : Uint32;
begin
with pPrimeDecomp^ do
Begin
pfR := pfRemain;
idx := pfMaxIdx;//always>0 for pfDivcnt>4
//check squarefree. i primefactors -> count of diviors = 2^i
p := 1 shl idx;
if pfR <> 1 then
begin
if (p+p <> pfdivcnt) then
EXIT(false);
if sqr(pfR)>=n then
EXIT(false);
//squarefree
result := (((n DIV pfR)-1)MOD pfR) = 0;
if result then
begin
idx -= 1;
repeat
p := SmallPrimes[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
if (p <> pfdivcnt) then
EXIT(false);
result := true;
repeat
dec(idx);
p := SmallPrimes[pfpotPrimIdx[idx]];
result := (((n DIV p)-1)MOD p) = 0;
if not(result) then
EXIT(false);
until idx<=0;
end;
end;
end;

const
const
LMT = 24423128562;//70*1000;//
LMT = 24423128562;//24423128562;//2214408306;//
var
var
pPrimeDecomp :tpPrimeFac;
pPrimeDecomp :tpPrimeFac;
T0:Int64;
T0:Int64;
n : NativeUInt;
n : NativeUInt;
idx,p,cnt : Int32;
cnt:Uint32;
Begin
chk: boolean;
BEGIN
InitSmallPrimes;
InitSmallPrimes;


T0 := GetTickCount64;
T0 := GetTickCount64;
n := 1;
Init_Sieve(1);
n := 0;
Init_Sieve(n);
pPrimeDecomp:= GetNextPrimeDecomp;
cnt := 0;
cnt := 0;
repeat
repeat
//only even numbers
inc(n);
inc(n,2);
GetNextPrimeDecomp;
pPrimeDecomp:= GetNextPrimeDecomp;
pPrimeDecomp:= GetNextPrimeDecomp;
//if not prime/semiprime

with pPrimeDecomp^ do
if pPrimeDecomp^.pfDivCnt >4 then
if ChkGuiga(n,pPrimeDecomp) then
begin
//if prime/semi-prime
if pfDivCnt <= 4 then continue;
//not even
if pfpotPrimIdx[0] <> 0 then continue;
idx := pfMaxIdx;//always>0 for pfDivcnt>2

if pfRemain = 1 then
begin
//not squarefree or biggest prime factor to big
if ((1 shl idx) <> pfdivcnt) or (sqr(SmallPrimes[pfpotPrimIdx[idx-1]])>=n) then
continue;
//check for Giuga number
chk := true;
For idx := idx-1 downto 0 do
begin
p := SmallPrimes[pfpotPrimIdx[idx]];
chk := (((n DIV p)-1)MOD p) = 0;
if not(chk) then
BREAK;
end;
end
else
begin
if (1 shl (idx+1) <> pfdivcnt) or (sqr(pfRemain)>=n) then
continue;
chk := (((n DIV pfRemain)-1)MOD pfRemain) = 0;
if chk then
For idx := idx-1 downto 0 do
begin
p := SmallPrimes[pfpotPrimIdx[idx]];
chk := (((n DIV p)-1)MOD p) = 0;
if not(chk) then
BREAK;
end;
end;
if chk then
begin
begin
inc(cnt);
inc(cnt);
writeln(cnt:3,'..',OutPots(pPrimeDecomp,n),' ',(GettickCount64-T0)/1000:6:3,' s');
writeln(cnt:3,'..',OutPots(pPrimeDecomp,n),' ',(GettickCount64-T0)/1000:6:3,' s');
end;
end;
end;
until n >= LMT;
until n >= LMT;
T0 := GetTickCount64-T0;
T0 := GetTickCount64-T0;
writeln('runtime ',T0/1000:0:3,' s');
writeln('Found ',cnt);
writeln('Count ',cnt);
writeln('Tested til ',n,' runtime ',T0/1000:0:3,' s');
writeln;
writeln;
writeln(OutPots(pPrimeDecomp,n));
END.</lang>
end.</lang>
{{out|@home AMD 5600G}}
{{out|@home AMD 5600G fpc3.2.2 -O4 -Xs}}
<pre>
<pre>
1..30 : 8 : 2*3*5_chk_30<_SoD_72< 0.001 s
1..30 : 8 : 2*3*5_chk_30<_SoD_72< 0.001 s
2..858 : 16 : 2*3*11*13_chk_858<_SoD_2016< 0.001 s
2..858 : 16 : 2*3*11*13_chk_858<_SoD_2016< 0.001 s
3..1722 : 16 : 2*3*7*41_chk_1722<_SoD_4032< 0.001 s
3..1722 : 16 : 2*3*7*41_chk_1722<_SoD_4032< 0.001 s
4..66198 : 32 : 2*3*11*17*59_chk_66198<_SoD_155520< 0.003 s
4..66198 : 32 : 2*3*11*17*59_chk_66198<_SoD_155520< 0.002 s
5..2214408306 : 64 : 2*3*11*23*31*47057_chk_2214408306<_SoD_5204238336< 52.020 s
5..2214408306 : 64 : 2*3*11*23*31*47057_chk_2214408306<_SoD_5204238336< 41.656 s
6..24423128562 : 64 : 2*3*7*43*3041*4447_chk_24423128562<_SoD_57154166784<
6..24423128562 : 64 : 2*3*7*43*3041*4447_chk_24423128562<_SoD_57154166784< 533.047 s
Found 6
runtime 673.643 s
Tested til 24423128562 runtime 533.047 s</pre>
Count 6
real 11m13,645s user 11m12,962s sys 0m0,247s</pre>


=={{header|Perl}}==
=={{header|Perl}}==

Revision as of 14:12, 4 July 2022

Giuga numbers is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
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

Translation of: Wren
Library: Go-rcu

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 around 2 hours 40 minutes. <lang go>package main

import (

   "fmt"
   "rcu"

)

func main() {

   const limit = 5
   var giuga []int
   for n := 4; len(giuga) < limit; n += 2 {
       factors := rcu.PrimeFactors(n)
       if len(factors) > 2 {
           isSquareFree := true
           for i := 1; i < len(factors); i++ {
               if factors[i] == factors[i-1] {
                   isSquareFree = false
                   break
               }
           }
           if isSquareFree {
               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

Pascal

Free Pascal

changing main routine at the end of Factors_of_an_integer#using_Prime_decomposition
Mostly, about 70%, the highest primefactor remains in pfRemain.
<lang pascal> function ChkGuiga(n:Uint64;pPrimeDecomp :tpPrimeFac):boolean;inline; var

 pFr : Uint64;
 idx: NativeInt;
 p : Uint32;

begin

with pPrimeDecomp^ do
Begin
 pfR := pfRemain;
 idx := pfMaxIdx;//always>0 for pfDivcnt>4
 //check squarefree. i primefactors -> count of diviors = 2^i 
 p := 1 shl idx;
 if pfR <> 1 then
 begin
   if (p+p <> pfdivcnt) then
     EXIT(false);
   if sqr(pfR)>=n then
     EXIT(false);
   //squarefree
   result := (((n DIV pfR)-1)MOD pfR) = 0;
   if result then
   begin
     idx -= 1;
     repeat
       p := SmallPrimes[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
   if (p <> pfdivcnt) then
     EXIT(false);
   result := true;
   repeat
     dec(idx);
     p := SmallPrimes[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;//24423128562;//2214408306;//

var

 pPrimeDecomp :tpPrimeFac;
 T0:Int64;
 n : NativeUInt;
 cnt:Uint32;

Begin

 InitSmallPrimes;
 T0 := GetTickCount64;
 Init_Sieve(1);
 n := 0;
 cnt := 0;
 repeat
   //only even numbers
   inc(n,2);
   GetNextPrimeDecomp;
   pPrimeDecomp:= GetNextPrimeDecomp;
   //if not prime/semiprime
   if pPrimeDecomp^.pfDivCnt >4 then
     if ChkGuiga(n,pPrimeDecomp) then
     begin
       inc(cnt);
       writeln(cnt:3,'..',OutPots(pPrimeDecomp,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(pPrimeDecomp,n));

end.</lang>

@home AMD 5600G fpc3.2.2 -O4 -Xs:
  1..30 :  8 : 2*3*5_chk_30<_SoD_72<   0.001 s
  2..858 : 16 : 2*3*11*13_chk_858<_SoD_2016<   0.001 s
  3..1722 : 16 : 2*3*7*41_chk_1722<_SoD_4032<   0.001 s
  4..66198 : 32 : 2*3*11*17*59_chk_66198<_SoD_155520<   0.002 s
  5..2214408306 : 64 : 2*3*11*23*31*47057_chk_2214408306<_SoD_5204238336<  41.656 s
  6..24423128562 : 64 : 2*3*7*43*3041*4447_chk_24423128562<_SoD_57154166784<  533.047 s
Found 6
Tested til 24423128562 runtime 533.047 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

Translation of: FreeBASIC

<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

Library: Wren-math
Library: Wren-seq

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.1 seconds to find the first four Giuga numbers but finding the fifth would take many hours using this approach, so I haven't bothered. Hopefully, someone can come up with a better method. <lang ecmascript>import "./math" for Int import "./seq" for Lst

var limit = 4 var giuga = [] var n = 4 while (giuga.count < limit) {

   var factors = Int.primeFactors(n)
   var c = factors.count
   if (c > 2) {
       var factors2 = Lst.prune(factors)
       if (c == factors2.count && factors2.all { |f| (n/f - 1) % f == 0 }) {
           giuga.add(n)
       }
   }
   n = n + 2

} 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