Goldbach's comet

Goldbach's comet is the name given to a plot of the function g(E), the so-called Goldbach function.

Task
Goldbach's comet
You are encouraged to solve this task according to the task description, using any language you may know.
This page uses content from Wikipedia. The original article was at Goldbach's comet. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)

The Goldbach function is studied in relation to Goldbach's conjecture. The function g(E) is defined for all even integers E>2 to be the number of different ways in which E can be expressed as the sum of two primes.

Examples
  • G(4) = 1, since 4 can only be expressed as the sum of one distinct pair of primes (4 = 2 + 2)
  • G(22) = 3, since 22 can be expressed as the sum of 3 distinct pairs of primes (22 = 11 + 11 = 5 + 17 = 3 + 19)


Task
  • Find and show (preferably, in a neat 10x10 table) the first 100 G numbers (that is: the result of the G function described above, for the first 100 even numbers >= 4)
  • Find and display the value of G(1000000)


Stretch
  • Calculate the values of G up to 2000 (inclusive) and display the results in a scatter 2d-chart, aka the Goldbach's Comet


See also



ALGOL 68Edit

Generates an ASCII-Art scatter plot - the vertical axis is n/10 and the hotizontal is G(n).

BEGIN # calculate values of the Goldbach function G where G(n) is the number #
      # of prime pairs that sum to n, n even and > 2                         #
      # generates an ASCII scatter plot of G(n) up to G(2000)                #
      # (Goldbach's Comet)                                                   #
    PR read "primes.incl.a68" PR          # include prime utilities          #
    INT max prime = 1 000 000;            # maximum number we will consider  #
    INT max plot  =     2 000;            # maximum G value for the comet    #
    []BOOL prime = PRIMESIEVE max prime;  # sieve of primes to max prime     #
    [ 0 : max plot ]INT g2;               # table of G values: g2[n] = G(2n) #
    # construct the table of G values #
    FOR n FROM LWB g2 TO UPB g2 DO g2[ n ] := 0 OD;
    g2[ 4 ] := 1;                     # 4 is the only sum of two even primes #
    FOR p FROM 3 BY 2 TO max plot OVER 2 DO
        IF prime[ p ] THEN
            g2[ p + p ] +:= 1;
            FOR q FROM p + 2 BY 2 TO max plot - p DO
                IF prime[ q ] THEN
                    g2[ p + q ] +:= 1
                FI
            OD
        FI
    OD;
    # show the first hundred G values #
    INT c := 0;
    FOR n FROM 4 BY 2 TO 202 DO
        print( ( whole( g2[ n ], -4 ) ) );
        IF ( c +:= 1 ) = 10 THEN print( ( newline ) ); c := 0 FI
    OD;
    # show G( 1 000 000 ) #
    INT gm := 0;
    FOR p FROM 3 TO max prime OVER 2 DO
        IF prime[ p ] THEN
            IF prime[ max prime - p ] THEN
                gm +:= 1
            FI
        FI
    OD;
    print( ( "G(", whole( max prime, 0 ), "): ", whole( gm, 0 ), newline ) );
    # find the maximum value of G up to the maximum plot size #
    INT max g := 0;
    FOR n FROM 2 BY 2 TO max plot DO
        IF g2[ n ] > max g THEN max g := g2[ n ] FI
    OD;
    # draw an ASCII scatter plot of G, each position represents 5 G values #
    # the vertical axis is n/10, the horizontal axis is G(n) #
    INT plot step = 10;
    STRING plot value = " .-+=*%$&#@";
    FOR g FROM 0 BY plot step TO max plot - plot step DO
        [ 0 : max g ]INT values;
        FOR v pos FROM LWB values TO UPB values DO values[ v pos ] := 0 OD;
        INT max v := 0;
        FOR g element FROM g BY 2 TO g + ( plot step - 1 ) DO
            INT g2 value = g2[ g element ];
            values[ g2 value ] +:= 1;
            IF g2 value > max v THEN max v := g2 value FI
        OD;
        print( ( IF g MOD 100 = 90 THEN "+" ELSE "|" FI ) );
        FOR v pos FROM 1 TO max v DO # exclude 0 values from the plot #
            print( ( plot value[ values[ v pos ] + 1 ] ) )
        OD;
        print( ( newline ) )
    OD
END
Output:
   1   1   1   2   1   2   2   2   2   3
   3   3   2   3   2   4   4   2   3   4
   3   4   5   4   3   5   3   4   6   3
   5   6   2   5   6   5   5   7   4   5
   8   5   4   9   4   5   7   3   6   8
   5   6   8   6   7  10   6   6  12   4
   5  10   3   7   9   6   5   8   7   8
  11   6   5  12   4   8  11   5   8  10
   5   6  13   9   6  11   7   7  14   6
   8  13   5   8  11   7   9  13   8   9
G(1000000): 5402

