11
$\begingroup$

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.

$\endgroup$

4 Answers 4

16
$\begingroup$

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

$\endgroup$
4
  • 1
    $\begingroup$ Yes, very, very fast!! :) $\endgroup$
    – martin
    Commented Aug 30, 2014 at 18:46
  • 2
    $\begingroup$ @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. ;^) $\endgroup$
    – Mr.Wizard
    Commented Aug 30, 2014 at 18:52
  • $\begingroup$ I would usually wait, but I couldn't see anything being much faster, as you say!! $\endgroup$
    – martin
    Commented Aug 30, 2014 at 18:52
  • 1
    $\begingroup$ @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. $\endgroup$
    – Mr.Wizard
    Commented Aug 30, 2014 at 20:17
16
$\begingroup$

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} *)
$\endgroup$
2
  • $\begingroup$ +1. It becomes about twice as fast as the list grows! $\endgroup$
    – RunnyKine
    Commented Aug 30, 2014 at 20:45
  • $\begingroup$ 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. $\endgroup$
    – Mr.Wizard
    Commented Aug 31, 2014 at 9:23
5
$\begingroup$
ss[x_, n_] :=  Flatten@Position[CoprimeQ[#, Sequence @@ Prime[Range@n]] & /@ Range@x, True]
$\endgroup$
4
$\begingroup$

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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