CRUS-like sieving challenge
 

The challenge: extend [url=http://oeis.org/A181356]A181356[/url]. This will require sieving a large range and then checking for probable primes.

I'm curious as to why the existing numbers are so large. (Or maybe it's just late at night and they are the expected size but I don't realize it.)

They are not quite large. 2^2^n-k is "tremendous". Primes get scarce at that size, especially when you looking for safe or sophie germain primes (that is, you have to find not only one, but two primes of that size).

One liner in pari, no optimization (school grade method):
[CODE]
for(n=1,20,a=2^(2^n); forstep(k=1,a,2,if(ispseudoprime(b=a-k) && ispseudoprime(c=((b-1)/2)), print(n", "k); write("A181356.txt",n", "k", "a", "b", "c); break)))
2, 5
3, 29
4, 269
5, 209
6, 1469
7, 15449
8, 36113
9, 38117
[/CODE]Remark that if I use "isprime()" then at n=8 the waiting time is already over 15 minutes. With "ispseudoprime()" the "bottleneck" is 10 (20 minutes passed, no output yet). For 9 it takes about 1 minute.

I was initially printing all on screen, but the screen got messy so I used the external file. The contents of the output file:


[CODE]2, 5, 16, 11, 5
3, 29, 256, 227, 113
4, 269, 65536, 65267, 32633
5, 209, 4294967296, 4294967087, 2147483543
6, 1469, 18446744073709551616, 18446744073709550147, 9223372036854775073
7, 15449, 340282366920938463463374607431768211456, 340282366920938463463374607431768196007, 170141183460469231731687303715884098003
8, 36113, 115792089237316195423570985008687907853269984665640564039457584007913129639936, 115792089237316195423570985008687907853269984665640564039457584007913129603823, 57896044618658097711785492504343953926634992332820282019728792003956564801911
9, 38117, 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084096, 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006045979, 6703903964971298549787012499102923063739682910296196688861780721860882015036773488400937149083451713845015929093243025426876941405973284973216824503022989

[/CODE]So, these numbers are tremendous...

[QUOTE=CRGreathouse;274075]The challenge: extend [url=http://oeis.org/A181356]A181356[/url]. This will require sieving a large range and then checking for probable primes.

I'm curious as to why the existing numbers are so large. (Or maybe it's just late at night and they are the expected size but I don't realize it.)[/QUOTE]

They're safe-primes, which are roughly probability 1/(log N)^2; so you'd expect the Nth entry to be of order 4^N, and it is.

Implying you'd need to check somewhere around a billion numbers of size 2^32768 to find the next one, which is really quite a lot of work (probably six CPU-months assuming you get 99.9% removal by sieving, which is realistic for safeprimes)

[QUOTE=fivemack;274104]They're safe-primes, which are roughly probability 1/(log N)^2; so you'd expect the Nth entry to be of order 4^N, and it is.[/QUOTE]

...and that's why I should leave math for the pre-midnight hours. :smile:

I believe the answer is k=718982153. :smile:

2^32768-718982153 is PRP
2^32767-359491077 is PRP

Attempted primality test using the -tc option in PFGW:
[code]
n=32768:
[Oddly several different recent versions of PFGW errored out when
attempting the -tc option on the n=32768 exponent.] Details:

Primality testing 2^32768-718982153 [N-1/N+1, Brillhart-Lehmer-Selfridge]
Running N-1 test using base 5
Error 1002 initializing FFT code

n=32767:
Primality testing 2^32767-359491077 [N-1/N+1, Brillhart-Lehmer-Selfridge]
Running N-1 test using base 2
Running N+1 test using discriminant 7, base 2+sqrt(7)
Calling N+1 BLS with factored part 0.07% and helper 0.04% (0.25% proof)
2^32767-359491077 is Fermat and Lucas PRP! (61.1761s+0.0023s)
[/code]

Some details:
This took ~4-6 CPU weeks using unmodified versions of srsieve for sieving and PFGW for PRP testing. All k<1.2G was tested. I sieved n=32768, converted the file, then sieved n=32767, both n-values to an approximate optimal depth of P=50M. There were ~798,000 tests to run after sieving both sides with an estimated ~1/720 chance of each test being PRP for either side. There were 1138 PRPs on the n=32768 side with also a ~1/720 chance of those being also PRP for n=32767...and one was found. The effective chance of a safe prime (or Sophie-Germain prime) for this range was ~75-80% so I very easily could have found nothing or even 2 or more safe primes.

Entries for the PRPs in the factoring DB are at:
[URL]http://factorization.ath.cx/index.php?query=2%5E32768-718982153[/URL]
[URL]http://factorization.ath.cx/index.php?query=2%5E32767-359491077[/URL]


Does the OEIS site need primality proof? The chance of either one of the PRPs being composite is very tiny. I don't care to run Primo tests on both a 9864 and a 9865-digit PRP. The CPU time needed would be more than finding the two PRPs. If they need proof, I'd be glad to share credit with anyone who would run the tests.

BTW, this sieving was nothing like the sieving that we do at CRUS, although the title of the thread is what caught my attention.

Now...does anyone care to try for the next term at n=65535/65536? :smile:


Gary

For the next one, PRPtop gives us a headstart [COLOR=darkred]...not[/COLOR]. :smile:

[URL]http://www.primenumbers.net/prptop/searchform.php?form=2%5E65536-x[/URL] - these are known PRPs and their counteparts are not PRPs

[QUOTE=Batalov;275843]For the next one, PRPtop gives us a headstart [COLOR=darkred]...not[/COLOR]. :smile:

[URL]http://www.primenumbers.net/prptop/searchform.php?form=2%5E65536-x[/URL] - these are known PRPs and their counteparts are not PRPs[/QUOTE]

Yeah, unfortunately their counterparts all have factors of 2 or 3 so they don't help us any, i.e.:
2^65535-2814 has factors: 2
2^65535-49634 has factors: 2
2^65535-61782 has factors: 2
2^65535-86009 has factors: 3
2^65535-134142 has factors: 2
2^65535-212768 has factors: 2^5
2^65535-241025 has factors: 3
2^65535-254177 has factors: 3^2

As implied by the title, the main challenge of this effort is sieving both sides to a reasonable limit. It's easy enough to find many PRPs on one side. The challenge is combining the two sides.

[QUOTE=CRGreathouse;274075]The challenge: extend [url=http://oeis.org/A181356]A181356[/url]. This will require sieving a large range and then checking for probable primes.

I'm curious as to why the existing numbers are so large. (Or maybe it's just late at night and they are the expected size but I don't realize it.)[/QUOTE]

I'm trying to figure it out a bit ( reading a bit):

[url]http://en.wikipedia.org/wiki/Safe_prime[/url]

[QUOTE]Combining both forms using lcm(6,4) we determine that a safe prime q > 7 also must be of the form 12k−1 or, equivalently, q ≡ 11 (mod 12).[/QUOTE]

PARI

[QUOTE](17:33)>for(n=2,20,print((2^(2^n))%12))
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4[/QUOTE]

so 4-(-1) = 5 so [TEX]k \equiv 5 \text { mod } 12[/TEX] that will help me personally.

now if only I could find something else to speed it up.

[QUOTE=science_man_88;275851]I'm trying to figure it out a bit ( reading a bit):

[url]http://en.wikipedia.org/wiki/Safe_prime[/url]



PARI



so 4-(-1) = 5 so [TEX]k \equiv 5 \text { mod } 12[/TEX] that will help me personally.

now if only I could find something else to speed it up.[/QUOTE]

and with that info I find different values:

[QUOTE](17:34)>for(n=2,15,forstep(k=5,2^(2^n-1)-1,12,if(isprime(((2^(2^n)-k)-1)/2),print(k);break())))
5
29
101
137
329
77
1469
677[/QUOTE]

[QUOTE=science_man_88;275852]and with that info I find different values:[/QUOTE]

nevermind found my flaw forgot to check if the first part was prime.

[QUOTE=science_man_88;275854]nevermind found my flaw forgot to check if the first part was prime.[/QUOTE]

[QUOTE]for(n=2,15,forstep(k=5,2^(2^n-1)-1,12,if(isprime(2^(2^n)-k) && isprime(((2^(2^n)-k)-1)/2),print(k);break())))[/QUOTE]

I'm guessing using stuff from the other pari code given might help.

[QUOTE=LaurV;274089]They are not quite large. 2^2^n-k is "tremendous". Primes get scarce at that size, especially when you looking for safe or sophie germain primes (that is, you have to find not only one, but two primes of that size).

One liner in pari, no optimization (school grade method):
[CODE]
for(n=1,20,a=2^(2^n); forstep(k=1,a,2,if(ispseudoprime(b=a-k) && ispseudoprime(c=((b-1)/2)), print(n", "k); write("A181356.txt",n", "k", "a", "b", "c); break)))
2, 5
3, 29
4, 269
5, 209
6, 1469
7, 15449
8, 36113
9, 38117
[/CODE]Remark that if I use "isprime()" then at n=8 the waiting time is already over 15 minutes. With "ispseudoprime()" the "bottleneck" is 10 (20 minutes passed, no output yet). For 9 it takes about 1 minute.

I was initially printing all on screen, but the screen got messy so I used the external file. The contents of the output file:


[CODE]2, 5, 16, 11, 5
3, 29, 256, 227, 113
4, 269, 65536, 65267, 32633
5, 209, 4294967296, 4294967087, 2147483543
6, 1469, 18446744073709551616, 18446744073709550147, 9223372036854775073
7, 15449, 340282366920938463463374607431768211456, 340282366920938463463374607431768196007, 170141183460469231731687303715884098003
8, 36113, 115792089237316195423570985008687907853269984665640564039457584007913129639936, 115792089237316195423570985008687907853269984665640564039457584007913129603823, 57896044618658097711785492504343953926634992332820282019728792003956564801911
9, 38117, 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084096, 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006045979, 6703903964971298549787012499102923063739682910296196688861780721860882015036773488400937149083451713845015929093243025426876941405973284973216824503022989

[/CODE]So, these numbers are tremendous...[/QUOTE]

using a version of this with a step of 12 and start of 5 I get to try line 10 in under 17 seconds.

SM88, please stop responding to yourself. It's annoying. Take your time, figure everything out ahead of time, and then make one post.

I don't know if this will help you any but for the n=32768/32767 challenge, for both sides to be prime, for n=32768, the condition k==(17, 29, or 53 mod 60) had to be true. What that condition did was eliminate all factors of 2, 3, and 5 from both n=32767 and 32768. That was the starting point that I used for sieving n=32768. For sieving n=65536, although the modulos will be different, the idea will be the same.

[QUOTE=gd_barnes;275860]SM88, please stop responding to yourself. It's annoying. Take your time, figure everything out ahead of time, and then make one post.

I don't know if this will help you any but for the n=32768/32767 challenge, for both sides to be prime, for n=32768, the condition k==(17, 29, or 53 mod 60) had to be true. What that condition did was eliminate all factors of 2, 3, and 5 from both n=32767 and 32768. That was the starting point that I used for sieving n=32768. For sieving n=65536, although the modulos will be different, the idea will be the same.[/QUOTE]

okay,sorry, it's probably just that I never catch everything so I'd never post under those conditions.

[QUOTE=gd_barnes;275860]...for n=32768, the condition k==(17, 29, or 53 mod 60) had to be true...[/QUOTE]
For sieving n=65536, the condition is the same.

For an estimate how large the k[SUB]16[/SUB] might be, here's the first PRP from 65536 side after the 2-side sieve (but not the other): 2^65536-2723249.

So, that's one PRP down, with a few thousand more to find and test on the other side until both are PRPs.

P.S. oh, here came another: 2^65536-2912969...

[QUOTE=Batalov;275903]For an estimate how large the k[SUB]16[/SUB] might be, here's the first PRP from 65536 side after the 2-side sieve (but not the other): 2^65536-2723249.

So, that's one PRP down, with a few thousand more to find and test on the other side until both are PRPs.

P.S. oh, here came another: 2^65536-2912969...[/QUOTE]

Using my time of 4-6 CPU weeks for the n=32768 testing, it would take 16 times as long to do the n=65536 testing, as follows:

1. The tests take 4 times as long.
2. Half as much chance of n=65536 (vs. n=32768) being prime so twice as many tests.
3. Half as much chance of n=65535 (vs. n=32767) being prime so twice as many tests.

4x2x2=16 times as long and 4 times the k-range.

So you can assume that it would take 64-96 CPU weeks or 1-2 CPU years with current sieving/testing technology to find an n=65536 safe prime using srsieve/PFGW like I did.

To have a 75-80% chance of safe prime at n=65536, you would want to sieve all k<~5G.

The key is for someone to come up with a better sieving program. P=50M was the approximate and very low optimum sieving limit for n=32768 because srsieve has to handle each k as a separate sequence for a single n and so is very slow. As a general rule, it can only effectively sieve ~100,000 sequences reasonably efficiently. Beyond that and it loses a lot of efficiency.

Personally I had to start out sieving in k=2M groups starting with all k==(17, 29 or 53 mod 60). 3 out of every 60 k's meant 100,000 k's in each k=2M group. I had to do that 600 times for the entire k<=1.2G range, i.e. 2M*600=1.2G. Once I had the # of sequences reduced after sieving the first go around, I could then sieve the other side in 120 k=10M groups. (Actually there were more sieving steps than that to increase CPU sieving efficiency but they are a moot point for this explanation. Also, this is definitely the limit of manual effort that anyone would want to do for such an excercise. It is very impractical for n=65536; requiring 4 times the # of sieving groups.)


Gary

[QUOTE=gd_barnes;275839]
Does the OEIS site need primality proof? The chance of either one of the PRPs being composite is very tiny. I don't care to run Primo tests on both a 9864 and a 9865-digit PRP. The CPU time needed would be more than finding the two PRPs. If they need proof, I'd be glad to share credit with anyone who would run the tests.[/QUOTE]

Can anyone give me an answer to the above? If PRPs are good enough for them, I'll go ahead and submit them to the OEIS site.

Usually, yes (unlike UTM Primes). With a clear note that they are PRP.

[QUOTE=Batalov;275912]Usually, yes (unlike UTM Primes). With a clear note that they are PRP.[/QUOTE]

I have submitted the proposed change. This is my first submission to OEIS. We'll see if it is accepted.

That was kind of fun. :smile:

[QUOTE=gd_barnes;275839]I don't care to run Primo tests on both a 9864 and a 9865-digit PRP.[/QUOTE]

You only need Primo on the smaller one. You can then use it in an N-1 proof for the larger one.

[QUOTE=gd_barnes;275910]The key is for someone to come up with a better sieving program. P=50M was the approximate and very low optimum sieving limit for n=32768 because srsieve has to handle each k as a separate sequence for a single n and so is very slow. As a general rule, it can only effectively sieve ~100,000 sequences reasonably efficiently. Beyond that and it loses a lot of efficiency...[/QUOTE]
Well, 2^(2^n)-k sieve in k should be an appallingly easy sieve to write from scratch and even without GMP (raising 2^(2^n) mod p is very easy and then just zero some bits), but will it do better than a "sieve" like this?

[CODE]# gp -p 80000000
write("kcan","ABC 2^65536-$a")
forstep(k=17,99999999,12, \
if(k%60==17||k%60==29||k%60==53, \
m=1;forprime(p=7,80000000,if(lift(Mod(2,p)^65536-k)<2,m=0;break)); \
if(m>0,write("kcan",k))) \
)
[/CODE]
and after a bit of time start "pfgw -f0 -llog kcan" in parallel?

[QUOTE=Batalov;276008]Well, 2^(2^n)-k sieve in k should be an appallingly easy sieve to write from scratch and even without GMP (raising 2^(2^n) mod p is very easy and then just zero some bits), but will it do better than a "sieve" like this?

[CODE]# gp -p 80000000
write("kcan","ABC 2^65536-$a")
forstep(k=17,99999999,12, \
if(k%60==17||k%60==29||k%60==53, \
m=1;forprime(p=7,80000000,if(lift(Mod(2,p)^65536-k)<2,m=0;break)); \
if(m>0,write("kcan",k))) \
)
[/CODE]
and after a bit of time start "pfgw -f0 -llog kcan" in parallel?[/QUOTE]

Is this PARI code? I don't have PARI and am not familiar with it.

Yes, it is. It is an invaluable tool and is easy to install (not as easy to use, but that's a different story; however, there are experts here on forum who will help). But this is only for prototyping.

I am sure that a plain C 50-liner siever will take care of the problem much easier and faster. I'll probably give it a crack over weekend. One could use the k%60==17||k%60==29||k%60==53 packing in it (or even tighter), or rather I am inclined to just keep only 1 (mod 12) 2Gb byte array (for k<24G) and then clear proper, equally spaced bytes relentlessly for each p until 2^32 or until too bored. :truck:

Quick'n'dirty... and lazy (no bit packing or any fancy stuff) and limited to sieving to 2B. Attached. Can be vastly improved (but of course Geoff could write a totally superior sieve with his left foot in 15 minutes).