little Explanation
Using second order lucas sequences $$U_{n + 2} = P\cdot{U_{n -1}} - Q\cdot{U_{n}}\qquad U_0=0, U_1=1$$ $$V_{n + 2} = P\cdot{V_{n -1}} - Q\cdot{V_{n}},\qquad V_0=2, V_1=P$$ Now our current primality tests use the divisibility property of These sequences to test the primality of number, $n$. Divisibility property only, is not enough to tell the primality of $n$ since some composite number also satisfy this property. Here is the solution.
Lucas Parameter Selection Tips
Let $D = P^2 - 4{Q}$, and Jacobi symbol $\big({\frac{D}{n}}\big)= -1$. $-1$ for Jacobi symbol ensure that the periods of all prime numbers, $p$, is greater than $p$ and the period of all composite numbers equals the multiple of phi of its lagest prime factor. The only exception is when $P$ and $Q$ $=\pm1$
if $\big({\frac{D}{n}}\big)= -1$, $(U_n, U_{n + 1}) = (n - 1, 0)$, this is true for any prime $p = n$ and some composite numbers. And let us call $(P,Q)$ is a solution if $(U_n, U_{n + 1}) = (n - 1, 0)$, and this happens only when $\big({\frac{D}{n}}\big)= -1$.
CONJECTURE:
if $P= 1$ and Let $Q < -1$, $(P, Q)$ are not solutions for composite numbers $n$ when $\left|{Q}\right|< \sqrt{n}$, but all are solution for prime numbers
The above conjecture is the base of this primality test. So now $(P,Q) = (1, -Q)$ and $D =4Q + 1$ The number of tries needed to abtain $\big({\frac{D}{n}}\big)= -1$, is always $<40$ an avarage of $17$, since we are using Jacob symbol, no factorization, its is efficient still
The Test
- Choose $Q$ using above criterion and set $P$ and $Q$
- compute $(U_n, U_{n + 1})$
- check if they satisfy $(U_n, U_{n + 1}) = (n - 1, 0)$,
- if that is True then $n$ is prime
I have tested numbers up to $10^{12}$ and no composite number passed the test
Question
- is this test a deterministic test?, It is True if the above conjecture is true
Here are codes for one want to test
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, 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]
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(p)[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, k = P**2 - 4 * Q, 1, 0
'''
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
k = 2 * k
#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**2) % p
Qk = (Qk * Q) % p
k = k + 1
else:
# doubling only
Uk = V * U
Vk = V**2 - 2 * Qk
U, V = Uk % p, Vk % p
Qk = (Qk**2) % p
k = 2 * k
#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
t = time.time()
for p in range(10000001, 15000000, 2):
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)
y = UV[1][1]
p1 = [[p - 1, 0], [P, y]]
prime = (UV == p1)
if prime:
strAp = str(UV).rjust(60)
strQ = str(-Q).zfill(4)
strA = str([P, strQ, D]).rjust(25)
strP = str(p).rjust(15)
isP = str(is_prime(p)).rjust(6)
print(strAp, strP, isP, strA)
list1.append(-Q)
list2.append(D)
count2 += 1
if not isP:
pseudo.append(p)
print(' ')
t2 = time.time()
print(t - t2)
print(count2, count)
print(list(set(list1)))
print(" ")
l = list(set(list2))
l.sort()
print(l)
print(pseudo)
