

 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)

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

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
 Wikipedia: Goldbach's conjecture
 OEIS: A045917 - From Goldbach problem: number of decompositions of 2n into unordered sums of two primes




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