I am trying to check for counterexamples to the conjecture stated in this MSE question, using the Pari-GP interpreter of Sage Cell Server.
I reproduce the statement of the conjecture here: If N > 8 is an even deficient-perfect number and Q = N/(2N - sigma(N)), then Q is prime.
Here, sigma(N) is the classical sum of divisors of N.
I am using the following code:
for(x=9, 1000, if(((Mod(x,(2*x - sigma(x))) == 0)) && ((fromdigits(Vecrev(digits(x / (2*x - sigma(x)))))) == (x / (2*x - sigma(x)))) && !(isprime((x / (2*x - sigma(x))))), print(x,factor(x))))
However, the Pari-GP interpreter of Sage Cell Server would not accept it, and instead gives the following error message:
*** at top-level: for(x=9,1000,if(((Mod(x,(2*x-sigma(x)))==0))&&
*** ^----------------------------
*** Mod: impossible inverse in %: 0.
What am I doing wrong?
Here's a better implementation of your algorithm
{
forfactored(X = 9, 10^7,
my (s = sigma(X), t = 2*X[1] - s);
if (t <= 0, next);
my ([q, r] = divrem(X[1], t));
if (r == 0 && fromdigits(Vecrev(digits(q))) == q && !ispseudoprime(q),
print(X)))
}
It's a bit more readable but most importantly it avoids factoring the same x
over and over again: each time you write sigma(x)
, we need to factor x
(the interpreter is not clever enough to compute subexpressions once). In fact, it doesn't perform a single factorization, through the use of forfactored
which performs a sieve instead (and the variable X
contains [x, factor(x)]
). This is about 3 times faster than the original implementation in this range.
I let it run to 10^9
(about 10 minutes), there was no further counterexample.