Has anyone worked with Giuga numbers?

A positive integer x > 1 with prime factors p1p2p3...pi that satifies the relationship 1/p1 + 1/p2 + 1/p3 +...+ 1/pi - 1/x = k, where k

is a positive integer. The first few Giuga numbers are 30, 858, 1772, and 66,198. For example, for x = 30 the prime factors are

2, 3, 5, so that we have 1/2 + 1/3 + 1/5 - 1/30 = 31/30 - 1/30 = 1 = k.

[Walter Roberson]

The code is not difficult to write if you use the Symbolic toolbox to add the fractions. But it is not all that fast, either, only processes about 20 per second.

[Ken Bannister]

Good tip. The issue at hand is: What are the first 10 unique combinations of primes, and their reciprocals, such that a modified RHS, namely k/(n+1), where n = 6, and 1 <= k <= n, can be found? No problem generating the primes and their unique combinations. And thanks to your tip, I can see where the symbolic toolbox can be used, with some suitable looping and testing, to do the rest of the arithmetic to sift out the 10 unique combinations.

Answer 1

[John D'Errico]

What does this have to do with MATLAB?

Anyway, the answer is easy. Don't use the symbolic toolbox, at least not unless the results become too large to handle using say INT64 arithmetic. And DON'T use fractions. That only makes things problematic. Now you need to worry about floating point arithmetic, or you NEED to use the symbolic toolbox, making things slow. For example, choose some subset of primes. Then multiply by x, in the governing relation. For example...

Plist = primes(20);
testP([2 3 5])

SUCCESS! 30

testP([2 3 7])

Failure! 42

testP([2 3 11 17 59])

SUCCESS! 66198

In your comment, you stated that you wanted to modify the governing relation, with the right hand side as k/(n+1). Again, this would seem trivial. But I won't write that code, since it looks like this must be some sort of assignment, as there is no reason to need to find the first 10 such solutions, UNLESS it is an assignment of some sort. (A possibility is a Project Euler problem, or something like that.) But even then, what happens if you multiply by x? Now look at the expression:

sum(X./P) - 1

For example, with

P = [2 3 7];
X = prod(P)

X = 42

sum(X./P) - 1

ans = 40

We need that result to be of the form X*k/(n+1), for some value of k between 1 and n+1, and apparently n==6. We can determine if that is true by looking carefully at the factors of (sum(X./P)-1). Again, if you stay in integer arithmetic, then things become fast, and you hve no need to play with floating point arithmetic.

The code I wrote in testP is trivial. And because I used simple double precision to compute the results, it will be exceedingly fast. Just loop over possible sets of small primes.

function testP(P)
  X = prod(P);
  if mod(sum(X./P) - 1,X) == 0 % think about why this works
    disp("SUCCESS! " + num2str(X))
  else
    disp("Failure! " + num2str(X))
  end
end

If it gets to the point where you start to exceed 2^53-1, so doubles will fail you to do the arithmetic, you can use int64, or even uint64 as long as you are careful in the code.

[Ken Bannister]

Sir: Thank you so much for showing a path to how this problem can be solved. I am still new

to MATLAB and also learning about primes. I would say that nothing about the primes is "easy"

or necessarily "trivial." Thank goodness we have MATLAB tools that help with laborious

calculations on the primes! I was very surprised to find that, apparently, no one in the brilliant and

diverse MATLAB community seems to have dealt with G-numbers. I was unaware of useful MATLAB functions such as "prod" and "sum," and the "./" notation. So, once again, thanks for pointing them out to me.

By the way, further info about Giuga Numbers can be found in the online

Encyclopedia of Integer Sequences: https://oeis.org/A007850

This Encyclopedia is a powerful assist on many aspects of number theory.

There one can see, as you have rightly pointed out, G-numbers grow very large, very quickly.

IDK if MATLAB, as received, was really meant to deal with so many digits. Perhaps some clever

programming with it has been done to cater for multiple-precision arithmetic. Alternatively, a computer algebra

system, such as Waterloo Maple or Wolfram's Mathematica, might be more appropriate. I have in

my collection somewhere an old 1970s FORTRAN-based multiple-precision arithmetic program by Dr. Fred Crary that could also be used in a pinch.

• KAB

[Paul]

As I'm sure @John D'Errico would tell you, the Symbolic Math Toolbox can be used for the bigger numbers. For example

% from OEIS
X = sym(split("30, 858, 1722, 66198, 2214408306, 24423128562, 432749205173838, 14737133470010574, 550843391309130318, 244197000982499715087866346, 554079914617070801288578559178, 1910667181420507984555759916338506",","))

X = 

for ii = 1:numel(X);
P = factor(X(ii));
success(ii) = isAlways(mod(sum(X(ii)./P) - 1,X(ii)) == 0);
if success(ii)
k(ii) = simplify(sum(1./P) - 1./X(ii));
end
end
k

k = 

success

success = 1×12 logical array
   1   1   1   1   1   1   1   1   1   1   1   1

[Ken Bannister]

Thank you for your insights and interesting programming tips. Perhaps now G-numbers

have finally been added to the community's portfolio !

[Ken Bannister]

I wanted to report on where I stand on computing the 6th cousins of the Giuga numbers. I have attached a table of the numbers found so far. I found right away that 7 must always be included in the list of prime

factors. I found out fairly soon that using the symbolic toolbox approach was in fact too slow, necessitating use of floating point arithmetic.

I have independently checked the cousins 1-8 I have found to verify they do yield the RHS fractions shpown in the table.

I used your "kernel": sum(1./P)-1/X, and just set up a brute force "for" loop in which I ranged over the integers 1-to-n using an index of 7. For a given integer n, I tested whether it had a unique list of prime factors, If it did, I used the "kernel" to test whether it was close to 6/7, 5/7, 4/7, 3/7, 2/7, or 1/7. Of these, I found so far that only the fractions 5/7, 4/7, 3/7, and 2/7 were of interest, moreover, the fractions 5/7 and 4/7 seem to show up most often. The tolerance I used was 1.0E-10. I thought perhaps a tighter tolerance of 1.0E-20 might improve the search for #'s 9 and 10, but that did not help, so I went back to 1.0E-10.

This approach works for cousins 1 through 8, but now #9 and #10 may be so huge I doubt I will ever

find them. Absent something obvious I have missed about "process orlogic" or using MATLAB, I am sure there is some clever approach based on prime number theory, e.g., the various discoveries about the distribution of primes, which I am not aware of that could be used to at least bound these two numbers. If they could be bounded, then I could search therein and find them.