|+ |.= | -+ | -.- | --. | --. | .. .- | +.. | -- . + ... . . | .- - | +. . | ... . . | ..... | .... . | .. . .. | .. . . . | .- . . | .. . .. + ... . . | -.. . | ... . . | - . .. | . + . | . .. . . | .- . . | .-. . | - . . . | . . . .. + - .. . | ... . . | .. . . . | . .- . | . . . . . | .. . .. | .- . . | - . . . | .. . . . | .. .. . + . - . . | . .. . . | .- . . | . .. . . | . .. . . | .-. . | . - . . | .. . . . | . ... . | . . . . . + - . .. | ..- . | - . . . | . . . .. | .- . . | .. . - | - . - | -. . . | .. . . . | .. . . . + . -. . | ... . . | . .. - | . ... . | . .. . . | ... . . | . . . . . | . . . . . | . . . . . | + . . + . . . . . | .. - . | .... . | .- . . | . . . . . | -- . | .. . . . | - . .. | .. . . . | .. . . . + . . . . . | . - . . | .. . . . | . .. . . | .. - . | . .. . . | .. . . . | -. . . | . .. . . | . . . . . + -. . . | - . . . | . . - . | - . . . | .. . . . | . . . .. | . .. . . | ... - | . . . . . | . . . . . + . . . . . | - . . . | . . . . . | - . . . | . . . . . | - . . . | .. . . . | - . . . | .. . . . | . .. . . + . . . . . | .- . . | . .. . . | . . . . . | . .. . . | . . . . . | . . . . . | . .. . . | .. . . . | - . . . + .- . . | .. . . . | .. . . . | - . . . | . . . . . | . . . .. | .. . . . | .. . . . | ... . . | . .. . . + . . . . . | - . . . | . .. . . | .. . . . | - . . . | . . - . | .. . . . | - . . . | . -. . | - . . . + - . .. | . . . . . | . . . . . | .. . . . | .. . . . | . .. . . | . . . .. | .. . . . | . . . . . | . . . - + .. .. . | . - . . | . . . . . | -. . . | . . . . . | . . - . | .. . . . | . .. . . | . . . . . | . . . . . + .. . . . | . . . . . | . . . . . | . . . . . | - . . . | ... . . | .. . . . | . . . .. | . . . . . | . . . . . + . . . . . | - . . . | . .. . . | . . . . . | . . .. . | - . . . | . . . . . | . .. . . | .. . . . | . . . . . + . . - . | . - . . | .. . - | . . . . . | . . . . . | . . . . . | .. . . . | . . . . . | . . . . . | . . . . . + .. . . . | . - . . | . . .. . | - . . . | . . . . . | .. . . . | . . . . . | . . . . . | -. . . | - . . . + .. . . .

ArturoEdit

G: function [n][
    size select 2..n/2 'x ->
        and? [prime? x][prime? n-x]
]

print "The first 100 G values:"
loop split.every: 10 map select 4..202 => even? => G 'row [
    print map to [:string] row 'item -> pad item 3
]

print ["\nG(1000000) =" G 1000000]

csv: join.with:",\n" map select 4..2000 => even? 'x ->
    ~"|x|, |G x|"

; write the CSV data to a file which we can then visualize
; via our preferred spreadsheet app
write "comet.csv" csv
Output:
The first 100 G values:
  1   1   1   2   1   2   2   2   2   3 
  3   3   2   3   2   4   4   2   3   4 
  3   4   5   4   3   5   3   4   6   3 
  5   6   2   5   6   5   5   7   4   5 
  8   5   4   9   4   5   7   3   6   8 
  5   6   8   6   7  10   6   6  12   4 
  5  10   3   7   9   6   5   8   7   8 
 11   6   5  12   4   8  11   5   8  10 
  5   6  13   9   6  11   7   7  14   6 
  8  13   5   8  11   7   9  13   8   9 

G(1000000) = 5402

Here, you can find the result of the visualization - or the "Goldbach's comet": 2D-chart of all G values up to 2000

AutoHotkeyEdit

c := 0
while (c<100)
    if (x := g(A_Index))
        c++, result .= x (Mod(c, 10) ? "`t" : "`n")
MsgBox % result "`ng(1000000) : " g(1000000)
return

g(n, i:=1) {
    if Mod(n, 2)
        return false
    while (++i <= n/2)
        if (is_prime(i) && is_prime(n-i))
            count++
    return count
}

is_prime(N) {
    Loop, % Floor(Sqrt(N))
        if (A_Index > 1 && !Mod(N, A_Index))
            Return false
    Return true
}
Output:
1	1	1	2	1	2	2	2	2	3
3	3	2	3	2	4	4	2	3	4
3	4	5	4	3	5	3	4	6	3
5	6	2	5	6	5	5	7	4	5
8	5	4	9	4	5	7	3	6	8
5	6	8	6	7	10	6	6	12	4
5	10	3	7	9	6	5	8	7	8
11	6	5	12	4	8	11	5	8	10
5	6	13	9	6	11	7	7	14	6
8	13	5	8	11	7	9	13	8	9
g(1000000) : 5402

