Large Primes<Text-field layout="_pstyle4" style="_pstyle4">Prime listing, primality testing, and prime factorizatrion</Text-field><Text-field layout="_pstyle5" style="_pstyle5">Prime listing</Text-field>primelist:=proc(n)
local candidatelist, i, propermultiple;
candidatelist:=array([seq(i,i=1..n)]);
candidatelist[1]:=0;
for i from 2 to floor(sqrt(n)) do
if candidatelist[i] <> 0 then
for propermultiple from 2*i to n by i do
candidatelist[propermultiple]:=0
end do;
end if;
end do;
{seq(candidatelist[i],i=1..n)} minus {0};
end proc:primelist(100);
map(x->[nops(primelist(x))],[seq(10^k,k=1..6)]);newprimelist:=proc(n)
local m, x,j,i;
m:=floor((n-1)/2);
x:=array([seq(1,i=1..m)]);
for j from 1 to floor(evalf((sqrt(n)-1)/2)) do
if x[j]=1 then
for i from 2*j*(j+1) to m by 2*j+1 do
x[i]:=0;
end do;
end if;
end do;
{seq(x[i]*(2*i+1),i=1..m)} minus {0} union {2};
end proc:newprimelist(100);ourlist:=convert(newprimelist(10^7),`list`):ourlist[500000];<Text-field layout="_pstyle5" style="_pstyle5">Primality testing</Text-field>primetest:=proc(n)
local upperlimit,trialdivisor;
upperlimit:=floor(evalf(sqrt(n)));
for trialdivisor from 2 to upperlimit do
if modp(n,trialdivisor)=0 then
return(False);
end if;
end do;
return(True);
end proc: starttime:=time():
primetest(2^31-1);time()-starttime;primetest(ithprime(500000)*ithprime(500001));time()-starttime;off;newprimetest:=proc(n)
local upperlimit, i;
upperlimit:=floor(evalf(sqrt(n)));
for i from 1 while ourlist[i]<=upperlimit do
if modp(n,ourlist[i])= 0 then
return(False);
end if;
end do;
return(True);
end proc:starttime:=time():newprimetest(2^31-1);
time()-starttime;newprimetest(ithprime(500000)*ithprime(500001));
time()-starttime;Prime factorizationprimefactorization:=proc(n)
local primefactors, nstar, upperlimit, trialdivisor;
primefactors:=[];
nstar:=n;
upperlimit:=floor(evalf(sqrt(nstar)));
trialdivisor:=2;
while trialdivisor <= upperlimit do
if modp(nstar,trialdivisor)=0
then
primefactors:=[op(primefactors),trialdivisor];
nstar:=nstar/trialdivisor;
upperlimit:=floor(evalf(sqrt(nstar)));
else
trialdivisor:=trialdivisor+1;
end if;
end do;
primefactors:=[op(primefactors),nstar];
{seq([x,ListTools[Occurrences](x,primefactors)],x in primefactors)};
end proc:starttime:=time():primefactorization(2973087990183122315334);
time()-starttime;newprimefactorization:=proc(n)
local primefactors, nstar, upperlimit, i, trialdivisor;
primefactors:=[];
nstar:=n;
upperlimit:=floor(evalf(sqrt(n)));
i:=1;
trialdivisor:=ourlist[i];
while trialdivisor <= upperlimit do
if modp(nstar,trialdivisor)=0
then
primefactors:=[op(primefactors),trialdivisor];
nstar:=nstar/trialdivisor;
upperlimit:=floor(evalf(sqrt(nstar)));
else
i:=i+1;
trialdivisor:=ourlist[i];
end if;
end do;
primefactors:=[op(primefactors),nstar];
{seq([x,ListTools[Occurrences](x,primefactors)],x in primefactors)};
end proc:starttime:=time():newprimefactorization(2973087990183122315334);
time()-starttime;newprimefactorization(9199838343019266106321698);off;ifactor(255678904536136571239);<Text-field layout="_pstyle4" style="_pstyle4">Fermat numbers</Text-field>fermatnumber:=n->2^(2^n)+1:seq(fermatnumber(n),n=0..4);seq(primetest(fermatnumber(n)),n=0..4);fermatnumber(5);primetest(%);starttime:=time():newprimefactorization(fermatnumber(5));time()-starttime;fermatnumber(6);starttime:=time():
newprimefactorization(fermatnumber(6));
time()-starttime;startime:=time():
f:=fermatnumber(10):
trialdivisor:=4097:
while modp(f,trialdivisor) <> 0 do
trialdivisor:=trialdivisor+4096:
end do:
trialdivisor;
time() - starttime;fprime:=f/trialdivisor:
while modp(fprime,trialdivisor) <> 0 do
trialdivisor:=trialdivisor+4096:
end do:
trialdivisor;
time() - starttime;st:=time():
f:=fermatnumber(7):
trialdivisor:=513:
while modp(f,trialdivisor) <> 0 do
trialdivisor:=trialdivisor+512:
end do:
trialdivisor;
time() -st;pepintest:= proc(n)
local m, a;
m:=fermatnumber(n);
a:=3;
to 2^n-1 do
a:=modp(a^2,m);
end do;
if modp(a,m) = m-1
then
return("prime")
else
return("composite");
end if;
end proc:starttime:=time():pepintest(7);time()-starttime;<Text-field layout="_pstyle4" style="_pstyle4">Mersenne numbers</Text-field>mersennenumber:= n->2^n-1:seq(mersennenumber(n),n=1..10);select(n->isprime(mersennenumber(n)),select(isprime,[seq(i,i=1..257)]));ifactor(mersennenumber(67));pseudoprimetest:=proc(n)
if modp(3&^(n-1),n)=1
then return("test fails")
else return("composite")
end if;
end proc:pseudoprimetest(mersennenumber(23));pseudoprimetest(mersennenumber(67));lucaslehmertest:=proc(p)
local m, r;
m:=mersennenumber(p);
r:=4;
to p-2 do
r:=modp(r^2-2,m)
end do;
if r=0
then return("prime")
else return("composite")
end if;
end proc:lucaslehmertest(521);select(n->evalb(lucaslehmertest(n)="prime"),select(isprime,[seq(i,i=1..5000)]));
st:=time():
primefactorization(mersennenumber(67));
time()-st;<Text-field layout="_pstyle4" style="_pstyle4">Prime Certificates</Text-field>primecertificate:=proc(n)
local primedivisors, exponents, a;
primedivisors:={seq(ifactors(n-1)[2][i][1],i=1..nops(ifactors(n-1)[2]))};
exponents:=map(x->(n-1)/x,primedivisors);
a:=2;
while member(1,map(x->modp(a&^x,n),exponents)) do
a:=a+1
end do;
return[n,a,primedivisors];
end proc:n:=ithprime(10^4);primecertificate(n);newprimecertificate:=proc(n)
local primedivisors, exponents, a;
if n<10^7
then return(n)
else
primedivisors:={seq(ifactors(n-1)[2][i][1],i=1..nops(ifactors(n-1)[2]))};
exponents:=map(x->(n-1)/x,primedivisors);
a:=2;
while member(1,map(x->modp(a&^x,n),exponents)) do
a:=a+1
end do;
return([n,a,map(newprimecertificate,primedivisors)])
end if;
end proc:
newprimecertificate(mersennenumber(521));<Text-field layout="_pstyle4" style="_pstyle4">Finding large primes</Text-field>ournextprime:=proc(n)
local p;
p:=n;
while primetest(p)=False
do
p:=p+1
end do:
return(p);
end proc:
ournextprime(24);newpseudoprimetest:=proc(n)
if modp(2&^(n-1),n)=1 and
modp(3&^(n-1),n)=1 and
modp(5&^(n-1),n)=1
then
return("test fails")
else
return("composite")
end if;
end proc:nextprimecandidate:=proc(n)
local p;
p:=n;
while newpseudoprimetest(p)="composite" do
p:=p+1
end do;
return(p);
end proc:nextprimecandidate(10^100);starttime:=time():
ournextprime(10^12);
time()-starttime;starttime:=time():
nextprimecandidate(10^12);
time()-starttime;primecertificatetest:=proc(n,p)
local primedivisors, exponents, a;
primedivisors:={seq(ifactors((n-1)/p)[2][i][1],i=1..nops(ifactors((n-1)/p)[2]))} union {p};
exponents:=map(x->(n-1)/x,primedivisors);
a:=2;
while modp(a&^(n-1),n)=1 do
if member(1,map(x->modp(a&^x,n),exponents))
then a:=a+1
else return("prime")
end if;
end do;
return("composite");
end proc:primecertificatetest(36,5);primecertificatetest(41,5);largeprime:=proc(d)
local p, pprime,i;
i:=1;
p:=ithprime(RandomTools[Generate](integer(range=1..10^4)));
while p<10^d do
pprime:=RandomTools[Generate](integer(range=10^6..10^7))*p+1;
while primecertificatetest(pprime,p)="composite" do
pprime:=RandomTools[Generate](integer(range=10^6..10^7))*p+1;
end do;
p:=pprime;
end do;
return(p);
end proc:starttime:=time():
largeprime(100);
time()-starttime;