Does Mathematica have a twin prime equivalent of `PrimePi`?
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 *)
The OP asked to "plot"
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] +
1.320323631693739147855624220029111556865246*LogIntegral[x]
The log integral approximation is much faster than a full count.
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 ->
Placed[{"PrimePi2[x] ", "LogIntegral2[x]"}, Below]
]
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}