Even numbers which cannot be expressed as the sum of two twin primes

Even numbers which cannot be expressed as the sum of two twin primes is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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

• Wikipedia: Dubner's conjecture
• Wikipedia: Goldbach's conjecture
• OEIS:A007534 - Even numbers that are not the sum of a pair of twin primes.
• Twin primes

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

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

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

Library: Wren-math

Library: Wren-fmt

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