Fast Sieve Implementation

Question

I'm going to write a code to compete in performance with NewPGen. The goal is to sieve numbers that have small factors. Mine is very slow in compare with NewPGen.

Here's my code:

random = RandomPrime[10^1000]; 
list = Table[k random + 1, {k, 1, 2000}];

counter = 2000;

Do[
 Do[
  If[list[[index]] != 0 && Divisible[list[[index]], Prime[p]], 
   list[[index]] = 0; counter--],
  {index, 1, 2000}],
 {p, 1, 10000}]

How can I improve it to sieve more in less time?

Updated

Timing for initial code:

31.125

then I've noticed if I remove list processing codes a whopping 15x speedup would occur!!!

Timing[
 Do[
  Do[
   If[Divisible[n, Prime[p]], list = Delete[list, Position[list, n]];
    counter--],
   {n, list}],
  {p, 1, 10000}]
 ]

1.969 seconds for Sieve up to 10,000th prime

38.11 seconds for Sieve up to 200,000th prime (20 fold increase)

The achieved speedup was for small 2K digits numbers so the only problem with my second implementation is that if my numbers have 100K or more digits then I guess iterating in first loop would makes the code so slow.

Any further optimize?

Updated 2

Timings for 1000 digits numbers and sieving up to 1,000,000th prime

{3.23, Daniel Lichtblau second solution}
{6.87, Daniel Lichtblau first solution}
{48.85, Simoon Woods}
{81.96 , b.gatessucks}
{212.81, myself!}

---

Timings for 100,000 digits numbers and sieving up to 10,000th prime !!!

{7.65, Daniel Lichtblau second solution}
{24.01, Daniel Lichtblau first solution}    
{176.87, Simoon Woods}
{187.92 , b.gatessucks}

---

Timings for 100,000 digits numbers and sieving up to 100,000th prime !

 {5.00, NewPGen} (Does anybody know its sieve method?)
 {148.46, Daniel Lichtblau second solution}

--- edit by DL ---

I hope you will pardon an edit from me (and feel free to change or delete it if you like). I am wondering if your comparison is to what you think. The numbers you generate as k*random+1 are not of the form actually used by NewPGen. The code below looks like it might more plausibly get what you want. Here I generate 2000 values each of a bit more than 100,000 digits.

SeedRandom[111111];
k = 5;
random = RandomInteger[{10^k, 10^k + 10^(k - 1)}];
list = Table[k 10^random + 1, {k, 1, 2000}];

Size check:

Length[IntegerDigits[First[list]]]

(* Out[222]= 103210 *)

I changed the product set to simply prepend 2 rather than multiply an element by 2. No real difference here, just wanted to try it out.

npower = 4;
pprods2 = 
  Apply[Times, 
   Table[Prime[Range[10^j + 1, 10^(j + 1)]], {j, 0, npower}], 2];

pprods2 = Prepend[pprods2, 2];

Here is the test.

Timing[
 res2 = Last[
    Reap[Do[Do[
       If[GCD[n, prod] != 1, Sow[n]; Break[]], {prod, pprods2}], {n, 
       list}]]];]

Length[res2[[1]]]

(* Out[220]= {26.170000, Null}

Out[221]= 1795 *)

It's not 5 seconds but also it's not 150 seconds. Somewhere around the geometric mean, actually. I do not know what NewPGen does but it is possible that it involves aggregating primes into products in a way that is more effective than what I did.

--- end DL edit ---

Answer 1

Might be faster to do batches by checking for nontrivial gcds. Below I group 100 primes into each batch, take their products, and do a gcd computation until finding that there is a divisor. The example sieves a million primes.

random = RandomPrime[10^1000];
list = Table[k random + 1, {k, 1, 2000}];

pprods = Apply[Times, Partition[Prime[Range[10^6]], 10^3], 1];

Timing[
 res = Last[
    Reap[Do[Do[
       If[GCD[n, prod] != 1, Sow[n]; Break[]], {prod, pprods}], {n, 
       list}]]];]

(* Out[47]= {11.216472, Null} *)

Length[res[[1]]]

(* Out[48]= 1926 *)

--- edit ---

We can improve speed another factor of two or so by using fewer blocks and fewer factors in the smaller blocks. Here is a way to create them in "decades" (power of 10 blocks).

pprods2 = 
  Apply[Times, Table[Prime[Range[10^j + 1, 10^(j + 1)]], {j, 0, 5}], 
   2];

In this construction the first product misses the fist prime, 2, so we throw that in.

pprods2[[1]] *= 2;

Here is the computation. I used a different random start, hence the slight difference from above.

Timing[
 res2 = Last[
    Reap[Do[Do[
       If[GCD[n, prod] != 1, Sow[n]; Break[]], {prod, pprods2}], {n, 
       list}]]];]

Length[res2[[1]]]

(* Out[23]= {5.397635, Null} *)

(* Out[24]= 1930 *)

--- end edit ---

Answer 2

Divisible is Listable, so a slight speedup (about a factor of 3) can be gained using something like:

result = Fold[Pick[#, Divisible[#, #2], False] &, list, Prime@Range@200000]

Answer 3

Just a small change and you get a free speedup : rather than looping over 10000 and then calculating Prime[p] each time, do a Prime[Range[10000]] once and for all.

SeedRandom[6];
random = RandomPrime[10^1000];
initialList = Table[k random + 1, {k, 1, 2000}];

list = initialList;
counter = 2000;
{elapsedOP, resultOP} = Do[Do[If[Divisible[n, Prime[p]], list = Delete[list, Position[list, n]];
     counter--], {n, list}], {p, 1, 200000}] // AbsoluteTiming;

elapsedOP
list // Length
(* 61.933193 *)
(* 79 *)

list1 = initialList;
counter = 2000;
{elapsed1, result1} = Do[Do[If[Divisible[n, p], list1 = Delete[list1, Position[list1, n]];
   counter--], {n, list1}], {p, Prime[Range[200000]]}] // AbsoluteTiming;

elapsed1
(* 27.499750 *)

list1 == list
(* True *)