Sz\303\241m\303\255t\303\263g\303\251pes sz\303\241melm\303\251letJ\303\241rai AntalEzek a programok csak szeml\303\251ltet\303\251sre szolg\303\241lnak1. A pr\303\255mek eloszl\303\241sa, szit\303\241l\303\241s2. Egyszer\305\261 faktoriz\303\241l\303\241si m\303\263dszerek3. Egyszer\305\261 pr\303\255mtesztel\303\251si m\303\263dszerekLUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0Yn4. Lucas-sorozatok5. Alkalmaz\303\241sok restart; with(numtheory);5.1. Fermat-sz\303\241mok.5.2. Feladat.#
# This procedure try to find factors i*2^(n+2)+1, 1<=i<=k
# of the Fermat number 2^(2^n)+1.
#
fermatfactors:=proc(n::posint,k::posint) local i,f;
f:=[];
for i from 1+2^(n+2) by 2^(n+2) to k*2^(n+2)+1 do
if 2&^(2^n)+1 mod i=0 then f:=[op(f),i] fi
od; f end;fermatfactors(5,10); fermatfactors(5,100000);5.3. Feladat.interface(verboseproc=2);print(fermat);5.4. Mersenne-sz\303\241mok.mersennes:=[2,3,5,7,13,17,19,31,61,89,107,127,521,607,1279,2203,2281,3217,4253,4423,9689,9941,11213,19937,21701,23209,44497,86243,110503,132049,216091,756839,859833,1257787,1398269,2976221,3021377,6972593,13466917,20996011,24036583,25964951,30402457,32582657,37156667,42643801,43112609];map(p->p mod 4, mersennes);logm:=evalf([[i,log[2.](mersennes[i])]$i=1..nops(mersennes)]);with(plots); plotm:=plot(logm):
plotline:=plot(evalf(exp(-gamma))*x,x=1..nops(mersennes)):
pointm:=plot(logm,style=POINT):
display({plotline,plotm,pointm});5.5. Feladat.#
# This procedure do a pretest for the Mersenne numbers
# M[n]=2^n-1. The exponents are taken from
# the sequence L. The result is the sequence of sequences,
# which contains all exponents n which is prime and the
# prime divisors having the form 2*k*n+1 for k until K, where
# k congruent to 0 or -n mod 4.
#
mersennepretest:=proc(L::list(posint),K::posint)
local nL,D,M,n,k,d1,d2; nL:=[];
for n in L do
if isprime(n) then
D:=[n]; d1:=modp(n,4); d2:=modp(-n,4);
for k from d2 while k<=K do
M:=2*k*n+1;
if isprime(M) then
if pmod(2&^n-1,M)=0 then D:=[op(D),M]; fi;
fi;
k:=k+d1; M:=2*k*n+1;
if isprime(M) then
if 2&^n-1 mod M=0 then D:=[op(D),M] fi;
fi;
k:=k+d2;
od;
nL:=[op(nL),D];
fi;
od; nL end;L:=[i$i=20..1000];LL:=mersennepretest(L,10000);nops(LL); select(x->evalb(nops(x)=1),LL); nops(%);5.6. Feladat.interface(verboseproc=2);print(mersenne);5.7. h*2^m+-1 alak\303\272 pr\303\255mek keres\303\251se.#
# A pretest will be done for (h0+c*h)*b^n+d, where h runs
# from 0 to H-1 and d in D. All prime divisors from P0 to
# P1 will be tried. P0 must be odd, and b,c must have prime
# factors smaller then P0. The values h, for which
# (h0+c*h)*b^n+d not divisible with primes from P0 to P1
# for all d in D are given back.
#
pretest:=proc(h0::nonnegint,H::posint,c::posint,b::posint,n::posint,
D::set(integer),P0::posint,P1::posint)
local T,h,p,d,hh,L,nn,bb,cc,bT,cT,i,j;
T:=array(0..H-1);
for h from 0 to H-1 do T[h]:=1 od;
bT:=array(0..b-1); # table of inverses for b
for i from 0 to b-1 do bT[i]:=0 od;
for i from 0 to b-1 do
for j from 0 to b-1 do
if irem(i*j+1,b)=0 then bT[i]:=j fi
od
od;
cT:=array(0..c-1); # table of inverses for c
for i from 0 to c-1 do cT[i]:=0 od;
for i from 0 to c-1 do
for j from 0 to c-1 do
if irem(i*j+1,c)=0 then cT[i]:=j fi
od
od;
for p from P0 to P1 by 2 do
if isprime(p) then
nn:=modp(n,p-1);
bb:=(bT[modp(p,b)]*p+1)/b;
bb:=modp(bb&^nn,p);
cc:=(cT[modp(p,c)]*p+1)/c;
for d in D do
hh:=modp((-d*bb-h0)*cc,p);
for h from hh to H-1 by p do T[h]:=0; od
od
fi
od;
L:=[];
for h from 0 to H-1 do
if T[h]=1 then L:=[op(L),h] fi
od; L end;#
# This procedure do the Proth test for the number n=h*2^m+1.
#
proth:=proc(h::posint,m::posint) local a,b,n; n:=h*2^m+1;
if h>=2^m then error "first parmeter is too large",h fi;
if not type(h,odd) then error "first parameter must be odd",h fi;
if m=1 or m=2 then return true fi;
if m=3 then if h=5 then return true else return false fi fi;
# 2 is not a quadratic residue mod n, we start with 3
# (2*a|n)=(a|n), hence enough to try odd bases a
for a from 3 by 2 do
b:=2&^m mod a; b:=h*b+1 mod a;
# by quadratic reciprocity, (a|n)=(n|a)=(b|a)
b:=modp(2&^m,a); b:=modp(h*b+1,a);
if jacobi(b,a)=0 then return false fi;
if jacobi(b,a)=-1 then
if modp(a&^((n-1)/2),n)=n-1 then return true else return false fi;
fi;
od end;#
# After a pretest the Proth test is used to find primes
# (h0+c*h)*2^n+1, where h runs from 0 to H-1. All prime
# divisors from interval P will be tried during the pretest.
#
prothprimes:=proc(h0::posint,c::posint,H::posint,n::posint,P::range(posint))
local h,L,LL;
L:=pretest(h0,H,c,2,n,{1},op(1,P),op(2,P)); LL:=[];
for h in L do
if proth(h0+c*h,n) then LL:=[op(LL),h] fi
od; LL end;prothprimes(3,30,3000,2000,7..1000);5.8. Feladat.5.9. Feladat.#
# This procedure calculate the Lucas sequence terms
# [V[n],V[n+1]]. If the parameter N also given
# then the calculations are done mod N. Here V[n]=a^n+b^n,
# where a,b are the roots of x^2-Px+Q=0.
#
lucasV:=proc(n,P,Q,N) local L,B,i,Qk;
B:=convert(n,base,2);
Qk:=1;
if nargs=3 then L:=[2,P];
for i from nops(B) to 1 by -1 do
if B[i]=0 then
L:=[L[1]^2-2*Qk,L[2]*L[1]-P*Qk]; Qk:=Qk^2;
else
L:=[L[2]*L[1]-P*Qk,L[2]^2-2*Q*Qk]; Qk:=Qk^2*Q;
fi;
od;
else L:=[2 mod N,P mod N];
for i from nops(B) to 1 by -1 do
if B[i]=0 then
L:=[L[1]&^2-2*Qk mod N,L[2]*L[1]-P*Qk mod N]; Qk:=Qk&^2 mod N;
else
L:=[L[2]*L[1]-P*Qk mod N,L[2]&^2-2*Q*Qk mod N]; Qk:=Qk&^2*Q mod N;
fi;
od;
fi; L end;#
# This procedure calculate a^h+a^(-h) for a unit a=r+s*sqrt(D)
# with norm 1. Here r,s rational numbers, D must be a square
# free integer. All computations are done mod N.
#
rieselsetup:=proc(r,s,D,h,N) local P,a;
a:=r+s*sqrt(D);
P:=r+s*sqrt(D)+(r-s*sqrt(D))/(r^2-s^2*D);
if not type(P,integer) then ERROR(a+1/a, `not an integer`) fi;
lucasV(h,P,1,N)[1]; end;#
# This procedure use iteration v<-v^2-2 starting with v=a^h+a^-h
# where the quadratic unit a=r+s*sqrt(D) is given for the
# primality testing of h*2^n-1. Here n>1 and h<2^n is odd.
#
riesel:=proc(h,n,r,s,D) local N,v,i;
if h>=2^n then error "first parameter is too large",h fi;
if type(h,even) then error "first parameter have to be odd",h fi;
if n<2 then error "second parameter is too small",h fi;
N:=h*2^n-1;
v:=rieselsetup(r,s,D,h,N);
for i to n-2 do v:=v^2-2 mod N od;
evalb(v=0); end;#
# After a pretest the Riesel test with D=5 and unit
# (1+sqrt(5))^2/4 is used to find primes (h0+30*h)*2^n-1,
# where h runs from 0 to H-1. Only the pairs [3,0],[9,0],
# [23,0],[29,0],[7,1],[9,1],[19,1],[27,1],[11,2],[17,2],
# [21,2],[27,2],[1,3],[3,3],[13,3],[21,3] for
# [h0 mod 30,n mod 4] gives combinations for which the
# test may used and N is not divisible by 2,3,5, hence
# this procedure should be used only with these pairs.
# All prime divisors until P will be tried during the
# pretest.
#
rieselsqrt5:=proc(h0,H,n,P) local h,L,LL;
L:=pretest(h0,H,30,2,n,{-1},7,P); LL:=[];
for h in L do
if riesel(h0+30*h,n,3/2,1/2,5) then LL:=[op(LL),h] fi
od; LL end;rieselsqrt5(3,3000,2000,1000);5.10. Ikerpr\303\255mek.#
# This filter left only those h's from the list L in the
# resulting list for which (h0+c*h)*2^n+d is base-2 pseudoprime.
#
prob:=proc(h0::nonnegint,c::posint,n::posint,d::integer,L::list(nonnegint)) local LL,h,N; LL:=[];
for h in L do
N:=(h0+c*h)*2^n+d;
if modp(3&^(N-1),N)=1 then LL:=[op(LL),h] fi
od; LL end;#
# This filter left only those h's from the list L in the
# resulting list for which (h0+c*h)*2^n+1 is prime. It use
# the Proth test.
#
exactplus:=proc(h0,c,n,L) local LL,h,a,S;
LL:=[];
for h in L do
if proth(h0+c*h,n) then LL:=[op(LL),h] fi
od; LL end;5.11. Feladat.#
# This filter left only those h's from the list L in the
# resulting list for which (h0+30*h)*2^n-1 is prime.
#
exactminus:=proc(h0,n,L) local LL,h;
LL:=[];
for h in L do
if riesel(h0+30*h,n,3/2,1/2,5) then LL:=[op(LL),h] fi
od; LL end;5.12. Feladat.#
# After a pretest the probabilistic test for N-1,
# the Riesel test with D=5 and unit (1+sqrt(5))^2/4 for N-1,
# the probabilistic test for N+1, and the Proth test for N+1
# is used to find twin primes N+-1 with N=(h0+30*h)*2^n,
# where h runs from 0 to H-1. Only the pairs
# [3,0],[9,1],[21,3],[27,2], for [h0 mod 30,n mod 4] gives
# combinations for which the Riesel test may used and
# N+-1 is not divisible by 2,3,5, hence this procedure
# should be used only with these pairs. All prime divisors
# until P will be tried during the pretest.
#
twin:=proc(h0,H,n,P) local LL;
LL:=pretest(h0,H,30,2,n,{1,-1},7,P);
print(nops(LL));
LL:=prob(h0,30,n,-1,LL);
print(nops(LL));
LL:=exactminus(h0,n,LL);
print(nops(LL));
LL:=prob(h0,30,n,1,LL);
print(nops(LL));
LL:=exactplus(h0,30,n,LL)
end;twin(3,20000,400,1000);#
# This is a test for measurement densities and times.
# First a pretests is done for N+-1 where N=(h0+30*h)*2^n
# and h runs from 0 to H-1. The pretest is repeated for
# prime limit P=2^e where e is integer from the interval
# eint. The time and the remaining number after the
# pretest is reported, and p*ln(P)^2 also calculated,
# where p is the relative frequency of numbers
# passed the pretest. The result of the last pretest
# became the input of a probabilistic test for N-1 and
# the result of this will be the input of the Riesel test
# with D=5 and unit (1+sqrt(5))^2/4 for N-1. The time of
# the probabilistic test and the result of booth tests is
# reported. The probabilistic test for N+1 and the
# Proth test for N+1 also done. Only the pairs
# [3,0],[9,1],[21,3],[27,2] for [h0 mod 30,n mod 4] gives
# combinations for which the Riesel test
# may used and N+-1 is not divisible
# by 2,3,5, hence this procedure
# should be used only with these pairs.
#
twintest:=proc(h0,H,n,eint) local e,p,P,L,t,N;
print(`Twintest; h0=`,h0,`H=`,H,`n=`,n);
N:=(h0+15*H)*2^n;
print(`Pretest for N+-1`);
for e from op(1,eint) to op(2,eint) do
P:=2^e;
t:=time();
L:=pretest(h0,H,30,2,n,{1,-1},7,P);
t:=time()-t;
p:=evalf(nops(L)/H);
print(`P=`,2^e,`p=`,p,`p*ln(P)^2=`,evalf((ln(P))^2*p),`time=`,t)
od;
print(`Test for N+-1`);
t:=time();
L:=prob(h0,30,n,-1,L);
L:=exactminus(h0,n,L);
L:=prob(h0,30,n,1,L);
L:=exactplus(h0,30,n,L);
t:=time()-t;
p:=evalf(nops(L)/H);
print(`p=`,p,`p*ln(N)^2=`,evalf((ln(N))^2*p),`time=`,t);
L
end;twintest(3,20000,400,3..15);5.13. Sophie Germain pr\303\255mek.#
# A pretest will be done for Sophie Germain prime candidates,
# to be more exact, for (h0+c*h)*b^n-1 and for (h0+c*h)*b^(n+1)-1,
# where h runs from 0 to H-1. All prime divisors from P0 to P1
# will be tried. P0 must be odd, and b,c must have prime factors
# smaller then P0. The values h0+c*h, for which (h0+c*h)*b^n-1 and
# (h0+c*h)*b^(n+1)-1 are not divisible with primes from P0 to P1
# are given back.
#
presophie:=proc(h0::nonnegint,H::posint,c::posint,b::posint,
n::posint,P0::posint,P1::posint)
local T,h,p,hh,L,nn,bb,bbb,cc,bT,cT,i,j;
T:=array(0..H-1);
for h from 0 to H-1 do T[h]:=1 od;
bT:=array(0..b-1); # table of inverses for b
for i from 0 to b-1 do bT[i]:=0 od;
for i from 0 to b-1 do
for j from 0 to b-1 do if irem(i*j+1,b)=0 then bT[i]:=j fi od
od;
cT:=array(0..c-1); # table of inverses for c
for i from 0 to c-1 do cT[i]:=0 od;
for i from 0 to c-1 do
for j from 0 to c-1 do if irem(i*j+1,c)=0 then cT[i]:=j fi od
od;
for p from P0 to P1 by 2 do
if isprime(p) then
cc:=(cT[p mod c]*p+1)/c;
nn:=modp(n,p-1);
bb:=(bT[modp(p,b)]*p+1)/b;
bbb:=modp(bb&^nn,p);
hh:=modp((bbb-h0)*cc,p);
for h from hh to H-1 by p do T[h]:=0; od;
bbb:=modp(bbb*bb,p);
hh:=modp((bbb-h0)*cc,p);
for h from hh to H-1 by p do T[h]:=0 od
fi
od; L:=[];
for h from 0 to H-1 do
if T[h]=1 then L:=[op(L),h] fi
od; L end;#
# After a pretest the probabilistic test for N-1, the Riesel
# test with D=5 and unit (1+sqrt(5))^2/4 for N-1, the
# probabilistic test for 2*N-1, and the Riesel test for 2*N-1
# with D=5 and unit (1+sqrt(5))^2/4 is used to find
# Sophie Germain primes with N=(h0+30*h)*2^n, where h runs
# from 0 to H-1. Only h0=9 and n mod 4=0 used; in this case
# the Riesel test can be used and N-1, 2*N-1 is not divisible
# by 2,3,5 hence this procedure should be used only with this
# pair. All prime divisors until P will be tried during the
# pretest.
#
sophiegermain:=proc(h0,H,n,P) local LL;
LL:=presophie(h0,H,30,2,n,7,P);
print(nops(LL));
LL:=prob(h0,30,n,-1,LL);
print(nops(LL));
LL:=exactminus(h0,n,LL);
print(nops(LL));
LL:=prob(h0,30,n+1,-1,LL);
print(nops(LL));
LL:=exactminus(h0,n+1,LL);
end;sophiegermain(9,20000,400,1000);5.14. Ikerpr\303\255m, amely Sophie Germain pr\303\255m is.#
# This part contains some special calculations
# which are useful for search of primes N
# for which N+2 and 2*N+1 are prime too.
# We try to find p in the form N=h*2^e-1.
# The test given by Riesel works with
# discriminnant D if the there are integers
# a,b,r such that (a+b*sqrt(D))^2/r is a unit
# such that r=|a^2-b^2*D|, (D|N)=-1,
# (r|D)(a^2-b^2*D)/r=-1. We are here interested
# for the choose D=13, a=3, b=1, r=4.
# The test is then applicable for those
# N, for which (13|N)=(N|13)=-1 and
# (13|2*N+1)=(2*N+1|13)=-1. This gives
# the condition h*7^e mod 13 is in {3,6,8}.
# In these cases N+2 mod 13 also not 0.
# We plan to use the modulus 30030=2*3*5*7*11*13.
# Because N, N+2 and 2*N+1 are certainly
# not divisible by 2, 3, 5, 7 and 11 if
# h divisible by 3, 5, 7 and 11, independently
# of the choose of n, we obtain a lot of cases
# in which the discriminant D=13 may be used,
# as follows:
#
# e mod 12 =0 : h mod 30030 = 5775,26565,10395
# e mod 12 =1 : h mod 30030 = 10395,5775,12705
# e mod 12 =2 : h mod 30030 = 12705,10395,28875
# e mod 12 =3 : h mod 30030 = 28875,12705,21945
# e mod 12 =4 : h mod 30030 = 21945,28875,3465
# e mod 12 =5 : h mod 30030 = 3465,21945,24255
# e mod 12 =6 : h mod 30030 = 24255,3465,19635
# e mod 12 =7 : h mod 30030 = 19635,24255,17325
# e mod 12 =8 : h mod 30030 = 17325,19635,1155
# e mod 12 =9 : h mod 30030 = 1155,17325,8085
# e mod 12 =10 : h mod 30030 = 8085,1155,26565
# e mod 12 =11 : h mod 30030 = 265775,8085,5775
#
#
# These where calculated by the following
# simple program:
#
calcsophietwinclasses:=proc() local n,L; L:=[];
for n from 0 to 11 do
[n,chrem([1,0,0,0,0,3*7^n mod 13],[2,3,5,7,11,13]),
chrem([1,0,0,0,0,6*7^n mod 13],[2,3,5,7,11,13]),
chrem([1,0,0,0,0,8*7^n mod 13],[2,3,5,7,11,13])];
L:=[op(L),%];
od; L end;
#
# More generally, such calculations can be done
# by the following program, which find all congruence classes
# for h for a given discriminant D=5,13,17,29,53 and primes from
# the list P, for which h*2^(n+e)+d<>0 mod p and all cases
# where d=-1 are appropriate for the Riesel test by D.
# List L contains the pairs [e,d].
#
classes:=proc(D::posint,n::integer,P::list(posint),L::list)
local l,i,ll,r,p,LL,q,x;
ll:=ilcm(phi(D),op(map(i->i-1,P)));
r:=[i$i=1..D];
for l in L do
if l[2]=-1 then
r:=remove(x->jacobi(D,x*2^(n+l[1])+l[2])<>-1,r);
else
r:=select(x->jacobi(D,x*2^(n+l[1])+l[2])<>0,r);
fi;
od; LL:=[[D,r]];
for p in P do
r:=[i$i=0..p-1];
for l in L do
r:=select(x->modp(x*2^(n+l[1])+l[2],p)<>0,r)
od; LL:=[op(LL),[p,r]];
od; ll,LL end;classes(5,60000,[2,3,7,11,13],[[0,-1],[1,-1]]);chrem([1,0,4],[2,3,5]);classes(13,1983*128,[2,3,5,7,11],[[0,-1],[1,-1],[0,+1]]);chrem([1,0,0,0,0,3],[2,3,5,7,11,13]);classes(13,1983*128,[2,3,5,7,11],[[0,-1],[1,+1],[0,+1]]);chrem([1,0,0,0,0,9],[2,3,5,7,11,13]);classes(13,21*2^14-128,[2,3,5,7,11],[[0,-1],[1,-1],[0,+1]]);chrem([0,0,0,0,0,1],[2,3,5,7,11,13]);classes(13,247*128,[2,3,5,7,11,13],[[0,-1],[1,-1],[2,-1]]);classes(13,273*128,[2,3,5,7,11,13],[[0,-1],[1,-1],[2,-1]]);273*128;chrem([1,0,0,0,0,3],[2,3,5,7,11,13]);ee:=273; log[10.](2.)*128*ee; (ee/247)^4.;a:=expand((3+1*sqrt(13))^2/4);riesel(5110664609396115,34944,11/2,3/2,13);riesel(5110664609396115,34945,11/2,3/2,13);riesel(5110664609396115,34946,11/2,3/2,13);2; log[2.](%); 2*3; log[2.](%); 2*3*5;log[2.](%); 2*3*5*7; log[2.](%);
2*3*5*7*11; log[2.](%); 2*3*5*7*11*13; log[2.](%);
2*3*5*7*11*13*17; log[2.](%); 2*3*5*7*11*13*17*19; log[2.](%);
2*3*5*7*11*13*17*19*23; log[2.](%); 2*3*5*7*11*13*17*19*23*29; log[2.](%);
2*3*5*7*11*13*17*19*23*29*31; log[2.](%);
2*3*5*7*11*13*17*19*23*29*31*37; log[2.](%);
2*3*5*7*11*13*17*19*23*29*31*37*41; log[2.](%);
2*3*5*7*11*13*17*19*23*29*31*37*41*47; log[2.](%);
2*3*5*7*11*13*17*19*23*29*31*37*41*47*53; log[2.](%);
2*3*5*7*11*13*17*19*23*29*31*37*41*47*53*59; log[2.](%);
2*3*5*7*11*13*17*19*23*29*31*37*41*47*53*59*61; log[2.](%);
2*3*5*7*11*13*17*19*23*29*31*37*41*47*53*59*61*67; log[2.](%);chrem([1,3],[2,13]);classes(13,61*128,[2,3,5,7,11,17],[[0,-1],[1,-1],[2,-1],[3,-1]]);chrem([1,0,0,11],[2,3,5,13]);classes(13,30*128,[2,3,5,7,11,17],[[0,-1],[1,-1],[2,-1],[3,-1],[4,-1]]);classes(13,30*128,[2,3,5,7,11,17],[[0,-1],[1,-1],[2,-1],[3,-1]]);classes(17,30*128,[2,3,5,7,11,13],[[4,-1]]);chrem([1,0,0,8],[2,3,5,13]);classes(13,29*64,[2,3,5,7,11,17,19,23,29,31],[[0,-1],[1,-1],[2,-1],[3,-1]]); cl13:=%[2];classes(17,29*64,[2,3,5,7,11,13,19,23,29,31],[[4,-1],[5,-1]]); cl17:=%[2];{op(cl13[6][2])}; {op(cl17[6][2])}; % intersect %%; nops(%);{op(cl13[8][2])}; {op(cl17[8][2])}; % intersect %%; nops(%);{op(cl13[9][2])}; {op(cl17[9][2])}; % intersect %%; nops(%);{op(cl13[10][2])}; {op(cl17[10][2])}; % intersect %%; nops(%);{op(cl13[11][2])}; {op(cl17[11][2])}; % intersect %%; nops(%);classes(17,1983*128,[2,3,5,7,11,13],[[4,-1],[5,-1],[6,-1]]);classes(5,1983*128,[2,3,7,11,13],[[0,-1],[0,+1]]);classes(5,1983*128,[2,3,7,11,13],[[0,-1],[0,+1],[1,+1]]);chrem([1,0,3,0,0,0],[2,3,5,7,11,13]);classes(5,247*128,[2,3,7,11,13],[[0,+1],[1,+1],[2,+1]]);classes(5,61*128,[2,3,7,11,13,17],[[0,+1],[1,+1],[2,+1],[3,+1]]);chrem([1,0,0],[2,3,5]);classes(5,30*128,[2,3,7,11,13,17,19],[[0,+1],[1,+1],[2,+1],[3,+1],[4,+1]]);classes(5,29*64,[2,3,7,11,13,17,19,23,29,31],[[0,+1],[1,+1],[2,+1],[3,+1],[4,+1],[5,+1]]);classes(17,0,[2,3,5,7,11],[[0,-1]]);classes(17,0,[2,3,5,7,11],[[0,-1],[1,-1]]);classes(17,0,[2,3,5,7,11],[[0,-1],[1,-1],[2,-1]]);classes(17,0,[2,3,5,7,11],[[0,-1],[1,-1],[2,-1],[3,-1]]);classes(29,0,[2,3,5,7,11],[[0,-1],[1,-1],[2,-1],[3,-1]]);classes(29,0,[2,3,5,7,11],[[0,-1],[1,-1],[2,-1],[3,-1],[4,-1]]);classes(53,0,[2,3,5,7,11],[[0,-1],[1,-1],[2,-1],[3,-1]]);classes(53,0,[2,3,5,7,11],[[0,-1],[1,-1],[2,-1],[3,-1],[4,-1]]);5.15. n^2+1 es n^4+1 alak\303\272 pr\303\255mek.#
# A pretest will be done for prime candidates having the
# form n^2+1, to be more exact, for (h0+c*h)^2*b^(2*n)+1,
# where h runs from 0 to H-1. All prime divisors from P0
# to P1 will be tried. P0 must be odd, and b,c must
# have prime factors smaller then P0. The values h, for
# which (h0+c*h)^2*b^(2*n)+1 is not divisible with primes
# from P0 to P1 are given back.
#
presquareplus:=proc(h0,H,c,b,n,P0,P1)
local T,h,p,hh,L,nn,bb,bbb,cc,bT,cT,i,j,R,r;
T:=array(0..H-1);
for h from 0 to H-1 do T[h]:=1 od;
bT:=array(0..b-1); # table of inverses for b
for i from 0 to b-1 do bT[i]:=0 od;
for i from 0 to b-1 do
for j from 0 to b-1 do
if irem(i*j+1,b)=0 then bT[i]:=j fi
od
od;
cT:=array(0..c-1); # table of inverses for c
for i from 0 to c-1 do cT[i]:=0 od;
for i from 0 to c-1 do
for j from 0 to c-1 do
if irem(i*j+1,c)=0 then cT[i]:=j fi
od
od;
for p from P0 to P1 by 2 do
if isprime(p) then
r:=msqrt(p-1,p);
if r<>FAIL then
R:=[r,p-r];
cc:=(cT[p mod c]*p+1)/c;
nn:=n mod (p-1);
bb:=(bT[p mod b]*p+1)/b;
bbb:=bb&^nn mod p;
for r in R do
hh:=(r*bbb-h0)*cc mod p;
for h from hh to H-1 by p do
T[h]:=0;
od;
od;
fi;
fi;
od;
L:=[];
for h from 0 to H-1 do
if T[h]=1 then L:=[op(L),h0+c*h] fi
od; L end;#
# This filter left only those h's from the list L in the
# resulting list for which (h0+c*h)^2*2^(2*n)+1 is
# probably prime.
#
squareplusprob:=proc(h0,c,n,L) local LL,h,N;
LL:=[];
for h in L do
N:=(h0+c*h)^2*2^(2*n)+1;
if modp(2&^(N-1),N)=1 then LL:=[op(LL),h] fi
od; LL end;#
# This filter left only those h's from the list L in the
# resulting list for which (h0+c*h)^2*2^(2*n)+1 is
# prime.
#
squareplusproth:=proc(h0,c,n,L) local LL,h,N;
LL:=[];
for h in L do
N:=(h0+c*h)^2*2^(2*n)+1;
if proth((h0+c*h)^2,2*n) then LL:=[op(LL),h] fi
od; LL end;#
# This is a test to measure densities and times.
# First pretests are done for N where
# N=(h0+2*h)^2*2^(2*n)+1 and h runs from 0 to H-1.
# The pretest is repeated for prime limit P=2^e where
# e is integer from the interval eint. The time
# and the number of remaining candidates after the pretest
# is reported, and p*ln(P) also calculated, where p is the
# relative frequency of numbers passed the pretest. The
# result of the last pretest became the input of the
# probabilistic test for N. The time and the result
# of this test is also reported.
#
squareplustest:=proc(h0,H,n,eint) local e,p,P,L,t,N;
print(`Squareplustest; h0=`,h0,`H=`,H,`n=`,n);
N:=(h0+H)^2*2^(2*n);
print(`Pretest for N`);
for e from op(1,eint) to op(2,eint) do
P:=2^e;
t:=time();
L:=presquareplus(h0,H,2,2,n,3,P);
t:=time()-t;
p:=evalf(nops(L)/H);
print(`P=`,2^e,`p=`,p,`p*ln(P)=`,evalf((ln(P))*p),`time=`,t)
od;
print(`Probabilistic test for N`);
t:=time();
L:=squareplusprob(h0,2,n,L);
t:=time()-t;
p:=evalf(nops(L)/H);
print(`p=`,p,`p*ln(N)=`,evalf((ln(N))*p),`time=`,t);
print(`Proth test for N`);
t:=time();
L:=squareplusproth(h0,2,n,L);
t:=time()-t;
p:=evalf(nops(L)/H);
print(`p=`,p,`p*ln(N)=`,evalf((ln(N))*p),`time=`,t);
L end;squareplustest(1,20000,1000,2..15);#
# Estimate of the powplusone density Cpowplusone for
# polynomial (X*2^n)^(2^e)+1 (n fix) by calculating the product for
# primes below x.
#
Cpowplusone:=proc(e::posint,x::posint) local P,p;
P:=2.; p:=nextprime(2);
while p<x do
if modp(p,2^(e+1))=1 then
P:=P*(1-2^e/p)/(1-1/p)
else
P:=P/(1-1/p)
fi;
p:=nextprime(p)
od; P end;Cpowplusone(1,10); Cpowplusone(1,100); Cpowplusone(1,1000);
Cpowplusone(1,10000); Cpowplusone(1,100000); Cpowplusone(1,1000000);Cpowplusone(2,10); Cpowplusone(2,100); Cpowplusone(2,1000);
Cpowplusone(2,10000); Cpowplusone(2,100000); Cpowplusone(2,1000000);#
# This procedure calculate the factor qeAB.
#
qeAB:=proc(e::posint,A::posint,B::posint) local P,p;
P:=1.; p:=nextprime(A-1);
while p<B do
if modp(p,2^(e+1))=1 then P:=P*(1-2^e/p) fi;
p:=nextprime(p)
od; P end;# square plus one primef:=h->((2*h+1)*2.^(340000.*4))^2+1;
g:=h->1/ln(f(h));H:=2.^21; evalf(H/6*(g(0)+4*g(H/2)+g(H)));%*2.745621017; # expected primesqeAB(1,3,1000000);%*(ln(1000000.)/(40*ln(2.)));%*H*160*20/3.054010268; # SPU time in sec%/24/3600/16; # time for one Cell blade in days# quad plus one primef:=h->((2*h+1)*2.^(340000.*2))^4+1;
g:=h->1/ln(f(h));H:=2.^15; evalf(H/6*(g(0)+4*g(H/2)+g(H)));%*5.357012604; # expected primesqeAB(2,3,1000000);%*(ln(1000000.)/(40*ln(2.)));%*H*160*20/0.09310438771; # SPU time in sec%/24/3600/16; # time for one Cell blade in days5.16. Egy\303\251b speci\303\241lis alak\303\272 pr\303\255mek.#
# Cullen and Woodall numbers: n*2^n+-1; all known below n<=8000000.
#
x:=20000000; evalf(int(y->2/y/ln(2.)/ln(ln(y)),8000000..x));#
# This procedure calculate the factor qsAB.
#
qsAB:=proc(s::posint,A::posint,B::posint) local P,p;
P:=1.; p:=nextprime(A-1);
while p<B do P:=P*(1-s/p); p:=nextprime(p) od;
P end;qsAB(1,3,1000000);%*(ln(1000000.)/(50*ln(2.))); # decreasing factor of sieve%*12000000*160*16*20; # SPU time in sec%/24/3600/16; # time for one Cell blade in days#
# Primorial numbers: p#+-1; p#><exp(p); record (02.2011): 843301#-1.
# Expected number of such primes for +1 and -1, respectively,
# up to x is exp(gamma)*ln(x), hence the next expected by
#
x:=evalf(exp(exp(-gamma)+ln(843301.))); # around this p
log[10.](exp(x)); # so many decimal digits
x/ln(x)-843301/ln(843301.); # so many to test, without sieve
%*300*128*3; # SPU time in sec%/24/3600/420; # time for all PS3 in days#
# Factorial numbers: n!+-1; n!><n^n/exp(n); record (02.2011): 103040!-1.
# Expected number of such primes for +1 and -1, respectively,
# up to x is exp(gamma)*ln(x), hence the next expected by
#
x:=evalf(exp(exp(-gamma)+ln(102030.))); # around this n
log[10.](x^x/exp(x)); # so many decimal digits
x-103040; # so many to test, without sieve
%*300*200; # SPU time in sec%/24/3600/420; # time for all PS3 in days#
# Primes in arithmetic sequences, 3 primes; record (02.2011):
# 1185*2^529445-1539*2^263359+1, d=1185*2^529444-1539*2^263359
# (159382 decimal digits).
#
evalf(2^64/3/5/7/11/13/17/19/23/29/31); # so many to sieve out
(2/1)*(3/2)*(5/4)*(7/6)*(11/10)*(13/12)*(17/16)*(19/18)*(23/22)*(29/28)*(31/30); evalf(%); # so many times more denseqsAB(1,32,1000000);%*(ln(1000000.)/(50*ln(2.))); # decreasing factor of sieve1/200000/ln(10.); # natural density of primes having 200000 digits%*6.542269709/0.1059825800; 1/%; # increased density and expected tests# after 30000000 test we have 400 primes, 80000 pairs, seems to be enough3000000*160*5; # SPU time in sec%/24/3600/16; evalf(%); # time for one Cell blade in days1/160000/ln(10.); # natural density of primes having 160000 digits%*6.542269709/0.1059825800; 1/%; # increased density and expected tests# after 1800000 test we have 300 primes, 45000 pairs, seems to be enough1800000*300*12; # SPU time in sec%/24/3600/420; evalf(%); # time for all PS3 in days# A ten-year jump:
1/310000/ln(10.); # natural density of primes having 310000 digits%*6.542269709/0.1059825800; 1/%; # increased density and expected tests# after 5800000 test we have 500 primes, 125000 pairs, seems to be enough
5800000*300*24; # SPU time in sec%/24/3600/420; evalf(%); # time for all PS3 in daysrepunit:=proc(n) local i;
for i to n do if isprime((10^i-1)/9) then print(i) fi od;
end;repunit(1000);5.17. K\303\241tai egy probl\303\251m\303\241ja.#
# This is a simple procedure to get all the prime pairs (p,q)
# for which the quotient (p+1)/(q+1)=n is an integer,
# 1<n<=N and 2<q<=Q. The limit Q is supposed to be small
# and N is supposed to be large, hence primality testing is used.
# The result is a list of pairs; each pair consist of n
# and a list of all pairs [p,q] corresponding to this n.
# The pairs are in increasing order.
#
prime_plus_one_quotient:=proc(Q,N)
local q,n,p,QQ,LL,LLL; QQ:=[];
for q from 3 by 2 while q<=Q do
if isprime(q) then QQ:=[op(QQ),q] fi;
od;
LL:=[];
for n from 2 while n<=N do
LLL:=[];
for q in QQ do
p:=n*q+n-1;
if isprime(p) then LLL:=[op(LLL),[p,q]] fi;
od;
LL:=[op(LL),[n,LLL]];
od; LL end;L:=prime_plus_one_quotient(10,100);#
# This routine collect statistics from the results of the
# above routine. For all number 0<=m<=M, where M is
# the number of odd prime less then Q, the number of those
# numbers 1<n<=N given back, for which exactly m solutions
# of n=(p+1)/(q+1) is found. L is the list calculated by the
# above program.
#
stat_prime_plus:=proc(Q,L) local q,M,MM,LL; M:=0;
for q from 3 by 2 while q<Q do
if isprime(q) then M:=M+1 fi;
od;
MM:=array(0..M);
for q from 0 to M do MM[q]:=0; od;
for LL in L do MM[nops(LL[2])]:=MM[nops(LL[2])]+1; od;
convert(MM,listlist) end;stat_prime_plus(10,L);5.18. P\303\251lda.5.19. Feladat.#
# This is a simple factorization
# procedure using trial division.
# The result is a sequence of pairs
# [p,e] where the p's are the prime
# factors and the e's are the exponents.
# The factors are anyway in increasing order.
# Only primes <= P are tried, hence the
# last "factor" may composite, if
# it is greater then P^2;
#
trialdiv:=proc(n::posint,P::posint) local L,p,i,d,nn;
L:=[]; nn:=n;
if type(nn,even) and 2<=P then
for i from 0 while type(nn,even) do nn:=nn/2; od;
L:=[[2,i]];
fi;
if nn mod 3=0 and 3<=P then
for i from 0 while nn mod 3=0 do nn:=nn/3; od;
L:=[op(L),[3,i]];
fi;
d:=2; p:=5;
while p<=P and nn>=p^2 do
if nn mod p=0 then
for i from 0 while nn mod p=0 do nn:=nn/p; od;
L:=[op(L),[p,i]];
fi;
p:=p+d; d:=6-d;
od;
if nn>1 then L:=[op(L),[nn,1]] fi;
L;
end;trialdiv(2^107-2^54+1,1000);n0:=%[3][1];modp(3&^(n0-1),n0);trialdiv(n0-1,1000);n1:=%[5][1];modp(3&^(n1-1),n1);trialdiv(n1,100000);n2:=%[2][1];modp(3&^(n2-1),n2);trialdiv(n2-1,1000);n3:=%[4][1];modp(3&^(n3-1),n3);trialdiv(n3,10000);n4:=%[2][1];modp(3&^(n4-1),n4);trialdiv(n4-1,1000);trialdiv(2^107+2^54+1,1000);n5:=%[1][1];modp(3&^(n5-1),n5);trialdiv(2^107+2^54+1,1000000);n6:=%[2][1];modp(3&^(n6-1),n6);5.20. Pr\303\255msz\303\241mk\303\263dol\303\241s.5.21. Feladat.p:=safeprime(1563456788814256178886765661555261342987645321345665432123456788772091275121098761232122333321233434123432123344454321234320948725467845467788859812342365); log[2.](p);
q:=safeprime(29841524475159001676561453467890987651234254321456490998767626788182514325678909987236514234154232396587874778993377722004988376667882767156363888377626677728888); log[2.](q);
n:=p*q;
e:=2876354132453678909987653432123409887635423125;
igcdex(e,(p-1)*(q-1),'d'); d; d*e mod (p-1)*(q-1);5.22. Feladat.M:="Mint v\303\255z alatti, elmer\303\274lt harangok
hint\303\241znak-e hajnalonk\303\251nt \303\241gyadn\303\241l
a tizennyolc \303\251ves iskol\303\241sok
kiket felakasztatt\303\241l";
convert(M,'bytes');
m:=sum(%[i]*256^(i-1),i=1..nops(%));
c:=m&^e mod n;c&^d mod n; convert(%,base,256); convert(%,'bytes');#
# This simple procedure generates a pair of safe primes p and q
# congruent 3 mod 4, and calculates the product n=p*q,
# a Blum integer. The output is a sequence of these numbers
# in this order. The input is two large numbers, these are the
# starting points for the search for p and q, respectively.
# For example, 100 and 104 digit starting points may be used.
# The best to choose them randomly and independently.
#
setBlum:=proc(p0,q0) local p,q,n,e,d,k;
p:=safeprime(p0);
while p mod 4 <> 3 do p:=safeprime(p) od;
q:=safeprime(q0);
while q mod 4 <> 3 do q:=safeprime(q) od;
n:=p*q; p,q,n; end;#
# This inner routine make the encryption and decryption
# for an arbitrary texts. Here n is the Blum integer,
# L is the text, x0 is the first quadratic residue,
# k is the number of bits used in one step.
#
Blumi:=proc(n,L,x0,k) local xi,b,i,j,nL;
nL:=[]; xi:=x0;
for i from 0 to nops(L)/k-1 do
b:=convert(xi mod 2^k,base,2);
for j to k do
if nops(b)>j then nL:=[op(nL),L[i*k+j]+b[j] mod 2]
else nL:=[op(nL),L[i*k+j]] fi
od;
xi:=xi &^ 2 mod n;
od; xi,nL end;#
# This make the encryption for arbitrary texts.
# Here n is the Blum integer, L is the text, x is the
# random number to generate the first quadratic residue,
# k is the number of bits used in one step.
#
Blume:=proc(n,x,L,k) local pt,i,j,nL;
if igcd(x,n)>1 then ERROR(`x is not relative prime to n`) fi;
convert(L,'bytes');
pt:=sum(%[i]*256^(i-1),i=1..nops(%));
pt:=convert(pt,base,2);
if k>1 then pmod(nops(pt),k);
if %<>0 then pt:=[op(pt),0$i=1..k-%] fi;
fi;
Blumi(n,pt,x &^ 2 mod n,k) end;#
# This make the decryption for arbitrary texts.
# Here p and q are the primes, x is the last quadratic residue,
# L is the string, and k is the number of bits used in one step.
#
Blumd:=proc(p,q,x,L,k) local a,b,u,v,i,n,t,nL;
igcdex(p,q,'a','b'); t:=nops(L);
u:=x &^ (((p+1)/4) &^ t mod (p-1)) mod p;
v:=x &^ (((q+1)/4) &^ t mod (q-1)) mod q;
nL:=Blumi(p*q,L,a*u+b*v mod p*q,k)[2]; n:=0;
for i from t to 1 by -1 do n:=n+n+nL[i]; od;
convert (convert(n,base,256),'bytes'); end;p,q,n:=setBlum(10^100,10^104);x,L:=Blume(n,55555,M,1);Blumd(p,q,x,L,1);#
# This simple procedure generates a safe prime q, and a generator
# g of the reduced residue classes. For the given exponent d,
# calculates e=g^d mod q. The output is a sequence
# of the numbers [q,g,e,d] in this order. The input is
# three large numbers, these are the starting points for the search
# for q and e, respectively, and the exponent d.
# For example, 200 digit starting points may be used.
# The best to choose them randomly and independently.
#
setElGamal:=proc(q0,g0,d) local f,p,q,g,e,P;
q:=safeprime(q0); P:=factorset(q-1); f:=false; g:=g0;
while not f do f:=true;
for p in P while f=true do
if modp(g &^ ((q-1)/p),q) = 1 then f:=false; g:=g+1 fi
od;
od; e:=modp(g &^ d,q); [q,g,e,d]; end;#
# This make the encryption for short texts. Here q is the modulus,
# g is the generator, e the encryption exponent, k the random
# exponent, L the text.
#
ElGamale:=proc(q,g,e,k,L) convert(L,'bytes');
modp(g &^ k,q),modp(sum(%[i]*256^(i-1),i=1..nops(%))*(e &^ k),q) end;#
# This make the decryption for short texts. Here q is the modulus,
# d the decryption exponent, x the number.
#
ElGamald:=proc(q,g,d,x,y)
convert(convert(modp((x &^ (q-d-1))*y,q),base,256),'bytes'); end;L:=setElGamal(10^200,2*10^199,3*10^199); q:=L[1]; g:=L[2]; e:=L[3]; d:=L[4];LL:=ElGamale(q,g,e,5555,M[1..36]);ElGamald(q,g,d,LL[1],LL[2]);debug(ElGamal);6. Sz\303\241mok \303\251s polinomok7. Gyors Fourier-transzform\303\241ci\303\2638. Elliptikus f\303\274ggv\303\251nyek9. Sz\303\241mol\303\241s elliptikus g\303\266rb\303\251ken10. Faktoriz\303\241l\303\241s elliptikus g\303\266rb\303\251kkel11. Pr\303\255mteszt elliptikus g\303\266rb\303\251kkel12. Polinomfaktoriz\303\241l\303\241s13. Az AKS-teszt14. A szita m\303\263dszerek alapjai15. Sz\303\241mtest szita16. Vegyes probl\303\251m\303\241kLUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0Yn