codepad

Language:

# lucas pseudoprimality test

def iSqrt(n):
    if type(n) not in (int, long): raise TypeError, "must be integer"
    if n < 1: raise ValueError, "must be positive"
    x = n
    while 1:
        y = (x + n // x) // 2
        if x <= y: return x
        x = y

def isSquare(n):
    n2 = iSqrt(n)
    n == n2 * n2

def gcd(a, b):
    while b != 0:
        a, b = b, a % b
    return a

def jacobi(a, m):
    a, t = a % m, 1
    while a != 0:
        while a % 2 == 0:
            a = a / 2
            if m % 8 in (3, 5):
                t = -t
        a, m = m, a
        if a % 4 == 3 and m % 4 == 3:
            t = -t
        a = a % m
    if m == 1:
        return t
    return 0

def isStrongPseudoprime(n, a):
    d, s = n - 1, 0
    while d % 2 == 0:
        d, s = d / 2, s + 1
    t = pow(a, d, n)
    if t == 1: return True
    while s > 0:
        if t == n - 1: return True
        t, s = (t * t) % n, s - 1
    return False

def selfridge(n):
    dAbs, sign, d = 5, 1, 5
    while 1:
        if 1 < gcd(d, n): return d, 0, 0
        if jacobi(d, n) == -1: break
        dAbs, sign = dAbs + 2, -sign
        d = dAbs * sign
    return d, 1, (1 - d) / 4

def chain(n, u, v, u2, v2, d, q, m):
    qkd = q
    while m > 0:
        u2, v2 = (u2 * v2) % n, (v2 * v2 - 2 * q) % n
        q = (q * q) % n
        if m % 2 == 1:
            t1, t2, t3, t4 = u2 * v, u * v2, v2 * v, u2 * u * d
            u = t1 + t2; u = u + n if u % 2 == 1 else u; u = (u / 2) % n
            v = t3 + t4; v = v + n if v % 2 == 1 else v; v = (v / 2) % n
            qkd = (qkd * q) % n
        m = m // 2
    return u, v, qkd

def isStandardLucasPrime(n):
    d, p, q = selfridge(n)
    if p == 0: return n == d
    u, v, qkd = chain(n, 0, 2, 1, p, d, q, (n + 1) / 2)
    return u == 0

def isStrongLucasPrime(n):
    d, p, q = selfridge(n)
    if p == 0: return n == d
    s, t = 0, n + 1
    while t % 2 == 0:
        s, t = s + 1, t / 2
    u, v, qkd = chain(n, 1, p, 1, p, d, q, t // 2)
    if u == 0 or v == 0: return True
    for r in xrange(1, s):
        v = (v * v - 2 * qkd) % n
        if v == 0: return True
        qkd = (qkd * qkd) % n
    return False

def isPrime(n):
    if type(n) not in (int, long): raise TypeError, "must be integer"
    if n < 2: return False
    if isSquare(n): return False
    for p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, \
              43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]:
        if n % p == 0: return n == p
    return isStrongPseudoprime(n, 2) \
       and isStrongPseudoprime(n, 3) \
       and isStrongLucasPrime(n)

for n in [79, 83, 89, 111, 323, 5777, 3825123056546413051, 2**89-1]:
    print n, isStrongPseudoprime(n, 2), \
             isStrongPseudoprime(n, 3), \
             isStandardLucasPrime(n), \
             isStrongLucasPrime(n), \
             isPrime(n)

Private [?]

Run code