Even numbers which cannot be expressed as the sum of two twin primes
In 1742, Goldbach and Euler corresponded about what is now known as the Goldbach Conjecture: that all even integers greater than 2 can be written as the sum of two primes.
In 2000, Henry Dubner proposed a stronger conjecture: that every even integer greater than 4208 can be written as the sum of two twin primes.
At the time the Goldbach conjecture was made, 1 was considered a prime number - this is not so now. So, the Goldbach conjecture was originally that all even natural numbers could be written as the sum of two primes.
- Task
- Find and display the positive even integers that cannot be expressed as the sum of two twin primes, up to a limit of 5,000.
- E.g.: The first 3 twin primes are 3, 5 and 7, so 6 (3+3), 8 (5+3) and 10 (5+5 or 7+3) can be formed but 2 and 4 cannot.
- Show the numbers that cannot be formed by summing two twin primes when 1 is treated as a prime (and so a twin prime).
- Note
- For the purposes of this task, twin prime refers to a prime that is 2 more or less than another prime, not the pair of primes.
- Stretch
- Verify that there no more such numbers up to 10,000 or more, as you like (note it has been verified up to at least 10^9).
- See also
ALGOL 68
BEGIN # find even numbers that are not the sum of two twin primes #
# prints the non-twin prime sum numbers #
PROC show non twin prime sums = ( STRING legend, []BOOL p2sums )VOID:
BEGIN
print( ( legend, ":", newline ) );
INT count := 0;
INT per line = 10;
FOR i FROM 2 BY 2 TO max number DO
IF NOT p2sum[ i ] THEN
print( ( whole( i, - 6 ) ) );
IF ( count +:= 1 ) MOD per line = 0 THEN print( ( newline ) ) FI
FI
OD;
IF count MOD per line /= 0 THEN print( ( newline ) ) FI;
print( ( "Found ", whole( count, 0 ), newline, newline ) )
END # print non twin primes sums # ;
INT max number = 100 000; # maximum number we will consider #
[ 0 : max number ]BOOL prime; # sieve the primes to max number #
prime[ 0 ] := prime[ 1 ] := FALSE;
prime[ 2 ] := TRUE;
FOR i FROM 3 BY 2 TO UPB prime DO prime[ i ] := TRUE OD;
FOR i FROM 4 BY 2 TO UPB prime DO prime[ i ] := FALSE OD;
FOR i FROM 3 BY 2 TO ENTIER sqrt( UPB prime ) DO
IF prime[ i ] THEN
FOR s FROM i * i BY i + i TO UPB prime DO prime[ s ] := FALSE OD
FI
OD;
prime[ 2 ] := FALSE; # restrict the sieve to twin primes only #
FOR i FROM 3 BY 2 TO max number - 2 DO
IF prime[ i ] AND NOT prime[ i - 2 ] AND NOT prime[ i + 2 ] THEN
prime[ i ] := FALSE
FI
OD;
# construct a table of the even numbers that can be formed by summing #
# two twin primes #
[ 1 : max number ]BOOL p2sum; FOR i TO UPB p2sum DO p2sum[ i ] := FALSE OD;
FOR i FROM 3 BY 2 TO max number DO
IF prime[ i ] THEN
FOR j FROM i TO ( max number + 1 ) - i DO
IF prime[ j ] THEN p2sum[ i + j ] := TRUE FI
OD
FI
OD;
# print the even numbers which aren't the sum of 2 twin primes #
show non twin prime sums( "Non twin prime sums", p2sum );
# adjust the table as if 1 was a twin prime #
p2sum[ 2 ] := TRUE;
FOR i FROM 3 BY 2 TO max number DO
IF prime[ i ] THEN p2sum[ i + 1 ] := TRUE FI
OD;
# print the revised even numbers which aren't the sum of 2 twin primes #
show non twin prime sums( "Non twin prime sums (including 1)", p2sum )
END
- Output:
Non twin prime sums: 2 4 94 96 98 400 402 404 514 516 518 784 786 788 904 906 908 1114 1116 1118 1144 1146 1148 1264 1266 1268 1354 1356 1358 3244 3246 3248 4204 4206 4208 Found 35 Non twin prime sums (including 1): 94 96 98 400 402 404 514 516 518 784 786 788 904 906 908 1114 1116 1118 1144 1146 1148 1264 1266 1268 1354 1356 1358 3244 3246 3248 4204 4206 4208 Found 33
J
For this task, it's convenient to have a routine to generate a sequence of twin primes up to a given limit:
twp=: [:(#~ ] e. (+,-)&2)i.&.(p:inv)
With this, the task becomes:
(2*1+i.2500)-.,+/~twp 5000
2 4 94 96 98 400 402 404 514 516 518 784 786 788 904 906 908 1114 1116 1118 1144 1146 1148 1264 1266 1268 1354 1356 1358 3244 3246 3248 4204 4206 4208
(2*1+i.2500)-.,+/~1,twp 5000 NB. including 1 in list of twin primes
94 96 98 400 402 404 514 516 518 784 786 788 904 906 908 1114 1116 1118 1144 1146 1148 1264 1266 1268 1354 1356 1358 3244 3246 3248 4204 4206 4208
For the stretch goal, we count primes and compare the count with the count we get with a raised limit:
#(2*1+i.2500)-.,+/~twp 5000
35
#(2*1+i.2500)-.,+/~twp 1000000
35
Julia
using Primes
""" Even numbers which cannot be expressed as the sum of two twin primes """
function nonpairsums(;include1=false, limit=20_000)
tpri = primesmask(limit + 2)
for i in 1:limit
tpri[i] && (i > 2 && !tpri[i - 2]) && !tpri[i + 2] && (tpri[i] = false)
end
tpri[2] = false
include1 && (tpri[1] = true)
twinsums = falses(limit * 2)
for i in 1:limit-2, j in 1:limit-2
if tpri[i] && tpri[j]
twinsums[i + j] = true
end
end
return [i for i in 2:2:limit if !twinsums[i]]
end
println("Non twin prime sums:")
foreach(p -> print(lpad(p[2], 6), p[1] % 10 == 0 ? "\n" : ""), pairs(nonpairsums()))
println("\n\nNon twin prime sums (including 1):")
foreach(p -> print(lpad(p[2], 6), p[1] % 10 == 0 ? "\n" : ""), pairs(nonpairsums(include1 = true)))
- Output:
Non twin prime sums: 2 4 94 96 98 400 402 404 514 516 518 784 786 788 904 906 908 1114 1116 1118 1144 1146 1148 1264 1266 1268 1354 1356 1358 3244 3246 3248 4204 4206 4208 Non twin prime sums (including 1): 94 96 98 400 402 404 514 516 518 784 786 788 904 906 908 1114 1116 1118 1144 1146 1148 1264 1266 1268 1354 1356 1358 3244 3246 3248 4204 4206 4208
Phix
with javascript_semantics constant lim = 10000 sequence prime = repeat(false,lim), p2sum = repeat(false,lim) for p in get_primes_le(lim)[2..$] do prime[p] = true end for -- restrict to twin primes only for i=3 to lim-2 by 2 do if prime[i] and not prime[i-2] and not prime[i+2] then prime[i] := false end if end for -- flag all summables for i=3 to lim by 2 do if prime[i] then for j=i to lim+1-i do if prime[j] then p2sum[i+j] = true end if end for end if end for sequence npts = {}, npts1 = {} for i=2 to lim by 2 do if not p2sum[i] then npts &= i end if end for -- adjust table as if 1 was a twin prime p2sum[2] = true for i=3 to lim by 2 do if prime[i] then p2sum[i+1] = true end if end for for i=2 to lim by 2 do if not p2sum[i] then npts1 &= i end if end for string fmt = "Non twin prime sums%s:\n%sFound %d\n\n", ntps = join_by(npts,1,10,"",fmt:="%6d"), ntps1 = join_by(npts1,1,10,"",fmt:="%6d") printf(1,fmt,{"",ntps,length(npts)}) printf(1,fmt,{" (including 1)",ntps1,length(npts1)})
- Output:
Non twin prime sums: 2 4 94 96 98 400 402 404 514 516 518 784 786 788 904 906 908 1114 1116 1118 1144 1146 1148 1264 1266 1268 1354 1356 1358 3244 3246 3248 4204 4206 4208 Found 35 Non twin prime sums (including 1): 94 96 98 400 402 404 514 516 518 784 786 788 904 906 908 1114 1116 1118 1144 1146 1148 1264 1266 1268 1354 1356 1358 3244 3246 3248 4204 4206 4208 Found 33
PL/M
... under CP/M (or an emulator)
100H: /* FIND EVEN NUMBERS THAT ARE NOT THE SUM OF TWO TWIN PRIMES */
/* CP/M SYSTEM CALL AND I/O ROUTINES */
BDOS: PROCEDURE( FN, ARG ); DECLARE FN BYTE, ARG ADDRESS; GOTO 5; END;
PR$CHAR: PROCEDURE( C ); DECLARE C BYTE; CALL BDOS( 2, C ); END;
PR$STRING: PROCEDURE( S ); DECLARE S ADDRESS; CALL BDOS( 9, S ); END;
PR$NL: PROCEDURE; CALL PR$CHAR( 0DH ); CALL PR$CHAR( 0AH ); END;
PR$NUMBER: PROCEDURE( N ); /* PRINTS A NUMBER IN THE MINIMUN FIELD WIDTH */
DECLARE N ADDRESS;
DECLARE V ADDRESS, N$STR ( 6 )BYTE, W BYTE;
V = N;
W = LAST( N$STR );
N$STR( W ) = '$';
N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
DO WHILE( ( V := V / 10 ) > 0 );
N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
END;
CALL PR$STRING( .N$STR( W ) );
END PR$NUMBER;
/* TASK */
DECLARE MAX$NUMBER LITERALLY '5001'
, FALSE LITERALLY '0'
, TRUE LITERALLY '0FFH'
;
/* PRINT THE NUMBERS THAT ARE NOT THE SUM OF TWO TWIN PRIMES */
SHOW$NON$TWIN$PRIME$SUMS: PROCEDURE( LEGEND, P2$PTR );
DECLARE ( LEGEND, P2$PTR ) ADDRESS;
DECLARE P2SUM BASED P2$PTR ( 0 )BYTE;
DECLARE ( COUNT, I )ADDRESS;
DECLARE PER$LINE LITERALLY '10';
CALL PR$STRING( LEGEND );CALL PR$CHAR( ':' );CALL PR$NL;
COUNT = 0;
DO I = 2 TO MAX$NUMBER - 1 BY 2;
IF NOT P2SUM( I ) THEN DO;
CALL PR$CHAR( ' ' );
IF I < 1000 THEN CALL PR$CHAR( ' ' );
IF I < 100 THEN CALL PR$CHAR( ' ' );
IF I < 10 THEN CALL PR$CHAR( ' ' );
CALL PR$CHAR( ' ' );
CALL PR$NUMBER( I );
IF ( COUNT := COUNT + 1 ) MOD PER$LINE = 0 THEN CALL PR$NL;
END;
END;
IF COUNT MOD PER$LINE <> 0 THEN CALL PR$NL;
CALL PR$STRING( .'FOUND $' );CALL PR$NUMBER( COUNT );CALL PR$NL;
END SHOW$NON$TWIN$PRIME$SUMS ;
DECLARE PRIME ( MAX$NUMBER )BYTE; /* SIEVE THE PRIMES TO MAX$NUMBER - 1 */
DECLARE I ADDRESS;
PRIME( 0 ), PRIME( 1 ) = FALSE;
PRIME( 2 ) = TRUE;
DO I = 3 TO LAST( PRIME ) BY 2; PRIME( I ) = TRUE; END;
DO I = 4 TO LAST( PRIME ) BY 2; PRIME( I ) = FALSE; END;
DO I = 3 TO LAST( PRIME ) BY 2;
IF PRIME( I ) THEN DO;
DECLARE S ADDRESS;
DO S = I + I TO LAST( PRIME ) BY I;
PRIME( S ) = FALSE;
END;
END;
END;
PRIME( 2 ) = FALSE; /* RESTRICT THE SIEVE TO TWIN PRIMES ONLY */
DO I = 3 TO LAST( PRIME ) - 2 BY 2;
IF PRIME( I ) AND NOT PRIME( I - 2 ) AND NOT PRIME( I + 2 ) THEN DO;
PRIME( I ) = FALSE;
END;
END;
/* CONSTRUCT A TABLE OF EVEN NUMBERS THAT CAN BE FORMED BY SUMMING TWO */
/* TWIN PRIMES */
DECLARE P2SUM ( MAX$NUMBER )BYTE;
DO I = 0 TO LAST( P2SUM ); P2SUM( I ) = FALSE; END;
DO I = 3 TO LAST( PRIME ) BY 2;
IF PRIME( I ) THEN DO;
DECLARE J ADDRESS;
DO J = I TO LAST( P2SUM ) - I;
IF PRIME( J ) THEN P2SUM( I + J ) = TRUE;
END;
END;
END;
/* PRINT THE EVEN NUMBERS THAT AREN'T THE SUM OF TWO TWIN PRIMES */
CALL SHOW$NON$TWIN$PRIME$SUMS( .'NON TWIN PRIME SUMS$', .P2SUM );
CALL PR$NL;
/* ADJUST THE TABLE AS IF 1 WAS A TWIN PRIME */
P2SUM( 2 ) = TRUE;
DO I = 3 TO LAST( P2SUM ) - 1 BY 2;
IF PRIME( I ) THEN P2SUM( I + 1 ) = TRUE;
END;
/* PRINT THE REVISED EVEN NUMBERS THAT AREN'T THE SUM OF TWO TWIN PRIMES */
CALL SHOW$NON$TWIN$PRIME$SUMS( .'NON TWIN PRIME SUMS (WITH 1 AS PRIME)$'
, .P2SUM
);
EOF
- Output:
NON TWIN PRIME SUMS: 2 4 94 96 98 400 402 404 514 516 518 784 786 788 904 906 908 1114 1116 1118 1144 1146 1148 1264 1266 1268 1354 1356 1358 3244 3246 3248 4204 4206 4208 FOUND 35 NON TWIN PRIME SUMS (WITH 1 AS PRIME): 94 96 98 400 402 404 514 516 518 784 786 788 904 906 908 1114 1116 1118 1144 1146 1148 1264 1266 1268 1354 1356 1358 3244 3246 3248 4204 4206 4208 FOUND 33
Python
''' rosettacode.org/wiki/Even_numbers_which_cannot_be_expressed_as_the_sum_of_two_twin_primes '''
from sympy import sieve
def nonpairsums(include1=False, limit=20_000):
""" Even numbers which cannot be expressed as the sum of two twin primes """
tpri = [i in sieve and (i - 2 in sieve or i + 2 in sieve)
for i in range(limit+2)]
if include1:
tpri[1] = True
twinsums = [False] * (limit * 2)
for i in range(limit):
for j in range(limit-i+1):
if tpri[i] and tpri[j]:
twinsums[i + j] = True
return [i for i in range(2, limit+1, 2) if not twinsums[i]]
print('Non twin prime sums:')
for k, p in enumerate(nonpairsums()):
print(f'{p:6}', end='\n' if (k + 1) % 10 == 0 else '')
print('\n\nNon twin prime sums (including 1):')
for k, p in enumerate(nonpairsums(include1=True)):
print(f'{p:6}', end='\n' if (k + 1) % 10 == 0 else '')
- Output:
Non twin prime sums: 2 4 94 96 98 400 402 404 514 516 518 784 786 788 904 906 908 1114 1116 1118 1144 1146 1148 1264 1266 1268 1354 1356 1358 3244 3246 3248 4204 4206 4208 Non twin prime sums (including 1): 94 96 98 400 402 404 514 516 518 784 786 788 904 906 908 1114 1116 1118 1144 1146 1148 1264 1266 1268 1354 1356 1358 3244 3246 3248 4204 4206 4208
Raku
my $threshold = 10000;
my @twins = unique flat (3..$threshold).grep(&is-prime).map: { $_, $_+2 if ($_+2).is-prime };
my @sums;
map {++@sums[$_]}, (@twins X+ @twins);
display 'Non twin prime sums:', @sums;
map {++@sums[$_]}, flat 2, (1 X+ @twins);
display "\nNon twin prime sums (including 1):", @sums;
sub display ($msg, @p) {
put "$msg\n" ~ @p[^$threshold].kv.map(-> $k, $v { $k if ($k %% 2) && !$v })\
.skip(1).batch(10)».fmt("%4d").join: "\n"
}
- Output:
Non twin prime sums: 2 4 94 96 98 400 402 404 514 516 518 784 786 788 904 906 908 1114 1116 1118 1144 1146 1148 1264 1266 1268 1354 1356 1358 3244 3246 3248 4204 4206 4208 Non twin prime sums (including 1): 94 96 98 400 402 404 514 516 518 784 786 788 904 906 908 1114 1116 1118 1144 1146 1148 1264 1266 1268 1354 1356 1358 3244 3246 3248 4204 4206 4208
Wren
import "./math" for Int
import "./fmt" for Fmt
var limit = 100000 // say
var primes = Int.primeSieve(limit)[2..-1] // exclude 2 & 3
var twins = [3]
for (i in 0...primes.count-1) {
if (primes[i+1] - primes[i] == 2) {
if (twins[-1] != primes[i]) twins.add(primes[i])
twins.add(primes[i+1])
}
}
var nonTwinSums = Fn.new {
var sieve = List.filled(limit+1, false)
for (i in 0...twins.count) {
for (j in i...twins.count) {
var sum = twins[i] + twins[j]
if (sum > limit) break
sieve[sum] = true
}
}
var res = []
var i = 2
while (i < limit) {
if (!sieve[i]) res.add(i)
i = i + 2
}
return res
}
System.print("Non twin prime sums:")
var ntps = nonTwinSums.call()
Fmt.tprint("$4d", ntps, 10)
System.print("Found %(ntps.count)")
System.print("\nNon twin prime sums (including 1):")
twins.insert(0, 1)
ntps = nonTwinSums.call()
Fmt.tprint("$4d", ntps, 10)
System.print("Found %(ntps.count)")
- Output:
Non twin prime sums: 2 4 94 96 98 400 402 404 514 516 518 784 786 788 904 906 908 1114 1116 1118 1144 1146 1148 1264 1266 1268 1354 1356 1358 3244 3246 3248 4204 4206 4208 Found 35 Non twin prime sums (including 1): 94 96 98 400 402 404 514 516 518 784 786 788 904 906 908 1114 1116 1118 1144 1146 1148 1264 1266 1268 1354 1356 1358 3244 3246 3248 4204 4206 4208 Found 33