CoCalc Shared Filessupport / 2015-12-19-number-theory-book / Chapter 16 - Large Primes.sagews
Authors: Harald Schilly, ℏal Snyder, William A. Stein
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))
%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