Of course, since n

_{0} < f/2, the cofactor (2*n

_{0}^2 - 1)/f < f/2. So if the cofactor is prime, it is certainly less than f.

Since I could think of no reason why the cofactor should always be prime, as a programming exercise I wrote a mindless script to look for counterexamples for 2^p - 1 with small prime exponent p.

I told it to print out the exponent p, the smallest factor f of 2^p - 1, the value n

_{0}, and the factorization of the cofactor (2*n

_{0}^2 - 1)/f when this was composite.

The smallest exponent giving a counterexample is p = 47. To my surprise, with the smallest example 47 (f = 2351, n

_{0} = 240) and 113 (f = 3391, n

_{0} = 700) the cofactor (2*n

_{0}^2 - 1)/f is the square of a single prime. The smallest prime exponent for which the cofactor (2*n

_{0}^2 - 1)/f has at least two distinct prime factors is 59.

Code:

? {
forprime(p=3,200,
M=factor(2^p-1);
f=M[1,1];
m=factormod(2*x^2-1,f);
n=lift(polcoeff(m[1,1],0,x));
if(n+n>f,n=f-n);
N=factor((2*n^2-1)/f);
if(#N[,1]>1||(#N[,1]==1&&N[1,2]>1),print(p" "f" "n" "N))
)
}
47 2351 240 Mat([7, 2])
59 179951 77079 [7, 1; 9433, 1]
67 193707721 66794868 [191, 1; 241177, 1]
71 228479 76047 [23, 1; 31, 1; 71, 1]
97 11447 5670 [41, 1; 137, 1]
101 7432339208719 3616686326055 [7, 1; 17, 1; 23, 1; 1583, 1; 812401, 1]
103 2550183799 270087243 [23, 1; 241, 1; 10321, 1]
109 745988807 298036466 [17, 1; 14008369, 1]
113 3391 700 Mat([17, 2])
137 32032215596496435569 6857964810884905735 [2503, 1; 358079, 1; 3276376633, 1]
151 18121 2513 [17, 1; 41, 1]
163 150287 31486 [79, 1; 167, 1]
173 730753 162850 [7, 1; 10369, 1]
179 359 170 [7, 1; 23, 1]
193 13821503 2664653 [7, 1; 146777, 1]
199 164504919713 50650852663 [7, 1; 4455809207, 1]
?