Twin primes
You are encouraged to solve this task according to the task description, using any language you may know.
Twin primes are pairs of natural numbers (P1 and P2) that satisfy the following:
- P1 and P2 are primes
- P1 + 2 = P2
- Task
Write a program that displays the number of pairs of twin primes that can be found under a user-specified number
(P1 < user-specified number & P2 < user-specified number).
- Extension
-
- Find all twin prime pairs under 100000, 10000000 and 1000000000.
- What is the time complexity of the program? Are there ways to reduce computation time?
- Examples
> Search Size: 100 > 8 twin prime pairs.
> Search Size: 1000 > 35 twin prime pairs.
- Also see
- The OEIS entry: A001097: Twin primes.
- The OEIS entry: A167874: The number of distinct primes < 10^n which are members of twin-prime pairs.
- The OEIS entry: A077800: List of twin primes {p, p+2}, with repetition.
- The OEIS entry: A007508: Number of twin prime pairs below 10^n.
ALGOL 68
Simplifies array bound checking by using the equivalent definition of twin primes: p and p - 2. <lang algol68>BEGIN
# count twin primes (where p and p - 2 are prime) #
PR heap=128M PR # set heap memory size for Algol 68G #
# sieve of Eratosthenes: sets s[i] to TRUE if i is a prime, FALSE otherwise #
PROC sieve = ( REF[]BOOL s )VOID:
BEGIN
FOR i TO UPB s DO s[ i ] := TRUE OD;
s[ 1 ] := FALSE;
FOR i FROM 2 TO ENTIER sqrt( UPB s ) DO
IF s[ i ] THEN FOR p FROM i * i BY i TO UPB s DO s[ p ] := FALSE OD FI
OD
END # sieve # ;
# find the maximum number to search for twin primes #
INT max;
print( ( "Maximum: " ) );
read( ( max, newline ) );
INT max number = max;
# construct a sieve of primes up to the maximum number #
[ 1 : max number ]BOOL primes;
sieve( primes );
# count the twin primes #
# note 2 cannot be one of the primes in a twin prime pair, so we start at 3 #
INT twin count := 0;
FOR p FROM 3 BY 2 TO max number - 1 DO IF primes[ p ] AND primes[ p - 2 ] THEN twin count +:= 1 FI OD;
print( ( "twin prime pairs below ", whole( max number, 0 ), ": ", whole( twin count, 0 ), newline ) )
END</lang>
- Output:
Maximum: 10 twin prime pairs below 10: 2
Maximum: 100 twin prime pairs below 100: 8
Maximum: 1000 twin prime pairs below 1000: 35
Maximum: 10000 twin prime pairs below 10000: 205
Maximum: 100000 twin prime pairs below 100000: 1224
Maximum: 1000000 twin prime pairs below 1000000: 8169
Maximum: 10000000 twin prime pairs below 10000000: 58980
Arturo
<lang rebol>pairsOfPrimes: function [upperLim][
count: 0
j: 0
k: 1
i: 0
while [i=<upperLim][
i: (6 * k) - 1
j: i + 2
if and? [prime? i] [prime? j] [
count: count + 1
]
k: k + 1
]
return count + 1
]
ToNum: 10 while [ToNum =< 1000000][
x: pairsOfPrimes ToNum print ["From 2 to" ToNum ": there are" x "pairs of twin primes"] ToNum: ToNum * 10
]</lang>
- Output:
From 2 to 10 : there are 3 pairs of twin primes From 2 to 100 : there are 9 pairs of twin primes From 2 to 1000 : there are 35 pairs of twin primes From 2 to 10000 : there are 205 pairs of twin primes From 2 to 100000 : there are 1224 pairs of twin primes From 2 to 1000000 : there are 8169 pairs of twin primes
AWK
<lang AWK>
- syntax: GAWK -f TWIN_PRIMES.AWK
BEGIN {
n = 1
for (i=1; i<=6; i++) {
n *= 10
printf("twin prime pairs < %8s : %d\n",n,count_twin_primes(n))
}
exit(0)
} function count_twin_primes(limit, count,i,p1,p2,p3) {
p1 = 0
p2 = p3 = 1
for (i=5; i<=limit; i++) {
p3 = p2
p2 = p1
p1 = is_prime(i)
if (p3 && p1) {
count++
}
}
return(count)
} function is_prime(x, i) {
if (x <= 1) {
return(0)
}
for (i=2; i<=int(sqrt(x)); i++) {
if (x % i == 0) {
return(0)
}
}
return(1)
} </lang>
- Output:
twin prime pairs < 10 : 2 twin prime pairs < 100 : 8 twin prime pairs < 1000 : 35 twin prime pairs < 10000 : 205 twin prime pairs < 100000 : 1224 twin prime pairs < 1000000 : 8169
BASIC256
<lang BASIC256> function isPrime(v)
if v < 2 then return False if v mod 2 = 0 then return v = 2 if v mod 3 = 0 then return v = 3 d = 5 while d * d <= v if v mod d = 0 then return False else d += 2 end while return True
end function
function paresDePrimos(limite)
p1 = 0 p2 = 1 p3 = 1 cont = 0 for i = 5 to limite p3 = p2 p2 = p1 p1 = isPrime(i) if (p3 and p1) then cont += 1 next i return cont
end function
n = 1 for i = 1 to 6
n = n * 10 print "pares de primos gemelos por debajo de < "; n; " : "; paresDePrimos(n)
next i end </lang>
- Output:
Similar a la entrada de FreeBASIC.
C
<lang c>#include <stdbool.h>
- include <stdint.h>
- include <stdio.h>
bool isPrime(int64_t n) {
int64_t i;
if (n < 2) return false; if (n % 2 == 0) return n == 2; if (n % 3 == 0) return n == 3; if (n % 5 == 0) return n == 5; if (n % 7 == 0) return n == 7; if (n % 11 == 0) return n == 11; if (n % 13 == 0) return n == 13; if (n % 17 == 0) return n == 17; if (n % 19 == 0) return n == 19;
for (i = 23; i * i <= n; i += 2) {
if (n % i == 0) return false;
}
return true;
}
int countTwinPrimes(int limit) {
int count = 0;
// 2 3 4 int64_t p3 = true, p2 = true, p1 = false; int64_t i;
for (i = 5; i <= limit; i++) {
p3 = p2;
p2 = p1;
p1 = isPrime(i);
if (p3 && p1) {
count++;
}
}
return count;
}
void test(int limit) {
int count = countTwinPrimes(limit);
printf("Number of twin prime pairs less than %d is %d\n", limit, count);
}
int main() {
test(10); test(100); test(1000); test(10000); test(100000); test(1000000); test(10000000); test(100000000); return 0;
}</lang>
- Output:
Number of twin prime pairs less than 10 is 2 Number of twin prime pairs less than 100 is 8 Number of twin prime pairs less than 1000 is 35 Number of twin prime pairs less than 10000 is 205 Number of twin prime pairs less than 100000 is 1224 Number of twin prime pairs less than 1000000 is 8169 Number of twin prime pairs less than 10000000 is 58980 Number of twin prime pairs less than 100000000 is 440312
C++
The primesieve library includes the functionality required for this task. RC already has plenty of C++ examples showing how to generate prime numbers from scratch. (The module Math::Primesieve, which is used by the Raku example on this page, is implemented on top of this library.) <lang cpp>#include <cstdint>
- include <iostream>
- include <string>
- include <primesieve.hpp>
void print_twin_prime_count(long long limit) {
std::cout << "Number of twin prime pairs less than " << limit
<< " is " << (limit > 0 ? primesieve::count_twins(0, limit - 1) : 0) << '\n';
}
int main(int argc, char** argv) {
std::cout.imbue(std::locale(""));
if (argc > 1) {
// print number of twin prime pairs less than limits specified
// on the command line
for (int i = 1; i < argc; ++i) {
try {
print_twin_prime_count(std::stoll(argv[i]));
} catch (const std::exception& ex) {
std::cerr << "Cannot parse limit from '" << argv[i] << "'\n";
}
}
} else {
// if no limit was specified then show the number of twin prime
// pairs less than powers of 10 up to 100 billion
uint64_t limit = 10;
for (int power = 1; power < 12; ++power, limit *= 10)
print_twin_prime_count(limit);
}
return 0;
}</lang>
- Output:
Number of twin prime pairs less than 10 is 2 Number of twin prime pairs less than 100 is 8 Number of twin prime pairs less than 1,000 is 35 Number of twin prime pairs less than 10,000 is 205 Number of twin prime pairs less than 100,000 is 1,224 Number of twin prime pairs less than 1,000,000 is 8,169 Number of twin prime pairs less than 10,000,000 is 58,980 Number of twin prime pairs less than 100,000,000 is 440,312 Number of twin prime pairs less than 1,000,000,000 is 3,424,506 Number of twin prime pairs less than 10,000,000,000 is 27,412,679 Number of twin prime pairs less than 100,000,000,000 is 224,376,048
C#
Runs in about 18 seconds. <lang csharp>using System; using System.Linq; using System.Collections; using System.Collections.Generic;
public static class TwinPrimes {
public static void Main()
{
CountTwinPrimes(Enumerable.Range(1, 9).Select(i => (int)Math.Pow(10, i)).ToArray());
}
private static void CountTwinPrimes(params int[] bounds)
{
Array.Sort(bounds);
int b = 0;
int count = 0;
string format = "There are {0:N0} twin primes below {1:N0}";
foreach (var twin in FindTwinPrimes(bounds[^1])) {
if (twin.p2 >= bounds[b]) {
Console.WriteLine(format, count, bounds[b]);
b++;
}
count++;
}
Console.WriteLine(format, count, bounds[b]);
}
private static IEnumerable<(int p1, int p2)> FindTwinPrimes(int bound) =>
PrimeSieve(bound).Pairwise().Where(pair => pair.p1 + 2 == pair.p2);
private static IEnumerable<int> PrimeSieve(int bound)
{
if (bound < 2) yield break;
yield return 2;
var composite = new BitArray((bound - 1) / 2);
int limit = (int)(Math.Sqrt(bound) - 1) / 2;
for (int i = 0; i < limit; i++) {
if (composite[i]) continue;
int prime = 2 * i + 3;
yield return prime;
for (int j = (prime * prime - 2) / 2; j < composite.Count; j += prime) {
composite[j] = true;
}
}
for (int i = limit; i < composite.Count; i++) {
if (!composite[i]) yield return 2 * i + 3;
}
}
private static IEnumerable<(T p1, T p2)> Pairwise<T>(this IEnumerable<T> source)
{
using var e = numbers.GetEnumerator();
if (!e.MoveNext()) yield break;
T p1 = e.Current;
while (e.MoveNext()) {
T p2 = e.Current;
yield return (p1, p2);
p1 = p2;
}
}
}</lang>
- Output:
There are 2 twin primes below 10 There are 8 twin primes below 100 There are 35 twin primes below 1,000 There are 205 twin primes below 10,000 There are 1,224 twin primes below 100,000 There are 8,169 twin primes below 1,000,000 There are 58,980 twin primes below 10,000,000 There are 440,312 twin primes below 100,000,000 There are 3,424,506 twin primes below 1,000,000,000
Delphi
<lang Delphi> program Primes;
{$APPTYPE CONSOLE}
{$R *.res}
uses
System.SysUtils;
function IsPrime(a: UInt64): Boolean; var
d: UInt64;
begin
if (a < 2) then exit(False);
if (a mod 2) = 0 then exit(a = 2);
if (a mod 3) = 0 then exit(a = 3);
d := 5;
while (d * d <= a) do
begin
if (a mod d = 0) then
Exit(false);
inc(d, 2);
if (a mod d = 0) then
Exit(false);
inc(d, 4);
end;
Result := True;
end;
function Sieve(limit: UInt64): TArray<Boolean>;
var
p, p2, i: UInt64;
begin
inc(limit); SetLength(Result, limit); FillChar(Result[2], sizeof(Boolean) * limit - 2, 0); // all false except 1,2 FillChar(Result[0], sizeof(Boolean) * 2, 1); // 1,2 are true
p := 3;
while true do
begin
p2 := p * p;
if p2 >= limit then
break;
i := p2;
while i < limit do
begin
Result[i] := true;
inc(i, 2 * p);
end;
while true do
begin
inc(p, 2);
if not Result[p] then
Break;
end;
end;
end;
function Commatize(const n: UInt64): string; var
str: string; digits: Integer; i: Integer;
begin
Result := ; str := n.ToString; digits := str.Length;
for i := 1 to digits do
begin
if ((i > 1) and (((i - 1) mod 3) = (digits mod 3))) then
Result := Result + ',';
Result := Result + str[i];
end;
end;
var
limit, start, twins: UInt64; c: TArray<Boolean>; i, j: UInt64;
begin
c := Sieve(Trunc(1e9 - 1));
limit := 10;
start := 3;
twins := 0;
for i := 1 to 9 do
begin
j := start;
while j < limit do
begin
if (not c[j]) and (not c[j - 2]) then
inc(twins);
inc(j, 2);
end;
Writeln(Format('Under %14s there are %10s pairs of twin primes.', [commatize
(limit), commatize(twins)]));
start := limit + 1; limit := 10 * limit; end;
readln;
end.
</lang>
- Output:
Under 10 there are 2 pairs of twin primes. Under 100 there are 8 pairs of twin primes. Under 1,000 there are 35 pairs of twin primes. Under 10,000 there are 205 pairs of twin primes. Under 100,000 there are 1,224 pairs of twin primes. Under 1,000,000 there are 8,169 pairs of twin primes. Under 10,000,000 there are 58,980 pairs of twin primes. Under 100,000,000 there are 440,312 pairs of twin primes. Under 1,000,000,000 there are 3,424,506 pairs of twin primes.
F#
This task uses Extensible Prime Generator (F#) <lang fsharp> printfn "twin primes below 100000: %d" (primes64()|>Seq.takeWhile(fun n->n<=100000L)|>Seq.pairwise|>Seq.filter(fun(n,g)->g=n+2L)|>Seq.length) printfn "twin primes below 1000000: %d" (primes64()|>Seq.takeWhile(fun n->n<=1000000L)|>Seq.pairwise|>Seq.filter(fun(n,g)->g=n+2L)|>Seq.length) printfn "twin primes below 10000000: %d" (primes64()|>Seq.takeWhile(fun n->n<=10000000L)|>Seq.pairwise|>Seq.filter(fun(n,g)->g=n+2L)|>Seq.length) printfn "twin primes below 100000000: %d" (primes64()|>Seq.takeWhile(fun n->n<=100000000L)|>Seq.pairwise|>Seq.filter(fun(n,g)->g=n+2L)|>Seq.length) printfn "twin primes below 1000000000: %d" (primes64()|>Seq.takeWhile(fun n->n<=1000000000L)|>Seq.pairwise|>Seq.filter(fun(n,g)->g=n+2L)|>Seq.length) printfn "twin primes below 10000000000: %d" (primes64()|>Seq.takeWhile(fun n->n<=10000000000L)|>Seq.pairwise|>Seq.filter(fun(n,g)->g=n+2L)|>Seq.length) printfn "twin primes below 100000000000: %d" (primes64()|>Seq.takeWhile(fun n->n<=100000000000L)|>Seq.pairwise|>Seq.filter(fun(n,g)->g=n+2L)|>Seq.length) </lang>
- Output:
twin primes below 100000: 1224 Real: 00:00:00.003, CPU: 00:00:00.015, GC gen0: 0, gen1: 0, gen2: 0 twin primes below 1000000: 8169 Real: 00:00:00.021, CPU: 00:00:00.031, GC gen0: 3, gen1: 3, gen2: 0 twin primes below 10000000: 58980 Real: 00:00:00.154, CPU: 00:00:00.171, GC gen0: 19, gen1: 19, gen2: 0 twin primes below 100000000: 440312 Real: 00:00:01.400, CPU: 00:00:01.406, GC gen0: 162, gen1: 162, gen2: 0 twin primes below 1000000000: 3424506 Real: 00:00:12.682, CPU: 00:00:12.671, GC gen0: 1428, gen1: 1426, gen2: 1 twin primes below 10000000000: 27412679 Real: 00:02:04.441, CPU: 00:02:04.406, GC gen0: 12771, gen1: 12768, gen2: 2 twin primes below 100000000000: 224376048 Real: 00:23:00.853, CPU: 00:23:00.125, GC gen0: 115562, gen1: 115536, gen2: 14
Factor
<lang factor>USING: io kernel math math.parser math.primes.erato math.ranges sequences tools.memory.private ;
- twin-pair-count ( n -- count )
[ 5 swap 2 <range> ] [ sieve ] bi [ over 2 - over [ marked-prime? ] 2bi@ and ] curry count ;
"Search size: " write flush readln string>number twin-pair-count commas write " twin prime pairs." print</lang>
- Output:
Search size: 100,000 1,224 twin prime pairs.
Search size: 10,000,000 58,980 twin prime pairs.
Search size: 1,000,000,000 3,424,506 twin prime pairs.
FreeBASIC
<lang freebasic> Function isPrime(Byval ValorEval As Integer) As Boolean
If ValorEval <=1 Then Return False
For i As Integer = 2 To Int(Sqr(ValorEval))
If ValorEval Mod i = 0 Then Return False
Next i
Return True
End Function
Function paresDePrimos(limite As Uinteger) As Uinteger
Dim As Uinteger p1 = 0, p2 = 1, p3 = 1, count = 0
For i As Uinteger = 5 To limite
p3 = p2
p2 = p1
p1 = isPrime(i)
If (p3 And p1) Then count += 1
Next i
Return count
End Function
Dim As Uinteger n = 1 For i As Byte = 1 To 6
n *= 10 Print Using "pares de primos gemelos por debajo de < ####### : ####"; n; paresDePrimos(n)
Next i Print !"\n--- terminado, pulsa RETURN---" Sleep </lang>
- Output:
pares de primos gemelos por debajo de < 10 : 2 pares de primos gemelos por debajo de < 100 : 8 pares de primos gemelos por debajo de < 1000 : 35 pares de primos gemelos por debajo de < 10000 : 205 pares de primos gemelos por debajo de < 100000 : 1224 pares de primos gemelos por debajo de < 1000000 : 8169
Go
<lang go>package main
import "fmt"
func sieve(limit uint64) []bool {
limit++
// True denotes composite, false denotes prime.
c := make([]bool, limit) // all false by default
c[0] = true
c[1] = true
// no need to bother with even numbers over 2 for this task
p := uint64(3) // Start from 3.
for {
p2 := p * p
if p2 >= limit {
break
}
for i := p2; i < limit; i += 2 * p {
c[i] = true
}
for {
p += 2
if !c[p] {
break
}
}
}
return c
}
func commatize(n int) string {
s := fmt.Sprintf("%d", n)
if n < 0 {
s = s[1:]
}
le := len(s)
for i := le - 3; i >= 1; i -= 3 {
s = s[0:i] + "," + s[i:]
}
if n >= 0 {
return s
}
return "-" + s
}
func main() {
c := sieve(1e10 - 1)
limit := 10
start := 3
twins := 0
for i := 1; i < 11; i++ {
for i := start; i < limit; i += 2 {
if !c[i] && !c[i-2] {
twins++
}
}
fmt.Printf("Under %14s there are %10s pairs of twin primes.\n", commatize(limit), commatize(twins))
start = limit + 1
limit *= 10
}
}</lang>
- Output:
Under 10 there are 2 pairs of twin primes. Under 100 there are 8 pairs of twin primes. Under 1,000 there are 35 pairs of twin primes. Under 10,000 there are 205 pairs of twin primes. Under 100,000 there are 1,224 pairs of twin primes. Under 1,000,000 there are 8,169 pairs of twin primes. Under 10,000,000 there are 58,980 pairs of twin primes. Under 100,000,000 there are 440,312 pairs of twin primes. Under 1,000,000,000 there are 3,424,506 pairs of twin primes. Under 10,000,000,000 there are 27,412,679 pairs of twin primes.
Alternative using primegen package
See Extensible prime generator#Go. <lang go>package main
import (
"fmt" "github.com/jbarham/primegen.go"
)
func main() {
p := primegen.New()
count := 0
previous := uint64(0)
power := 1
limit := uint64(10)
for {
prime := p.Next()
if prime >= limit {
fmt.Printf("Number of twin prime pairs less than %d: %d\n", limit, count)
power++
if power > 10 {
break
}
limit *= 10
}
if previous > 0 && prime == previous + 2 {
count++
}
previous = prime
}
}</lang>
- Output:
Number of twin prime pairs less than 10: 2 Number of twin prime pairs less than 100: 8 Number of twin prime pairs less than 1000: 35 Number of twin prime pairs less than 10000: 205 Number of twin prime pairs less than 100000: 1224 Number of twin prime pairs less than 1000000: 8169 Number of twin prime pairs less than 10000000: 58980 Number of twin prime pairs less than 100000000: 440312 Number of twin prime pairs less than 1000000000: 3424506 Number of twin prime pairs less than 10000000000: 27412679
Java
BigInteger Implementation: <lang Java> import java.math.BigInteger; import java.util.Scanner;
public class twinPrimes {
public static void main(String[] args) {
Scanner input = new Scanner(System.in);
System.out.println("Search Size: ");
BigInteger max = input.nextBigInteger();
int counter = 0;
for(BigInteger x = new BigInteger("3"); x.compareTo(max) <= 0; x = x.add(BigInteger.ONE)){
BigInteger sqrtNum = x.sqrt().add(BigInteger.ONE);
if(x.add(BigInteger.TWO).compareTo(max) <= 0) {
counter += findPrime(x.add(BigInteger.TWO), x.add(BigInteger.TWO).sqrt().add(BigInteger.ONE)) && findPrime(x, sqrtNum) ? 1 : 0;
}
}
System.out.println(counter + " twin prime pairs.");
}
public static boolean findPrime(BigInteger x, BigInteger sqrtNum){
for(BigInteger divisor = BigInteger.TWO; divisor.compareTo(sqrtNum) <= 0; divisor = divisor.add(BigInteger.ONE)){
if(x.remainder(divisor).compareTo(BigInteger.ZERO) == 0){
return false;
}
}
return true;
}
} </lang>
- Output:
> Search Size: > 100 > 8 twin prime pairs.
> Search Size: > 1000 > 35 twin prime pairs.
jq
Slightly modified from the C entry
<lang jq>def odd_gt2_is_prime:
. as $n
| if ($n % 3 == 0) then $n == 3
elif ($n % 5 == 0) then $n == 5
elif ($n % 7 == 0) then $n == 7
elif ($n % 11 == 0) then $n == 11
elif ($n % 13 == 0) then $n == 13
elif ($n % 17 == 0) then $n == 17
elif ($n % 19 == 0) then $n == 19
else {i:23}
| until( (.i * .i) > $n or ($n % .i == 0); .i += 2)
| .i * .i > $n
end;
def twin_primes($max):
{count:0, i:3, isprime:true}
| until(.i >= $max;
.i += 2
| if .isprime
then if .i|odd_gt2_is_prime then .count+=1 else .isprime = false end
else .isprime = (.i|odd_gt2_is_prime)
end )
| .count;
pow(10; range(1;8)) | "Number of twin primes less than \(.) is \(twin_primes(.))."</lang>
- Output:
Number of twin primes less than 10 is 2. Number of twin primes less than 100 is 8. Number of twin primes less than 1000 is 35. Number of twin primes less than 10000 is 205. Number of twin primes less than 100000 is 1224. Number of twin primes less than 1000000 is 8169. Number of twin primes less than 10000000 is 58980.
Julia
<lang julia>using Formatting, Primes
function counttwinprimepairsbetween(n1, n2)
npairs, t = 0, nextprime(n1)
while t < n2
p = nextprime(t + 1)
if p - t == 2
npairs += 1
end
t = p
end
return npairs
end
for t2 in (10).^collect(2:8)
paircount = counttwinprimepairsbetween(1, t2)
println("Under", lpad(format(t2, commas=true), 12), " there are",
lpad(format(paircount, commas=true), 8), " pairs of twin primes.")
end
</lang>
- Output:
Under 100 there are 8 pairs of twin primes. Under 1,000 there are 35 pairs of twin primes. Under 10,000 there are 205 pairs of twin primes. Under 100,000 there are 1,224 pairs of twin primes. Under 1,000,000 there are 8,169 pairs of twin primes. Under 10,000,000 there are 58,980 pairs of twin primes. Under 100,000,000 there are 440,312 pairs of twin primes.
Extension to large n and other tuples
Task Extension: to get primes up to a billion, it becomes important to cache the results so that large numbers do not need to be factored more than once. This trades memory for speed. The time complexity is dominated by the prime sieve time used to create the primes mask, which is n log(log n).
We can generalize pairs to reflect any tuple of integer differences between the first prime and the successive primes: see http://www.rosettacode.org/wiki/Successive_prime_differences.
If we ignore the first difference from the index prime with itself (always 0), we can express a prime pair as a difference tuple of (2,), and a prime quadruplet such as [11, 13, 17, 19] as the tuple starting with 11 of type (2, 6, 8). <lang julia>using Formatting, Primes
const PMAX = 1_000_000_000 const pmb = primesmask(PMAX) const primestoabillion = [i for i in 2:PMAX if pmb[i]]
tuplefitsat(k, tup, arr) = all(i -> arr[k + i] - arr[k] == tup[i], 1:length(tup))
function countprimetuples(tup, n)
arr = filter(i -> i <= n, primestoabillion) return count(k -> tuplefitsat(k, tup, arr), 1:length(arr) - length(tup))
end
println("Count of prime pairs from 1 to 1 billion: ",
format(countprimetuples((2,), 1000000000), commas=true))
println("Count of a form of prime quads from 1 to 1 million: ",
format(countprimetuples((2, 6, 8), 1000000), commas=true))
println("Count of a form of prime octets from 1 to 1 million: ",
format(countprimetuples((2, 6, 12, 14, 20, 24, 26), 1000000), commas=true))
</lang>
- Output:
Count of prime pairs from 1 to 1 billion: 3,424,506 Count of a form of prime quads from 1 to 1 million: 166 Count of a form of prime octets from 1 to 1 million: 3
Kotlin
<lang scala>import java.math.BigInteger import java.util.*
fun main() {
val input = Scanner(System.`in`)
println("Search Size: ")
val max = input.nextBigInteger()
var counter = 0
var x = BigInteger("3")
while (x <= max) {
val sqrtNum = x.sqrt().add(BigInteger.ONE)
if (x.add(BigInteger.TWO) <= max) {
counter += if (findPrime(
x.add(BigInteger.TWO),
x.add(BigInteger.TWO).sqrt().add(BigInteger.ONE)
) && findPrime(x, sqrtNum)
) 1 else 0
}
x = x.add(BigInteger.ONE)
}
println("$counter twin prime pairs.")
}
fun findPrime(x: BigInteger, sqrtNum: BigInteger?): Boolean {
var divisor = BigInteger.TWO
while (divisor <= sqrtNum) {
if (x.remainder(divisor).compareTo(BigInteger.ZERO) == 0) {
return false
}
divisor = divisor.add(BigInteger.ONE)
}
return true
}</lang>
Mathematica/Wolfram Language
<lang Mathematica>ClearAll[TwinPrimeCount] TwinPrimeCount[mx_] := Module[{pmax, min, max, total},
pmax = PrimePi[mx];
total = 0;
Do[
min = 10^6 i;
min = Max[min, 1];
max = 10^6 (i + 1);
max = Min[max, pmax];
total += Count[Differences[Prime[Range[min, max]]], 2]
,
{i, 0, Ceiling[pmax/10^6]}
];
total
]
Do[Print[{10^i, TwinPrimeCount[10^i]}], {i, 9}]</lang>
- Output:
{10,2}
{100,8}
{1000,35}
{10000,205}
{100000,1224}
{1000000,8169}
{10000000,58980}
{100000000,440312}
{1000000000,3424506}
Nim
We use a sieve of Erathostenes which needs a lot of memory. It is possible to reduce memory usage by using bit strings for the sieve (one bit per boolean instead of eight bits), but the price is a significant loss of performance.
As, except for the pair (3, 5), all twins pairs are composed of a number congruent to 2 modulo 3 and a number congruent to 1 modulo 3, we can save some time by using a step of 6. Unfortunately, this is the filling of the sieve which is the most time consuming, so the gain is not very important (on our computer, half a second on a total time of 8.3 s).
<lang Nim>import math, strformat, strutils
const N = 1_000_000_000
proc sieve(n: Positive): seq[bool] =
## Build and fill a sieve of Erathosthenes.
result.setLen(n + 1) # Default to false which means prime.
result[0] = true
result[1] = true
for n in countup(3, sqrt(N.toFloat).int, 2):
if not result[n]:
for k in countup(n * n, N, 2 * n):
result[k] = true
let composite = sieve(N)
proc findTwins(composite: openArray[bool]) =
var
lim = 10
count = 1 # Start with 3, 5 which is a special case.
n = 7 # First prime congruent to 1 modulo 3.
while true:
if not composite[n] and not composite[n - 2]: inc count
inc n, 6 # Next odd number congruent to 1 modulo 3.
if n > lim:
echo &"There are {insertSep($count)} pairs of twin primes under {insertSep($lim)}."
lim *= 10
if lim > N: break
composite.findTwins()</lang>
- Output:
There are 2 pairs of twin primes under 10. There are 8 pairs of twin primes under 100. There are 35 pairs of twin primes under 1_000. There are 205 pairs of twin primes under 10_000. There are 1_224 pairs of twin primes under 100_000. There are 8_169 pairs of twin primes under 1_000_000. There are 58_980 pairs of twin primes under 10_000_000. There are 440_312 pairs of twin primes under 100_000_000. There are 3_424_506 pairs of twin primes under 1_000_000_000.
Perl
<lang perl>use strict; use warnings;
use Primesieve;
sub comma { reverse ((reverse shift) =~ s/(.{3})/$1,/gr) =~ s/^,//r }
printf "Twin prime pairs less than %14s: %s\n", comma(10**$_), comma count_twins(1, 10**$_) for 1..10;</lang>
- Output:
Twin prime pairs less than 10: 2 Twin prime pairs less than 100: 8 Twin prime pairs less than 1,000: 35 Twin prime pairs less than 10,000: 205 Twin prime pairs less than 100,000: 1,224 Twin prime pairs less than 1,000,000: 8,169 Twin prime pairs less than 10,000,000: 58,980 Twin prime pairs less than 100,000,000: 440,312 Twin prime pairs less than 1,000,000,000: 3,424,506 Twin prime pairs less than 10,000,000,000: 27,412,679
Phix
Added both parameter to reflect the recent task specification changes, as shown for a limit of 6 you can count {3,5} and {5,7} as one pair (the default, matching task description) or two. Obviously delete the "6 --" if you actually want a prompt.
The time complexity here is all about building a table of primes. It turns out that using the builtin get_prime() is actually faster than using an explicit sieve (as per Delphi/Go/Wren) due to retaining all the intermediate 0s, not that I particularly expect this to win any performance trophies.
with javascript_semantics atom t0 = time() function twin_primes(integer maxp, bool both=true) integer n = 0, -- result pn = 2, -- next prime index p, -- a prime, <= maxp prev_p = 2 while true do p = get_prime(pn) if both and p>=maxp then exit end if n += (prev_p = p-2) if (not both) and p>=maxp then exit end if prev_p = p pn += 1 end while return n end function integer mp = 6 -- prompt_number("Enter limit:") printf(1,"Twin prime pairs less than %,d: %,d\n",{mp,twin_primes(mp)}) printf(1,"Twin prime pairs less than %,d: %,d\n",{mp,twin_primes(mp,false)}) for p=1 to 9 do integer p10 = power(10,p) printf(1,"Twin prime pairs less than %,d: %,d\n",{p10,twin_primes(p10)}) end for ?elapsed(time()-t0)
- Output:
Twin prime pairs less than 6: 1 Twin prime pairs less than 6: 2 Twin prime pairs less than 10: 2 Twin prime pairs less than 100: 8 Twin prime pairs less than 1,000: 35 Twin prime pairs less than 10,000: 205 Twin prime pairs less than 100,000: 1,224 Twin prime pairs less than 1,000,000: 8,169 Twin prime pairs less than 10,000,000: 58,980 Twin prime pairs less than 100,000,000: 440,312 Twin prime pairs less than 1,000,000,000: 3,424,506 "16.2s"
PureBasic
<lang PureBasic> Procedure isPrime(v.i)
If v <= 1 : ProcedureReturn #False
ElseIf v < 4 : ProcedureReturn #True
ElseIf v % 2 = 0 : ProcedureReturn #False
ElseIf v < 9 : ProcedureReturn #True
ElseIf v % 3 = 0 : ProcedureReturn #False
Else
Protected r = Round(Sqr(v), #PB_Round_Down)
Protected f = 5
While f <= r
If v % f = 0 Or v % (f + 2) = 0
ProcedureReturn #False
EndIf
f + 6
Wend
EndIf
ProcedureReturn #True
EndProcedure
Procedure paresDePrimos(limite.d)
p1.i = 0
p2.i = 1
p3.i = 1
count.i = 0
For i.i = 5 To limite
p3 = p2
p2 = p1
p1 = isPrime(i)
If p3 And p1
count + 1
EndIf
Next i
ProcedureReturn count
EndProcedure
OpenConsole() n.i = 1 For i.i = 1 To 6
n = n * 10
PrintN("pares de primos gemelos por debajo de < " + Str(n) + " : " + Str(paresDePrimos(n)))
Next i PrintN(#CRLF$ + "--- terminado, pulsa RETURN---"): Input() CloseConsole() End </lang>
- Output:
Similar a la entrada de FreeBASIC.
Python
<lang python>prime_dict = {}
def is_prime(n: int) -> bool:
global prime_dict
if n not in prime_dict:
prime_dict[n] = is_prime_(n)
if n >= 23:
for x in range(n * 2, n * 50, n):
prime_dict[x] = False # Imitation of sieve of Eratosthenes
return prime_dict[n]
def is_prime_(n: int) -> bool:
for x in [2, 3, 5, 7, 11, 13, 17, 19]:
if n % x == 0:
return n == x
i = 23
while i ** 2 <= n:
if n % i == 0 or n % (i + 2) == 0:
return False
i += 6
return True
def count_twin_primes(limit: int) -> int:
p2 = True
p1 = False
count = 0
for i in range(5, limit):
p3 = p2
p2 = p1
p1 = is_prime(i)
if p3 and p1:
count += 1
return count
def test(limit: int):
count = count_twin_primes(limit)
print(f"Number of twin prime pairs less than {limit} is {count}")
if __name__ == '__main__':
test(10) test(100) test(1000) test(10000) test(100000) test(1000000) test(10000000) test(100000000)</lang>
- Output:
Number of twin prime pairs less than 10 is 2 Number of twin prime pairs less than 100 is 8 Number of twin prime pairs less than 1000 is 35 Number of twin prime pairs less than 10000 is 205 Number of twin prime pairs less than 100000 is 1224 Number of twin prime pairs less than 1000000 is 8169 Number of twin prime pairs less than 10000000 is 58980 Number of twin prime pairs less than 100000000 is 440312
Raku
<lang perl6>use Lingua::EN::Numbers;
use Math::Primesieve;
my $p = Math::Primesieve.new;
printf "Twin prime pairs less than %14s: %s\n", comma(10**$_), comma $p.count(10**$_, :twins) for 1 .. 10;</lang>
- Output:
Twin prime pairs less than 10: 2 Twin prime pairs less than 100: 8 Twin prime pairs less than 1,000: 35 Twin prime pairs less than 10,000: 205 Twin prime pairs less than 100,000: 1,224 Twin prime pairs less than 1,000,000: 8,169 Twin prime pairs less than 10,000,000: 58,980 Twin prime pairs less than 100,000,000: 440,312 Twin prime pairs less than 1,000,000,000: 3,424,506 Twin prime pairs less than 10,000,000,000: 27,412,679
REXX
straight-forward prime generator
The genP function could be optimized for higher specifications of the limit(s). <lang rexx>/*REXX pgm counts the number of twin prime pairs under a specified number N (or a list).*/ parse arg $ . /*get optional number of primes to find*/ if $= | $="," then $= 10 100 1000 10000 100000 1000000 10000000 /*No $? Use default.*/ w= length( commas( word($, words($) ) ) ) /*get length of the last number in list*/ @found= ' twin prime pairs found under ' /*literal used in the showing of output*/
do i=1 for words($); x= word($, i) /*process each N─limit in the $ list.*/
say right( commas(genP(x)), 20) @found right(commas(x), max(length(x), w) )
end /*i*/
exit 0 /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ commas: parse arg _; do ?=length(_)-3 to 1 by -3; _=insert(',', _, ?); end; return _ /*──────────────────────────────────────────────────────────────────────────────────────*/ genP: parse arg y; @.1=2; @.2=3; @.3=5; @.4=7; @.5=11; @.6=13; #= 6; tp= 2; sq.6= 169
if y>10 then tp= tp+1
do j=@.#+2 by 2 for max(0, y%2-@.#%2-1) /*find odd primes from here on. */
parse var j -1 _ /*obtain the last digit of the J var.*/
if _==5 then iterate; if j// 3==0 then iterate /*J ÷ by 5? J ÷ by 3? */
if j//7==0 then iterate; if j//11==0 then iterate /*" " " 7? " " " 11? */
/* [↓] divide by the primes. ___ */
do k=6 to # while sq.k<=j /*divide J by other primes ≤ √ J */
if j//@.k == 0 then iterate j /*÷ by prev. prime? ¬prime ___ */
end /*k*/ /* [↑] only divide up to √ J */
prev= @.#; #= #+1; sq.#= j*j; @.#= j /*save prev. P; bump # primes; assign P*/
if j-2==prev then tp= tp + 1 /*This & previous prime twins? Bump TP.*/
end /*j*/; return tp</lang>
- output when using the default inputs:
2 twin prime pairs found under 10
8 twin prime pairs found under 100
35 twin prime pairs found under 1,000
205 twin prime pairs found under 10,000
1,224 twin prime pairs found under 100,000
8,169 twin prime pairs found under 1,000,000
58,980 twin prime pairs found under 10,000,000
optimized prime number generator
This REXX version has some optimization for prime generation.
This version won't return a correct value (for the number of twin pairs) for a limit < 73 (because of the manner in
which low primes are generated from a list), however, the primes are returned from the function.
<lang rexx>/*REXX pgm counts the number of twin prime pairs under a specified number N (or a list).*/
parse arg $ . /*get optional number of primes to find*/
if $= | $="," then $= 100 1000 10000 100000 1000000 10000000 /*No $? Use default.*/
w= length( commas( word($, words($) ) ) ) /*get length of the last number in list*/
@found= ' twin prime pairs found under ' /*literal used in the showing of output*/
do i=1 for words($); x= word($, i) /*process each N─limit in the $ list.*/
say right( commas(genP(x)), 20) @found right(commas(x), max(length(x), w) )
end /*i*/
exit 0 /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ commas: parse arg _; do ?=length(_)-3 to 1 by -3; _=insert(',', _, ?); end; return _ /*──────────────────────────────────────────────────────────────────────────────────────*/ genP: arg y; _= 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101
tp=8; #= words(_); sq.103=103*103 /*#: number of prims; TP: # twin pairs.*/
do aa=1 for #; @.aa= word(_, aa) /*assign some low primes for quick ÷'s.*/
end /*aa*/
do j=@.#+2 by 2 while j<y /*continue with the next prime past 101*/
parse var j -1 _ /*obtain the last digit of the J var.*/
if _ ==5 then iterate /*is this integer a multiple of five? */
if j//3 ==0 then iterate /* " " " " " " three? */
do a=4 for 23 /*divide low primes starting with seven*/
if j//@.a ==0 then iterate j /*is integer a multiple of a low prime?*/
end /*a*/
/* [↓] divide by the primes. ___ */
do k=27 to # while sq.k<= j /*divide J by other primes ≤ √ J */
if j//@.k ==0 then iterate j /*÷ by prev. prime? ¬prime ___ */
end /*k*/ /* [↑] only divide up to √ J */
prev= @.#; #= #+1; sq.#= j*j; @.#= j /*save prev. P; bump # primes; assign P*/
if j-2==prev then tp= tp + 1 /*This & previous prime twins? Bump TP.*/
end /*j*/; return tp</lang>
Ring
<lang ring> load "stdlib.ring"
limit = list(7) for n = 1 to 7
limit[n] = pow(10,n)
next
TwinPrimes = []
for n = 1 to limit[7]-2
bool1 = isprime(n)
bool2 = isprime(n+2)
bool = bool1 and bool2
if bool =1
add(TwinPrimes,[n,n+2])
ok
next
numTwin = list(7) len = len(TwinPrimes)
for n = 1 to len
for p = 1 to 6
if TwinPrimes[n][2] < pow(10,p) and TwinPrimes[n+1][1] > pow(10,p)-2
numTwin[p] = n
ok
next
next
numTwin[7] = len
for n = 1 to 7
see "Maximum: " + pow(10,n) + nl see "twin prime pairs below " + pow(10,n) + ": " + numTwin[n] + nl + nl
next </lang> Output:
Maximum: 10 twin prime pairs below 10: 2 Maximum: 100 twin prime pairs below 100: 8 Maximum: 1000 twin prime pairs below 1000: 35 Maximum: 10000 twin prime pairs below 10000: 205 Maximum: 100000 twin prime pairs below 100000: 1224 Maximum: 1000000 twin prime pairs below 1000000: 8169 Maximum: 10000000 twin prime pairs below 10000000: 58980
Ruby
<lang ruby>require 'prime'
(1..8).each do |n|
count = Prime.each(10**n).each_cons(2).count{|p1, p2| p2-p1 == 2}
puts "Twin primes below 10**#{n}: #{count}"
end </lang>
- Output:
Twin primes below 10**1: 2 Twin primes below 10**2: 8 Twin primes below 10**3: 35 Twin primes below 10**4: 205 Twin primes below 10**5: 1224 Twin primes below 10**6: 8169 Twin primes below 10**7: 58980 Twin primes below 10**8: 440312
Rust
Limits can be specified on the command line, otherwise the twin prime counts for powers of ten from 1 to 10 are shown. <lang rust>// [dependencies] // primal = "0.3" // num-format = "0.4"
use num_format::{Locale, ToFormattedString};
fn twin_prime_count_for_powers_of_ten(max_power: u32) {
let mut count = 0;
let mut previous = 0;
let mut power = 1;
let mut limit = 10;
for prime in primal::Primes::all() {
if prime > limit {
println!(
"Number of twin prime pairs less than {} is {}",
limit.to_formatted_string(&Locale::en),
count.to_formatted_string(&Locale::en)
);
limit *= 10;
power += 1;
if power > max_power {
break;
}
}
if previous > 0 && prime == previous + 2 {
count += 1;
}
previous = prime;
}
}
fn twin_prime_count(limit: usize) {
let mut count = 0;
let mut previous = 0;
for prime in primal::Primes::all().take_while(|x| *x < limit) {
if previous > 0 && prime == previous + 2 {
count += 1;
}
previous = prime;
}
println!(
"Number of twin prime pairs less than {} is {}",
limit.to_formatted_string(&Locale::en),
count.to_formatted_string(&Locale::en)
);
}
fn main() {
let args: Vec<String> = std::env::args().collect();
if args.len() > 1 {
for i in 1..args.len() {
if let Ok(limit) = args[i].parse::<usize>() {
twin_prime_count(limit);
} else {
eprintln!("Cannot parse limit from string {}", args[i]);
}
}
} else {
twin_prime_count_for_powers_of_ten(10);
}
}</lang>
- Output:
Number of twin prime pairs less than 10 is 2 Number of twin prime pairs less than 100 is 8 Number of twin prime pairs less than 1,000 is 35 Number of twin prime pairs less than 10,000 is 205 Number of twin prime pairs less than 100,000 is 1,224 Number of twin prime pairs less than 1,000,000 is 8,169 Number of twin prime pairs less than 10,000,000 is 58,980 Number of twin prime pairs less than 100,000,000 is 440,312 Number of twin prime pairs less than 1,000,000,000 is 3,424,506 Number of twin prime pairs less than 10,000,000,000 is 27,412,679
Sidef
<lang ruby>func twin_primes_count(upto) {
var count = 0
var p1 = 2
each_prime(3, upto, {|p2|
if (p2 - p1 == 2) {
++count
}
p1 = p2
})
return count
}
for n in (1..9) {
var count = twin_primes_count(10**n)
say "There are #{count} twin primes <= 10^#{n}"
}</lang>
- Output:
There are 2 twin primes <= 10^1 There are 8 twin primes <= 10^2 There are 35 twin primes <= 10^3 There are 205 twin primes <= 10^4 There are 1224 twin primes <= 10^5 There are 8169 twin primes <= 10^6 There are 58980 twin primes <= 10^7 There are 440312 twin primes <= 10^8 There are 3424506 twin primes <= 10^9
Visual Basic
<lang vb>Function IsPrime(x As Long) As Boolean
Dim i As Long
If x Mod 2 = 0 Then
Exit Function
Else
For i = 3 To Int(Sqr(x)) Step 2
If x Mod i = 0 Then Exit Function
Next i
End If
IsPrime = True
End Function
Function TwinPrimePairs(max As Long) As Long
Dim p1 As Boolean, p2 As Boolean, count As Long, i As Long
p2 = True
For i = 5 To max Step 2
p1 = p2
p2 = IsPrime(i)
If p1 And p2 Then count = count + 1
Next i
TwinPrimePairs = count
End Function
Sub Test(x As Long)
Debug.Print "Twin prime pairs below" + Str(x) + ":" + Str(TwinPrimePairs(x))
End Sub
Sub Main()
Test 10 Test 100 Test 1000 Test 10000 Test 100000 Test 1000000 Test 10000000
End Sub</lang>
- Output:
Twin prime pairs below 10: 2 Twin prime pairs below 100: 8 Twin prime pairs below 1000: 35 Twin prime pairs below 10000: 205 Twin prime pairs below 100000: 1224 Twin prime pairs below 1000000: 8169 Twin prime pairs below 10000000: 58980
Wren
<lang ecmascript>import "/math" for Int import "/fmt" for Fmt
var c = Int.primeSieve(1e8-1, false) var limit = 10 var start = 3 var twins = 0 for (i in 1..8) {
var j = start
while (j < limit) {
if (!c[j] && !c[j-2]) twins = twins + 1
j = j + 2
}
Fmt.print("Under $,11d there are $,7d pairs of twin primes.", limit, twins)
start = limit + 1
limit = limit * 10
}</lang>
- Output:
Under 100 there are 8 pairs of twin primes. Under 1,000 there are 35 pairs of twin primes. Under 10,000 there are 205 pairs of twin primes. Under 100,000 there are 1,224 pairs of twin primes. Under 1,000,000 there are 8,169 pairs of twin primes. Under 10,000,000 there are 58,980 pairs of twin primes. Under 100,000,000 there are 440,312 pairs of twin primes.
Yabasic
<lang yabasic> sub isPrime(v)
if v < 2 then return False : fi
if mod(v, 2) = 0 then return v = 2 : fi
if mod(v, 3) = 0 then return v = 3 : fi
d = 5
while d * d <= v
if mod(v, d) = 0 then return False else d = d + 2 : fi
wend
return True
end sub
sub paresDePrimos(limite)
p1 = 0 : p2 = 1 : p3 = 1 : count = 0
for i = 5 to limite
p3 = p2
p2 = p1
p1 = isPrime(i)
if (p3 and p1) then count = count + 1 : fi
next i
return count
end sub
n = 1 for i = 1 to 6
n = n * 10 print "pares de primos gemelos por debajo de < ", n, " : ", paresDePrimos(n)
next i end </lang>
- Output:
Igual que la entrada de FreeBASIC.