Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question

4 Answers 4

up vote 15 down vote accepted

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]

enter image description here

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:

enter image description here

share|improve this answer
1  
Yes, very, very fast!! :) –  martin Aug 30 at 18:46
2  
@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 at 18:52
    
I would usually wait, but I couldn't see anything being much faster, as you say!! –  martin Aug 30 at 18:52
1  
@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 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} *)
share|improve this answer
    
+1. It becomes about twice as fast as the list grows! –  RunnyKine Aug 30 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 at 9:23
ss[x_, n_] :=  Flatten@Position[CoprimeQ[#, Sequence @@ Prime[Range@n]] & /@ Range@x, True]
share|improve this answer

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.

share|improve this answer
6  
What is fastCompile? –  RunnyKine Aug 31 at 1:37
1  
I think there must be an n/logn factor to account for the number of primes in the sieving process. –  Daniel Lichtblau Aug 31 at 13:48
    
Sure, you're right. I guess I was considering n to be constant. –  Chip Hurst Aug 31 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 at 16:25

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.