# Fast Sieve Implementation

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

-
 @DanielLichtblau, the 100,000 digits number I've used for test was of primorial type k*p#+1 not the random one I have used in smaller tests. Actually, my numbers are of this type. – Mohsen Afshin Jan 3 at 18:14

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

-
 I was exactly looking for this solution that checks bunch of primes in a single check.If you look at my Updated Timings, you'll see how fast your solution is. Thanks a lot @DanielLichtblau. – Mohsen Afshin Jan 2 at 17:15

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]

-
 Another Impressive answer, timing is now (9.625 seconds). Mathematica is really the knowledge of KNOWING .... – Mohsen Afshin Dec 31 '12 at 12:55 Can you give a hint on this question (mathematica.stackexchange.com/questions/15400/…), too? – Mohsen Afshin Dec 31 '12 at 13:46

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 *)

-
 Impressive, my timings in now half (15.8 seconds). It seems successive search for primes does really matter. Thanks – Mohsen Afshin Dec 31 '12 at 12:48