Does Mathematica have a twin prime equivalent of `PrimePi`?

Question

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?

Answer 1

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}

Answer 2

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"
 ]

Answer 3

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]
   ]

Answer 4

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]&]]], 
  PadRight[FactorInteger[P][[All, 1]], EulerPhi[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}

Answer 5

    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]