# Set of integers not divisible by smaller set of primes

Let $p_n$ be the sequence of prime numbers, and $s(x,n)=$ the set of integers less or equal than $x$ that are not divisible by $p_1,\dots,p_n.$ I can define it as follows:

s[x_,n_]:=DeleteCases[Map[If[Total[Table[Map[If[CoprimeQ[Prime[Range[n]],a][[#]]==True,1,0]&,Range[n]],{a,1,x}],{2}][[#]]==n,Range[x][[#]],0]&,Range[x]],0]


But it is inefficient. I am sure there is a more efficient way of doing this.

I believe this is correct, and very fast:

fn[x_Integer, n_Integer] :=
Complement[Range @ x, Join @@ Range[#, x, #]] & @ Prime @ Range @ n


Test:

fn[10000, 1223]

{1, 9929, 9931, 9941, 9949, 9967, 9973}


It seems I am a bit late to return to this problem and Simon Woods already provided a memory optimized approach. His sieve is comparatively fast when the return list is long, that is to say when n is small relative to x, but there is a much more efficient approach when the return list is short, that is when n is large.

The core of my method is this:

fx[x_Integer, n_Integer] /; x < Prime[n + 1]^2 := Prime@Range[n + 1, PrimePi@x] ~Prepend~ 1


Note the condition; this method is not valid for small n values, but that is exactly where Simon's code is superior anyway. However as n increases this method becomes faster, ultimately being instantaneous when the output is {1} which is where Simon's code is slowest:

tbl = Table[f[1*^7, n] // Length // Timing, {n, 5*^4, 7*^5, 5*^4}, {f, {fn2, fx}}];

ListPlot[tbl\[Transpose], PlotLegends -> {"SparseArray", "Prime"},
AxesLabel -> {"Seconds", "Output Length"}, ImageSize -> 500]


Clearly these are complementary methods! Therefore I propose this:

fnHybrid[x_Integer, n_Integer] :=
With[{pp = PrimePi @ x},
If[pp - n < 2 n,
Prime @ Range[n + 1, pp] ~Prepend~ 1,
Module[{y = Range @ x},
(y[[# ;; x ;; #]] = 0) & /@ Prime @ Range @ Min[n, pp];
SparseArray[y]["NonzeroValues"]]]]


The crossover point may need to be tuned for other x values but this is surely the best of both worlds:

• Yes, very, very fast!! :) – martin Aug 30 '14 at 18:46
• @martin Thanks for the Accept. I usually recommend waiting 24 hours first, to let everyone around the world have a chance to answer. However in this case I doubt there is a generally faster method available as this already makes use of highly optimized operations. Compilation to C could be faster still of course, but I don't do that. ;^) – Mr.Wizard Aug 30 '14 at 18:52
• I would usually wait, but I couldn't see anything being much faster, as you say!! – martin Aug 30 '14 at 18:52
• @martin Although I haven't found anything faster than the code above there are memory optimization to be made especially when n is large relative to x, e.g. in the case where the result is {1}. I'll add these later today if I have time. – Mr.Wizard Aug 30 '14 at 20:17

This is competitive with Mr Wizards code and seems faster in some cases:

fn2[x_Integer, n_Integer] := Module[{y = Range @ x},
(y[[# ;; x ;; #]] = 0) & /@ Prime[Range @ Min[n, PrimePi @ x]];
SparseArray[y]["NonzeroValues"]]

AbsoluteTiming[fn[10000, 1223];]
(* {0.004000, Null} *)

AbsoluteTiming[fn2[10000, 1223];]
(* {0.010001, Null} *)

AbsoluteTiming[fn[2000000, 100000];]
(* {0.828047, Null} *)

AbsoluteTiming[fn2[2000000, 100000];]
(* {0.412023, Null} *)

• +1. It becomes about twice as fast as the list grows! – RunnyKine Aug 30 '14 at 20:45
• This is related to the memory optimization I had in mind but you executed it better than I would have. (And Min[n, PrimePi @ x] didn't occur to me.) :-) I found a method complementary to this; see my updated answer. – Mr.Wizard Aug 31 '14 at 9:23
ss[x_, n_] :=  Flatten@Position[CoprimeQ[#, Sequence @@ Prime[Range@n]] & /@ Range@x, True]


We can use a simple sieve to find these numbers in $O(x \log \log x)$ time. I went ahead and compiled my solution to make it as fast as possible.

PrimesUpTo = Compile[{{n, _Integer}},
Block[{S = Range[2, n]},
Do[
If[S[[i]] != 0,
S[[2i+1 ;; -1 ;; i+1]] *= 0;
],
{i, Sqrt[n]}
];
Select[S, Positive]
],
CompilationTarget -> "C",
Parallelization -> True,
RuntimeOptions -> "Speed",
CompilationOptions -> {"InlineCompiledFunctions" -> True, "InlineExternalDefinitions" -> True}
];

F = Compile[{{x, _Integer}, {n, _Integer}},
Block[{S = Range[x], primes = PrimesUpTo[Prime[1223]]},
Do[
If[S[[p]] != 0,
S[[p ;; -1 ;; p]] *= 0;
],
{p, primes}
];
Select[S, Positive]
],
CompilationTarget -> "C",
Parallelization -> True,
RuntimeOptions -> "Speed",
CompilationOptions -> {"InlineCompiledFunctions" -> True, "InlineExternalDefinitions" -> True}
];


Here's timings of all functions so far:

F[10000, 1223] // AbsoluteTiming

(* {0.001991, {1, 9929, 9931, 9941, 9949, 9967, 9973}} *)

fn[10000, 1223] // AbsoluteTiming

(* {0.007860, {1, 9929, 9931, 9941, 9949, 9967, 9973}} *)

fn2[10000, 1223] // AbsoluteTiming

(* {0.004625, {1, 9929, 9931, 9941, 9949, 9967, 9973}} *)


## Edit

I just realized Simon Wood's method is the same as mine, but he uses sparse arrays.

• What is fastCompile? – RunnyKine Aug 31 '14 at 1:37
• I think there must be an n/logn factor to account for the number of primes in the sieving process. – Daniel Lichtblau Aug 31 '14 at 13:48
• Sure, you're right. I guess I was considering n to be constant. – Chip Hurst Aug 31 '14 at 16:24
• @RunnyKine, whoops I copy and pasted PrimesUpTo from my init file and that's where fastCompile is defined. I'll change it. – Chip Hurst Aug 31 '14 at 16:25