Chapter 16 - Large Primes.sagews

Chapter 16 - Large Primes

Fermat numbers

We define a function to output Fermat numbers:

def fermat(n):
    return(2^(2^n)+1)

fermat(3)

257

︠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 F₅:

factor(fermat(5))

641 * 6700417

The procedure below will compute appropriate powers of 3 mod F₃ 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 F₃ is prime. Now we try it for F₇:

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 F₇ 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 F₁₀ 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 F₁₀ 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 M₇ is composite.

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

95591506202441271281

Thus M₆₇ 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 M₇ 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