Prime Gap Theory

[robert44444uk]

This thread captures theory about prime gaps, links to key mathematical papers etc

[mart_r]

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

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[1]=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]
        )
    )
}

I just noticed, for V1.1, I really should make the sieve more compact. Maybe tomorrow.

[science_man_88]

Originally Posted by mart_r

(though I can't really promise I'm able to implement all of them:)

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.

 if(d%2,
        d=d+1;
        print("set d to "d)
    );

can be turned into:

d=d+d%2;
print("set d to "d);

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

d+=d%2;
print("set d to "d);

a few things like:

i=i*e[s,2]

or

floor(2^z/3)

can be shortened in length ( number of typed keys)

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

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

for(s=1,#e[,1],i*=e[s,2]);if(i>1

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]

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]

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

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]
        )
    )
}

[Antonio]

Originally Posted by mart_r;425557
[CODE

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

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]

Originally Posted by mart_r

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

the d=d+d%2 thing is that if d%2 is non-zero (what if looks for at last check, in this case equivalent to if d%2==1) then it adds 1 just like your program does, and adds 0 otherwise ( effectively what your program does by doing nothing, mine doesn't miss that print either way but that's about the difference). incorporating it into the print is just my way of doing two things at once. one thought I had is if checking wasn't potentially going to waste time with squarefree numbers is at very least you could cut the loop short once e[s,2] isn't 1 therefore stopping the checks early in more than 50% of cases as the resulting d that got factored was even as odd+odd ==even+even == even and half of even numbers are divisible by 4 so half of the numbers tested could be stopped at the first check and have i>1 when it's checked by the if statement. I see a potential way to work parts of the two if statements about d together but it may be slower than what you already have. the for loop search through e could also be made into a parallel for loop. in the version I made ( though I'd have to remake it I didn't save it I sped things up by assigning log(r) it's own variable ( okay technically constant or name) and using that from then on. in fact anything not repeatedly changing and used enough times can get that treatment and speed things up. anyways I'm probably repeating myself on some things. edit: oh and dividing by two things is like dividing by the two multiplied. etc. it might be best to take the code conversations to the PARI commands thread there's a lot I could randomly say about things. in fact this brings up the fact that if d=3 mod 4 adding one will automatically make it fail your check.

[mart_r]

Originally Posted by Antonio

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#

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

[Antonio]

Originally Posted by mart_r

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

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 .

[Antonio]

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]

@ 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?

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]
        )
    )
}