consider lucas parameters $(P, Q)$ and $D = P^2 - 4Q$. Let $n>0$,$\big(\frac{D}{n}\big)= - 1$ then $U_{n + 1}\equiv{0 \pmod{n}}$ and $n$ is a Lucas probable prime. This test base only on the divisibility property of second order recurrence relation which is not sufficient for primality testing.
After studying alot of third order recurrences, i realized that it is possible to come up with third order deterministic primality test. if we can shift this knowledge to second order the product is the following test
Few Points
- In second order recurrence relation, two terms should be used to determine the primality of $n$, $(U_n, U_{n+1})$
- Let $P$ be fixed and equals any odd number $<\sqrt{n}$ and $Q<-1$. IN my test i used $(P, Q) = (1, Q< - 1)$
- $n$ is prime if and only if $(U_n, U_{n +1}) = (-1, 0)\pmod{n}$
there is no pseudoprime up to $6.5\cdot 10^{14}$, it takes me 3 weeks to test using core !5 hp Pavilion Laptop
here is the codes that you can test
from math import sqrt
def mrange(start, stop, step):
while start < stop:
yield start
start += step
def is_prime(num):
if num == 2:
return True
if (num < 2) or (num % 2 == 0):
return False
return all(num % i for i in mrange(3, int(sqrt(num)) + 1, 2))
def matrixDotMod(a, b, m):
'''
this function multiply two mutrices
number of column in a should be equal to number of rows in b
'''
abDot = []
a_rows, a_columns = len(a), len(a[0])
b_rows, b_columns = len(b), len(b[0])
# check if all internal list length are the same
for row in range(len(a)):
ab_row = []
a_row = a[row]
psum = 0
for column in range(a_rows * b_columns):
a_row_el = a_row[column % a_columns]
b_column_el = b[column % b_rows][column//b_rows]
product = a_row_el * b_column_el
psum = (psum + product) % m
if column % a_columns == a_columns - 1:
ab_row.append(psum)
psum = 0
abDot.append(ab_row)
ab_row = []
return abDot
def matrixExpMod(matrix, exp, prime):
result = []
exp = exp - 1
row, column = len(matrix), len(matrix[0])
for x in range(row):
ab = []
for y in range(column):
if x == y:
ab.append(1)
else:
ab.append(0)
result.append(ab)
ab = []
while(exp > 0):
if exp % 2 == 1:
result = matrixDotMod(result, matrix, prime)
exp = exp - 1
else:
matrix = matrixDotMod(matrix, matrix, prime)
exp = exp//2
return result
def calculateGcd(a, b):
if a > b:
a, b = b, a
while a != 0:
t = a
a = b % a
b = t
return b
def calculateJacob(a, n):
a = a % n
if calculateGcd(a, n) > 1:
return 0
if a == 1:
return 1
else:
if a % 2 == 0:
y = (-1) ** ((n**2 - 1) // 8)
return calculateJacob(a // 2, n) * y
else:
y = ((-1) ** ((a - 1) * (n - 1) // 4))
return calculateJacob(n % a, a) * y
def LucasUVthTerm(A, t, p):
'''
This function calculate nth U and V term in Lucas sequences
by manipulating bit string of number being tested for primality
Keyword arguments:
A -- list, carry lucas parameters, P and Q, [P, Q]
t -- nth term to be calculated
p -- int, number to be tested for primality
it returns Lucas Terms in a list, [[U(p), U(p + 1)], [V(p), V(p + 1)]]
'''
#string of bits of number p and its length
bitString = str(bin(t)[2:])
length = len(bitString)
#set Lucas parameters, initial terms, Uo, Vo, D, P and Q
#from lucas identity, V(2k) = V^2 - 2 * Q^k, need to compute Q^k eerytime we double k
#For comutational effeciency, we use modulo exponentions
#instead of computing Q^k, set Qk, if we are doubling, set Qk = (Qk**2) % p
#if we are adding 1 to get V(2k + 1), set Qk = (Qk * Q) % p
U, V, P, Q = 0, 2, A[0], A[1]
D, Qk = P**2 - 4 * Q, 1
'''
BIT MANIPULATION OF P TO CREATE LUCAS CHAIN.
Lucas chain help us to know when to double and when to add one to 0btain next terms
example:let nth term be 24 and 112, below are its lucas chains
24 <- 12 <- 6 <- 3 <- 2 <- 1 <- 0
112 <- 61 <- 60 <- 30 <- 15 <- 14 <- 7 <- 6 <- 3 <- 2 <- 1 <- 0
We need to read this sequence from right to left.
here is how we do it without pre calculation of chain
1. we double the value add one when nth bit in the bitstring of p is 1,
so one carry two processes
2. we only double when nth bit in bitstring is zero
Note our terms start from zero index
'''
for y in range(length):
bit = bitString[y]
# when nth bit is !, double and add one
if bit == "1":
#doubling from k to 2k
Uk = V * U
Vk = V**2 - 2 * Qk
Qk = (Qk**2) % p
#adding one from 2k to 2k + 1
Uk1 = (P * Uk + Vk) * (p + 1) // 2
Vk1 = (D * Uk + P * Vk) * (p + 1) // 2
#cahnging alues of our Lucas parameters
U, V = Uk1 % p, Vk1 % p
Qk = (Qk * Q) % p
else:
# doubling only
Uk = V * U
Vk = V**2 - 2 * Qk
U, V = Uk % p, Vk % p
Qk = (Qk**2) % p
#since we need to term for primality testing, U(p) and U(p + 1)
#add one to the U(p) terms to get U(p + 1)
U1 = ((P * U + V) * (p + 1) // 2) % p
V1 = ((D * U + P * V) * (p + 1) // 2) % p
return [[U, U1], [V, V1]]
list1 = []
list2 = []
pseudo = []
count, count2 = 0, 0
import time
t1 = time.time()
f = open("carmichael.txt", 'r')
for x in range(600000000000001, 650000000000000, 2):
p = 0
if 0:
t = f.readline()
if t == '':
break
p = int(t[len(str(x)):])
else:
p = x
D, Jacob = 0, 0
P, Q = 1, 0
import math
R = math.floor(math.sqrt(p))
#check if a number is a square
if R**2 < p:
for n in range(2, p):
D = P**2 - 4 * (-n)
Jacob = calculateJacob(D, p)
if Jacob == -1:
Q = -n
break
A = [P, Q]
UV = LucasUVthTerm(A, p, p)
y = UV[1][1]
p1 = [[p - 1, 0], [1, p + 2 * Q]]
prime = (UV == p1)
if prime:
strAp = str(UV).rjust(55)
strQ = str(-Q).zfill(4)
strA = str([P, strQ, D]).rjust(25)
strP = str(p).rjust(25)
#isP = str(is_prime(p)).rjust(4)
print(strAp, strP, strA)
list1.append(-Q)
list2.append(D)
count2 += 1
print(' ')
t2 = time.time()
print(t2 - t1)
print(count2, count)
st = list(set(list1))
st.sort()
print(st, len(st))
print(" ")
l = list(set(list2))
l.sort()
print(l)
p = 319889369713946602502766595032347
P, Q = 1, 0
for n in range(2, p):
D = P**2 - 4 * (-n)
Jacob = calculateJacob(D, p)
if Jacob == -1:
Q = -n
break
A = [P, Q]
print(A)
print(LucasUVthTerm(A, p, p))
Trying to summarize, take the standard Baillie/Wagstaff (1980) Lucas pseudoprime definition and test from page 1, use their recommendation of $(D|n)=-1$ (page 2), and use P = 1, Q < -1 as the parameter selection. Then $U_{n+1} = 0 \pmod n$ is tested just like the paper says.
Add the requirement that $U_n = -1 \pmod n$. This property is noted on the Wikipedia Lucas Sequence page where the list of divisibility properties includes: "if p is an odd prime then $U_p = (D|p) \pmod p$ and $V_p = P \pmod p$", where our parameter selection makes the right hand sides $-1$ and $1$ respectively.
We can intelligently implement the two conditions on $U_n$ and $U_{n+1}$ as a single Lucas sequence for $U_n$ and a simple step to get $U_{n+1}$.
It is interesting in that there aren't counterexamples for smallish numbers (e.g. $n < 10^{12}$), and known lists of PSPs and various Lucas-type pseudoprimes also do not show counterexamples. Surprising for such a simple test. We could add some additional conditions that are likely to be nearly computationally free, such as $V_n = 1 \pmod n$ and $V_{n+1} = 2Q \pmod n$, as the paper mentions in the remarks of section 5. I just noticed your Python code is testing both of these in the carmichael section.
In terms of speed, it doesn't seem to be as fast as the extra strong test due to both optimization with Q=1 as well as fewer steps due to removal of factors of 2 from n+1. But of course the ES test has some pseudoprimes so needs additional conditions or a combined test. Using the SPSP-2 test first, like BPSW, is of course faster for almost all composites though could lengthen the time for primes and base 2 pseudoprimes. In my C implementation, BPSW was over 2x faster for odd composites and just slightly faster for primes. I cannot guarantee the code is optimal however.
Unfortunately I think this is in the same state as similar tests like the Frobenius-Khashin and Frobenius-Underwood tests. They have been tested to fairly high values but there is no proof that they are deterministic, nor honestly do we have a strong belief that they are. In practice the BPSW test gives faster and more trustworthy results -- it is known deterministic for all 64-bit inputs, with no larger counterexamples found after extensive study and use.