# This files contains data to compute function h(n) # # for n small (say n <= 10^12) or n = 10^a (a <= 35) or n <= 2^a (a<= 117) # # See the paper : "Minimal product of primes whose sum is bounded" # # jointly written by Marc Deleglise and Jean-Louis Nicolas ########################################################## ########################################################## ## Section 1 : Computation of G(p,m) ########################################################## ########################################################## ########################################################## ## Combinatorial algorithm for G(p,m) ########################################################## cherche:=proc(liste,m,ini,ifi) local a,b,c; # "cherche" means look for; "liste" is a list of numbers in ascending order # If m >= liste[ifi], returns ifi; if m < liste[ini] returns "error"; # if not, returns c such that ini <= c < ifi and liste[c] <=m < liste[ifi]. # If the list is not-decreasing, returns the maximal such c if nops(liste)=0 or liste[ini] > m or ini <= 0 or ifi > nops(liste) then error "bad choice of m=%1, ini=%2, ifi=%3 or nopsliste=%4",m,ini,ifi,nops(liste); fi; if m >= liste[ifi] then return ifi; fi; a:=ini; b:=ifi; # here, liste[ini] <= m and liste[ifi] < m while b-a > 1 do c:=iquo(a+b,2); if liste[c] <= m then a:=c; else b:=c; fi; # liste[b] > m holds od; a; end; mj:=proc(j,K,r) local a; global S; # This short procedure calculates m_j(P_r) defined in the article, (7.17). if j <= r then a:=S[K+j]-S[K]-S[r]+S[r-j]; else a:=infinity; fi; a; end; Gcomb:=proc(pk,Mmax,imp) # imp is a printing parameter. If imp < 5 nothing is printed. # This procedure computes G(pk,m) for m even, 0 <= m <= Mmax <= pkp1-3 # where pk is a prime and pkp1 is the prime following pk. # The result is returned in a table GG[m]. Mmax is supposed not too large. # First, we compute the primes P[1],...,P[K]=pk,...,P[Rhat] # This procedure uses the combinatorial method explained in the article, § 8. # By induction on r and j, H(j,P[r],m) is saved in the memory HH[j,m]. local pkp1,P,r,Rhat,K,i,j,j1,mj1,m,bou,maxP,maxPgen,pkjr, GG,HH,pro,rmes,mmes,nu, tem,ti; # variables used for the time global S; tem:=time(); ti:=tem; # initialization of the time if not isprime(pk) then error "pk=%1 is not prime",pk;fi; pkp1:=nextprime(pk); if Mmax > pkp1-3 then error "Mmax=%1 is too big",Mmax; fi; if Mmax < pkp1-pk then for m from 0 by 2 to Mmax do GG[m]:=1; od; return GG; fi; # Determination of P[1] (cf. (8.7)) P[1]:=nextprime(pkp1-Mmax-1); r:=1; # Note that Mmax >= pkp1-pk implies P[1] <= pk. # Determination of the primes P[r], r=1..K while P[r] < pk do r:=r+1; P[r]:=nextprime(P[r-1]); od; K:=r; # here, we have P[K]=pk # Determination of large primes. # These primes are limited to P[Rhat] (cf. § 8.4) Rhat:=K+10; # 10 seems a reasonnable value. for r from K+1 to Rhat+1 do P[r]:=nextprime(P[r-1]); od; # P[Rhat+1] will only be used to compute maxPgen S[0]:=0; #S is the sommatory function of P[r]. for r from 1 to Rhat do S[r]:=S[r-1]+P[r]; od; if imp >= 6 then print(Function Gcomb, pk=,pk,Mmax=,Mmax,Rhat=,Rhat, K=,K,P[K]=,P[K],small primes=, [P[r1] $r1=1..K-1], large primes=, [P[r1]$ r1=K+1..Rhat], time of computation of the primes=,time()-ti); ti:=time(); fi; pro:=1; # computation of the product pro:=P[K+1]...P[Rhat] for r from K+1 to Rhat do pro:=pro*P[r]; od; # The boolean variable "bou" will be used to eventually start again the # "while" loop, if maxPgen (= rhs of (8.23)) is not large enough. bou:=true; while bou do for m from 0 by 2 to Mmax do HH[0,m]:=1; od; # formula (8.22) for r from 1 to Rhat do # computation by induction of H(j,p[r],m)=HH[j,m] for j from min(Rhat-K,r) by -1 to max(1,r-K) do j1:=j-1; mj1:=mj(j,K,r-1); pkjr:=P[K+j]-P[r]; for m from mj(j,K,r) by 2 to min(Mmax,mj1-2) do HH[j,m]:= P[r]*HH[j1,m-pkjr]; # formula (8.21) od; for m from mj1 by 2 to Mmax do HH[j,m]:= min(HH[j,m],P[r]*HH[j1,m-pkjr]); # formula (8.20) od; od; # (j) if imp >= 8 then print(r=,r,time=,time()-ti,time()-tem); ti:=time(); fi; od; # (r) # Determination of maxPgen. This parameter is an upper bound # for the large primes (cf. formula (8.23)). maxPgen:=0; for m from 0 by 2 to Mmax do GG[m]:=pro/HH[Rhat-K,m]; # computation of G(pk,m) by formula (8.16) if GG[m] > 1 then maxP:=min(m+pk,floor(GG[m]*m/(GG[m]-1))); # maximum for m else maxP:=m+pk; # maximum for m fi; maxPgen:=max(maxPgen,maxP); # maximum for all m's od; # (m) # If maxPgen >= P[Rhat+1], P[Rhat+1] can occur, we have to continue. if maxPgen < P[Rhat+1] then bou:=false; # the computation is over if imp >= 5 then print(P[Rhat+1]=,P[Rhat+1],is > maxPgen=,maxPgen, Rhat-K=,Rhat-K,real increment=, cherche([P[r1] $r1=K..Rhat+1],maxPgen,1,Rhat+2-K)-1, the computation is over,time=,time()-ti); ti:=time(); fi; else if imp >= 5 then print(P[Rhat+1]=,P[Rhat+1],is <= maxPgen=,maxPgen, we start again with more primes,time=,time()-ti); ti:=time(); fi; r:=Rhat+1; # We have to add primes while P[r] <= maxPgen do pro:=pro*P[r]; S[r]:=S[r-1]+P[r]; r:=r+1; P[r]:=nextprime(P[r-1]); od; # here, P[r] > maxPgen holds Rhat:=r-1; fi; # If maxPgen was too small, we start again the loop on "bou" od; # (bou) if imp >= 5 then # we measure the largest prime which has been used rmes:=K+1; for m from 0 by 2 to Mmax do nu:=numer(GG[m]); for r from rmes+1 to Rhat do if irem(nu,P[r])=0 then rmes:=r; mmes:=m; fi; od; # (r) od; # (m) print(pmaxreel=,P[rmes],increm=,rmes-K,mmes=,mmes,total time=, time()-tem); fi; GG; end; ########################################################## ## Computation of G(p,m) by delta-method ########################################################## Gpkmh:=proc(pk,m,maxdinit,imp) local mm,G1,G2,delta1max,maxd,maxdelta,pkp1,pkp2,q,d,oldd, Gprov,Gpr,Gout,bou,qhat,delta,pkp1mmm,dmax,lisq,lisd,comp,compqmax,ecart,qmax,typ, qdelta,compdelta,maxiq,maxid,mot; # This procedure computes G(pk,m) as explained in the article [DNZ] # Deleglise, Nicolas, Zimmermann, J. de Bordeaux, # "Landau's function for one million billions", § 9. # Or in the paper Deleglise, Nicolas [DN] # "Maximal product of primes whose sum is bounded" # imp is a printing parameter. if imp=0, nothing is written # if imp <= 4 no remarks about the combinatorial algorithm is written # The parameter "maxdinit" is the initial value for "maxd" which is # the upper bound for the variable "d" # "maxd" can possibly increase during the execution # A good choice for maxdinit is 500 # The proc returns a list Gout whose three first elements are # Gout = [G(p_k,m),ecart,typ,...] # ecart is the greatest integer st G(p_k,m-ecart) = G(p_k,m) # typ is a number in {0,1,2,3,4} which indicates the way of computation. # if typ = 0, 1 or 2, Gout = [G(p_k,m),ecart,typ] has only 3 terms. # if typ = 0, then m satisfies p_{k+1}-3 < m < p_{k+1}. # The value of G(p_k,m) = p_{k+1}/2 follows from [DN, Prop. 5.3]. # If typ = 1 then m satisfies 0 \leq m < p_{k+1}-p_k # The value G(p_k,m) = 1 follows from [DNZ, (89)]. # If typ = 2 then m satisfies p_{k+1} - m is prime. # The value G(p_k,m) = p_{k+1}/(p_{k+1}-m) follows from [DNZ, Prop. 8]. # if typ = 3, the list Gout has 6 terms : # Gout = [G(p_k,m),ecart,typ,delta,maxid,nops(lisq)] # delta is the value given by [DNZ, Prop. 10] # maxid is the greatest value of d for which we have to compute # G(p_{k+1},d) in [DNZ, (113)] # lisq is the list of those q's which occurs as indices in [DNZ, (113)]. # if typ = 4, the combinatorial algo is used to compute G(p_k,m) # Gout = [G(p_k,m),ecart,typ,mm = m-irem(m,2)] ########## ########## if not isprime(pk) then error "pk=%1 is not prime",pk;fi; pkp1:=nextprime(pk); if m > pkp1-1 or m < 0 or not type(m,integer) then error "m=%1 is wrong",m; fi; if imp >= 4 then print(&&& p_k=,pk,n'=,m,&&&); fi; if m >= pkp1-2 then # typ = 0, m is greater than p_{k+1}-3; typ:=0; ecart:=m-pkp1+2; Gprov:=pkp1/2; Gout:=[Gprov,ecart,typ]; if imp >= 2 then print("nprime is greater than p_{k+1}-3"); fi; return Gout; fi; if m < pkp1-pk then # typ = 1, m is smaller than p_{k+1}-p_k; typ:=1; ecart:=m; if imp >= 2 then print("nprime is smaller than p_{k+1}-p_k"); fi; Gout:= [1,ecart,typ]; return Gout; fi; mm:=m-irem(m,2); # mm is the greatest even number <= m if irem(m,2)=0 then mot:="mm=nprime" else mot:="mm=nprime - 1"; fi; # We have G(pk,m)=G(pk,mm). (see [DNZ, § 8.1]) pkp1mmm:=pkp1-mm; if isprime(pkp1mmm) then # typ = 2, see [DNZ, Prop. 8] typ:= 2; ecart:=irem(m,2); if imp >= 2 then print(mot," p_{k+1}-mm is prime"); fi; return [pkp1/(pkp1mmm),ecart,typ]; fi; pkp2:=nextprime(pkp1); delta1max:=2*trunc(1.275*log(evalf(pk))^2); # delta1max is an empirical upper bound for # delta_1(pk) defined in [DNZ, (118)] if irem(mm,9,'qq')=0 then maxdelta:=2*qq-2 else maxdelta:=2*qq; fi; # here, maxdelta is the greatest even integer < 2m/9, cf. [DNZ, (112)] maxdelta:=min(maxdelta,delta1max); # it is our upper bound to look for delta maxd:=min(pkp2-3,maxdinit-irem(maxdinit,2)); # we choose "maxd" even if imp >= 3 then if delta1max <> maxdelta then print(maxdinit=,maxdinit,maxd=,maxd,delta1max=,delta1max,maxdelta=,maxdelta); else print(maxdinit=,maxdinit,maxd=,maxd,maxdelta=,maxdelta); fi; fi; G1:=Gcomb(pkp1,maxd,imp); # Here, we use the combinatorial algorithm q:=nextprime(pkp1mmm); # here, we are sure that pkp1-mm is not prime # q is the smallest prime which may occur in the denominator of G(pk,m) lisq:=[]; # It will be the list of the indices "q" in [DNZ, (113)] lisd:=[]; # It will be the list of the corresponding "d = mm-p_{k+1}+q" comp:=0; # it is a counter for the primes "q" # Now, we are looking for delta and, simultaneously, compute the max [DNZ, (113)] d:=q-pkp1mmm; # d is the smallest possible value for delta bou:=false; # bou will turn true when delta is found Gprov:=1; while d <= maxdelta and not bou do lisq:=[op(lisq),q]; comp:=comp+1; lisd:=[op(lisd),d]; if d > maxd then while d > maxd do maxd:=min(2*maxd,pkp2-3); od; print(new value of maxd=,maxd); G1:=Gcomb(pkp1,maxd,imp); # Again, the combinatorial algorithm fi; Gpr:=pkp1/q*G1[d]; # formula [DNZ, (113)] if Gpr > Gprov then Gprov:=Gpr; dmax:=d; qmax:=q; compqmax:=comp; fi; bou:=evalb(G1[d] >= 1+d/pkp1); # formula [DNZ, (111)] q:=nextprime(q); oldd:=d; d:=q-pkp1mmm; od; if bou then # here, we have found delta delta:=oldd; qdelta:=q; compdelta:=comp; qhat:=iquo(pkp1*pkp2*(pkp1mmm+delta),(pkp1+delta)*(pkp1-3*delta/2)); maxiq:=prevprime(qhat+1); # this is the largest q used in [DNZ, (113)] maxid:=maxiq-pkp1mmm; # the largest value of d = mm-p_{k+1}+maxiq in [DNZ, (113)] if maxid > maxd then maxd:=maxid; # It can be shown that this value is < pkp2-3 print(new value of maxd = ,maxd, (= maxid)); G1:=Gcomb(pkp1,maxd,imp); # Again, the combinatorial algorithm fi; while q <= qhat do lisq:=[op(lisq),q]; comp:=comp+1; lisd:=[op(lisd),d]; Gpr:=pkp1/q*G1[d]; # formula (9.4) if Gpr > Gprov then Gprov:=Gpr; dmax:=d; qmax:=q; compqmax:=comp; fi; q:=nextprime(q); d:=q-pkp1mmm; od; # here the delta-method is finished if imp >= 2 then print("[number of q's, compdelta, compqmax]=",[comp,compdelta,compqmax]); fi; if imp >= 2 and comp < 100 then print("list of trialed q's=",lisq); print("list of corresponding d's=",lisd, "delta=",delta); fi; d:=dmax-2; # computation of "ecart" while G1[d] = G1[dmax] do d:=d-2; od; ecart:=dmax-d-2+irem(m,2); typ:=3; Gout:=[Gprov,ecart,typ,delta,maxid,comp]; else # here, the delta-method has failed if imp >= 1 then if mm > 9*delta1max/2 then print("m=",m,"The delta-method has failed, we should have found delta"); else print("The delta-method has failed, but nprime is small"); fi; fi; G2:=Gcomb(pk,mm,imp); # we use the comb algo to directly compute G(pk,m) Gprov:=G2[mm]; d:=mm-2; # computation of "ecart" while G2[d] = G2[mm] do d:=d-2; od; ecart:=mm-d-2+irem(m,2); typ:=4; Gout:=[Gprov,ecart,typ,mm]; fi; Gout; end; fact:=proc(GG,pp) local i,inu,lisinu,lispnu,lisiden,lispden,qq,p,com; # factos the value of GG=G(p_k,m) with p=p_k obtained from the proc Gpkmh if GG=1 then return [[],[],[],[]]; fi; qq:=numer(GG); lisinu:=[]; lispnu:=[]; inu:=1; p:=nextprime(pp); while qq > 1 do while irem(qq,p) <> 0 do inu:=inu+1; p:=nextprime(p); od; qq:=qq/p; lisinu:=[op(lisinu),inu]; lispnu:=[op(lispnu),p]; od; # ihere, lisinu = [i1,i2,...] and lispnu = [p_{k+i1},p_{k+i2},...] qq:=denom(GG); lisiden:=[]; lispden:=[]; i:=0; p:=pp; com:=0; while com < nops(lisinu)-1 do while irem(qq,p) <> 0 do i:=i+1; p:=prevprime(p); od; qq:=qq/p; com:=com+1; lisiden:=[op(lisiden),-i]; lispden:=[op(lispden),p]; od; # here, com = nops(lisinu)-1, nops(lisiden) = com - 1, # lisidem = [-i0,-i1,...] and lispnu = [p_{k-i0},p_{k-i1},...,qq] if not isprime(qq) then error "qq=%1 should be prime", qq; fi; lispden:=[op(lispden),qq]; [lisinu,lispnu,lisiden,lispden]; end; ############################################################### ############################################################### # Section 2 : computation of h(n) for n small (say <= 3 10^12) ############################################################### ############################################################### # The procedure "hdennext" computes the sum of consecutive primes # by generating primes with the procedure "nextprime". # It takes a couple of minutes for n around 10^12. hdennext:=proc(n) local sig,pk,nprime,A,lis; # computes first p_k, sigma_k and nprime = n-sigma-k # and then G(p_k,nprime) sig:=0; pk:=0; while sig <= n do pk:=nextprime(pk); sig :=sig+pk; od; sig:=sig -pk; pk:=prevprime(pk); nprime:=n-sig; print("***** n=",n,"p_k=",pk,"nprime=",nprime,"*****"); print(); A:=Gpkmh(pk,nprime,500,2); lis:=fact(A[1],pk); print("factors of G(p_k,nprime)=",op(lis)); print(A[1],"ecart=",A[2],"type=",A[3]); end; ############################################################### # A sieving procedure to calculate h(n) ############################################################### # Before actuvating the procedure "hdensieve" to calculate h(n), we # generate the primes and the sum of consecutive primes # by Eratosthenes's sieve. raci:=proc(n) local rac; # returns the largest odd integer rac s.t. rac^2 <= n rac:=isqrt(n); # rac is the closest integer of the square root of n if rac*rac > n then rac:=rac-1; fi; # rac is the largest integer s.t. rac^2 <= n rac-1+irem(rac,2); end; crible:=proc(nmax) local tab,n,k,a,q,rac,sigma; global p,sigmalis; # generates the primes p_k up to nmax and the sums # sigma_k = p_1+p_2+...+p_k # Writes the maximal value of n for which h(n) can be computed p:='p'; for n from 3 by 2 to nmax do tab[n]:=true; od; rac:=raci(nmax); k:=1; p[k]:=2; sigma[k]:=2; for q from 3 by 2 to rac do; if tab[q] then k:=k+1; p[k]:=q; sigma[k]:=sigma[k-1]+q; a:=q+q; for n from q*q by a to nmax do tab[n]:=false; od; fi; od; for q from rac + 2 by 2 to nmax do if tab[q] then k:=k+1; p[k]:=q; sigma[k]:=sigma[k-1]+q; fi; od; print("You are allowed to compute hden(n) by proc hdensieve up to n=",sigma[k]); sigmalis:=[sigma[kk]$kk=1..k]; end; hdensieve:=proc(n) local sig,pk,nprime,k,nosig,A,lis; global p,sigmalis; # First, be sure that the procedure "crible" has been activated # with nmax >= n # first finds in the tables "p" and "sigma", k, p_k, sigma_k, # nprime = n-sigma_k and further, G(p_k,nprime) nosig:=nops(sigmalis); if n > sigmalis[nosig] then error "the sieving is not long enough" fi; k:=cherche(sigmalis,n,1,nosig); sig:=sigmalis[k]; pk:=p[k]; nprime:=n-sig; print("***** n=",n,"p_k=",pk,"nprime=",nprime,"*****"); print(); A:=Gpkmh(pk,nprime,500,2); lis:=fact(A[1],pk); print("factors of G(p_k,nprime)=",op(lis)); print(A[1],"ecart=",A[2],"type=",A[3]); end; ########################################################## ########################################################## ## Section 3 : computation of h(n) for n = 10^a ########################################################## ########################################################## # Below are the data # which allows the computation of h(n) for n=10^a. # They have been obtained by Marc Deleglise # in a way described in the paper # "Minimal product of primes whose sum is bounded" # jointly written by Marc Deleglise and Jean-Louis Nicolas # More precisely, it contains the variables # pp[a] = largest prime s.t. pi_{id}(p) <= 10^a # pi_{id}(x) is the sum of primes <= x # qq[a] = nextprime(pp[a]) # Del[a] = 10^a-pi_{id}(pp[a]) # for 1 <= a <= 35. pp[1]:=5 ; qq[1]:=7 ; Del[1]:=0; pp[2]:=23 ; qq[2]:=29 ; Del[2]:=0; pp[3]:=89 ; qq[3]:=97 ; Del[3]:=37; pp[4]:=331 ; qq[4]:=337 ; Del[4]:=146; pp[5]:=1151 ; qq[5]:=1153 ; Del[5]:=315; pp[6]:=3931 ; qq[6]:=3943 ; Del[6]:=2339; pp[7]:=13421 ; qq[7]:=13441 ; Del[7]:=8555; pp[8]:= 45191 ; qq[8]:=45197 ; Del[8]:=44820; pp[9]:=151057 ; qq[9]:=151091 ; Del[9]:=133215; pp[10]:=502247 ; qq[10]:=502259 ; Del[10]:=78132; pp[11]:=1662361 ; qq[11]:=1662377 ; Del[11]:=679220; pp[12]:=5477081 ; qq[12]:=5477083 ; Del[12]:=4935150; pp[13]:=17993471 ; qq[13]:=17993477 ; Del[13]:=17687910; pp[14]:=58954409 ; qq[14]:=58954411 ; Del[14]:=18206505; pp[15]:=192682397 ; qq[15]:=192682537 ; Del[15]:=79297359; pp[16]:=628420087 ; qq[16]:=628420129 ; Del[16]:=468817588; pp[17]:=2045812129 ; qq[17]:=2045812157 ; Del[17]:=530434512; pp[18]:=6649205947 ; qq[18]:=6649205951 ; Del[18]:=4975017306; pp[19]:=21579536887 ; qq[19]:=21579536893 ; Del[19]:=7802551562; pp[20]:=69943284083 ; qq[20]:=69943284097 ; Del[20]:=37512814187; pp[21]:=226431374401 ; qq[21]:=226431374423 ; Del[21]:=142613249827; pp[22]:=732253012133 ; qq[22]:=732253012171 ; Del[22]:=451575818329; pp[23]:=2365707119569 ; qq[23]:=2365707119579 ; Del[23]:=2123854653548; pp[24]:=7636102747781 ; qq[24]:=7636102747789 ; Del[24]:=3265879922023; pp[25]:=24627776204923 ; qq[25]:=24627776205037 ; Del[25]:=19762898291976; pp[26]:=79368662592301 ; qq[26]:=79368662592367 ; Del[26]:=41809634354534; pp[27]:=255604186937189 ; qq[27]:=255604186937213 ; Del[27]:=6459083731725; pp[28]:=822628809813037 ; qq[28]:=822628809813047 ; Del[28]:=815508130025257; pp[29]:=2645919057430751 ; qq[29]:=2645919057430769 ; Del[29]:=2501502719639271; pp[30]:=8505572989358131 ; qq[30]:=8505572989358167 ; Del[30]:=7343449508390437; pp[31]:=27327564443962931 ; qq[31]:=27327564443962997 ; Del[31]:=13877499683805723; pp[32]:=87757258184451287 ; qq[32]:=87757258184451299 ; Del[32]:=25633740738769901; pp[33]:=281684563123416169 ; qq[33]:=281684563123416247 ; Del[33]:=227862450865627006; pp[34]:=903759496776593377 ; qq[34]:=903759496776593423 ; Del[34]:=679498720026090893; pp[35]:=2898434150644708999 ; qq[35]:=2898434150644709023 ; Del[35]:=1886081812111845520; # pp[36]:=9291904456147945403 ; qq[36]:=9291904456147945441 ; # Del[36]:=991398840103841912; ####################################################################### # Computation of h(10^a) ####################################################################### hdedixpuisa:=proc(amin,amax) local a,nprime,pk,A,n,lis; # computes h(10^a) for amin <= a <=amax; for a from amin to amax do nprime:=Del[a]; pk:=pp[a]; print("*** n=10^",a,"p_k=",pk,"nprime=",nprime,"***"); print(); A:=Gpkmh(pk,nprime,500,2); lis:=fact(A[1],pk); # fact is a factoring procedure print("factors of G(p_k,nprime)=",op(lis)); print(A[1],"ecart=",A[2],"type=",A[3]); print();print(); od; end; ########################################################## ########################################################## ## Section 4 : computation of h(n) for n = 2^a ########################################################## ########################################################## # Below are the data # which allows the computation of h(n) for n = 2^a. # They have been obtained by Marc Deleglise # in a way described in the paper # "Minimal product of primes whose sum is bounded" # jointly written by Marc Deleglise and Jean-Louis Nicolas # More precisely, it contains the variables # ppp[a] = largest prime s.t. pi_{id}(p) <= 10^a # pi_{id}(x) is the sum of primes <= x # qqq[a] = nextprime(ppp[a]) # Delta[a] = 2^a-pi_{id}(ppp[a]) # for 1 <= a <= 117. a:= 1 ;nnn[a]:= 2 ;ppp[a]:= 2 ;qqq[a]:= 3 ;Delta[a]:= 0; a:= 2 ;nnn[a]:= 4 ;ppp[a]:= 2 ;qqq[a]:= 3 ;Delta[a]:= 2; a:= 3 ;nnn[a]:= 8 ;ppp[a]:= 3 ;qqq[a]:= 5 ;Delta[a]:= 3; a:= 4 ;nnn[a]:= 16 ;ppp[a]:= 5 ;qqq[a]:= 7 ;Delta[a]:= 6; a:= 5 ;nnn[a]:= 32 ;ppp[a]:= 11 ;qqq[a]:= 13 ;Delta[a]:= 4; a:= 6 ;nnn[a]:= 64 ;ppp[a]:= 17 ;qqq[a]:= 19 ;Delta[a]:= 6; a:= 7 ;nnn[a]:= 128 ;ppp[a]:= 23 ;qqq[a]:= 43 ;Delta[a]:= 28; a:= 8 ;nnn[a]:= 256 ;ppp[a]:= 41 ;qqq[a]:= 43 ;Delta[a]:= 18 ;a:= 9 ;nnn[a]:= 512 ;ppp[a]:= 61 ;qqq[a]:= 67 ;Delta[a]:= 11 ;a:= 10 ;nnn[a]:= 1024 ;ppp[a]:= 89 ;qqq[a]:= 97 ;Delta[a]:= 61 ;a:= 11 ;nnn[a]:= 2048 ;ppp[a]:= 137 ;qqq[a]:= 139 ;Delta[a]:= 60 ;a:= 12 ;nnn[a]:= 4096 ;ppp[a]:= 197 ;qqq[a]:= 199 ;Delta[a]:= 68 ;a:= 13 ;nnn[a]:= 8192 ;ppp[a]:= 283 ;qqq[a]:= 293 ;Delta[a]:= 210 ;a:= 14 ;nnn[a]:= 16384 ;ppp[a]:= 431 ;qqq[a]:= 433 ;Delta[a]:= 416 ;a:= 15 ;nnn[a]:= 32768 ;ppp[a]:= 619 ;qqq[a]:= 631 ;Delta[a]:= 415 ;a:= 16 ;nnn[a]:= 65536 ;ppp[a]:= 919 ;qqq[a]:= 929 ;Delta[a]:= 2 ;a:= 17 ;nnn[a]:= 131072 ;ppp[a]:= 1319 ;qqq[a]:= 1321 ;Delta[a]:= 334 ;a:= 18 ;nnn[a]:= 262144 ;ppp[a]:= 1933 ;qqq[a]:= 1949 ;Delta[a]:= 922 ;a:= 19 ;nnn[a]:= 524288 ;ppp[a]:= 2789 ;qqq[a]:= 2791 ;Delta[a]:= 2628 ;a:= 20 ;nnn[a]:= 1048576 ;ppp[a]:= 4049 ;qqq[a]:= 4051 ;Delta[a]:= 2929 ;a:= 21 ;nnn[a]:= 2097152 ;ppp[a]:= 5851 ;qqq[a]:= 5857 ;Delta[a]:= 5142 ;a:= 22 ;nnn[a]:= 4194304 ;ppp[a]:= 8521 ;qqq[a]:= 8527 ;Delta[a]:= 1847 ;a:= 23 ;nnn[a]:= 8388608 ;ppp[a]:= 12269 ;qqq[a]:= 12277 ;Delta[a]:= 2068 ;a:= 24 ;nnn[a]:= 16777216 ;ppp[a]:= 17657 ;qqq[a]:= 17659 ;Delta[a]:= 12129 ;a:= 25 ;nnn[a]:= 33554432 ;ppp[a]:= 25457 ;qqq[a]:= 25463 ;Delta[a]:= 12848 ;a:= 26 ;nnn[a]:= 67108864 ;ppp[a]:= 36677 ;qqq[a]:= 36683 ;Delta[a]:= 11619 ;a:= 27 ;nnn[a]:= 134217728 ;ppp[a]:= 52807 ;qqq[a]:= 52813 ;Delta[a]:= 46278 ;a:= 28 ;nnn[a]:= 268435456 ;ppp[a]:= 75937 ;qqq[a]:= 75941 ;Delta[a]:= 22242 ;a:= 29 ;nnn[a]:= 536870912 ;ppp[a]:= 109159 ;qqq[a]:= 109169 ;Delta[a]:= 55375 ;a:= 30 ;nnn[a]:= 1073741824 ;ppp[a]:= 156799 ;qqq[a]:= 156817 ;Delta[a]:= 9189 ;a:= 31 ;nnn[a]:= 2147483648 ;ppp[a]:= 225263 ;qqq[a]:= 225287 ;Delta[a]:= 52318 ;a:= 32 ;nnn[a]:= 4294967296 ;ppp[a]:= 323377 ;qqq[a]:= 323381 ;Delta[a]:= 125320 ;a:= 33 ;nnn[a]:= 8589934592 ;ppp[a]:= 464437 ;qqq[a]:= 464447 ;Delta[a]:= 188704 ;a:= 34 ;nnn[a]:= 17179869184 ;ppp[a]:= 665789 ;qqq[a]:= 665801 ;Delta[a]:= 42158 ;a:= 35 ;nnn[a]:= 34359738368 ;ppp[a]:= 954469 ;qqq[a]:= 954491 ;Delta[a]:= 910863 ;a:= 36 ;nnn[a]:= 68719476736 ;ppp[a]:= 1368229 ;qqq[a]:= 1368233 ;Delta[a]:= 1313452 ;a:= 37 ;nnn[a]:= 137438953472 ;ppp[a]:= 1959739 ;qqq[a]:= 1959751 ;Delta[a]:= 56059 ;a:= 38 ;nnn[a]:= 274877906944 ;ppp[a]:= 2807477 ;qqq[a]:= 2807479 ;Delta[a]:= 340473 ;a:= 39 ;nnn[a]:= 549755813888 ;ppp[a]:= 4019809 ;qqq[a]:= 4019831 ;Delta[a]:= 1222039 ;a:= 40 ;nnn[a]:= 1099511627776 ;ppp[a]:= 5753453 ;qqq[a]:= 5753467 ;Delta[a]:= 4612647 ;a:= 41 ;nnn[a]:= 2199023255552 ;ppp[a]:= 8232613 ;qqq[a]:= 8232659 ;Delta[a]:= 5678746 ;a:= 42 ;nnn[a]:= 4398046511104 ;ppp[a]:= 11776829 ;qqq[a]:= 11776859 ;Delta[a]:= 10786702 ;a:= 43 ;nnn[a]:= 8796093022208 ;ppp[a]:= 16843297 ;qqq[a]:= 16843313 ;Delta[a]:= 2478963 ;a:= 44 ;nnn[a]:= 17592186044416 ;ppp[a]:= 24082589 ;qqq[a]:= 24082607 ;Delta[a]:= 2282340 ;a:= 45 ;nnn[a]:= 35184372088832 ;ppp[a]:= 34426313 ;qqq[a]:= 34426321 ;Delta[a]:= 21040805 ;a:= 46 ;nnn[a]:= 70368744177664 ;ppp[a]:= 49198559 ;qqq[a]:= 49198573 ;Delta[a]:= 45448874 ;a:= 47 ;nnn[a]:= 140737488355328 ;ppp[a]:= 70292291 ;qqq[a]:= 70292309 ;Delta[a]:= 29290709 ;a:= 48 ;nnn[a]:= 281474976710656 ;ppp[a]:= 100416713 ;qqq[a]:= 100416721 ;Delta[a]:= 88527836 ;a:= 49 ;nnn[a]:= 562949953421312 ;ppp[a]:= 143413759 ;qqq[a]:= 143413789 ;Delta[a]:= 57216853 ;a:= 50 ;nnn[a]:= 1125899906842624 ;ppp[a]:= 204789517 ;qqq[a]:= 204789523 ;Delta[a]:= 102090071 ;a:= 51 ;nnn[a]:= 2251799813685248 ;ppp[a]:= 292368871 ;qqq[a]:= 292368877 ;Delta[a]:= 267431603 ;a:= 52 ;nnn[a]:= 4503599627370496 ;ppp[a]:= 417329221 ;qqq[a]:= 417329249 ;Delta[a]:= 342115280 ;a:= 53 ;nnn[a]:= 9007199254740992 ;ppp[a]:= 595600669 ;qqq[a]:= 595600673 ;Delta[a]:= 413704712 ;a:= 54 ;nnn[a]:= 18014398509481984 ;ppp[a]:= 849875249 ;qqq[a]:= 849875251 ;Delta[a]:= 803469692 ;a:= 55 ;nnn[a]:= 36028797018963968 ;ppp[a]:= 1212532679 ;qqq[a]:= 1212532681 ;Delta[a]:= 649377128 ;a:= 56 ;nnn[a]:= 72057594037927936 ;ppp[a]:= 1729634189 ;qqq[a]:= 1729634233 ;Delta[a]:= 1133374585 ;a:= 57 ;nnn[a]:= 144115188075855872 ;ppp[a]:= 2466937619 ;qqq[a]:= 2466937621 ;Delta[a]:= 830682139 ;a:= 58 ;nnn[a]:= 288230376151711744 ;ppp[a]:= 3517968721 ;qqq[a]:= 3517968737 ;Delta[a]:= 148045078 ;a:= 59 ;nnn[a]:= 576460752303423488 ;ppp[a]:= 5016124667 ;qqq[a]:= 5016124727 ;Delta[a]:= 3838405139 ;a:= 60 ;nnn[a]:= 1152921504606846976 ;ppp[a]:= 7151293243 ;qqq[a]:= 7151293301 ;Delta[a]:= 2259746163 ;a:= 61 ;nnn[a]:= 2305843009213693952 ;ppp[a]:= 10193990659 ;qqq[a]:= 10193990663 ;Delta[a]:= 5198647657 ;a:= 62 ;nnn[a]:= 4611686018427387904 ;ppp[a]:= 14529489427 ;qqq[a]:= 14529489437 ;Delta[a]:= 10108557436 ;a:= 63 ;nnn[a]:= 9223372036854775808 ;ppp[a]:= 20706267137 ;qqq[a]:= 20706267151 ;Delta[a]:= 3577204465 ;a:= 64 ;nnn[a]:= 18446744073709551616 ;ppp[a]:= 29505444469 ;qqq[a]:= 29505444491 ;Delta[a]:= 16168326584 ;a:= 65 ;nnn[a]:= 36893488147419103232 ;ppp[a]:= 42038976929 ;qqq[a]:= 42038977031 ;Delta[a]:= 33826102762 ;a:= 66 ;nnn[a]:= 73786976294838206464 ;ppp[a]:= 59889909799 ;qqq[a]:= 59889909839 ;Delta[a]:= 21572996461 ;a:= 67 ;nnn[a]:= 147573952589676412928 ;ppp[a]:= 85311533303 ;qqq[a]:= 85311533311 ;Delta[a]:= 26304530075 ;a:= 68 ;nnn[a]:= 295147905179352825856 ;ppp[a]:= 121511451227 ;qqq[a]:= 121511451233 ;Delta[a]:= 91095475794 ;a:= 69 ;nnn[a]:= 590295810358705651712 ;ppp[a]:= 173053904629 ;qqq[a]:= 173053904657 ;Delta[a]:= 130748031024 ;a:= 70 ;nnn[a]:= 1180591620717411303424 ;ppp[a]:= 246435048317 ;qqq[a]:= 246435048341 ;Delta[a]:= 220789910868 ;a:= 71 ;nnn[a]:= 2361183241434822606848 ;ppp[a]:= 350898854197 ;qqq[a]:= 350898854209 ;Delta[a]:= 63149630811 ;a:= 72 ;nnn[a]:= 4722366482869645213696 ;ppp[a]:= 499598330159 ;qqq[a]:= 499598330237 ;Delta[a]:= 453777514684 ;a:= 73 ;nnn[a]:= 9444732965739290427392 ;ppp[a]:= 711246472639 ;qqq[a]:= 711246472667 ;Delta[a]:= 246093977943 ;a:= 74 ;nnn[a]:= 18889465931478580854784 ;ppp[a]:= 1012466288003 ;qqq[a]:= 1012466288017 ;Delta[a]:= 393105146011 ;a:= 75 ;nnn[a]:= 37778931862957161709568 ;ppp[a]:= 1441131375731 ;qqq[a]:= 1441131375737 ;Delta[a]:= 606145831006 ;a:= 76 ;nnn[a]:= 75557863725914323419136 ;ppp[a]:= 2051115029477 ;qqq[a]:= 2051115029497 ;Delta[a]:= 1673698573035 ;a:= 77 ;nnn[a]:= 151115727451828646838272 ;ppp[a]:= 2919047298757 ;qqq[a]:= 2919047298773 ;Delta[a]:= 2804477424618 ;a:= 78 ;nnn[a]:= 302231454903657293676544 ;ppp[a]:= 4153913403373 ;qqq[a]:= 4153913403383 ;Delta[a]:= 787915533962 ;a:= 79 ;nnn[a]:= 604462909807314587353088 ;ppp[a]:= 5910713670031 ;qqq[a]:= 5910713670053 ;Delta[a]:= 1759543655792 ;a:= 80 ;nnn[a]:= 1208925819614629174706176 ;ppp[a]:= 8409873290459 ;qqq[a]:= 8409873290509 ;Delta[a]:= 1976700506884 ;a:= 81 ;nnn[a]:= 2417851639229258349412352 ;ppp[a]:= 11964838507781 ;qqq[a]:= 11964838507793 ;Delta[a]:= 11223571927097 ;a:= 82 ;nnn[a]:= 4835703278458516698824704 ;ppp[a]:= 17021306080813 ;qqq[a]:= 17021306080817 ;Delta[a]:= 330323646390 ;a:= 83 ;nnn[a]:= 9671406556917033397649408 ;ppp[a]:= 24212985402113 ;qqq[a]:= 24212985402139 ;Delta[a]:= 12145573107408 ;a:= 84 ;nnn[a]:= 19342813113834066795298816 ;ppp[a]:= 34440851298773 ;qqq[a]:= 34440851298779 ;Delta[a]:= 11341905109330 ;a:= 85 ;nnn[a]:= 38685626227668133590597632 ;ppp[a]:= 48985800814999 ;qqq[a]:= 48985800815009 ;Delta[a]:= 11062144781712 ;a:= 86 ;nnn[a]:= 77371252455336267181195264 ;ppp[a]:= 69668758192043 ;qqq[a]:= 69668758192063 ;Delta[a]:= 29031345531597 ;a:= 87 ;nnn[a]:= 154742504910672534362390528 ;ppp[a]:= 99078196477573 ;qqq[a]:= 99078196477609 ;Delta[a]:= 37956250152266 ;a:= 88 ;nnn[a]:= 309485009821345068724781056 ;ppp[a]:= 140893471551149 ;qqq[a]:= 140893471551163 ;Delta[a]:= 112612550390062 ;a:= 89 ;nnn[a]:= 618970019642690137449562112 ;ppp[a]:= 200344335550459 ;qqq[a]:= 200344335550501 ;Delta[a]:= 52369382699130 ;a:= 90 ;nnn[a]:= 1237940039285380274899124224 ;ppp[a]:= 284863767904379 ;qqq[a]:= 284863767904427 ;Delta[a]:= 284252006676246 ;a:= 91 ;nnn[a]:= 2475880078570760549798248448 ;ppp[a]:= 405015743379671 ;qqq[a]:= 405015743379701 ;Delta[a]:= 273696384133842 ;a:= 92 ;nnn[a]:= 4951760157141521099596496896 ;ppp[a]:= 575813320142441 ;qqq[a]:= 575813320142443 ;Delta[a]:= 252513156113169 ;a:= 93 ;nnn[a]:= 9903520314283042199192993792 ;ppp[a]:= 818591311846669 ;qqq[a]:= 818591311846697 ;Delta[a]:= 542609494740313 ;a:= 94 ;nnn[a]:= 19807040628566084398385987584 ;ppp[a]:= 1163666945148943 ;qqq[a]:= 1163666945149103 ;Delta[a]:= 848850227750354 ;a:= 95 ;nnn[a]:= 39614081257132168796771975168 ;ppp[a]:= 1654119631980029 ;qqq[a]:= 1654119631980047 ;Delta[a]:= 1417695125364628 ;a:= 96 ;nnn[a]:= 79228162514264337593543950336 ;ppp[a]:= 2351160427227343 ;qqq[a]:= 2351160427227361 ;Delta[a]:= 764139511665826 ;a:= 97 ;nnn[a]:= 158456325028528675187087900672 ;ppp[a]:= 3341759674235299 ;qqq[a]:= 3341759674235317 ;Delta[a]:= 1909623792272991 ;a:= 98 ;nnn[a]:= 316912650057057350374175801344 ;ppp[a]:= 4749481664698139 ;qqq[a]:= 4749481664698163 ;Delta[a]:= 4557540044699731 ;a:= 99 ;nnn[a]:= 633825300114114700748351602688 ;ppp[a]:= 6749874763859599 ;qqq[a]:= 6749874763859603 ;Delta[a]:= 5415322522313957 ;a:= 100 ;nnn[a]:= 1267650600228229401496703205376 ;ppp[a]:= 9592330540937003 ;qqq[a]:= 9592330540937041 ;Delta[a]:= 2077789589606763 ;a:= 101 ;nnn[a]:= 2535301200456458802993406410752 ;ppp[a]:= 13631131151055013 ;qqq[a]:= 13631131151055029 ;Delta[a]:= 2714875794432079 ;a:= 102 ;nnn[a]:= 5070602400912917605986812821504 ;ppp[a]:= 19369544065478017 ;qqq[a]:= 19369544065478081 ;Delta[a]:= 10123908719389466 ;a:= 103 ;nnn[a]:= 10141204801825835211973625643008 ;ppp[a]:= 27522446269314197 ;qqq[a]:= 27522446269314209 ;Delta[a]:= 1548737012611438 ;a:= 104 ;nnn[a]:= 20282409603651670423947251286016 ;ppp[a]:= 39105260093248747 ;qqq[a]:= 39105260093248751 ;Delta[a]:= 36643908758087768 ;a:= 105 ;nnn[a]:= 40564819207303340847894502572032 ;ppp[a]:= 55560252023842141 ;qqq[a]:= 55560252023842147 ;Delta[a]:= 17911529241757605 ;a:= 106 ;nnn[a]:= 81129638414606681695789005144064 ;ppp[a]:= 78935883898354877 ;qqq[a]:= 78935883898354919 ;Delta[a]:= 28686274256363761 ;a:= 107 ;nnn[a]:= 162259276829213363391578010288128 ;ppp[a]:= 112141493087954879 ;qqq[a]:= 112141493087954893 ;Delta[a]:= 57107815856427089 ;a:= 108 ;nnn[a]:= 324518553658426726783156020576256 ;ppp[a]:= 159308929519397491 ;qqq[a]:= 159308929519397641 ;Delta[a]:= 75642919247661640 ;a:= 109 ;nnn[a]:= 649037107316853453566312041152512 ;ppp[a]:= 226306049041980031 ;qqq[a]:= 226306049041980107 ;Delta[a]:= 174406070796167550 ;a:= 110 ;nnn[a]:= 1298074214633706907132624082305024 ;ppp[a]:= 321465805729626037 ;qqq[a]:= 321465805729626079 ;Delta[a]:= 102284797271067708 ;a:= 111 ;nnn[a]:= 2596148429267413814265248164610048 ;ppp[a]:= 456621435656355677 ;qqq[a]:= 456621435656355689 ;Delta[a]:= 154368688220482762 ;a:= 112 ;nnn[a]:= 5192296858534827628530496329220096 ;ppp[a]:= 648576188739732977 ;qqq[a]:= 648576188739733043 ;Delta[a]:= 541008846215755201 ;a:= 113 ;nnn[a]:= 10384593717069655257060992658440192 ;ppp[a]:= 921189957002831419 ;qqq[a]:= 921189957002831447 ;Delta[a]:= 537657706675359459 ;a:= 114 ;nnn[a]:= 20769187434139310514121985316880384 ;ppp[a]:= 1308341646613338367 ;qqq[a]:= 1308341646613338413 ;Delta[a]:= 353615181767568420 ;a:= 115 ;nnn[a]:= 41538374868278621028243970633760768 ;ppp[a]:= 1858134705721126627 ;qqq[a]:= 1858134705721126697 ;Delta[a]:= 1392019486205659815 ;a:= 116 ;nnn[a]:= 83076749736557242056487941267521536 ;ppp[a]:= 2638867319610280727 ;qqq[a]:= 2638867319610280729 ;Delta[a]:= 1164786857674569691 ;a:= 117 ;nnn[a]:= 166153499473114484112975882535043072 ;ppp[a]:= 3747507505182025829 ;qqq[a]:= 3747507505182025961 ;Delta[a]:= 631716657140780641; ####################################################################### # Computation of h(2^a) ####################################################################### hdedeuxpuisa:=proc(amin,amax) local a,nprime,pk,A,n,lis; # computes h(2^a) for amin <= a <=amax; for a from amin to amax do nprime:=Delta[a]; pk:=ppp[a]; print("*** n=2^",a,"p_k=",pk,"nprime=",nprime,"***"); print(); A:=Gpkmh(pk,nprime,500,2); lis:=fact(A[1],pk); # fact is a factoring procedure print("factors of G(p_k,nprime)=",op(lis)); print(A[1],"ecart=",A[2],"type=",A[3]); print();print(); od; end;