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

8

2

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?

Richard Burke-Ward

Posted 2018-07-28T09:42:06.903

Reputation: 1 923

1Is 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). – Daniel Lichtblau – 2018-07-28T14:52:28.997

Wolfram Challenges – J. M.'s ennui – 2019-03-10T11:25:54.580

Answers

7

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}

Niki Estner

Posted 2018-07-28T09:42:06.903

Reputation: 34 978

13

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

rhermans

Posted 2018-07-28T09:42:06.903

Reputation: 26 859

8

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

KennyColnago

Posted 2018-07-28T09:42:06.903

Reputation: 14 269

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}

Chip Hurst

Posted 2018-07-28T09:42:06.903

Reputation: 29 735

Nice. +1, but what and where is fastCompile? – KennyColnago – 2018-08-01T20:28:08.653

@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. – Chip Hurst – 2018-08-01T20:30:53.947

2

    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]

David G. Stork

Posted 2018-07-28T09:42:06.903

Reputation: 31 784