import mrpt
import math
import random
def prime_census(threshold):
#The prime number theorem gives the probability that a number less
#than threshold is prime.
prob = 1/math.log(threshold)
#since we sample one hundred thousand values, the expected number of primes is
c2=100000*prob
c1=0
for j in range(100000):
candidate=random.randint(1,threshold)
if mrpt.is_probable_prime(candidate):
c1+=1
return (c1,c2)
Here is a trial run, with threshold taking on an an assortment of large values.
for j in [10,20,50,100]:
for k in range(3):
print j,' digits:',prime_census(10**j)
The approximation is pretty accurate. We should not expect exact accuracy: On one hand, there is sampling error. On the other hand, the prime number theorem is an asymptotic statement: It says that the ratio of N/ln(N) to the number of primes less than N approaches 1 as a limit as N approaches infinity. For a specific value of N this ratio will be different from 1.
The naive Fermat test will never give a false negative: If Fermat says the number is composite, then it is definitely composite. The code below counts the number of false positives for 10000 randomly sampled values less than the threshold.
def fermat_vs_mrpt(threshold):
count=0
for j in range(10000):
candidate=random.randint(1,threshold)
witness=random.randint(1,candidate)
fermat_prime=(pow(witness,candidate-1,candidate)==1)
probable_prime=mrpt.is_probable_prime(candidate)
if fermat_prime and not probable_prime:
count+=1
return count
Here is a test run with different threshold values.
for t in [10**6,10**10,10**20,10**50,10**100]:
print fermat_vs_mrpt(t)
This result never ceases to astound me. For large values--in particular for the large values encountered in cryptography---the simple test, which computes the p-1 power mod p of a single randomly chosen a less than p is in practice one hundred percent reliable!
Since the answer is different for every student, to demonstrate the solution, we'll do something a little different and crack one of the posted ciphertexts without the secret information. Let's work with the pair (1039954589329399 ,120651427100177). Since the modulus is about 1015, there will be a prime factor less the square root of this, which is less than 4 X 107, so we can do a brute-force factorization:
N=1039954589329399
p=3
while N%p !=0:
p+=2
print p
q=N/p
print q
Given the factors p and q, we generate the smallest possible encryption exponent by finding the smallest odd prime that does not divide (p-1)(q-1), and then derive the encryption exponent.
import mrpt
m=(p-1)*(q-1)
e=3
while m%e==0:
e+=2
print e
d=mrpt.modinverse(e,m)
print d
Now we'll use this information to decrypt the ciphertext and convert it to ASCII
import bytestuff
ciphertext=120651427100177
plaintext=pow(ciphertext,d,N)
print bytestuff.bytes_to_ascii(bytestuff.long_to_bytes(plaintext))