3
$\begingroup$

I've written the code that follows:

tails = Function[l, If[ l == {}, {}, Prepend[tails[Drop[l, 1]], l]]]

Te = Composition[AllTrue[PrimeQ], Map[FromDigits], tails, IntegerDigits]

Timing[Select[Range[2, 10^6], Te]]

This takes around 21seconds on my computer, but the equivalent Haskell code takes 1second.

Is there something I can do to improve the performance of the code in Mathematica?

$\endgroup$
  • 4
    $\begingroup$ What's the tl;dr of what this code is supposed to do? In case we don't feel like reverse engineering it... $\endgroup$ – user6014 Dec 11 '17 at 21:56
  • 2
    $\begingroup$ You can get somewhere by replacing tails with NestList[#[[2 ;;]] &, #, Length@# - 1] &. Prepend is generally pretty slow, and should be avoided in functions that get used repeatedly. $\endgroup$ – aardvark2012 Dec 11 '17 at 22:00
9
$\begingroup$

[Not really an answer but too long for a comment.]

Not sure if there is a direct Haskell-like way that can be made quite that fast (I guess it depends on exactly how equivalent is that equivalent Haskell code). Here is an approach that short-circuites when it hits a non-prime.

Te4[n_?PrimeQ] :=
 Module[{dig = Rest[IntegerDigits[n]]}, 
  While[Length[dig] > 0 && PrimeQ[FromDigits[dig]], 
   dig = dig[[2 ;; -1]]];
  Length[dig] == 0]
Te4[_Integer] = False;

Check speed.

Timing[t4 = Select[Range[2, 10^6], Te4];]

(* Out[199]= {1.055633, Null} *)

Compare in speed to original version.

tails = Function[l, If[l == {}, {}, Prepend[tails[Drop[l, 1]], l]]];
Te1 = Composition[AllTrue[PrimeQ], Map[FromDigits], tails, 
   IntegerDigits];
Timing[t1 = Select[Range[2, 10^6], Te1];]

(* Out[207]= {13.95006, Null} *)

Check correctness.

t4 === t1

(* Out[208]= True *)
$\endgroup$
7
$\begingroup$

Updated to use PrimePi per @MichaelE2's helpful comment!

Here's another approach. First, let's get a list of all the primes, and then define a prime vector where 1 indicates the index is a prime:

primeList = Prime[Range[PrimePi[10^7]]]; // AbsoluteTiming

primeVector = ConstantArray[0, 10^7]; //AbsoluteTiming
primeVector[[primeList]] = 1; //AbsoluteTiming

{1.21477, Null}

{0.022988, Null}

{0.018425, Null}

Now, rather than create a predicate to test each element, I will use a predicate that works on lists:

primeModQ[p_List, pow_] := primeVector[[Mod[p, 10^pow]]]

The function primeModQ determines which elements of p are primes modulo 10^pow, returning a list of the same length where 1 indicates the number is a prime, and 0 indicates the number is composite. Finally, I will use primeModQ inside a Fold:

r1 = Fold[
    Pick[#, primeModQ[#, #2], 1]&,
    Pick[Range[10^6], primeVector[[Range[10^6]]], 1],
    Range[5, 1, -1]
]; //AbsoluteTiming

Length[r1]

{0.014613, Null}

1672

The above is in agreement with @DanielLichtblau's answer. Let's try primes up to 10^7:

r7 = Fold[
    Pick[#, primeModQ[#, #2], 1]&,
    Pick[Range[10^7], primeVector[[Range[10^7]]], 1],
    Range[6, 1, -1]
]; //AbsoluteTiming

Length[r7]

{0.273536, Null}

3984

Addendum

A much faster algorithm is to start with single digits instead. Here is a version that does this:

Clear[p]
p[0] = Prime @ Range @ PrimePi[10];
p[n_] := p[n] = Join[
    p[n-1],
    q[p[n-1], n]
]

q[plist_, n_] := With[
    {next = Range[10^n, 9 10^n, 10^n]},
    With[
        {new = Tuples[Unevaluated[Plus[next, plist]]]},
        Pick[new,PrimeQ[new]]
    ]
]

Here is the number of primes less than $10^6$:

Length @ p[5] //AbsoluteTiming

{0.004048, 1672}

And the number of primes less than $10^7$ and $10^8$:

Length @ p[6] //AbsoluteTiming
Length @ p[7] //AbsoluteTiming

{0.012358, 3984}

{0.030985, 8779}

$\endgroup$
  • $\begingroup$ primeList = Table[Prime[n], {n, PrimePi[10^7]}] is a bit faster. $\endgroup$ – Michael E2 Dec 12 '17 at 2:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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