Sharedsupport / 2015-12-19-number-theory-book / Chapter 16 - Large Primes.sagewsOpen in CoCalc
Authors: Harald Schilly, ℏal Snyder, William A. Stein
License: GNU General Public License v3.0
Description: Examples for support purposes.

Chapter 16 - Large Primes

Fermat numbers

We define a function to output Fermat numbers:

def fermat(n): return(2^(2^n)+1)
fermat(5)
4294967297
89694ebf-b1c4-493b-8149-9732b9b03b29i %md A table of the first few Fermat numbers:

A table of the first few Fermat numbers:

c64d24e4-d0d6-4f76-9d84-a4a0aa77745e table([(n,fermat(n),is_prime(fermat(n))) for n in [0..5]], header_row=["n", "F_n","Prime?"], frame=True)
+---+------------+--------+ | n | F_n | Prime? | +===+============+========+ | 0 | 3 | True | +---+------------+--------+ | 1 | 5 | True | +---+------------+--------+ | 2 | 17 | True | +---+------------+--------+ | 3 | 257 | True | +---+------------+--------+ | 4 | 65537 | True | +---+------------+--------+ | 5 | 4294967297 | False | +---+------------+--------+

Factorization of F5:

factor(fermat(5))
641 * 6700417

The procedure below will compute appropriate powers of 3 mod F3 to allow us to apply Peppin's test:

a=3 for n in [1..(2^3)-1]: a=mod(a^2,fermat(3)) print(a,fermat(3))
(256, 257)

Thus, Peppin's test implies F3 is prime. Now we try it for F7:

a=3 for n in [1..2^7-1]: a=mod(a^2,fermat(7)) print(a,fermat(7))
(110780954395540516579111562860048860420, 340282366920938463463374607431768211457)

In this case, Peppin's test tell us that F7 is composite.

We now define a function to implement Peppin's test:

def peppin(n): m=fermat(n) a=3 for k in [1..(2^n-1)]: a=mod(a^2,m) if mod(a,m)==m-1: return("prime") else: return("composite")
peppin(4)
'prime'
peppin(5)
'composite'

Using the function we defined, we can produce a table summarizing what goes on for the first 8 Fermat numbers.

table([(n,peppin(n)) for n in [1..8]], header_row=["n", "Result"], frame=True)
+---+-----------+ | n | Result | +===+===========+ | 1 | prime | +---+-----------+ | 2 | prime | +---+-----------+ | 3 | prime | +---+-----------+ | 4 | prime | +---+-----------+ | 5 | composite | +---+-----------+ | 6 | composite | +---+-----------+ | 7 | composite | +---+-----------+ | 8 | composite | +---+-----------+

The built-in factor command can't factor F10 in a reasonable amount of time. Try it:

factor(fermat(10)) 1dd6302b-55a7-445b-a180-4115f9ad15ddi %md However, our Peppin's test function can tell us instantly that F<sub>10</sub> is composite:

However, our Peppin's test function can tell us instantly that F10 is composite:

peppin(10)
'composite'

Mersenne numbers

We define a function to output Mersenne numbers:

def mersenne(n): return(2^n-1)

The procedure below will find the Mersenne primes with exponent less than 258:

p=2 mp=[] while p<258: if is_prime(mersenne(p)): mp.append([p,mersenne(p)]) p=next_prime(p) print(table(mp))
2 3 3 7 5 31 7 127 13 8191 17 131071 19 524287 31 2147483647 61 2305843009213693951 89 618970019642690137449562111 107 162259276829213363391578010288127 127 170141183460469231731687303715884105727
0295196c-ef16-4908-92c8-db510954a66ci %md Using Fermat's Little Theorem we can sometimes show a Mersenne number is composite.

Using Fermat's Little Theorem we can sometimes show a Mersenne number is composite.

3.powermod(mersenne(11)-1,mersenne(11))
1013

Thus M7 is composite.

3.powermod(mersenne(67)-1,mersenne(67))
95591506202441271281

Thus M67 is composite.

factor((2^67-1))
193707721 * 761838257287

Lucas-Lehmer test

r=4 print(r) for n in [1..5]: r=mod(r^2-2,127) print(r)
4 14 67 42 111 0

Thus M7 is prime.

The following function implements the Lucas-Lehmer test:

def lltest(p): m=mersenne(p) r=4 for k in [2..p-1]: r=mod(r^2-2,m) if r==0: return("prime") else: return("composite")
lltest(521)
'prime'

We now use our function to find the Mersenee primes with exponent less than 500.

p=2 mp=[] while p<5000: if lltest(p)=="prime": mp.append(p) p=next_prime(p) print(mp)
[3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423]

Prime Certificates

e88d2833-df74-4bc0-8813-7b8d6d379782i %md The command **Primes()** will produce the set of primes, and if we apply the command **unrank(n)** to it, we obtain the n-th prime.

