Search code examples
mathsagepari-gp

Checking for a counterexample to a conjecture on even deficient-perfect numbers using Pari-GP


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?


Solution

  • 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.