8
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ 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). $\endgroup$ – Daniel Lichtblau Jul 28 '18 at 14:52
  • $\begingroup$ Wolfram Challenges $\endgroup$ – J. M. will be back soon Mar 10 at 11:25
7
$\begingroup$

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}

$\endgroup$
13
$\begingroup$

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

$\endgroup$
8
$\begingroup$

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

$\endgroup$
4
$\begingroup$

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}
$\endgroup$
  • $\begingroup$ Nice. +1, but what and where is fastCompile? $\endgroup$ – KennyColnago Aug 1 '18 at 20:28
  • $\begingroup$ @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. $\endgroup$ – Chip Hurst Aug 1 '18 at 20:30
2
$\begingroup$
    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]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.