# Has anyone worked with Giuga numbers?

3 views (last 30 days)
Ken Bannister on 31 Dec 2021
Commented: Ken Bannister on 17 Jan 2022
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.
##### 2 CommentsShowHide 1 older comment
Ken Bannister on 4 Jan 2022
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.

John D'Errico on 5 Jan 2022
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 on 17 Jan 2022
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.

R2021b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!