PrimalityTests

32 days ago by bcr49

Trial Divison:

S = [] for a in range(10^5): if is_prime(a): S.append(a) def SmallDivisors(n): divs = [] for p in S: if n%p ==0: divs.append(p) return divs print len(S) print float(10^5/log(10^5)) 
       

Fermat Test:

def FermatTest(n,a): R = Integers(n) a = R(a) return a^(n-1)==1 
       

This is the supposed "Socat prime"

Socatn = 143319364394905942617148968085785991039146683740268996579566827015580969124702493833109074343879894586653465192222251909074832038151585448034731101690454685781999248641772509287801359980318348021809541131200479989220793925941518568143721972993251823166164933334796625008174851430377966394594186901123322297453 print "Size of n is", ceil(log(Socatn,2)) SD = SmallDivisors(Socatn) print "Fermat test of n:", FermatTest(Socatn,2) print "Small prime factors of n:", SD N = Socatn/SD[0]/SD[1] print "Fermat test of N = n/smallprimefactors:", FermatTest(N,2) print "Size of N is", ceil(log(N,2)) #We can't factor further, but someone out there might be able to... 
       
def MRwitness(n,a): s = (n-1).valuation(2) d = (n-1)/2^s A = Integers(n)(a) ad = A^d #print ad if ad == 1: #print "a =",A,"is not a MR witness for n =",n, "(case i)" return false else: for r in range(s): if ad == -1: #print "a =",a,"is not a MR witness for n =",n, "(case ii)" return false ad = ad^2 #print ad return true def MRtest(n,k): #Test n with max k iterations of MR iterations = 0 while iterations < k: a = randrange(n-2) + 2 # randrange(N) selects from [0,N-1] inclusive # we want from [2,n-1] witnessfound = MRwitness(n,a) if witnessfound: #print "a =",a,"is an MR witness for n =",n return false iterations += 1 return true 
       

How Many MR witnesses does a given integer have?

n = randrange(2^10) if n%2 == 0: n += 1 Witnesses = [] for a in range(n-2): if MRwitness(n,a+1): Witnesses.append(a) m = len(Witnesses) print "n =",n print "number of witnesses =", m print "proportion of witnesses =", float(m)/float(n) is_prime(n) 
       

Generating a large prime:

Here we select a random odd integer of size k and use MR test to check for primality. In practice it is good to look for small prime factors using trial division. This implementation does not do this.

def randomprime(k): #k is the number of bits for our prime attempts = 0 foundprime = false while not foundprime: attempts += 1 n = randrange(2^(k-2))+1 + 2^(k-1) #a random integer of size k if n.valuation(2) > 0: n += 1 #make sure n is odd foundprime = MRtest(n,10) print "After", attempts, "attempts we found a probable prime: p =",n return n k = 1000 time n = randomprime(k) #time print is_prime(n) print ceil(float(log(n,2))) 
       

The prime number theorem can be used to estimate the probability that a random k-bit integer is prime as:

prop = 2*float(((2^k)/(k*log(2)) - (2^(k-1))/((k-1)*log(2)))/2^(k-1)) #the 2 in front is because we only consider odd integers print "probability that a random k-bit integer is prime =", prop print "that is 1 in", ceil(1/prop) print "probability of finding a prime after 100 attempts is", 1-(1-prop)^100 print "probability of finding a prime after 200 attempts is", 1-(1-prop)^200 print "probability of finding a prime after 300 attempts is", 1-(1-prop)^300 print "probability of finding a prime after 400 attempts is", 1-(1-prop)^400 print "probability of finding a prime after 500 attempts is", 1-(1-prop)^500 
       

Fermat Factoring:

def FermatFactoring(n,k): #stop after k attempts for y in range(k): if is_square(n+y^2): #print y x = sqrt(n+y^2) #print x^2 return [(x+y),(x-y)] return 0 print FermatFactoring(315,10) 21*15 == 315 
       

Fermat factoring works well when n = ab is a product with a ~ b. We can generate some such examples

a = randomprime(1024) b = next_prime(a) print "n =", n n = a*b print "MRtest of n:", MRtest(a*b,10) factors = FermatFactoring(a*b,1000) print "Found a factor", factors[0] print factors[0] * factors[1] == n print "difference between the two factors is", a-b ## this works because the difference of the factors is small 
       

Unfortunately, the difference between 2 random k-bit primes is usually around O(2^k)

a = randomprime(1000) b = randomprime(1000) print "Difference between the two factors is", a-b print "This is an integer of size", ceil(log(abs(a-b),2)) 
       

 

A MR witness that is not a Fermat witness allows us to factor n. These are however, are usually not very common. An exception to this is Carmichael Numbers.

n = randrange(2^15) if n%2 == 0: n += 1 n = 561 FactoringWitnesses = [] MRWitnesses = [] for b in range(n-1): a = b+1 if GCD(a,n) == 1: if MRwitness(n,a): MRWitnesses.append(a) if FermatTest(n,a): FactoringWitnesses.append(a) m1 = len(MRWitnesses) m2 = len(FactoringWitnesses) print "n =",n print "number of MR witnesses =", m1 print "number of factoring witnesses =", m2 print "proportion of factoring witnesses =", float(m2)/float(n) is_prime(n) print FactoringWitnesses prime_factors(n) 
       

Fermat Factoring when a multiple of phi(n) is known:

def Fermat2(n,m,a): s = valuation(m,2) d = m/2^s a = Integers(n)(a) ad = a^d adr = ad print "r =", 0, "has a^(d*2^r) =", adr if ad == 1: return 1 if ad == -1: return 1 for r in range(s-1): adr = ad^2 print "r =", r+1, "has a^(d*2^r) =", adr if adr == 1: return ad ad = adr p = randomprime(1000) q = randomprime(1000) #p = 89 #q = 97 m = (p-1)*(q-1) n = p*q print "n = ", n print "m = ", m s = valuation(m,2) d = m/2^s print "d = ", d print "s = ", s #a = 3 #a = 24 a = randrange(n-2)+2 print "a =", a x = Fermat2(n,m,a) print "x-1 =", x-1 factor = GCD(x-1,n) print "yields the factor t =", factor print "Did we succeed?", p == factor or q == factor 
       
print a^(33*2^2) % 89 print a^(33*2^2) % 97 
       

Largest known prime:

p = 2^77232917 -1 MRwitness(p,2)