Working with primes in Python -- Updated.

primes.py

from fractions import gcd
from random import randrange

PRIMES_LE_31 = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31)
PRIMONIAL_31 = 200560490130


def invmul(x, mod):
    if mod <= 0:
        raise ValueError('Modulus must be greater than zero.')
    a = abs(x)
    b = mod
    c, d, e, f = 1, 0, 0, 1
    while b > 0:
        q, r = divmod(a, b)
        g = c - q * d
        h = e - q * f
        a, b, c, d, e, f = b, r, d, g, f, h
    if a != 1:
        raise ValueError('{x} has no multiplicative '
                         'inverse modulo {m}'.format(x=x, m=mod))
    return c * -1 if x < 0 else c * 1


def isqrt(n):
    if n < 0:
        raise ValueError('Square root is not defined for negative numbers.')
    x = int(n)
    if x == 0:
        return 0
    a, b = divmod(x.bit_length(), 2)
    n = 2 ** (a + b)
    while True:
        y = (n + x // n) >> 1
        if y >= n:
            return n
        n = y


def is_square(n):
    s = isqrt(n)
    return s * s == n


def factor(n, p=2):
    s = 0
    d = n - 1
    q = p

    while not d & q - 1:
        s += 1
        q *= p

    return s, d // (q // p)


def jacobi(a, p):
    if (not p & 1) or (p < 0):
        raise ValueError('p must be a positive odd number.')

    if (a == 0) or (a == 1):
        return a

    a = a % p
    t = 1

    while a != 0:
        while not a & 1:
            a >>= 1
            if p & 7 in (3, 5):
                t = -t

        a, p = p, a
        if (a & 3 == 3) and (p & 3) == 3:
            t = -t

        a = a % p

    if p == 1:
        return t

    return 0


def selfridge(n):
    d = 5
    s = 1
    ds = d * s

    while True:
        if gcd(ds, n) > 1:
            return ds, 0, 0

        if jacobi(ds, n) == -1:
            return ds, 1, (1 - ds) / 4

        d += 2
        s *= -1
        ds = d * s


def chain(n, u1, v1, u2, v2, d, q, m):
    k = q
    while m > 0:
        u2 = (u2 * v2) % n
        v2 = (v2 * v2 - 2 * q) % n
        q = (q * q) % n

        if m & 1 == 1:
            t1, t2 = u2 * v1, u1 * v2
            t3, t4 = v2 * v1, u2 * u1 * d
            u1, v1 = t1 + t2, t3 + t4

            if u1 & 1 == 1:
                u1 = u1 + n

            if v1 & 1 == 1:
                v1 = v1 + n

            u1, v1 = (u1 / 2) % n, (v1 / 2) % n
            k = (q * k) % n

        m = m >> 1

    return u1, v1, k


def strong_pseudoprime(n, a, s=None, d=None):
    if not n & 1:
        return False

    if (s is None) or (d is None):
        s, d = factor(n, 2)

    x = pow(a, d, n)

    if x == 1:
        return True

    for i in xrange(s):
        if x == n - 1:
            return True

        x = pow(x, 2, n)

    return False


def lucas_pseudoprime(n):
    if not n & 1:
        return False

    d, p, q = selfridge(n)
    if p == 0:
        return n == d

    u, v, k = chain(n, 0, 2, 1, p, d, q, (n + 1) >> 1)
    return u == 0


def strong_lucas_pseudoprime(n):
    if not n & 1:
        return False

    d, p, q = selfridge(n)
    if p == 0:
        return n == d

    s, t = factor(n + 2)

    u, v, k = chain(n, 1, p, 1, p, d, q, t >> 1)

    if (u == 0) or (v == 0):
        return True

    for i in xrange(1, s):
        v = (v * v - 2 * k) % n
        k = (k * k) % n
        if v == 0:
            return True

    return False


def miller_rabin(n, k=10):
    if n == 2:
        return True

    if not n & 1:
        return False

    if n < 2 or is_square(n):
        return False

    if gcd(n, PRIMONIAL_31) > 1:
        return n in PRIMES_LE_31

    s, d = factor(n)

    for i in xrange(k):
        a = randrange(2, n - 1)
        if not strong_pseudoprime(n, a, s, d):
            return False

    return True


def baillie_psw(n, limit=100):
    if n == 2:
        return True

    if not n & 1:
        return False

    if n < 2 or is_square(n):
        return False

    if gcd(n, PRIMONIAL_31) > 1:
        return n in PRIMES_LE_31

    bound = min(limit, isqrt(n))
    for i in xrange(3, bound, 2):
        if not n % i:
            return False

    return strong_pseudoprime(n, 2) \
        and strong_pseudoprime(n, 3) \
        and strong_lucas_pseudoprime(n)


def next_prime(n):
    if n < 2:
        return 2

    if n < 5:
        return [3, 5, 5][n - 2]

    gap = [1, 6, 5, 4, 3, 2, 1, 4, 3, 2, 1, 2, 1, 4, 3, 2, 1, 2, 1, 4, 3, 2, 1,
           6, 5, 4, 3, 2, 1, 2]

    n += 1 if not n & 1 else 2

    while not baillie_psw(n):
        n += gap[n % 30]

    return n

In the strong Lucas checking code, you've used:

    for i in xrange(1, s):
        v = (v * v - 2 * k) % n
        k = (k * k) % n
        if v == 0:
            return True

You are implementing the subscript doubling recurrence relation there, which is correct.
However, the value of subscript k is not doubled, but squared in your code.
I think that should be

        k = (k << 1) % n            or of course            k = (2 * k) % n

Kind regards,
Maurice

Hello, where should I add an int(input()) to be able to input any number of my choice? It would be very important for me. Preferably an exact line, but just an "at the start" or "at the end" works too. Thank you in advance, whoever replies!

@mmakaay Any guesses or ideas? No offense

@pleasehugme line 5 is the primorial prime of the product of primes 2 to 31. It’s used in the baille_psw check gcd(n, PRIMORIAL_31)

https://en.wikipedia.org/wiki/Primorial_prime

The two primality checks here are baillie_psw and miller_rabin, so you use a user’s input it would be

n = int(input())

is_prime = baillie_psw(n)
# or
is_prime = miller_rabin(n)