1
$\begingroup$

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))
$\endgroup$
  • $\begingroup$ You describe using $U_n$ and $U_{n-1}$ but it looks like your code uses $U_n$ and $U_{n+1}$. The former has clear failures on primes so can't be right. $\endgroup$ – DanaJ Apr 13 '18 at 17:40
  • $\begingroup$ Sorry it is U(n+1) $\endgroup$ – Esdras E E Dansha Apr 14 '18 at 17:36
  • $\begingroup$ You exhaustively tested all integers up to $6.5 \times 10^{14}$? The Python code is very slow for me, taking something like 3 minutes just to test all odds up to $10^5$. In C+asm it is over 1000x faster but that's still a lot of work. $\endgroup$ – DanaJ Apr 16 '18 at 20:57
  • $\begingroup$ I think this slowness is caused by one of the function that i am using to confirm the primarity of number in the test. This function test all odd number up to sqrt(n). We can improve it by using BPSW test $\endgroup$ – Esdras E E Dansha Apr 17 '18 at 8:11
  • $\begingroup$ Your is_prime (and mrange, matrixDotMod, matrixDotExp) aren't called at all. Two changes that make it run about 4x faster are replacing the Jacobi function with a non-recursive version, and inserting a perfect square test in the Q-finding loop. $\endgroup$ – DanaJ Apr 17 '18 at 17:42
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ In PBSW test they use spsp2 test, but this test does not need any spsp2 test... i think this test is very strong and this may be deterministic. Let us wait for proof $\endgroup$ – Esdras E E Dansha Apr 14 '18 at 16:41
  • $\begingroup$ Can you give me one counterexample of composite that pass this test? $\endgroup$ – Esdras E E Dansha Apr 14 '18 at 16:52
  • 1
    $\begingroup$ It seems to be strong, but there is no evidence it is deterministic (i.e. that it is provably correct for all inputs, not in the sense of a single computation path). I have no counterexamples. Nor do we for BPSW after 38 years of daily use in many applications, and much more exhaustive testing of small numbers. While this doesn't need SPSP-2, my point was that BPSW ends up being faster overall. The same argument applies to the Khashin and Underwood Frobenius tests which are also quite fast and similarly have no known counterexamples. $\endgroup$ – DanaJ Apr 16 '18 at 3:07
  • $\begingroup$ PBSW cant be faster than this test, because of SPSP2 test. The only thing here is bad implementation of this test and we dont test for small divisors before running the test $\endgroup$ – Esdras E E Dansha Apr 16 '18 at 8:12
  • $\begingroup$ SPSP-2 test removes almost all composites in less than half the time of this test, so of course it is faster for those cases. The AES and ES Lucas test use Q=1 which means Q and Qk drop entirely, saving 2-3 mulmods per loop, which is substantial. $\endgroup$ – DanaJ Apr 16 '18 at 16:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.