How do you efficiently sieve for prime 3/4-tuples?

[Puzzle-Peter]

It took a while, but I finally got a running executable. Thanks!

[Puzzle-Peter]

Originally Posted by Puzzle-Peter

It took a while, but I finally got a running executable. Thanks!

After quite some time I wanted to start a real effort and compared polysieve to polysieve2. Something doesn't quite work out.
I am using 62601*2^36641-1 for s
P0=-6
P1=-3
Q0=0
Q1=0
Q2=2
Q3=1

c=-1, 1, 5

with these parameters, polysieve will produce a file while polysieve 2 exits and tells me that 18841 will divide every candidate for c=5. What is happening here?

[henryzz]

I can't find a version of polysieve2 on my pc that gives that behavior. Can you provide the exact parameters(bound_small_primes etc). Is this using the gmp version?

[Puzzle-Peter]

Originally Posted by henryzz

I can't find a version of polysieve2 on my pc that gives that behavior. Can you provide the exact parameters(bound_small_primes etc). Is this using the gmp version?

It's been a while since I built it, but I believe it is the gmp version. I tried to rebuild but didn't remember about the -lgmp. Anyway this is what I used and what I got:

Params in the declarations:
Denom 1024
SmallPrimes2 256
Sieve_len 16384
Hash_bits3 27


Input and output:

  Sieve P(s)+a*Q(s)+c for multiple c values, with fixed s=k*b^n+d; P,Q is polynom.
  Give k: 62601
  Give b: 2
  Give n: 36641
  Give d: -1
  Give the degree of the P polynom: 1
  Give the 0-th coefficient of P: -6
  Give the 1-th coefficient of P: -3
  Give the degree of the Q polynom: 3
  Give the 0-th coefficient of Q: 0
  Give the 1-th coefficient of Q: 0
  Give the 2-th coefficient of Q: 2
  Give the 3-th coefficient of Q: 1
  Give the number of c values for the sieve: 3
  0-th c value: -1
  1-th c value: 1
  2-th c value: 5
  Give start and end value for 'a' (in billions)! 0 1000
  Give the limit for the wheel sieve (bound_small_primes): 11
  Give the limit for sieving primes (smallPrimes): 1000000000
  Give the limit for sieving primes (maxp): 10000000000
  Give the name of the file to output the numbers! test.txt
  Using primes for wheelsieve up to 11
  initialization of all small primes done, primepi(1000000000)=50847534
  Prime divisors of Q(s) up to 1000000000: 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 
53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 
167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 
277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 
401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 
523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 
653 659 661 673 677 683 691 701 709 719 727 733 739 743 751 757 761 769 773 787 
797 809 811 821 823 827 829 839 853 857 859 863 877 881 883 887 907 911 919 929 
937 941 947 953 967 971 977 983 991 997 1009 1013 1019 1021 1031 1033 1039 1049 
1051 1061 1063 1069 1087 1091 1093 1097 1103 1109 1117 1123 1129 1151 1153 1163 
1171 1181 1187 1193 1201 1213 1217 1223 1229 1231 1237 1249 1259 1277 1279 


[snip, exceeding 10,000 characters]

13421 13441 13451 13457 13463 13469 13477 13487 13499 13513 13523 13537 13553 
13567 13577 13591 13597 13613 13619 13627 13633 13649 13669 13679 13681 13687 
13691 13693 13697 13709 13711 13721 13723 13729 13751 13757 13759 13763 13781 
13789 13799 13807 13829 13831 13841
  Warning: 13841 | P(s)+a*Q(s)+(5) for every a value. Exit.
  25.282u 1.173s 1:16.00 34.8%

[henryzz]

I still can't reproduce this. Could you attach your source code so I can do a comparison. Also what platform/compiler are you using?

I have attached a binary that works for me.

[henryzz]

I managed a fair speedup today by replacing the mulmod with a faster version.
This reduced one run from 130 seconds to 100 seconds.

lli mulmod(lli a, lli b, lli n) {
    lli d, dummy; /* d will get a*b mod c */
    asm ("imulq %3\n\t" /* mul a*b -> rdx:rax */
         "idivq %4\n\t" /* (a*b)/c -> quot in rax remainder in rdx */
         :"=a"(dummy), "=&d"(d) /* output */
         :"a"(a), "rm"(b), "rm"(n) /* input */
         :"cc" /* imulq and idivq can set conditions */
        );
    return d;
  }

Instead of:

lli mulmod(lli a,lli b,lli p){// return (a*b)%p
    lli y=(lli)((double)a*(double)b/p+0.5);
    lli r=a*b-y*p;
    if(r<0)r+=p;
    return r;
}

[Puzzle-Peter]

This is confusing me big time.

I seem to be unable to produce a new executable for polysieve2-gmp. I try

 gcc polysieve2-gmp.c -lgmp

and I get the message

   /usr/bin/ld: cannot find –lgmp

I did successfully compile the non-gmp version however and this time it's the number 53 that supposedly divides all of the +5 candidates. Comparing the code to polysieve.c shows only some very minor differences that should not be responsible for the difference in behavior. I do not understand this at all. For this test I used the code versions from post #119.

Platform is Linux

Redhat enterprise 2.6.18-371.3.1.e15

Compiler is gcc 4.1.2-54

[henryzz]

Try modifying the order of the commands in your compile string. -lgmp might want to be first on your system. It might be worth seeing if adding optimization(-O3) changes the behavior. It would be amusing if there was a bug without optimization which is fixed with optimization.

Your gcc version is 7 years old. I am using 4.8.2. Might be worth updating if you can.

Robert do you have any ideas?

[Puzzle-Peter]

It's getting more confusing every day. I managed to compile the gmp version, had to tell the linker the path to libgmp.a. The resulting executable showed the same behavior as the old one, regardless of any -O flags.

I then compared polysieve to polysieve2 and the only real difference seemed to be that bound_small_primes and smallPrimes are set in one version and user input in the other. So I edited polysieve2 to have them set rather than input. Now it worked, even when I set these variables to exactly the same values that I always entered during test runs. Setting only one of them and having the user input the other was no good, I tried it both ways round.

So I tried that with polysieve2-gmp and it did not give me that "all candidates have a divisor" message, but it produced a completely different list of candidates. I'll plug them into PFGW and see what happens using -l. Polysieve and polysieve2-non-gmp give identical output.

I can't say I understand anything that's happening.

[R. Gerbicz]

Originally Posted by Puzzle-Peter

I
So I tried that with polysieve2-gmp and it did not give me that "all candidates have a divisor" message, but it produced a completely different list of candidates. I'll plug them into PFGW and see what happens using -l. Polysieve and polysieve2-non-gmp give identical output.

I have written the first code (polysieve.c), what was my first impression is that even this line is wrong: "Prime divisors of Q(s) up to 1000000000: 2 3 5 7 (...)"
as Q(s)=s^3+2*s^2 is odd (in fact Q(s)=s^2*(s+2) is the prime factorization, so there is no small prime divisor of Q(s)).

[Puzzle-Peter]

OK, I tried to feed the output of the gmp executable into PFGW and many of the candidates had small factors. Something is wrong here, I just don't know what it is.

Edit:
Will these codes do the trick mentioned in post #113 and following? If so, how do I tell it to do it - just set bound_small_primes to 0?