# Performance Problem

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?

• What's the tl;dr of what this code is supposed to do? In case we don't feel like reverse engineering it... – user6014 Dec 11 '17 at 21:56
• 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. – aardvark2012 Dec 11 '17 at 22:00

[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= {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= {13.95006, Null} *)


Check correctness.

t4 === t1

(* Out= True *)


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

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

Clear[p]
p = Prime @ Range @ PrimePi;
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 //AbsoluteTiming


{0.004048, 1672}

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

Length @ p //AbsoluteTiming
Length @ p //AbsoluteTiming


{0.012358, 3984}

{0.030985, 8779}

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