The command Primes() will produce the set of primes, and if we apply the command unrank(n) to it, we obtain the n-th prime.

Primes().unrank(10^4-1)
104729
factor(104729-1)
2^3 * 13 * 19 * 53

We look for a witness to the primality of 104729:

[2.powermod(104728/k,104729) for k in [2,13,19,53]]
[1, 2522, 48162, 1]
[3.powermod(104728/k,104729) for k in [2,13,19,53]]
[104728, 1, 785, 78661]
[12.powermod(104728/k,104729) for k in [2,13,19,53]]
[104728, 76744, 48162, 78661]
12.powermod(104728,104729)
1

Thus 12 is a witness to the primality of 104729.

The following function will return the list of prime divisors of a number:

def prime_divisors(n): primelist=[] L=list(factor(n)) for m in [0..len(L)-1]: primelist.append(L[m][0]) return(primelist)
prime_divisors(60)
[2, 3, 5]
prime_divisors(2^107-2)
[2, 3, 107, 6361, 69431, 20394401, 28059810762433]

The following function will produce a prime certificate when a prime number is given as input.

def prime_certificate(n): pd = prime_divisors(n-1) exponents= [(n-1)/k for k in pd] a=2 while 1 in [a.powermod(k,n) for k in exponents]: a=a+1 if a.powermod(n-1,n)==1: return([n,a,pd]) else: return("composite")
prime_certificate(12)
'composite'
prime_certificate(104729)
[104729, 12, [2, 13, 19, 53]]
prime_certificate(mersenne(107))
[162259276829213363391578010288127, 3, [2, 3, 107, 6361, 69431, 20394401, 28059810762433]]
prime_certificate(20394401)
[20394401, 3, [2, 5, 13, 37, 53]]
prime_certificate(28059810762433)
[28059810762433, 5, [2, 3, 41, 53, 67254877]]
prime_certificate(mersenne(521))
[6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151, 3, [2, 3, 5, 11, 17, 31, 41, 53, 131, 157, 521, 1613, 2731, 8191, 42641, 51481, 61681, 409891, 858001, 5746001, 7623851, 34110701, 308761441, 2400573761, 65427463921, 108140989558681, 145295143558111, 173308343918874810521923841]]

Example 16.11 - Finding a large prime

p=Primes().unrank(int(10^4*random())) p
47947
k=int(10^6*random()) k
646766
pprime=47947*646766+1 pprime
31010489403
2.powermod(pprime-1,pprime)
4
pprime=706020*47947+1 pprime
33851540941
6.powermod(pprime-1,pprime)
1
prime_certificate(pprime)
[33851540941, 6, [2, 3, 5, 7, 41, 47947]]
pprime=509248*33851540941+1 pprime
17238829521122369
3.powermod(pprime-1,pprime)
1
prime_certificate(pprime)
[17238829521122369, 3, [2, 73, 109, 33851540941]]

Procedure for finding a large prime

We begin by modifying out prime_certificate command to make it into a prime testing procedure. The input is a number n to be testes and prime divisor p of n-1.

def prime_certificate_test(n,p): pd = prime_divisors((n-1)/p) pd.append(p) exponents= [(n-1)/k for k in pd] a=2 while a.powermod(n-1,n) ==1: if 1 in [a.powermod(k,n) for k in exponents]: a=a+1 else: return("prime") return("composite")
prime_certificate_test(41,5)
'prime'
prime_certificate_test(46,5)
'composite'

We use this function to test the primality of numbers of the form kp+1 for random values of k.

k=int(10^6*random()) pprime=k*p+1 while prime_certificate_test(pprime,p)=="composite": k=int(10^6*random()) pprime=k*p+1 p=pprime (p,k)
(33851540941, 706020)
47947*706020+1
33851540941
prime_certificate(33851540941)
[33851540941, 6, [2, 3, 5, 7, 41, 47947]]
k=int(10^6*random()) pprime=k*p+1 while prime_certificate_test(pprime,p)=="composite": k=int(10^6*random()) pprime=k*p+1 p=pprime (p,k)
(17238829521122369, 509248)
509248 * 33851540941+1
17238829521122369
prime_certificate(17238829521122369)
[17238829521122369, 3, [2, 73, 109, 33851540941]]

We now build a function which automates the steps above to produce a prime with more than d digits.

def large_prime(d): p=Primes().unrank(int(10^4*random())) while p<10^d: pprime=int(10^6*random())*p+1 while prime_certificate_test(pprime,p)=="composite": pprime=int(10^6*random())*p+1 p=pprime return(p)
large_prime(10)
3124353799825117
large_prime(100)
101146885292260685134952367990338838962219308894528819916336354941587210297210635989873113194068712643749