AWKEdit

# syntax: GAWK -f GOLDBACHS_COMET.AWK
BEGIN {
    print("The first 100 G numbers:")
    for (n=4; n<=202; n+=2) {
      printf("%4d%1s",g(n),++count%10?"":"\n")
    }
    n = 1000000
    printf("\nG(%d): %d\n",n,g(n))
    n = 4
    printf("G(%d): %d\n",n,g(n))
    n = 22
    printf("G(%d): %d\n",n,g(n))
    exit(0)
}
function g(n,  count,i) {
    if (n % 2 == 0) { # n must be even
      for (i=2; i<=(1/2)*n; i++) {
        if (is_prime(i) && is_prime(n-i)) {
          count++
        }
      }
    }
    return(count)
}
function is_prime(n,  d) {
    d = 5
    if (n < 2) { return(0) }
    if (n % 2 == 0) { return(n == 2) }
    if (n % 3 == 0) { return(n == 3) }
    while (d*d <= n) {
      if (n % d == 0) { return(0) }
      d += 2
      if (n % d == 0) { return(0) }
      d += 4
    }
    return(1)
}
Output:
The first 100 G numbers:
   1    1    1    2    1    2    2    2    2    3
   3    3    2    3    2    4    4    2    3    4
   3    4    5    4    3    5    3    4    6    3
   5    6    2    5    6    5    5    7    4    5
   8    5    4    9    4    5    7    3    6    8
   5    6    8    6    7   10    6    6   12    4
   5   10    3    7    9    6    5    8    7    8
  11    6    5   12    4    8   11    5    8   10
   5    6   13    9    6   11    7    7   14    6
   8   13    5    8   11    7    9   13    8    9

G(1000000): 5402
G(4): 1
G(22): 3

BASICEdit

FreeBASICEdit

Function isPrime(Byval ValorEval As Uinteger) 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 g(n As Uinteger) As Uinteger
    Dim As Uinteger i, count = 0
    If (n Mod 2 = 0) Then     'n in goldbach function g(n) must be even
        For i = 2 To (1/2) * n
            If isPrime(i) And isPrime(n - i) Then count += 1
        Next i
    End If
    Return count
End Function

Print "The first 100 G numbers are:"

Dim As Uinteger col = 1
For n As Uinteger = 4 To 202 Step 2
    Print Using "####"; g(n);
    If (col Mod 10 = 0) Then Print
    col += 1
Next n

Print !"\nThe value of G(1000000) is "; g(1000000)
Sleep
Output:
The first 100 G numbers are:
   1   1   1   2   1   2   2   2   2   3
   3   3   2   3   2   4   4   2   3   4
   3   4   5   4   3   5   3   4   6   3
   5   6   2   5   6   5   5   7   4   5
   8   5   4   9   4   5   7   3   6   8
   5   6   8   6   7  10   6   6  12   4
   5  10   3   7   9   6   5   8   7   8
  11   6   5  12   4   8  11   5   8  10
   5   6  13   9   6  11   7   7  14   6
   8  13   5   8  11   7   9  13   8   9

The value of G(1000000) is 5402

SmileBASICEdit

(click the image to view full-size)

OPTION STRICT: OPTION DEFINT
VAR MAX_G = 4000, MAX_P = 1000000
VAR ROOT_MAX_P = FLOOR(SQR(MAX_P))
VAR HALF_MAX_G = MAX_G DIV 2
VAR G[MAX_G + 1], P[MAX_P + 1]
VAR I, J
CLS: GCLS
P[0] = FALSE: P[1] = 0: P[2] = TRUE
FOR I = 4 TO MAX_P STEP 2 P[I] = FALSE: NEXT
FOR I = 3 TO MAX_P STEP 2 P[I] = TRUE: NEXT
FOR I = 3 TO ROOT_MAX_P STEP 2
  IF P[I] THEN
    FOR J = I * I TO MAX_P STEP I
      P[J] = FALSE
    NEXT
  ENDIF
NEXT
FOR I = 1 TO MAX_G G[I] = 0: NEXT
G[4] = 1 ' 4 is the only sum of 2 even primes
FOR I = 3 TO HALF_MAX_G STEP 2
  IF P[I] THEN
    INC G[I + 1]
    FOR J = I + 2 TO MAX_G - 1
      IF P[J] THEN
        INC G[I + 1]
      ENDIF
    NEXT
  ENDIF
NEXT
VAR C = 0
FOR I = 4 TO 202 STEP 2
  PRINT FORMAT$("%3D", G[I]),
  INC C
  IF C == 10 THEN PRINT: C = 0: ENDIF
NEXT
VAR GM = 0
FOR I = 3 TO MAX_P DIV 2 STEP 2
  IF P[I] THEN
    IF P[MAX_P - I] THEN INC GM: ENDIF
  ENDIF
