PRIMALITY.sagews

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] )

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

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

7919

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

True

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

[1, 4, 6, 8, 9, 10, 12, 14, 15, 16]

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

True

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

179769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393357658789768814416622492847430639474124377767893424865485276302219601246094119453082952085005768838150682342462881473913110540827237163350510684586298239947245938479716304835356329624224137859

%time isprime( p )

True CPU time: 2.76 s, Wall time: 2.83 s

%time isprime( c )

False CPU time: 0.00 s, Wall time: 0.00 s

%time p in Primes( )

True CPU time: 5.70 s, Wall time: 5.80 s

%time c in Primes( )

False CPU time: 0.00 s, Wall time: 0.00 s