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

Mathematica graphics

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

Mathematica graphics


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

twin prime count


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}