20160117, 08:28  #1 
Jun 2003
Oxford, UK
2^{2}·3^{2}·53 Posts 
Prime Gap Theory
This thread captures theory about prime gaps, links to key mathematical papers etc

20160205, 22:47  #2 
Dec 2008
you know...around...
1001100000_{2} Posts 
How difficult is it to find a merit >= x?
While I'm pausing my search for prime gaps, I've put together a small Pari program that can help finding optimal parameters for prime gap searches using primorial numbers.
It compares the actually targetted merit to an "effective merit", meaning how hard it is to find from a statistical point of view. The average number of trials to find one merit > x symmetrical around the selected center number is then exp(corresponding effective merit). Thus, you can see what parameters have a positive (or negative) effect on the numbers you are searching, or the magnitude of the merit you want to find. I hope it is instructive enough. I've only taken about three hours to put it together. Ideas for improvements are welcome (though I can't really promise I'm able to implement all of them:). Code:
print("merit(p [,d [,m [,u ] ] ] ):") print("set center number p#") print("optional d as divisor (default = 2)") print("optional m as multiplier (default = 1)") print("optional u as upper bound for merit calculation (default = 20)") print() print("Example: merit(9973,2*3*5*7,,30) calculates 9973#/210 for merits up to 30") print() print("> Be sure to set allocatemem(x) with x>180*p if p>22000 (resp. x>9*u*p) ! <") print() merit(p,d,m,u)={ if(!m,m=1); if(!d, d=2; print("set d to 2 by default") ); if(d%2, d=d+1; print("set d to "d) ); e=factor(d); i=1; for(s=1,#e[,1], i=i*e[s,2] ); if(i>1  p<e[#e[,1],1], print(d" does not divide "p"#") , if(!u,u=20); print("sieve around "m"*"p"#/"d" ..."); r=1; forprime(q=2,p,r=r*q); r=m*r/d; z=floor(1+u/4*log(r))*2; a=binary(2^(z1)+floor(2^z/3));a[1]=0; b=a; w=.5; forprime(q=3,p, w=w*(q1)/q; t=r%q; v=qt; if(!t,t=t+q); forstep(s=t,z,q,a[s]=0); forstep(s=v,z,q,b[s]=0) ); c=0; n=1; forstep(s=2,z,2, if(2*s>n*log(r), print("merit "n" = effective merit "c/w/log(r)); n=n+1 ); c=c+a[s]+b[s] ) ) } 
20160205, 23:15  #3  
"Forget I exist"
Jul 2009
Dumbassville
2^{6}·131 Posts 
Quote:
Code:
if(d%2, d=d+1; print("set d to "d) ); Code:
d=d+d%2; print("set d to "d); Code:
d+=d%2; print("set d to "d); Code:
i=i*e[s,2] Code:
floor(2^z/3) Code:
forstep(s=t,z,q,a[s]=0); forstep(s=v,z,q,b[s]=0) other than shortening things ( admittedly which I don't do), or parallel code, I still haven't looked at what it's doing ( I bet I couldn't figure it out if I tried). edit: oh and at least one of the forprime loops can be replaced with factorback(primes([2,p])) the prints at the start can be combined together into one with \n in it to symbolize the new lines. anyways I'm not really helping optimization. edit: is Code:
for(s=1,#e[,1],i*=e[s,2]);if(i>1 Last fiddled with by science_man_88 on 20160205 at 23:51 

20160206, 03:20  #4 
"Dana Jacobsen"
Feb 2011
Bangkok, TH
2×3×151 Posts 
I think discussion of what one can do with the script would be more useful than its performance. It runs essentially instantly, so it would only matter if one were running it many thousands of times (which would mean modifications since it currently generates many lines of output per run).

20160207, 18:57  #5 
Dec 2008
you know...around...
2^{5}×19 Posts 
@ science_man_88:
Thanks for your comments. Strangely, last time I checked, "a=a+1" or "a=a*2" was slightly faster than "a+=1" or "a*=2", but I've done yet another check today with a different result. I adopted the shorter version now anyway. Your version for the d%2 case works a bit differently than mine, so I didn't change that. I haven't looked into the "issquarefree" function, but I agree, though it is surely faster than a normal factorisation, it is still slower than a quick check on the powers in the factorization which I already need to look for the biggest prime factor of the divisor d. Also, though there may be more recently implemented functions which are faster, I vote for downward compatibility :) I'm more content with the program now. It can keep track of a certain specified merit value during a predefined loop, so it essentially works like my former VBAprogram, albeit more userfriendly. Code:
print() print("merit(p [,d [,m [,u [,tr ] ] ] ] ):") print(" set center number p#") print(" optional d as divisor (default = 2)") print(" optional m as multiplier (default = 1)") print(" optional u as upper bound for merit calculation (default = 20)") print(" optional tr: track values for a specific merit in file 'merit [tr]'") print() print(" Example: merit(9973,2*3*5*7,,30) calculates 9973#/210 for merits up to 30") print() print("> Be sure to set allocatemem(x) with x>90*p if p>44000 (resp. x>5*u*p) ! <") print() merit(p,d,m,u,tr)={ if(!m,m=1); if(!d, d=2; print("set d to 2 by default") ); if(d%2, d+=1; print("set d to "d) ); e=factor(d); i=1; for(s=1,#e[,1], i*=e[s,2] ); if(i>1  p<e[#e[,1],1], print(d" does not divide "p"#") , if(!u,u=20); print("sieve around "m"*"p"#/"d" ..."); r=1; forprime(q=2,p,r*=q); r*=m/d; z=floor(1+u/4*log(r)); a=binary(2^z1); b=a; w=.5; forprime(q=3,p, w*=(q1)/q; t=r%q; v=qt; t=(t+q*(t%2))/2; if(!t,t=q); v=(v+q*(v%2))/2; forstep(s=t,z,q,a[s]=0); forstep(s=v,z,q,b[s]=0) ); c=0; n=1; for(s=1,z, if(s>n/4*log(r), o=Str("merit "n" = effective merit "c/w/log(r)); print(o); if(tr==n, write("merit "tr".txt",m"*"p"#/"d" : "o,Strchr(13)) ); n+=1 ); c+=a[s]+b[s] ) ) } 
20160207, 19:45  #6  
"Antonio Key"
Sep 2011
UK
3^{2}·59 Posts 
Quote:
If I want the divisor to be, say, 3*5 then you force it to become 15+1 = 16 which will (should) then fail as not dividing p#. Any divisor, d, even 1, should be legitimate, as long as:  1) d has no repeated factors 2) the largest factor of d is less than or equal to p and 3) d < p# 

20160207, 19:57  #7  
"Forget I exist"
Jul 2009
Dumbassville
10000011000000_{2} Posts 
Quote:
Last fiddled with by science_man_88 on 20160207 at 20:50 

20160208, 20:04  #8  
Dec 2008
you know...around...
2^{5}×19 Posts 
Quote:
The sieve also only looks for even numbers (well, that might be easily changed, I guess). Admittedly, to be consistent, I should've just halted the program for an odd input of d and not try to change it automatically. I'm still learning... 

20160209, 18:42  #9  
"Antonio Key"
Sep 2011
UK
3^{2}·59 Posts 
Quote:
To explain I will define the following (just to be clear): for prime, p, and divisor, d, with multiplier, i, relative prime to d we compute C = i * p# / d then: P1 = prev_prime(C) P2 = next_prime(C) gap = P2  P1 As an example, I have a version of Dana's code which searches over a specified digit range of C. If I use this with parameters: i from 1e5 to 2e5 C from 20 to 100 digits d = 24090 (2*3*5*11*73) and record the number of values found with merit > 20 the program finds 216 (from 789096 possible) i.e. 3653 tries per result If I then use d = 12045 (3*5*11*73), all other parameters unchanged, then the program finds 272 (from 1578093 possible) i.e. 5802 tries per result. I will leave you to decide if that is ineffective or not . 

20160210, 07:49  #10 
"Antonio Key"
Sep 2011
UK
3^{2}×59 Posts 
Just to add (too late for edit) that the program chose p to be over the range p = 73 to p = 241 (inclusive).

20160210, 21:19  #11 
Dec 2008
you know...around...
260_{16} Posts 
@ sm88: You surely are concerned about my program :) (and nearly doubled the number of messages in my PM inbox ;)
I've changed a few things today and thereby made it quite a bit faster, especially "l=log(r)/4" outside the output loop helped a lot. Also, fixed the main concern about odd inputs of d. Apart from some cosmetics (maybe addhelp sometime later...), how's it now? Code:
print() print("merit(p [,d [,m [,u [,tv ] ] ] ] ):") print(" set center number p#") print(" optional d as divisor (default = 2)") print(" optional m as multiplier (default = 1)") print(" optional u as upper bound for merit calculation (default = 20)") print(" optional tv: track values for a specific merit in file 'merit [tv]'") print() print(" Example: merit(9973,2*3*5*7,,30) calculates 9973#/210 for merits up to 30") print() print("> Be sure to set allocatemem(x) with x>90*p if p>44000 (resp. x>5*u*p) ! <") print() merit(p,d,m,u,tv)={ if(!m,m=1); if(d<1, d=2; print("set d to 2 by default") ); k=d%2; e=factor(d); i=1; for(s=1,#e[,1], i*=e[s,2] ); if(i>1  p<e[#e[,1],1], print(d" does not divide "p"#") , if(!u,u=20); print("sieve around "m"*"p"#/"d" ..."); r=prod(g=1,primepi(p),prime(g)); r*=m/d; l=log(r)/4; z=floor(1+u*l); a=binary(2^z1); b=a; w=.5; forprime(q=3,p, w*=(11/q); t=r%q; v=qt; if(k, t=(t+q*((t+1)%2)+1)/2; v=(v+q*((v+1)%2)+1)/2 , if(t, t=(t+q*(t%2))/2 , t=q ); v=(v+q*(v%2))/2 ); forstep(s=t,z,q,a[s]=0); forstep(s=v,z,q,b[s]=0) ); c=0; n=1; for(s=1,z, if(s>n*l, o=Str("merit "n" = effective merit "c/(w*log(r))); print(o); if(tv==n, write("merit "tv".txt",m"*"p"#/"d" : "o,Strchr(13)) ); n+=1 ); c+=a[s]+b[s] ) ) } Last fiddled with by mart_r on 20160210 at 21:57 Reason: yet some quick modifications... 
Thread Tools  
Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Basic Number Theory 4: a first look at prime numbers  Nick  Number Theory Discussion Group  6  20161014 19:38 
Before you post your new theory about prime, remember  firejuggler  Math  0  20160711 23:09 
Mersene Prime and Number Theory  Ricie  Miscellaneous Math  24  20090814 15:31 
online tutoring in prime number theory  jasong  Math  3  20050515 04:01 
Prime Theory  clowns789  Miscellaneous Math  5  20040108 17:09 