from random import randint
import heapq

def PRIME ( k , p ) :

    """
        p is odd
        see ed. 3, p. 401
    """

    _1 = p - 1

    for i in range( k ) :

        a = randint( 1 , _1 )

        if power_mod( a , _1 , p ) != 1 : return False

        s = _1

        while not s & 1 :

            x = power_mod( a , s , p )

            if x != 1 :

                if x != _1 : return False

                break

            s //= 2

    return True

def isprime ( p ) :

    if p <= 1 : return False

    if p == 2 : return True

    if not p & 1 : return False

    d = 2

    e = 1 / 2**d

    a = -log( 4 * e * ( 1 - e ) , 2 )

    f = lambda n : n / 8

    n = ceil( log( p , 2 ) )

    k = int( ceil( f( n ) / a ) )

    return PRIME( d * k , p )

n = 1000
primes = primes_first_n( n )
print( primes[:10] )

last = primes[-1]
primes = set( primes )
last

all( map( isprime , sorted( primes ) ) )

composites = set( range( 1 , last ) ) - primes
print( heapq.nsmallest( 10 , composites ) )

not any( map( isprime , sorted( composites ) ) )

p = next_prime( 2**1024 )
c = p - 2
print( p )

%time isprime( p )

%time isprime( c )

%time p in Primes( )

%time c in Primes( )