''' 
    Fermat Primality test
    Exercise 13 (i)
    
    Melina Zikou (2021)
'''
import random
import math
import sys
import time

# b = base number
# e = exponential
# m = modulo 
#     
# calculate b^e mod m
def fast(b,e,m):
    x = b
    g = e
    d = 1

    while g > 0:
        if g % 2 == 0:
            x = (x * x) % m
            g = g // 2
        else:
            d = (x * d) % m
            g = g - 1
    return d

# Fermats Primality Test     
# a ^ (n - 1) == 1 -> prime
# a ^ (n - 1) != 1 -> composite
def fermat(n):
    k = 0
    while k < 128:
        a = random.randint(2, n - 2)

        if fast(a, n - 1, n) != 1:
            return False
        k += 2
    return True

# decrease the number of potential primes to be tested
def isPrime(n):
    # lowPrimes => all primes under 1000
    
    lowPrimes =   [3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97
                ,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179
                ,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269
                ,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367
                ,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461
                ,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571
                ,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661
                ,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773
                ,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883
                ,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997]

    # numbers greater than 3
    if (n >= 3):
        # n & 1 = 0 +> n is composite
        if (n & 1 != 0):
            for p in lowPrimes:
                if (n == p):
                    return True
                if (n % p == 0):
                    return False
            return fermat(n)
        
    return False

def generateLargePrime(k):
     # k is the desired bit length
     r = 100 * (math.log(k,2) + 1) # max number of attempts

     while r > 0:
         n = random.randrange(pow(2, k - 1), pow(2,k))
         r -= 1
         if isPrime(n) == True:
            return n
     return r

# fermat prime
start = time.time()
print(generateLargePrime(2048))
end = time.time()
print("\nFERMAT - 2048 bits:", end - start)
