 mersenneforum.org (https://www.mersenneforum.org/index.php)
-   Prime Gap Searches (https://www.mersenneforum.org/forumdisplay.php?f=131)
-   -   Prime Gap Theory (https://www.mersenneforum.org/showthread.php?t=20868)

 robert44444uk 2016-01-17 08:28

Prime Gap Theory

 mart_r 2016-02-05 22:47

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^(z-1)+floor(2^z/3));a=0;
b=a;
w=.5;
forprime(q=3,p,
w=w*(q-1)/q;
t=r%q;
v=q-t;
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]
)
)
}[/CODE]I just noticed, for V1.1, I really should make the sieve more compact. Maybe tomorrow.

 science_man_88 2016-02-05 23:15

[QUOTE=mart_r;425374] (though I can't really promise I'm able to implement all of them:)[/QUOTE]

my problem is although I know of ways in theory you could do the same I don't know what is the fastest a lot of the forsteps with 1 step in the step values could all be turned into for or parfor in newer versions of pari perhaps also the forprimes could be written as parforprimes depending on how many threads and cores the person running it has etc.

[CODE] if(d%2,
d=d+1;
print("set d to "d)
);[/CODE]

can be turned into:

[CODE]d=d+d%2;
print("set d to "d);[/CODE]

edit: which can then be turned into ( and probably even shortened more by printing and assigning value in the same statement):

[CODE]d+=d%2;
print("set d to "d);[/CODE]

a few things like:

[CODE]i=i*e[s,2][/CODE]

or [CODE]floor(2^z/3)[/CODE] can be shortened in length ( number of typed keys)

[CODE]forstep(s=t,z,q,a[s]=0);
forstep(s=v,z,q,b[s]=0)[/CODE]

could be turned into one loop with a few short calculations:

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[/CODE] a faster form of issquarefree ? if so e is not needed for that part and i may not be needed in that part either that may free up some memory. I'll see what else I can think of. doh I see why you did that loop and then checked i you used the factorization of d in the || part so you needed it so it was faster to check if all the values of e[,2] are 1 which can be done without the loop still admittedly but at least I realized the reason calling issquarefree would likely factor the number again costing more time.

 danaj 2016-02-06 03:20

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).

 mart_r 2016-02-07 18:57

@ 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 VBA-program, albeit more user-friendly.

[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^z-1);
b=a;
w=.5;
forprime(q=3,p,
w*=(q-1)/q;
t=r%q;
v=q-t;
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]
)
)
}
[/CODE]

 Antonio 2016-02-07 19:45

[QUOTE=mart_r;425557
[CODE]
if(d%2,
d+=1;
print("set d to "d)
);
[/CODE][/QUOTE]

I don't understand why you force the divisor to be even?
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#

 science_man_88 2016-02-07 19:57

[QUOTE=mart_r;425557]@ 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 VBA-program, albeit more user-friendly.

[/QUOTE]

 mart_r 2016-02-08 20:04

[QUOTE=Antonio;425561]I don't understand why you force the divisor to be even?
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#[/QUOTE]

Doesn't an odd divisor make the search for such gaps dramatically ineffective?
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...:wink:

 Antonio 2016-02-09 18:42

[QUOTE=mart_r;425657]Doesn't an odd divisor make the search for such gaps dramatically ineffective?
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...:wink:[/QUOTE]

It depends on your definition of 'ineffective'.

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 :smile:.

 Antonio 2016-02-10 07:49

Just to add (too late for edit) that the program chose p to be over the range p = 73 to p = 241 (inclusive).

 mart_r 2016-02-10 21:19

@ 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^z-1);
b=a;
w=.5;
forprime(q=3,p,
w*=(1-1/q);
t=r%q;
v=q-t;
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]
)
)
}
[/CODE]

All times are UTC. The time now is 07:38.