from math import sqrt, ceil
from itertools import compress

TRIAL_DIVISION_UNTIL=1000


def isPrime(n):
    """This primality check implements the composite Baillie-PSW test.
       See: https://en.wikipedia.org/wiki/Baillie%E2%80%93PSW_primality_test

       Input:
           n: the number to check
       Output:
           a tuple containing two elements:
           - True if (probable) prime, False otherwise
           - The name of the test that lead to that conclusion
       """
    if not isInteger(n):
        return (False, "not an integer")
    if isBelowLowestPrime(n):
        return (False, "below lowest prime")
    if isLowPrime(n):
        return (True, "low prime")
    if isTrialDivisionCompound(n):
        return (False, "trial division")
    if not isRabinMillerProbablePrime(n):
        return (False, "rabin-miller")
    if isPerfectSquare(n):
        return (False, "perfect square")
    if not isLucasProbablePrime(n):
        return (False, "lucas")
    return (True, "passed all tests")


def isInteger(n):
    return int(n) == n
        
        
def isBelowLowestPrime(n):
    return n < 2


def isLowPrime(n):
    return (
        n <= TRIAL_DIVISION_UNTIL and
        n in lowPrimes
    )


def isTrialDivisionCompound(n):
    return any(
        n % prime == 0
        for prime in lowPrimes
        if prime < n
    )


def isRabinMillerProbablePrime(n):
    """Based on pseudo code for the Rabin-Miller algorithm.
       See: https://en.wikipedia.org/wiki/Miller-Rabin_primality_test
       This simple version implements only a single pass using base 2,
       as prescribed by the Baillie-PSW specification."""
    (factor, exponent) = factorizeBase2(n - 1)
    x = pow(2, factor, n)
    if x == 1 or x == n - 1:
        return True
    for _ in range(exponent - 1):
        x = pow(x, 2, n)
        if x == 1:
            return False
        elif x == n - 1:
            return True
    return False


def isPerfectSquare(n):
    """A quick check to make sure that no perfect squares are fed to the
       Lucas probable prime check."""
    return isqrt(n) ** 2 == n


def isqrt(n):
    """Integer square root, see:
      https://en.wikipedia.org/wiki/Integer_square_root"""
    x = n
    y = (x + 1) >> 1
    while y < x:
        x = y
        y = (x + n // x) >> 1
    return x


def isLucasProbablePrime(n):
    """Code based on "Implementing a Lucas probable prime test" to be found at:
       https://en.wikipedia.org/wiki/Lucas_pseudoprime"""

    def getLucasSequenceInputAsSuggestedBySelfridge():
        seq = ((-1 if i&1 else 1) * (i*2 + 5) for i in range(0, n>>1))
        D = next((D for D in seq if jacobi(D, n) == -1))
        Q = (1 - D) // 4
        P = 1
        return D, Q, P

    D, Q, P = getLucasSequenceInputAsSuggestedBySelfridge()

    def computeLucasValue(subscript):
        U = V = k = 0
        for digit in format(subscript, "b"):
            U, V, k = doubleK(U, V, k)
            if digit == "1":
                U, V, k = incrementK(U, V, k)
        return U, V, subscript

    def doubleK(U, V, subscript):
        return (U*V) % n,  (pow(V, 2, n) - 2 * pow(Q, subscript, n)) % n, subscript << 1

    def incrementK(U, V, subscript):
        U, V = (P*U + V), (D*U + P*V)
        makeEven = lambda x: x+n if x&1 else x
        div2modn = lambda x: (x >> 1) % n
        return (div2modn(makeEven(U)), div2modn(makeEven(V)), subscript+1)

    # Weak check.
    U, V, k = computeLucasValue(n+1)
    if U != 0:
        return False

    # Strong check.
    (factor, exponent) = factorizeBase2(n + 1)
    U, V, k = computeLucasValue(factor)
    if U == 0 or V == 0:
        return True
    for r in range(exponent-1):
        U, V, k = doubleK(U, V, k)
        if V == 0:
            return True


def factorizeBase2(factor, exponent=0):
    """Factor out 2 from a number. For example factorizeBase2(80) returns (5, 4),
       which indicates that 2 ** 5 + 4 == 80."""
    if factor&1:
        return (factor, exponent)
    return factorizeBase2(factor >> 1, exponent + 1)


def jacobi(a, n):
    """Compute the Jacobi symbol. The code uses ideas from:
       https://www.academia.edu/1012630/A_binary_algorithm_for_the_Jacobi_symbol"""
    t = 1
    while a:
        if a < 0:
            a = -a
            if n&3 == 3: t = -t
        while not a&1:
            a = a >> 1
            if n&7 == 3 or n&7 == 5: t = -t
        if (a < n):
            a, n = n, a
            if a&3 == 3 and n&3 == 3: t = -t
        a = (a - n) >> 1
        if n&7 == 3 or n&7 == 5: t = -t
    return t if n == 1 else 0


def sieve(until):
    """This is an optimized sieve-algorithm that I created for an exercism.io assignment.
       On my system, it is able to produce primes until 10,000,000 in half a second.
       See: http://exercism.io/submissions/a6085623b76a6747d300420e"""
    if until < 2:
        return []
    store_len = until >> 1
    is_prime = [True] * store_len
    for offset, number in enumerate(range(3, int(sqrt(until)) + 1, 2)):
        if is_prime[offset]:
            first_compound_at = (number*number-3) >> 1
            nr_of_compounds = int(ceil((store_len-first_compound_at) / number))
            is_prime[first_compound_at::number] = [False] * nr_of_compounds
    return [2] + list(compress(range(3, until + 1, 2), is_prime))


lowPrimes = set(sieve(TRIAL_DIVISION_UNTIL))