import random, decimal, time
from fractions import gcd

D = decimal.Decimal
decimal.getcontext().prec = 50

#v11 = 5689 71935 26942 02437 03269 72321

def tword( p ):
    """ Changes the number to the form 2^r * d.  Outputs
        r and d in a list. [r, d]"""

    #possible prime to be sorted to form p = 2^r * d
    r = 0
    d = 0

    while (p % pow(2, r+1)) == 0:
        r += 1
    d = p / pow(2, r)
    return [r, d]

def testA( p, a, r, d ):

    rd = tword( p - 1 )
    r = rd[0]
    d = rd[1]
    #print "r: ", r
    #print "d: ", d

    i = 1;

    test = pow(a, d, p)

    if test == 1 or test == p - 1:   #if (p^d != 1) mod p
        return 1

    #print "i:", i
    #print "r:", r
    while i < r:
        #print pow(a, pow(2, i)*d, p)
        if pow(a, pow(2, i)*d, p) == p-1:  #if ( a^(2^r *d) != (p - 1) ) mod p
            return 1
        i += 1
    return 0

def millerRabin( p, k ):
    i = 0;

    rd = tword( p - 1 )
    r = rd[0]
    d = rd[1]
    elapsed = 0;
    gcd_elapsed = 0;
    start = time.clock();

    while i < k:
        #print "Spot"
        a = random.randint(2, p-1)
        if a%2 == 0:
            a += 1
        start_segment = time.clock();
        if gcd( p, a ) != 1:
            print "This is composite."
            return 0
        end_segment = time.clock();
        gcd_elapsed += end_segment - start_segment


        if testA(p, a, r, d) == 1:    #if base returns prime
            i += 1              #try another base
        else:
            print "Not prime. a =", a;
            return 0
    end = time.clock()
    elapsed = end - start;
    #print "The gcd function took", (gcd_elapsed / elapsed)*100, "% of \
#the total time."
    print "It took", elapsed, " seconds."
    print "Probably prime, with", D(D(1) - D( str(pow(4, -k)) ))*D(100), " % certainty."
    return 1

def millerRabinBases( p, bases ):
    print "bases"
    i = 0;

    rd = tword( p - 1 )
    r = rd[0]
    d = rd[1]
    elapsed = 0;
    start = time.clock();

    for a in bases:
        #print "Spot"

        if testA(p, a, r, d) == 1:    #if base returns prime
            i += 1              #try another base
        else:
            #print "Not prime. a =", a;
            return 0
    end = time.clock()
    elapsed = end - start;
    #print "The gcd function took", (gcd_elapsed / elapsed)*100, "% of \
#the total time."
    print "It took", elapsed, " seconds."
    return 1

def millerRabinProof( p ):
    """Determines if number "p" is prime.  P must be a 336 digit
        number or less.  This is faster than any other algorithm
        I think."""
    if p < 1373653: #7 digits.  V2
        return millerRabinBases( p, [2, 3] );
    elif p < 4759123141:    #10 digits
        return millerRabinBases( p, [2, 7, 61] )
    elif p < 41550071728321:    #14 digits.  V7
        return millerRabinBases( p, [2, 3, 5, 7, 11, 13, 17] )
    elif p < 41234316135705689041:  #20 digits.  V9: smallest strong pseudoprime
        print "V9"
        return millerRabinBases( p, sieve.sieveSmall(9))#to first 9 prime bases
    elif p < 1955097530374556503981: #22 digits.  V10
        return millerRabinBases( p, sieve.sieveSmall(10))
    elif p < 7395010240794120709381: #22 digits.  V11
        return millerRabinBases( p, sieve.sieveSmall(11))
    elif p < 318665857834031151167461: #24 digits.  V12
        return millerRabinBases( p, sieve.sieveSmall(12))
    elif p < 8038374574536394912570796143419421081388376882875581458374889175222974273765333652186502336163960045457915042023603208766569966760987284043965408232928738791850869166857328267761771029389697739470167082304286871099974399765441448453411558724506334092790222752962294149842306881685404326457534018329786111298960644845216191652872597534901:
        #337 digits.  V46
        return millerRabinBases( p, sieve.sieveSmall(46))
    else:
        return "This number is too big.  It must be 336 digits or less."

def millerRabinGRH( p ):
    """Input: An integer to be tested for primality.
        Assuming the Generalized Riemann Hypothesis is true, this algorithm
        proves an integer is prime."""
    import math

    return millerRabinBases( p, [2] + range(3, math.ceil(2*pow(math.log(p), 2))))

def millerRabinGRHsieve( p ):
    """Input: An integer to be tested for primality.
        Assuming the Generalized Riemann Hypothesis is true, this algorithm
        proves an integer is prime."""
    import math

    return millerRabinBases( p, sieve.sieve(math.ceil((2*pow(math.log(p), 2)))))

def naivePrime( p ):
    import math
    i = 3;
    while i <= math.ceil(math.sqrt(p)):
        if p % i == 0:
            return 0
        i += 2
    return 1