mersenneforum.org  

Go Back   mersenneforum.org > Prime Search Projects > Prime Gap Searches

Reply
 
Thread Tools
Old 2016-01-17, 08:28   #1
robert44444uk
 
robert44444uk's Avatar
 
Jun 2003
Oxford, UK

2·13·73 Posts
Default Prime Gap Theory

This thread captures theory about prime gaps, links to key mathematical papers etc
robert44444uk is offline   Reply With Quote
Old 2016-02-05, 22:47   #2
mart_r
 
mart_r's Avatar
 
Dec 2008
you know...around...

25716 Posts
Default 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[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.
mart_r is offline   Reply With Quote
Old 2016-02-05, 23:15   #3
science_man_88
 
science_man_88's Avatar
 
"Forget I exist"
Jul 2009
Dumbassville

8,369 Posts
Default

Quote:
Originally Posted by mart_r View Post
(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.

Code:
 if(d%2,
        d=d+1;
        print("set d to "d)
    );
can be turned into:

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

Code:
d+=d%2;
print("set d to "d);
a few things like:

Code:
i=i*e[s,2]
or
Code:
floor(2^z/3)
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)
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
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.

Last fiddled with by science_man_88 on 2016-02-05 at 23:51
science_man_88 is offline   Reply With Quote
Old 2016-02-06, 03:20   #4
danaj
 
"Dana Jacobsen"
Feb 2011
Bangkok, TH

90610 Posts
Default

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).
danaj is offline   Reply With Quote
Old 2016-02-07, 18:57   #5
mart_r
 
mart_r's Avatar
 
Dec 2008
you know...around...

11278 Posts
Default

@ 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]
        )
    )
}
mart_r is offline   Reply With Quote
Old 2016-02-07, 19:45   #6
Antonio
 
Antonio's Avatar
 
"Antonio Key"
Sep 2011
UK

32×59 Posts
Default

Quote:
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#
Antonio is offline   Reply With Quote
Old 2016-02-07, 19:57   #7
science_man_88
 
science_man_88's Avatar
 
"Forget I exist"
Jul 2009
Dumbassville

8,369 Posts
Default

Quote:
Originally Posted by mart_r View Post
@ 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.

Last fiddled with by science_man_88 on 2016-02-07 at 20:50
science_man_88 is offline   Reply With Quote
Old 2016-02-08, 20:04   #8
mart_r
 
mart_r's Avatar
 
Dec 2008
you know...around...

11278 Posts
Default

Quote:
Originally Posted by Antonio View Post
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...
mart_r is offline   Reply With Quote
Old 2016-02-09, 18:42   #9
Antonio
 
Antonio's Avatar
 
"Antonio Key"
Sep 2011
UK

32·59 Posts
Default

Quote:
Originally Posted by mart_r View Post
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 is offline   Reply With Quote
Old 2016-02-10, 07:49   #10
Antonio
 
Antonio's Avatar
 
"Antonio Key"
Sep 2011
UK

10238 Posts
Default

Just to add (too late for edit) that the program chose p to be over the range p = 73 to p = 241 (inclusive).
Antonio is offline   Reply With Quote
Old 2016-02-10, 21:19   #11
mart_r
 
mart_r's Avatar
 
Dec 2008
you know...around...

599 Posts
Default

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

Last fiddled with by mart_r on 2016-02-10 at 21:57 Reason: yet some quick modifications...
mart_r is offline   Reply With Quote
Reply

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 2016-10-14 19:38
Before you post your new theory about prime, remember firejuggler Math 0 2016-07-11 23:09
Mersene Prime and Number Theory Ricie Miscellaneous Math 24 2009-08-14 15:31
online tutoring in prime number theory jasong Math 3 2005-05-15 04:01
Prime Theory clowns789 Miscellaneous Math 5 2004-01-08 17:09

All times are UTC. The time now is 12:26.

Wed Oct 28 12:26:55 UTC 2020 up 48 days, 9:37, 1 user, load averages: 2.13, 1.76, 1.67

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2020, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.