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?

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?

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}

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

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

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

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}`

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

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.997Wolfram Challenges – J. M.'s ennui – 2019-03-10T11:25:54.580