# N+1 Primality Proving is slow

I am trying to implement the N+1 method of proving primality. Here is my description, based on the Brillhart, Lehmer, Selfridge paper, Theorem 13 and Corollary 8:

Choose P and Q such that D = P^2 - 4*Q is not a square modulo N. Let N+1 = F*R with F > R, where R is odd and the prime factorization of F is known. If there exists a Lucas sequence of discriminant D with U(N+1) == 0 (mod N) and gcd(U((N+1)/q), N) = 1 for each prime q dividing F, then N is prime; if no such sequence exists for a given P and Q, a new P' and Q' with the same D can be computed as P' = P + 2 and Q' = P + Q + 1 (the same D must be used for all the factors q).

I compute D as the first number in the sequence 5, -7, 9, -11, ... for which the jacobi symbol (D/N) = -1, with P = 1 and Q = (1-D)/4, and find the primes q dividing F by trial division. Here is my code, in Python, which you can see and run at ideone.com/HVXIGR; function u(p, q, m, n) returns the n'th element of the Lucas sequence U(p,q) modulo m:

def provePrime(n):
f, fs, fp, r, d = 2, [], 1, n+1, 5
while fp*fp < n:
while r % f == 0:
r /= f; fp *= f
if f not in fs:
fs.append(f)
f += 1
while jacobi(d, n) <> -1:
d = (abs(d)+2) * signum(d) * -1
p, q = 1, (1-d)/4
print fs, d, p, q
if gcd(d, n) > 1: return False
if u(p, q, n, n+1) <> 0: return False
for x in fs:
while True:
if gcd(u(p,q,n,(n+1)/x), n) == 1: break
p, q = p + 2, p + q + 1
return True


The problem is that for some N this is very slow because the final while loop that searches for P and Q takes a very long time. For instance, it takes over two minutes to prove the primality of 10^12+61 with D=-7.

So my questions: 1) Is my description of the algorithm correct? 2) Does my code correctly implement the algorithm? 3) How can I make it faster? 4) Can I make it faster by switching to a new D and restarting the computation? 5) If so, how can I determine when to abandon the current D and start over?

EDIT WITH SOLUTION:

This code works; you can see it running at ideone.com/HVXIGR:

def provePrime(n):
f, fs, fp, r, d = 2, [], 1, n+1, 5
while fp*fp < n:
if f*f > r:
fs.append(r)
break
while r % f == 0:
r /= f; fp *= f
if f not in fs:
fs.append(f)
f += 1
while jacobi(d, n) <> -1:
d = (abs(d)+2) * signum(d) * -1
p, q = 1, (1-d)/4
if gcd(d, n) > 1: return False
if u(p, q, n, n+1) <> 0: return False
for x in fs:
while True:
if gcd(u(p,q,n,(n+1)/x), n) == 1:
print x, p, q
break
p, q = p + 2, p + q + 1
return True


The factorization is mildly tricky. R (initially N+1) is the number currently being factored, F is the current trial factor, FS is the current list of factors without duplicates, and FP is the current product of the factors including multiplicity (so FP in the code corresponds to F in the description). If F*F > R, then the current R is prime, so it is the last factor of N+1 and factoring is complete. Alternately, if FP*FP > N, then even though factorization isn't complete, it is sufficient to prove primality, so the factoring can be terminated.

• If you are interested in speed, why are you using Python? Commented Feb 5, 2014 at 0:58
• It is convenient to work out the algorithm in a language like Python, then implement it for real in C/GMP. Commented Feb 5, 2014 at 3:15
• Notice Remark 2 after Theorem 13: From one Lucas sequence with $P_1, Q_1, \text{and }D$, another with the same D can be obtained by setting $P_2= P_1 + 2\;$ and $Q_2 = P_1+ Q_1+ 1$. It is necessary to check that $(N,Q_i) = 1$. You do not check this gcd, but I do not know how much this would speed up your algorithm. Another obvious speed up would be to use only 2 and the odd f for trial in the first part. Commented Feb 5, 2014 at 12:17
• Your initial $p$ and $q$ are 1 and 2, respectively. What are the ending values? Commented Feb 5, 2014 at 12:38
• @gammatester: Well, it's not really necessary. If (N, Qi) > 1, then Qi | N and we have found a factorization. But assuming N is prime, that will never happen. Certainly before attempting to prove N is prime, you have done enough work to suggest that it has no small factors. Commented Feb 5, 2014 at 13:22