Does Mathematica have a twin prime equivalent of PrimePi?

Well, it's all there in the title! I'd like to be able to plot the number of twin primes =<x. Is there an inbuilt function that can do this?

• Is an exact value required? If not, one can take a "probabilistic" approach using log integral. If one wants to get a rough idea of how such a plot would look, this might suffice. Obviously it is too heuristic to give finer structure or anything like that. Main advantage is it is cheap and can be computed over arbitrarily large ranges (so not limited in the way that Prime-based methods would be). Jul 28 '18 at 14:52
• Wolfram Challenges Mar 10 '19 at 11:25

Here's a very short version:

primePairs[x_] := With[{primes = Prime[Range[PrimePi[x]]]},
Intersection[primes, primes - 2]]


It's returns the same numbers as @KennyColnago's TwinPrimeLesser, but it's a bit faster:

primePairs[10^8]; // AbsoluteTiming


{9.16803, Null}

TwinPrimeLesser[10^8]; // AbsoluteTiming


{11.8908, Null}

If you need all the primes that are twin primes up to n then.

TwinPrimesUpTo[n_] := Select[
Most@NestWhileList[
NextPrime,
Prime[1], # <= n &
]
, (PrimeQ[# + 2]||PrimeQ[# + 2]) &]


It was pointed out by @KennyColnago in other answer that it's customary to count the number of pairs of twin primes, in which case a list of the lesser of the twin pairs is given by

TwinPrimesLesserUpTo[n_] := Select[
Most@NestWhileList[
NextPrime,
Prime[1], # <= n &
]
, (PrimeQ[# + 2]) &]


This gives the same output as @KennyColnago's answer, but his performance is better (+1), as it takes advantage of PrimePi and Range.

TwinPrimesLesserUpTo[99999] == TwinPrimeLesser[99999]
(* True *)


ListPlot[
TwinPrimesUpTo[1000]
, Filling -> Axis
, PlotTheme -> "Scientific"
, FrameLabel -> {"Index", "Value"}
, PlotLabel -> "Twin Primes"
]


ListPlot[
Prepend[MapIndexed[Flatten[{##}] &, TwinPrimesUpTo[1000]], {0, 0}]
, InterpolationOrder -> 0
, Joined -> True
, PlotTheme -> "Scientific"
, FrameLabel -> {"n", "Number of Twin primes"}
, PlotLabel -> "Number of twin primes up to n"
]


One convention (the usual?) for counting twin primes is to count just 1 for each pair. The original question is ambiguous, but it seems that @rhermans (+1) has a different convention of counting 2 for each pair.

The following definition returns the lesser of all twin prime pairs $\{p,p+2\}$ such that $p \le x$. Make a list of primes up to $x$, then test for prime $x+2$.

TwinPrimeLesser[x_] := Pick[#, PrimeQ[# + 2]] &[Prime[Range[PrimePi[x]]]]


The function PrimePi2[x] counts 1 for each twin-prime pair $\{p,p+2\}$ with $p \le x$.

PrimePi2[x_] := Length[Pick[#, PrimeQ[# + 2]]] &[Prime[Range[PrimePi[x]]]]


The twin prime count is fairly fast.

AbsoluteTiming[PrimePi2[10^7]]


{1.38108, 58980}

As @DanielLichtblau mentioned in a comment, an approximate count is given by a log integral formula (equation 5) on MathWorld.

The integral has an explicit form for $x>2$.

LogIntegral2[x_] :=
2.42969427374664261373628276610088274856979843439 -
1.320323631693739147855624220029111556865246*x / Log[x] +


AbsoluteTiming[Round[LogIntegral2[10^7]]]


{0.000634, 58754}

As for plotting the number of twin primes, don't forget ListStepPlot.

TwinPrimePlot[x_] :=
ListStepPlot[
{Transpose[{#, Range[Length[#]]}] &[TwinPrimeLesser[x]],
Table[{n, LogIntegral2[n]}, {n, Range[3, x, x/100]}]},
Frame -> True, FrameLabel -> {"x", "Number of Twin Primes"},
PlotLabel -> "Number of Twin Primes {p,p+2} With p<=x",
ImageSize -> 500, BaseStyle -> {FontSize -> 14},
PlotRange -> {{0, Automatic}, {0, Automatic}},
PlotLegends ->
]


A way to make things faster is to build our own sieve. This still requires a large time and space complexity though.

Here's a compiled sieve of Eratosthenes:

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


Now write a function to pick out consecutive primes and one that counts them:

TwinPrimeList[n_] :=
With[{primes = PrimesSieve[n+1]},
Pick[Most[primes], Differences[primes], 2]
]

PrimePi2[n_] := Length[TwinPrimeList[n]]


Try a large value:

PrimePi2[10^8] // AbsoluteTiming

{4.361947, 440312}


We can get even faster with a wheel sieve

WheelSieve = Compile[{{n, _Integer}},
Module[{k, res, loc, pp, S, i = 1, j = 2, x, p = 0, sqrt = Floor[Sqrt[n]], max, a = 0, mod},
k = #1;
res = #2;
loc = #3;
pp = #4;

S = Table[If[j < 0, pp[[loc[[r]]]], With[{m = k j + r}, Boole[m <= n]m]], {j, -1, Ceiling[n/#]-1}, {r, res}];
max = Max[S];

x = Length[res];

While[p < sqrt,
If[(p = S[[j, i++]]) != 0,
a = 0;
While[a >= 0,
Do[
If[m > max, a = -1; Break[]];
mod = Mod[m, k, 2];
S[[Quotient[m, k]+1+Boole[mod < k], loc[[mod]]]] = 0,
{m, p*res + p*k*(a++)}
]
]
];
If[i > x, i = 1; j++]
];

Select[Flatten[S], Positive]
],
CompilationTarget -> "C",
RuntimeOptions -> "Speed",
CompilationOptions -> {"ExpressionOptimization" -> True}
]&[
P = Times @@ Prime[Range[6]],
Select[Range[2, P+1], CoprimeQ[#, P]&],
ReplacePart[ConstantArray[0, P+1], MapIndexed[#1 -> First[#2]&, Select[Range[2, P+1], CoprimeQ[#, P]&]]],
];

TwinPrimeListFast[n_] :=
With[{primes = WheelSieve[n+1]},
Pick[Most[primes], Differences[primes], 2]
]

PrimePi2Fast[n_] := Length[TwinPrimeListFast[n]]


This speeds things up by a factor of 3:

PrimePi2Fast[10^8] // AbsoluteTiming

{1.348864, 440312}

• Nice. +1, but what and where is fastCompile? Aug 1 '18 at 20:28
• @KennyColnago Sorry, fastCompile is my own version of Compile that has all of those fast options turned on by default. I forgot to change the name when I copied over to here. Aug 1 '18 at 20:30
    j = 100000;
myTwinPrime = {};
For[k = 1, Prime[k] < j, k++,
If[Prime[k + 1] == Prime[k] + 2,
myTwinPrime = Append[myTwinPrime, {Prime[k], Prime[k + 1]}]]];
myAccumulatedTwinPrimes = Accumulate[myTwinPrimes];


then

myAccumulatedTwinPrimes[[69]]


or

ListPlot[myAccumulatedTwinPrimes]

• There is a slight syntax problem. Accumulate[myTwinPrimes] should be Accumulate[myTwinPrime]. When I run it. Jun 19 '21 at 1:08