NEXT
PRINT FORMAT$("G(%D): ", MAX_P); GM
FOR I = 2 TO MAX_G - 10 STEP 10
  FOR J = 1 TO I + 8 STEP 2
    GPSET I DIV 10, 240-G[J],RGB(255,255,255)
  NEXT
NEXT

JEdit

Task implementation:

   10 10$#/.~4,/:~ 0-.~,(<:/~ * +/~) p:1+i.p:inv 202
 1  1  1  2  1  2  2  2  2  3
 3  3  2  3  2  4  4  2  3  4
 3  4  5  4  3  5  3  4  6  3
 5  6  2  5  6  5  5  7  4  5
 8  5  4  9  4  5  7  3  6  8
 5  6  8  6  7 10  6  6 12  4
 5 10  3  7  9  6  5  8  7  8
11  6  5 12  4  8 11  5  8 10
 5  6 13  9  6 11  7  7 14  6
 8 13  5  8 11  7  9 13  8  9

And, for G(1e6):

   -:+/1 p: 1e6-p:i.p:inv 1e6
5402

Explanation:

For the first part, the easiest approach seems to be to brute force it. The first 100 numbers starting with 4 end at 202, so we start by finding all primes less than 202 (and sum all pairs of these primes).

Also, we can eliminate an odd/even test on these sums of primes by excluding 2 from our list of primes and including 4 as an explicit constant.

Also, since we are only concerned about distinct sums, we eliminate all pairs where the first prime in the sum exceeds the second prime in the sum.

So.. we compute a couple thousand sums, sort them in numeric order (prepending that list with 4), count how many times each unique value appears, and order the first 100 of those frequencies in a 10x10 table.

---

For G(1e6) that brute force approach becomes inefficient -- instead of about 2000 sums, almost 600 of which are relevant, we would need to find over 6e9 sums, and about 6e9 of them would be irrelevant.

Instead, for G(1e6), we find all primes less than a million, subtract each from 1 million and count how many of the differences are prime and cut that in half. We cut that sum in half because this approach counts each pair twice (once with the smallest value first, again with the smallest value second -- since 1e6 is not the square of a prime we do not have a prime which appears twice in one of these sums).

JuliaEdit

Run in VS Code or REPL to view and save the plot.

using Combinatorics
using Plots
using Primes

g(n) = iseven(n) ? count(p -> all(isprime, p), partitions(n, 2)) : error("n must be even")

println("The first 100 G numbers are: ")

foreach(p -> print(lpad(p[2], 4), p[1] % 10 == 0 ? "\n" : ""), map(g, 4:2:202) |> enumerate)

println("\nThe value of G(1000000) is ", g(1_000_000))

x = collect(2:2002)
y = map(g, 2x)
scatter(x, y, markerstrokewidth = 0, color = ["red", "blue", "green"][mod1.(x, 3)])
Output:

The first 100 G numbers are:

  1   1   1   2   1   2   2   2   2   3
  3   3   2   3   2   4   4   2   3   4
  3   4   5   4   3   5   3   4   6   3
  5   6   2   5   6   5   5   7   4   5
  8   5   4   9   4   5   7   3   6   8
  5   6   8   6   7  10   6   6  12   4
  5  10   3   7   9   6   5   8   7   8
 11   6   5  12   4   8  11   5   8  10
  5   6  13   9   6  11   7   7  14   6
  8  13   5   8  11   7   9  13   8   9

The value of G(1000000) is 5402

 

LuaEdit

function T(t) return setmetatable(t, {__index=table}) end
table.range = function(t,n) local s=T{} for i=1,n do s[i]=i end return s end
table.map = function(t,f) local s=T{} for i=1,#t do s[i]=f(t[i]) end return s end
table.batch = function(t,n,f) for i=1,#t,n do local s=T{} for j=1,n do s[j]=t[i+j-1] end f(s) end return t end

function isprime(n)
  if n < 2 then return false end
  if n % 2 == 0 then return n==2 end
  if n % 3 == 0 then return n==3 end
  for f = 5, n^0.5, 6 do
    if n%f==0 or n%(f+2)==0 then return false end
  end
  return true
end

function goldbach(n)
  local count = 0
  for i = 1, n/2 do
    if isprime(i) and isprime(n-i) then
      count = count + 1
    end
  end
  return count
end

print("The first 100 G numbers:")
g = T{}:range(100):map(function(n) return goldbach(2+n*2) end)
g:map(function(n) return string.format("%2d ",n) end):batch(10,function(t) print(t:concat()) end)
print("G(1000000) = "..goldbach(1000000))
Output:
The first 100 G numbers:
 1  1  1  2  1  2  2  2  2  3
 3  3  2  3  2  4  4  2  3  4
 3  4  5  4  3  5  3  4  6  3
 5  6  2  5  6  5  5  7  4  5
 8  5  4  9  4  5  7  3  6  8
 5  6  8  6  7 10  6  6 12  4
 5 10  3  7  9  6  5  8  7  8
11  6  5 12  4  8 11  5  8 10
 5  6 13  9  6 11  7  7 14  6
 8 13  5  8 11  7  9 13  8  9
