0

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?

Charles
  • 11,269
  • 13
  • 67
  • 105
Jose Arnaldo Bebita Dris
  • 103
  • 3

2 Answers2

2

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.

K.B.
  • 861
  • 5
  • 14
0

I got it to work myself.

Here is the code that I used:

for(x=9, 10000000, if((2*x > sigma(x)) && ((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))))

The search returns the odd counterexample N = 9018009, which is expected.

It did not return any even counterexamples, in the specified range.

Jose Arnaldo Bebita Dris
  • 103
  • 3