Has anyone worked with Giuga numbers?
4 views (last 30 days)
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.
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])
testP([2 3 7])
testP([2 3 11 17 59])
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)
sum(X./P) - 1
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.
X = prod(P);
if mod(sum(X./P) - 1,X) == 0 % think about why this works
disp("SUCCESS! " + num2str(X))
disp("Failure! " + num2str(X))
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.