G(1000000) = 5402


Mathematica/Wolfram LanguageEdit

ClearAll[GoldbachFuncion]
GoldbachFuncion[e_Integer] := Module[{ps},
  ps = Prime[Range[PrimePi[e/2]]];
  Total[Boole[PrimeQ[e - ps]]]
]
Grid[Partition[GoldbachFuncion /@ Range[4, 220, 2], 10]]
GoldbachFuncion[10^6]
DiscretePlot[GoldbachFuncion[e], {e, 4, 2000}, Filling -> None]
Output:
1	1	1	2	1	2	2	2	2	3
3	3	2	3	2	4	4	2	3	4
3	4	5	4	3	5	3	4	6	3
5	6	2	5	6	5	5	7	4	5
8	5	4	9	4	5	7	3	6	8
5	6	8	6	7	10	6	6	12	4
5	10	3	7	9	6	5	8	7	8
11	6	5	12	4	8	11	5	8	10
5	6	13	9	6	11	7	7	14	6
8	13	5	8	11	7	9	13	8	9

5402

[graphical object showing Goldbach's comet]


PerlEdit

Library: ntheory
use strict;
use warnings;
use feature 'say';

use List::Util 'max';
use GD::Graph::bars;
use ntheory 'is_prime';

sub table { my $t = shift() * (my $c = 1 + max map {length} @_); ( sprintf( ('%'.$c.'s')x@_, @_) ) =~ s/.{1,$t}\K/\n/gr }

sub G {
    my($n) = @_;
    scalar grep { is_prime($_) and is_prime($n - $_) } 2 .. $n/2;
}

my @y;
push @y, G(2*$_ + 4) for my @x = 0..1999;

say $_ for table 10, @y;
printf "G $_: %d", G($_) for 1e6;

my @data = ( \@x, \@y);
my $graph = GD::Graph::bars->new(1200, 400);
$graph->set(
    title          => q/Goldbach's Comet/,
    y_max_value    => 170,
    x_tick_number  => 10,
    r_margin       => 10,
    dclrs          => [ 'blue' ],
) or die $graph->error;
my $gd = $graph->plot(\@data) or die $graph->error;

open my $fh, '>', 'goldbachs-comet.png';
binmode $fh;
print $fh $gd->png();
close $fh;
Output:
  1  1  1  2  1  2  2  2  2  3
  3  3  2  3  2  4  4  2  3  4
  3  4  5  4  3  5  3  4  6  3
  5  6  2  5  6  5  5  7  4  5
  8  5  4  9  4  5  7  3  6  8
  5  6  8  6  7 10  6  6 12  4
  5 10  3  7  9  6  5  8  7  8
 11  6  5 12  4  8 11  5  8 10
  5  6 13  9  6 11  7  7 14  6
  8 13  5  8 11  7  9 13  8  9

G 1000000: 5402

Stretch goal: (offsite image) goldbachs-comet.png

PhixEdit

Library: Phix/pGUI
Library: Phix/online

You can run this online here.

--
-- demo\rosetta\Goldbachs_comet.exw
-- ================================
--
-- Note: this plots n/2 vs G(n) for n=6 to 4000 by 2, matching wp and
--       Algol 68, Python, and Raku. However, while not wrong, Arturo 
--       and Wren apparently plot n vs G(n) for n=6 to 2000 by 2, so
--       should you spot any (very) minor differences, that'd be why.
--
with javascript_semantics
requires("1.0.2")

constant limit = 4000
sequence primes = get_primes_le(limit)[2..$],
       goldbach = reinstate(repeat(0,limit),{4},{1})
for i=1 to length(primes) do
    for j=i to length(primes) do
        integer s = primes[i] + primes[j]
        if s>limit then exit end if
        goldbach[s] += 1
    end for
end for
 
sequence fhg = extract(goldbach,tagstart(4,100,2))
string fhgs = join_by(fhg,1,10," ","\n","%2d")
integer gm = sum(apply(sq_sub(1e6,get_primes_le(499999)[2..$]),is_prime))
printf(1,"The first 100 G values:\n%s\n\nG(1,000,000) = %,d\n",{fhgs,gm})

include pGUI.e
include IupGraph.e
 
function get_data(Ihandle /*graph*/)
    return {{tagset(limit/2,3,3),extract(goldbach,tagset(limit,6,6)),CD_RED},
            {tagset(limit/2,4,3),extract(goldbach,tagset(limit,8,6)),CD_BLUE},
            {tagset(limit/2,5,3),extract(goldbach,tagset(limit,10,6)),CD_DARK_GREEN}}
end function

IupOpen()
Ihandle graph = IupGraph(get_data,"RASTERSIZE=640x440,MARKSTYLE=PLUS")
IupSetAttributes(graph,"XTICK=%d,XMIN=0,XMAX=%d",{limit/20,limit/2})
IupSetAttributes(graph,"YTICK=20,YMIN=0,YMAX=%d",{max(goldbach)+20})
Ihandle dlg = IupDialog(graph,`TITLE="Goldbach's comet",MINSIZE=400x300`)
IupShow(dlg)
if platform()!=JS then 
    IupMainLoop()
    IupClose()
end if
Output:
The first 100 G values:
 1  1  1  2  1  2  2  2  2  3
 3  3  2  3  2  4  4  2  3  4
 3  4  5  4  3  5  3  4  6  3
 5  6  2  5  6  5  5  7  4  5
 8  5  4  9  4  5  7  3  6  8
 5  6  8  6  7 10  6  6 12  4
 5 10  3  7  9  6  5  8  7  8
11  6  5 12  4  8 11  5  8 10
 5  6 13  9  6 11  7  7 14  6
 8 13  5  8 11  7  9 13  8  9


G(1,000,000) = 5,402

PicatEdit

main =>
  println("First 100 G numbers:"),
  foreach({G,I} in zip(take([G: T in 1..300, G=g(T),G>0],100),1..100))
    printf("%2d %s",G,cond(I mod 10 == 0,"\n",""))
  end,
  nl,
  printf("G(1_000_000): %d\n", g(1_000_000)).

g(N) = cond((N > 2, N mod 2 == 0),
              {1 : I in 1..N // 2,
                 prime(I),prime(N-I)}.len,
              0).
Output:
First 100 G numbers:
 1  1  1  2  1  2  2  2  2  3 
 3  3  2  3  2  4  4  2  3  4 
 3  4  5  4  3  5  3  4  6  3 
 5  6  2  5  6  5  5  7  4  5 
 8  5  4  9  4  5  7  3  6  8 
 5  6  8  6  7 10  6  6 12  4 
 5 10  3  7  9  6  5  8  7  8 
11  6  5 12  4  8 11  5  8 10 
 5  6 13  9  6 11  7  7 14  6 
 8 13  5  8 11  7  9 13  8  9 

G(1_000_000): 5402

PythonEdit

from matplotlib.pyplot import scatter, show
from sympy import isprime

def g(n):
    assert n > 2 and n % 2 == 0, 'n in goldbach function g(n) must be even'          
    count = 0
    for i in range(1, n//2 + 1):
        if isprime(i) and isprime(n - i):
            count += 1
    return count

print('The first 100 G numbers are:')
 
col = 1
for n in range(4, 204, 2):
    print(str(g(n)).ljust(4), end = '\n' if (col % 10 == 0) else '')
    col += 1

print('\nThe value of G(1000000) is', g(1_000_000))
 
x = range(4, 4002, 2)
y = [g(i) for i in x]
colors = [["red", "blue", "green"][(i // 2) % 3] for i in x]
scatter([i // 2 for i in x], y, marker='.', color = colors)
show()
Output:
The first 100 G numbers are:
1   1   1   2   1   2   2   2   2   3   
3   3   2   3   2   4   4   2   3   4   
3   4   5   4   3   5   3   4   6   3   
5   6   2   5   6   5   5   7   4   5   
8   5   4   9   4   5   7   3   6   8   
5   6   8   6   7   10  6   6   12  4   
5   10  3   7   9   6   5   8   7   8   
11  6   5   12  4   8   11  5   8   10  
5   6   13  9   6   11  7   7   14  6   
8   13  5   8   11  7   9   13  8   9   

The value of G(1000000) is 5402

 

RakuEdit

For the stretch, actually generates a plot, doesn't just calculate values to be plotted by a third party application. Deviates slightly from the stretch goal, plots the first two thousand defined values rather than the values up to two thousand that happen to be defined. (More closely matches the wikipedia example image.)

sub G (Int $n) { +(2..$n/2).grep: { .is-prime && ($n - $_).is-prime } }

# Task
put "The first 100 G values:\n", (^100).map({ G 2 × $_ + 4 }).batch(10)».fmt("%2d").join: "\n";

put "\nG 1_000_000 = ", G 1_000_000;

# Stretch
use SVG;
use SVG::Plot;

my @x = map 2 × * + 4, ^2000;
my @y = @x.map: &G;

'Goldbachs-Comet-Raku.svg'.IO.spurt: SVG.serialize: SVG::Plot.new(
    width       => 1000,
    height      => 500,
    background  => 'white',
    title       => "Goldbach's Comet",
    x           => @x,
    values      => [@y,],
).plot: :points;
Output:
The first 100 G values:
 1  1  1  2  1  2  2  2  2  3
 3  3  2  3  2  4  4  2  3  4
 3  4  5  4  3  5  3  4  6  3
 5  6  2  5  6  5  5  7  4  5
 8  5  4  9  4  5  7  3  6  8
 5  6  8  6  7 10  6  6 12  4
 5 10  3  7  9  6  5  8  7  8
11  6  5 12  4  8 11  5  8 10
 5  6 13  9  6 11  7  7 14  6
 8 13  5  8 11  7  9 13  8  9

G 1_000_000 = 5402

Stretch goal: (offsite SVG image) Goldbachs-Comet-Raku.svg

RustEdit

// [dependencies]
// primal = "0.3"
// plotters = "0.3.2"

use plotters::prelude::*;

fn goldbach(n: u64) -> u64 {
    let mut p = 2;
    let mut count = 0;
    loop {
        let q = n - p;
        if q < p {
            break;
        }
        if primal::is_prime(p) && primal::is_prime(q) {
            count += 1;
        }
        if p == 2 {
            p += 1;
        } else {
            p += 2;
        }
    }
    count
}

fn goldbach_plot(filename: &str) -> Result<(), Box<dyn std::error::Error>> {
    let gvalues : Vec<u64> = (1..=2000).map(|x| goldbach(2 * x + 2)).collect();
    let mut gmax = *gvalues.iter().max().unwrap();
    gmax = 10 * ((gmax + 9) / 10);

    let root = SVGBackend::new(filename, (1000, 500)).into_drawing_area();
    root.fill(&WHITE)?;

    let mut chart = ChartBuilder::on(&root)
        .x_label_area_size(20)
        .y_label_area_size(20)
        .margin(10)
        .caption("Goldbach's Comet", ("sans-serif", 24).into_font())
        .build_cartesian_2d(0usize..2000usize, 0u64..gmax)?;

    chart
        .configure_mesh()
        .disable_x_mesh()
        .disable_y_mesh()
        .draw()?;

    chart.draw_series(
        gvalues
            .iter()
            .cloned()
            .enumerate()
            .map(|p| Circle::new(p, 2, BLUE.filled())),
    )?;

    Ok(())
}

fn main() {
    println!("First 100 G numbers:");
    for i in 1..=100 {
        print!(
            "{:2}{}",
            goldbach(2 * i + 2),
            if i % 10 == 0 { "\n" } else { " " }
        );
    }

    println!("\nG(1000000) = {}", goldbach(1000000));

    match goldbach_plot("goldbach.svg") {
        Ok(()) => {}
        Err(error) => eprintln!("Error: {}", error),
    }
}
Output:
First 100 G numbers:
 1  1  1  2  1  2  2  2  2  3
 3  3  2  3  2  4  4  2  3  4
 3  4  5  4  3  5  3  4  6  3
 5  6  2  5  6  5  5  7  4  5
 8  5  4  9  4  5  7  3  6  8
 5  6  8  6  7 10  6  6 12  4
 5 10  3  7  9  6  5  8  7  8
11  6  5 12  4  8 11  5  8 10
 5  6 13  9  6 11  7  7 14  6
 8 13  5  8 11  7  9 13  8  9

G(1000000) = 5402

Media:Goldbach's comet rust.svg

uBasic/4tHEdit

Translation of: FreeBASIC

For performance reasons only g(100000) is calculated. The value "810" has been verified and is correct.

Print "The first 100 G numbers are:"
c = 1

For n = 4 To 202 Step 2
  Print Using "___#"; FUNC(_g(n));
  If (c % 10) = 0 Then Print
  c = c + 1
Next

Print "\nThe value of G(100000) is "; FUNC(_g(100000))

End

_isPrime
  Param (1)
  Local (1)
  
  If a@ < 2 Then Return (0)
  For b@ = 2 To FUNC(_Sqrt(a@))
    If (a@ % b@) = 0 Then Unloop: Return (0)
  Next
Return (1)
 
_g
  Param (1)
  Local (2)
  
  c@ = 0
  
  If (a@ % 2) = 0 Then     'n in goldbach function g(n) must be even
    For b@ = 2 To a@/2
      If FUNC(_isPrime(b@)) * FUNC(_isPrime(a@ - b@)) Then c@ = c@ + 1
    Next
  EndIf
Return (c@)

_Sqrt
  Param (1)
  Local (3)

  Let b@ = 1
  Let c@ = 0

  Do Until b@ > a@
    Let b@ = b@ * 4
  Loop

  Do While b@ > 1
    Let b@ = b@ / 4
    Let d@ = a@ - c@ - b@
    Let c@ = c@ / 2
    If d@ > -1 Then
      Let a@ = d@
      Let c@ = c@ + b@
    Endif
  Loop

Return (c@)
Output:
The first 100 G numbers are:
   1   1   1   2   1   2   2   2   2   3
   3   3   2   3   2   4   4   2   3   4
   3   4   5   4   3   5   3   4   6   3
   5   6   2   5   6   5   5   7   4   5
   8   5   4   9   4   5   7   3   6   8
   5   6   8   6   7  10   6   6  12   4
   5  10   3   7   9   6   5   8   7   8
  11   6   5  12   4   8  11   5   8  10
   5   6  13   9   6  11   7   7  14   6
   8  13   5   8  11   7   9  13   8   9

The value of G(100000) is 810

0 OK, 0:210 

WrenEdit

Library: DOME
Library: Wren-math
Library: Wren-trait
Library: Wren-fmt
Library: Wren-plot

This follows the Raku example in plotting the first two thousand G values rather than the values up to G(2000) in order to produce an image something like the image in the Wikipedia article.

import "dome" for Window
import "graphics" for Canvas, Color
import "./math2" for Int
import "./trait" for Stepped
import "./fmt" for Fmt
import "./plot" for Axes

var limit = 4002
var primes = Int.primeSieve(limit-1).skip(1).toList
var goldbach = {4: 1}
for (i in Stepped.new(6..limit, 2)) goldbach[i] = 0
for (i in 0...primes.count) {
    for (j in i...primes.count) {
        var s = primes[i] + primes[j]
        if (s > limit) break
        goldbach[s] = goldbach[s] + 1
    }
}

System.print("The first 100 G values:")
var count = 0
for (i in Stepped.new(4..202, 2)) {
    count = count + 1
    Fmt.write("$2d ", goldbach[i])
    if (count % 10 == 0) System.print()
}

primes = Int.primeSieve(499999).skip(1)
var gm = 0
for (p in primes) {
    if (Int.isPrime(1e6 - p)) gm = gm + 1
}
System.print("\nG(1000000) = %(gm)")

var Red   = []
var Blue  = []
var Green = []

// create lists for the first 2000 G values for plotting by DOME.
for(e in Stepped.new(4..limit, 2)) {
    var c = e % 6
    var n = e/2 - 1
    if (c == 0) {
        Red.add([n, goldbach[e]])
    } else if (c == 2) {
        Blue.add([n, goldbach[e]])
    } else {
        Green.add([n, goldbach[e]])
    }
}

class Main {
    construct new() {
        Window.title = "Goldbach's comet"
        Canvas.resize(1000, 600)
        Window.resize(1000, 600)
        Canvas.cls(Color.white)
        var axes = Axes.new(100, 500, 800, 400, 0..2000, 0..200)
        axes.draw(Color.black, 2)
        var xMarks = Stepped.new(0..2000, 200)
        var yMarks = Stepped.new(0..200, 20)
        axes.mark(xMarks, yMarks, Color.black, 2)
        var xMarks2 = Stepped.new(0..2000, 400)
        var yMarks2 = Stepped.new(0..200, 40)
        axes.label(xMarks2, yMarks2, Color.black, 2, Color.black)
        axes.plot(Red, Color.red, "+")
        axes.plot(Blue, Color.blue, "+")
        axes.plot(Green, Color.green, "+")
    }

    init() {}

    update() {}

    draw(alpha) {}
}

var Game = Main.new()
Output:
The first 100 G values:
 1  1  1  2  1  2  2  2  2  3 
 3  3  2  3  2  4  4  2  3  4 
 3  4  5  4  3  5  3  4  6  3 
 5  6  2  5  6  5  5  7  4  5 
 8  5  4  9  4  5  7  3  6  8 
 5  6  8  6  7 10  6  6 12  4 
 5 10  3  7  9  6  5  8  7  8 
11  6  5 12  4  8 11  5  8 10 
 5  6 13  9  6 11  7  7 14  6 
 8 13  5  8 11  7  9 13  8  9 

G(1000000) = 5402

XPL0Edit

func IsPrime(N); \Return 'true' if N is prime
int  N, I;
[if N <= 2 then return N = 2;
if (N&1) = 0 then \even >2\ return false;
for I:= 3 to sqrt(N) do
    [if rem(N/I) = 0 then return false;
    I:= I+1;
    ];
return true;
];

int  PT(1_000_000);

func G(E);      \Ways E can be expressed as sum of two primes
int  E, C, I, J, T;
[C:= 0;  I:= 0;
loop    [J:= I;
        if PT(J) + PT(I) > E then return C;
        loop    [T:= PT(J) + PT(I);
                if T = E then C:= C+1;
                if T > E then quit;
                J:= J+1;
                ];
        I:= I+1;
        ];
];

int I, N;
[I:= 0; \make prime table
for N:= 2 to 1_000_000 do
        if IsPrime(N) then
                [PT(I):= N;  I:= I+1];
I:= 4;  \show first 100 G numbers
Format(4, 0);
for N:= 1 to 100 do
        [RlOut(0, float(G(I)));
        if rem(N/10) = 0 then CrLf(0);
        I:= I+2;
        ];
CrLf(0);
Text(0, "G(1,000,000) = ");  IntOut(0, G(1_000_000));
CrLf(0);
]
Output:
   1   1   1   2   1   2   2   2   2   3
   3   3   2   3   2   4   4   2   3   4
   3   4   5   4   3   5   3   4   6   3
   5   6   2   5   6   5   5   7   4   5
   8   5   4   9   4   5   7   3   6   8
   5   6   8   6   7  10   6   6  12   4
   5  10   3   7   9   6   5   8   7   8
  11   6   5  12   4   8  11   5   8  10
   5   6  13   9   6  11   7   7  14   6
   8  13   5   8  11   7   9  13   8   9

G(1,000,000) = 5402