Goldbach's comet: Difference between revisions
added RPL |
|||
(33 intermediate revisions by 16 users not shown) | |||
Line 28: | Line 28: | ||
=={{header|11l}}== |
|||
{{trans|Python}} |
|||
<syntaxhighlight lang="11l"> |
|||
F is_prime(a) |
|||
I a == 2 |
|||
R 1B |
|||
I a < 2 | a % 2 == 0 |
|||
R 0B |
|||
L(i) (3 .. Int(sqrt(a))).step(2) |
|||
I a % i == 0 |
|||
R 0B |
|||
R 1B |
|||
F g(n) |
|||
assert(n > 2 & n % 2 == 0, ‘n in goldbach function g(n) must be even’) |
|||
V count = 0 |
|||
L(i) 1 .. n I/ 2 |
|||
I is_prime(i) & is_prime(n - i) |
|||
count++ |
|||
R count |
|||
print(‘The first 100 G numbers are:’) |
|||
V col = 1 |
|||
L(n) (4.<204).step(2) |
|||
print(String(g(n)).ljust(4), end' I (col % 10 == 0) {"\n"} E ‘’) |
|||
col++ |
|||
print("\nThe value of G(1000000) is "g(1'000'000)) |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|ALGOL 68}}== |
=={{header|ALGOL 68}}== |
||
Generates an ASCII-Art scatter plot - the vertical axis is n/10 and the hotizontal is G(n). |
Generates an ASCII-Art scatter plot - the vertical axis is n/10 and the hotizontal is G(n). |
||
{{libheader|ALGOL 68-primes}} |
{{libheader|ALGOL 68-primes}} |
||
< |
<syntaxhighlight lang="algol68">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 # |
# of prime pairs that sum to n, n even and > 2 # |
||
# generates an ASCII scatter plot of G(n) up to G(2000) # |
# generates an ASCII scatter plot of G(n) up to G(2000) # |
||
Line 94: | Line 143: | ||
print( ( newline ) ) |
print( ( newline ) ) |
||
OD |
OD |
||
END</ |
END</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 313: | Line 362: | ||
=={{header|Arturo}}== |
=={{header|Arturo}}== |
||
< |
<syntaxhighlight lang="rebol">G: function [n][ |
||
size select 2..n/2 'x -> |
size select 2..n/2 'x -> |
||
and? [prime? x][prime? n-x] |
and? [prime? x][prime? n-x] |
||
Line 330: | Line 379: | ||
; write the CSV data to a file which we can then visualize |
; write the CSV data to a file which we can then visualize |
||
; via our preferred spreadsheet app |
; via our preferred spreadsheet app |
||
write "comet.csv" csv</ |
write "comet.csv" csv</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 349: | Line 398: | ||
Here, you can find the result of the visualization - or the "Goldbach's comet": [https://i.ibb.co/JvLDRc0/Screenshot-2022-05-07-at-11-35-23.png 2D-chart of all G values up to 2000] |
Here, you can find the result of the visualization - or the "Goldbach's comet": [https://i.ibb.co/JvLDRc0/Screenshot-2022-05-07-at-11-35-23.png 2D-chart of all G values up to 2000] |
||
=={{header|AutoHotkey}}== |
|||
<syntaxhighlight lang="autohotkey">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 |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>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</pre> |
|||
=={{header|AWK}}== |
=={{header|AWK}}== |
||
<syntaxhighlight lang="awk"> |
|||
<lang AWK> |
|||
# syntax: GAWK -f GOLDBACHS_COMET.AWK |
# syntax: GAWK -f GOLDBACHS_COMET.AWK |
||
BEGIN { |
BEGIN { |
||
Line 388: | Line 473: | ||
return(1) |
return(1) |
||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 408: | Line 493: | ||
</pre> |
</pre> |
||
=={{ |
=={{Header|BASIC}}== |
||
==={{header|BASIC256}}=== |
|||
<lang freebasic>Function isPrime(Byval ValorEval As Uinteger) As Boolean |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="vb">#arraybase 1 |
|||
print "The first 100 G numbers are:" |
|||
col = 1 |
|||
for n = 4 to 202 step 2 |
|||
print rjust(string(g(n)), 4); |
|||
if col mod 10 = 0 then print |
|||
col += 1 |
|||
next n |
|||
print : print "G(1000000) = "; g(1000000) |
|||
end |
|||
function isPrime(v) |
|||
if v <= 1 then return False |
|||
for i = 2 to int(sqrt(v)) |
|||
if v mod i = 0 then return False |
|||
next i |
|||
return True |
|||
end function |
|||
function g(n) |
|||
cont = 0 |
|||
if n mod 2 = 0 then |
|||
for i = 2 to (1/2) * n |
|||
if isPrime(i) = 1 and isPrime(n - i) = 1 then cont += 1 |
|||
next i |
|||
end if |
|||
return cont |
|||
end function</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Similar to FreeBASIC entry.</pre> |
|||
==={{header|FreeBASIC}}=== |
|||
<syntaxhighlight lang="freebasic">Function isPrime(Byval ValorEval As Uinteger) As Boolean |
|||
If ValorEval <= 1 Then Return False |
If ValorEval <= 1 Then Return False |
||
For i As Integer = 2 To Int(Sqr(ValorEval)) |
For i As Integer = 2 To Int(Sqr(ValorEval)) |
||
Line 437: | Line 558: | ||
Print !"\nThe value of G(1000000) is "; g(1000000) |
Print !"\nThe value of G(1000000) is "; g(1000000) |
||
Sleep</ |
Sleep</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>The first 100 G numbers are: |
<pre>The first 100 G numbers are: |
||
Line 453: | Line 574: | ||
The value of G(1000000) is 5402</pre> |
The value of G(1000000) is 5402</pre> |
||
==={{header|Gambas}}=== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="vbnet">Use "isprime.bas" |
|||
Public Sub Main() |
|||
Print "The first 100 G numbers are:" |
|||
Dim n As Integer, col As Integer = 1 |
|||
For n = 4 To 202 Step 2 |
|||
Print Format$(Str(g(n)), "####"); |
|||
If col Mod 10 = 0 Then Print |
|||
col += 1 |
|||
Next |
|||
Print "\nG(1.000.000) = "; g(1000000) |
|||
End |
|||
Function g(n As Integer) As Integer |
|||
Dim i As Integer, count As Integer = 0 |
|||
If n Mod 2 = 0 Then |
|||
For i = 2 To n \ 2 '(1/2) * n |
|||
If isPrime(i) And isPrime(n - i) Then count += 1 |
|||
Next |
|||
End If |
|||
Return count |
|||
End Function</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Same as FreeBASIC entry.</pre> |
|||
==={{Header|SmileBASIC}}=== |
|||
(click the image to view full-size) |
|||
[[File:GoldbachsComet SmileBASIC 3DS screenshot.png|thumb]] |
|||
<syntaxhighlight lang="smilebasic">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</syntaxhighlight> |
|||
==={{header|uBasic/4tH}}=== |
|||
{{trans|FreeBASIC}} |
|||
For performance reasons only '''g(100000)''' is calculated. The value "810" has been verified and is correct. |
|||
<syntaxhighlight lang="text">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@)</syntaxhighlight> |
|||
{{out}} |
|||
<pre>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 </pre> |
|||
==={{header|Yabasic}}=== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="vb">import isprime |
|||
print "The first 100 G numbers are:" |
|||
col = 1 |
|||
for n = 4 to 202 step 2 |
|||
print g(n) using ("####"); |
|||
if mod(col, 10) = 0 print |
|||
col = col + 1 |
|||
next n |
|||
print "\nG(1000000) = ", g(1000000) |
|||
end |
|||
sub g(n) |
|||
count = 0 |
|||
if mod(n, 2) = 0 then |
|||
for i = 2 to (1/2) * n |
|||
if isPrime(i) and isPrime(n - i) count = count + 1 |
|||
next i |
|||
fi |
|||
return count |
|||
end sub</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Same as FreeBASIC entry.</pre> |
|||
=={{Header|C++}}== |
|||
<syntaxhighlight lang="c++"> |
|||
#include <cmath> |
|||
#include <cstdint> |
|||
#include <iomanip> |
|||
#include <iostream> |
|||
#include <stdexcept> |
|||
#include <vector> |
|||
std::vector<bool> primes; |
|||
void initialise_primes(const int32_t& limit) { |
|||
primes.resize(limit); |
|||
for ( int32_t i = 2; i < limit; ++i ) { |
|||
primes[i] = true; |
|||
} |
|||
for ( int32_t n = 2; n < sqrt(limit); ++n ) { |
|||
for ( int32_t k = n * n; k < limit; k += n ) { |
|||
primes[k] = false; |
|||
} |
|||
} |
|||
} |
|||
int32_t goldbach_function(const int32_t& number) { |
|||
if ( number <= 2 || number % 2 == 1 ) { |
|||
throw std::invalid_argument("Argument must be even and greater than 2: " + std::to_string(number)); |
|||
} |
|||
int32_t result = 0; |
|||
for ( int32_t i = 1; i <= number / 2; ++i ) { |
|||
if ( primes[i] && primes[number - i] ) { |
|||
result++; |
|||
} |
|||
} |
|||
return result; |
|||
} |
|||
int main() { |
|||
initialise_primes(2'000'000); |
|||
std::cout << "The first 100 Goldbach numbers:" << std::endl; |
|||
for ( int32_t n = 2; n < 102; ++n ) { |
|||
std::cout << std::setw(3) << goldbach_function(2 * n) << ( n % 10 == 1 ? "\n" : "" ); |
|||
} |
|||
std::cout << "\n" << "The 1,000,000th Goldbach number = " << goldbach_function(1'000'000) << std::endl; |
|||
} |
|||
</syntaxhighlight> |
|||
{{ out }} |
|||
<pre> |
|||
The first 100 Goldbach 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 |
|||
The 1,000,000th Goldbach number = 5402 |
|||
</pre> |
|||
=={{header|Delphi}}== |
|||
{{works with|Delphi|6.0}} |
|||
{{libheader|SysUtils,StdCtrls}} |
|||
<syntaxhighlight lang="Delphi"> |
|||
function GetGoldbachCount(N: integer): integer; |
|||
{Count number of prime number combinations add up to N } |
|||
var I: integer; |
|||
begin |
|||
Result:=0; |
|||
{Look at all number pairs that add up to N} |
|||
{And see if they are prime} |
|||
for I:=1 to N div 2 do |
|||
if IsPrime(I) and IsPrime(N-I) then Inc(Result); |
|||
end; |
|||
procedure ShowGoldbachComet(Memo: TMemo); |
|||
{Show first 100 Goldback numbers} |
|||
var I,N,Cnt,C: integer; |
|||
var S: string; |
|||
begin |
|||
Cnt:=0; N:=2; S:=''; |
|||
while true do |
|||
begin |
|||
C:=GetGoldbachCount(N); |
|||
if C>0 then |
|||
begin |
|||
Inc(Cnt); |
|||
S:=S+Format('%3d',[C]); |
|||
if (Cnt mod 10)=0 then S:=S+CRLF; |
|||
if Cnt>=100 then break; |
|||
end; |
|||
Inc(N,2); |
|||
end; |
|||
Memo.Lines.Add(S); |
|||
end; |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
Elapsed Time: 1.512 ms. |
|||
</pre> |
|||
=={{header|EasyLang}}== |
|||
{{trans|AWK}} |
|||
<syntaxhighlight lang="easylang"> |
|||
func isprim n . |
|||
if n mod 2 = 0 and n > 2 |
|||
return 0 |
|||
. |
|||
i = 3 |
|||
sq = sqrt n |
|||
while i <= sq |
|||
if n mod i = 0 |
|||
return 0 |
|||
. |
|||
i += 2 |
|||
. |
|||
return 1 |
|||
. |
|||
func goldbach n . |
|||
for i = 2 to n div 2 |
|||
if isprim i = 1 |
|||
cnt += isprim (n - i) |
|||
. |
|||
. |
|||
return cnt |
|||
. |
|||
numfmt 0 3 |
|||
for n = 4 step 2 to 202 |
|||
write goldbach n |
|||
if n mod 20 = 2 |
|||
print "" |
|||
. |
|||
. |
|||
print goldbach 1000000 |
|||
</syntaxhighlight> |
|||
=={{header|J}}== |
=={{header|J}}== |
||
Line 458: | Line 931: | ||
Task implementation: |
Task implementation: |
||
< |
<syntaxhighlight lang="j"> 10 10$#/.~4,/:~ 0-.~,(<:/~ * +/~) p:1+i.p:inv 202 |
||
1 1 1 2 1 2 2 2 2 3 |
1 1 1 2 1 2 2 2 2 3 |
||
3 3 2 3 2 4 4 2 3 4 |
3 3 2 3 2 4 4 2 3 4 |
||
Line 468: | Line 941: | ||
11 6 5 12 4 8 11 5 8 10 |
11 6 5 12 4 8 11 5 8 10 |
||
5 6 13 9 6 11 7 7 14 6 |
5 6 13 9 6 11 7 7 14 6 |
||
8 13 5 8 11 7 9 13 8 9</ |
8 13 5 8 11 7 9 13 8 9</syntaxhighlight> |
||
And, for G(1e6): |
And, for G(1e6): |
||
<syntaxhighlight lang="j"> |
|||
<lang j> |
|||
-:+/1 p: 1e6-p:i.p:inv 1e6 |
-:+/1 p: 1e6-p:i.p:inv 1e6 |
||
5402</ |
5402</syntaxhighlight> |
||
Explanation: |
Explanation: |
||
Line 490: | Line 963: | ||
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). |
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). |
||
=={{header|Java}}== |
|||
<syntaxhighlight lang="java"> |
|||
import java.awt.Color; |
|||
import java.awt.Graphics; |
|||
import java.awt.image.BufferedImage; |
|||
import java.io.File; |
|||
import java.io.IOException; |
|||
import java.util.List; |
|||
import javax.imageio.ImageIO; |
|||
public final class GoldbachsComet { |
|||
public static void main(String[] aArgs) { |
|||
initialisePrimes(2_000_000); |
|||
System.out.println("The first 100 Goldbach numbers:"); |
|||
for ( int n = 2; n < 102; n++ ) { |
|||
System.out.print(String.format("%3d%s", goldbachFunction(2 * n), ( n % 10 == 1 ? "\n" : "" ))); |
|||
} |
|||
System.out.println(); |
|||
System.out.println("The 1,000,000th Goldbach number = " + goldbachFunction(1_000_000)); |
|||
createImage(); |
|||
} |
|||
private static void createImage() { |
|||
final int width = 1040; |
|||
final int height = 860; |
|||
BufferedImage image = new BufferedImage(width, height, BufferedImage.TYPE_INT_RGB); |
|||
Graphics graphics = image.getGraphics(); |
|||
graphics.setColor(Color.WHITE); |
|||
graphics.fillRect(0, 0, width, height); |
|||
List<Color> colours = List.of( Color.BLUE, Color.GREEN, Color.RED ); |
|||
for ( int n = 2; n < 2002; n++ ) { |
|||
graphics.setColor(colours.get(n % 3)); |
|||
graphics.fillOval(n / 2, height - 5 * goldbachFunction(2 * n), 10, 10); |
|||
} |
|||
try { |
|||
ImageIO.write(image, "png", new File("GoldbachsCometJava.png")); |
|||
} catch (IOException ioe) { |
|||
ioe.printStackTrace(); |
|||
} |
|||
} |
|||
private static int goldbachFunction(int aNumber) { |
|||
if ( aNumber <= 2 || aNumber % 2 == 1 ) { |
|||
throw new AssertionError("Argument must be even and greater than 2: " + aNumber); |
|||
} |
|||
int result = 0; |
|||
for ( int i = 1; i <= aNumber / 2; i++ ) { |
|||
if ( primes[i] && primes[aNumber - i] ) { |
|||
result += 1; |
|||
} |
|||
} |
|||
return result; |
|||
} |
|||
private static void initialisePrimes(int aLimit) { |
|||
primes = new boolean[aLimit]; |
|||
for ( int i = 2; i < aLimit; i++ ) { |
|||
primes[i] = true; |
|||
} |
|||
for ( int n = 2; n < Math.sqrt(aLimit); n++ ) { |
|||
for ( int k = n * n; k < aLimit; k += n ) { |
|||
primes[k] = false; |
|||
} |
|||
} |
|||
} |
|||
private static boolean[] primes; |
|||
} |
|||
</syntaxhighlight> |
|||
{{ out }} |
|||
[[Media:GoldbachsCometJava.png]] |
|||
<pre> |
|||
The first 100 Goldbach 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 |
|||
The 1,000,000th Goldbach number = 5402 |
|||
</pre> |
|||
=={{header|jq}}== |
|||
'''Works with jq and gojq, the C and Go implementations of jq.''' |
|||
'''Preliminaries''' |
|||
<syntaxhighlight lang=jq> |
|||
def count(s): reduce s as $_ (0; .+1); |
|||
def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .; |
|||
def nwise($n): |
|||
def n: if length <= $n then . else .[0:$n] , (.[$n:] | n) end; |
|||
n; |
|||
def is_prime: |
|||
. as $n |
|||
| if ($n < 2) then false |
|||
elif ($n % 2 == 0) then $n == 2 |
|||
elif ($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 |
|||
($n | sqrt) as $rt |
|||
| 23 |
|||
| until( . > $rt or ($n % . == 0); .+2) |
|||
| . > $rt |
|||
end; |
|||
</syntaxhighlight> |
|||
'''The Tasks''' |
|||
<syntaxhighlight lang=jq> |
|||
# emit nothing if . is odd |
|||
def G: |
|||
select(. % 2 == 0) |
|||
| count( range(2; (./2)+1) as $i |
|||
| select(($i|is_prime) and ((.-$i)|is_prime)) ); |
|||
def task1: |
|||
"The first 100 G numbers:", |
|||
([range(4; 203; 2) | G] | nwise(10) | map(lpad(4)) | join(" ")); |
|||
def task($n): |
|||
$n, 4, 22 |
|||
|"G(\(.)): \(G)"; |
|||
task1, "", task(1000000) |
|||
</syntaxhighlight> |
|||
{{output}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
Run in VS Code or REPL to view and save the plot. |
Run in VS Code or REPL to view and save the plot. |
||
< |
<syntaxhighlight lang="julia">using Combinatorics |
||
using Plots |
using Plots |
||
using Primes |
using Primes |
||
Line 508: | Line 1,146: | ||
y = map(g, 2x) |
y = map(g, 2x) |
||
scatter(x, y, markerstrokewidth = 0, color = ["red", "blue", "green"][mod1.(x, 3)]) |
scatter(x, y, markerstrokewidth = 0, color = ["red", "blue", "green"][mod1.(x, 3)]) |
||
</ |
</syntaxhighlight>{{out}}<pre> |
||
The first 100 G numbers are: |
The first 100 G numbers are: |
||
1 1 1 2 1 2 2 2 2 3 |
1 1 1 2 1 2 2 2 2 3 |
||
Line 523: | Line 1,161: | ||
The value of G(1000000) is 5402 |
The value of G(1000000) is 5402 |
||
</pre> |
</pre> |
||
[[File:Jgoldbachplot.png]] |
|||
=={{header|Lua}}== |
|||
<syntaxhighlight lang="lua">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))</syntaxhighlight> |
|||
{{out}} |
|||
<pre>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</pre> |
|||
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
||
< |
<syntaxhighlight lang="mathematica">ClearAll[GoldbachFuncion] |
||
GoldbachFuncion[e_Integer] := Module[{ps}, |
GoldbachFuncion[e_Integer] := Module[{ps}, |
||
ps = Prime[Range[PrimePi[e/2]]]; |
ps = Prime[Range[PrimePi[e/2]]]; |
||
Line 532: | Line 1,216: | ||
Grid[Partition[GoldbachFuncion /@ Range[4, 220, 2], 10]] |
Grid[Partition[GoldbachFuncion /@ Range[4, 220, 2], 10]] |
||
GoldbachFuncion[10^6] |
GoldbachFuncion[10^6] |
||
DiscretePlot[GoldbachFuncion[e], {e, 4, 2000}, Filling -> None]</ |
DiscretePlot[GoldbachFuncion[e], {e, 4, 2000}, Filling -> None]</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>1 1 1 2 1 2 2 2 2 3 |
<pre>1 1 1 2 1 2 2 2 2 3 |
||
Line 549: | Line 1,233: | ||
[graphical object showing Goldbach's comet]</pre> |
[graphical object showing Goldbach's comet]</pre> |
||
=={{header|Nim}}== |
|||
{{libheader|nim-plotly}} |
|||
{{libheader|chroma}} |
|||
To display the Golbach’s comet, we use a library providing a Nim interface to “plotly”. The graph is displayed into a browser. |
|||
<syntaxhighlight lang="Nim">import std/[math, strformat, strutils, sugar] |
|||
import chroma, plotly |
|||
const |
|||
N1 = 100 # For part 1 of task. |
|||
N2 = 1_000_000 # For part 2 of task. |
|||
N3 = 2000 # For stretch part. |
|||
# Erathostenes sieve. |
|||
var isPrime: array[1..N1, bool] |
|||
for i in 2..N1: isPrime[i] = true |
|||
for n in 2..sqrt(N1.toFloat).int: |
|||
for k in countup(n * n, N1, n): |
|||
isPrime[k] = false |
|||
proc g(n: int): int = |
|||
## Goldbach function. |
|||
assert n > 2 and n mod 2 == 0, "“n” must be even and greater than 2." |
|||
for i in 1..(n div 2): |
|||
if isPrime[i] and isPrime[n - i]: |
|||
inc result |
|||
# Part 1. |
|||
echo &"First {N1} G numbers:" |
|||
var col = 1 |
|||
for n in 2..N1: |
|||
stdout.write align($g( 2 * n), 3) |
|||
stdout.write if col mod 10 == 0: '\n' else: ' ' |
|||
inc col |
|||
# Part 2. |
|||
echo &"\nG({N2}) = ", g(N2) |
|||
# Stretch part. |
|||
const Colors = collect(for name in ["red", "blue", "green"]: name.parseHtmlName()) |
|||
var x, y: seq[float] |
|||
var colors: seq[Color] |
|||
for n in 2..N3: |
|||
x.add n.toFloat |
|||
y.add g(2 * n).toFloat |
|||
colors.add Colors[n mod 3] |
|||
let trace = Trace[float](type: Scatter, mode: Markers, marker: Marker[float](color: colors), xs: x, ys: y) |
|||
let layout = Layout(title: "Goldbach’s comet", width: 1200, height: 400) |
|||
Plot[float64](layout: layout, traces: @[trace]).show(removeTempFile = true) |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre>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 |
|||
</pre> |
|||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
{{libheader|ntheory}} |
{{libheader|ntheory}} |
||
< |
<syntaxhighlight lang="perl">use strict; |
||
use warnings; |
use warnings; |
||
use feature 'say'; |
use feature 'say'; |
||
Line 587: | Line 1,337: | ||
binmode $fh; |
binmode $fh; |
||
print $fh $gd->png(); |
print $fh $gd->png(); |
||
close $fh;</ |
close $fh;</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> 1 1 1 2 1 2 2 2 2 3 |
<pre> 1 1 1 2 1 2 2 2 2 3 |
||
Line 608: | Line 1,358: | ||
{{libheader|Phix/online}} |
{{libheader|Phix/online}} |
||
You can run this online [http://phix.x10.mx/p2js/goldbach.htm here]. |
You can run this online [http://phix.x10.mx/p2js/goldbach.htm here]. |
||
<!--< |
<!--<syntaxhighlight lang="phix">(phixonline)--> |
||
<span style="color: #000080;font-style:italic;">-- |
<span style="color: #000080;font-style:italic;">-- |
||
-- demo\rosetta\Goldbachs_comet.exw |
-- demo\rosetta\Goldbachs_comet.exw |
||
Line 656: | Line 1,406: | ||
<span style="color: #7060A8;">IupClose</span><span style="color: #0000FF;">()</span> |
<span style="color: #7060A8;">IupClose</span><span style="color: #0000FF;">()</span> |
||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
||
<!--</ |
<!--</syntaxhighlight>--> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 676: | Line 1,426: | ||
=={{header|Picat}}== |
=={{header|Picat}}== |
||
< |
<syntaxhighlight lang="picat">main => |
||
println("First 100 G numbers:"), |
println("First 100 G numbers:"), |
||
foreach({G,I} in zip(take([G: T in 1..300, G=g(T),G>0],100),1..100)) |
foreach({G,I} in zip(take([G: T in 1..300, G=g(T),G>0],100),1..100)) |
||
Line 687: | Line 1,437: | ||
{1 : I in 1..N // 2, |
{1 : I in 1..N // 2, |
||
prime(I),prime(N-I)}.len, |
prime(I),prime(N-I)}.len, |
||
0).</ |
0).</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 705: | Line 1,455: | ||
=={{header|Python}}== |
=={{header|Python}}== |
||
< |
<syntaxhighlight lang="python">from matplotlib.pyplot import scatter, show |
||
from sympy import isprime |
from sympy import isprime |
||
Line 729: | Line 1,479: | ||
colors = [["red", "blue", "green"][(i // 2) % 3] 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) |
scatter([i // 2 for i in x], y, marker='.', color = colors) |
||
show()</ |
show()</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
The first 100 G numbers are: |
The first 100 G numbers are: |
||
Line 745: | Line 1,495: | ||
The value of G(1000000) is 5402 |
The value of G(1000000) is 5402 |
||
</pre> |
</pre> |
||
[[File:PGoldbachplot.png]] |
|||
=={{header|Raku}}== |
=={{header|Raku}}== |
||
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 [[wp:Goldbach's_comet|wikipedia example image]].) |
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 [[wp:Goldbach's_comet|wikipedia example image]].) |
||
<lang |
<syntaxhighlight lang="raku" line>sub G (Int $n) { +(2..$n/2).grep: { .is-prime && ($n - $_).is-prime } } |
||
# Task |
# Task |
||
Line 771: | Line 1,522: | ||
values => [@y,], |
values => [@y,], |
||
).plot: :points; |
).plot: :points; |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre>The first 100 G values: |
<pre>The first 100 G values: |
||
Line 788: | Line 1,539: | ||
'''Stretch goal:''' (offsite SVG image) [https://raw.githubusercontent.com/thundergnat/rc/master/img/Goldbachs-Comet-Raku.svg Goldbachs-Comet-Raku.svg] |
'''Stretch goal:''' (offsite SVG image) [https://raw.githubusercontent.com/thundergnat/rc/master/img/Goldbachs-Comet-Raku.svg Goldbachs-Comet-Raku.svg] |
||
=={{header|RPL}}== |
|||
{{works with|HP|49}} |
|||
« '''IF''' 2 MOD '''THEN''' |
|||
"GOLDB Error: Odd number" DOERR |
|||
'''ELSE''' |
|||
0 |
|||
2 PICK3 2 / CEIL '''FOR''' j |
|||
'''IF''' j ISPRIME? '''THEN''' |
|||
'''IF''' OVER j - ISPRIME? '''THEN''' 1 + '''END''' |
|||
'''END''' |
|||
'''NEXT''' NIP |
|||
'''END''' |
|||
» '<span style="color:blue">GOLDB</span>' STO |
|||
« « n <span style="color:blue">GOLDB</span> » 'n' 4 202 2 SEQ |
|||
OBJ→ DROP { 10 10 } →ARRY |
|||
1000000 <span style="color:blue">GOLDB</span> "G(1000000)" →TAG |
|||
» '<span style="color:blue">TASK</span>' STO |
|||
{{out}} |
|||
<pre> |
|||
2: [[ 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 ]] |
|||
1: G(1000000): 5402 |
|||
</pre> |
|||
=={{header|Ruby}}== |
|||
following the J comments, but only generating primes up-to half a million |
|||
<syntaxhighlight lang="ruby">require 'prime' |
|||
n = 100 |
|||
puts "The first #{n} Godbach numbers are: " |
|||
sums = Prime.each(n*2 + 2).to_a[1..].repeated_combination(2).map(&:sum) |
|||
sums << 4 |
|||
sums.sort.tally.values[...n].each_slice(10){|slice| puts "%4d"*slice.size % slice} |
|||
n = 1000000 |
|||
puts "\nThe value of G(#{n}) is #{Prime.each(n/2).count{|pr| (n-pr).prime?} }." |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre>The first 100 Godbach 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. |
|||
</pre> |
|||
=={{header|Rust}}== |
|||
<syntaxhighlight lang="rust">// [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), |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
[[Media:Goldbach's comet rust.svg]] |
|||
=={{header|Wren}}== |
=={{header|Wren}}== |
||
{{libheader|DOME}} |
|||
{{libheader|Wren-math}} |
{{libheader|Wren-math}} |
||
{{libheader|Wren- |
{{libheader|Wren-iterate}} |
||
{{libheader|Wren-fmt}} |
{{libheader|Wren-fmt}} |
||
{{libheader|Wren-plot}} |
|||
<lang ecmascript>import "./math" for Int |
|||
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 "./trait" for Stepped |
|||
<syntaxhighlight lang="wren">import "dome" for Window |
|||
import "graphics" for Canvas, Color |
|||
import "./math2" for Int |
|||
import "./iterate" for Stepped |
|||
import "./fmt" for Fmt |
import "./fmt" for Fmt |
||
import " |
import "./plot" for Axes |
||
var limit = |
var limit = 4002 |
||
var primes = Int.primeSieve(limit-1).skip(1).toList |
var primes = Int.primeSieve(limit-1).skip(1).toList |
||
var goldbach = {4: 1} |
var goldbach = {4: 1} |
||
Line 825: | Line 1,739: | ||
System.print("\nG(1000000) = %(gm)") |
System.print("\nG(1000000) = %(gm)") |
||
var Red = [] |
|||
// create .csv file for values up to 2000 for display by an external plotter |
|||
var Blue = [] |
|||
// the third field being the color (red = 0, blue = 1, green = 2) |
|||
var Green = [] |
|||
File.create("goldbachs_comet.csv") { |file| |
|||
for(i in Stepped.new(4..limit, 2)) { |
|||
// create lists for the first 2000 G values for plotting by DOME. |
|||
file.writeBytes("%(i), %(goldbach[i]), %(i/2 % 3)\n") |
|||
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]]) |
|||
} |
} |
||
} |
|||
}</lang> |
|||
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()</syntaxhighlight> |
|||
{{out}} |
{{out}} |
||
Line 851: | Line 1,802: | ||
=={{header|XPL0}}== |
=={{header|XPL0}}== |
||
<syntaxhighlight lang="xpl0"> |
|||
<lang XPL0> |
|||
func IsPrime(N); \Return 'true' if N is prime |
func IsPrime(N); \Return 'true' if N is prime |
||
int N, I; |
int N, I; |
||
Line 894: | Line 1,845: | ||
Text(0, "G(1,000,000) = "); IntOut(0, G(1_000_000)); |
Text(0, "G(1,000,000) = "); IntOut(0, G(1_000_000)); |
||
CrLf(0); |
CrLf(0); |
||
]</ |
]</syntaxhighlight> |
||
{{out}} |
{{out}} |
Latest revision as of 14:08, 2 January 2024
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
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) |
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
11l
F is_prime(a)
I a == 2
R 1B
I a < 2 | a % 2 == 0
R 0B
L(i) (3 .. Int(sqrt(a))).step(2)
I a % i == 0
R 0B
R 1B
F g(n)
assert(n > 2 & n % 2 == 0, ‘n in goldbach function g(n) must be even’)
V count = 0
L(i) 1 .. n I/ 2
I is_prime(i) & is_prime(n - i)
count++
R count
print(‘The first 100 G numbers are:’)
V col = 1
L(n) (4.<204).step(2)
print(String(g(n)).ljust(4), end' I (col % 10 == 0) {"\n"} E ‘’)
col++
print("\nThe value of G(1000000) is "g(1'000'000))
- 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
ALGOL 68
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
|+ |.= | -+ | -.- | --. | --. | .. .- | +.. | -- . + ... . . | .- - | +. . | ... . . | ..... | .... . | .. . .. | .. . . . | .- . . | .. . .. + ... . . | -.. . | ... . . | - . .. | . + . | . .. . . | .- . . | .-. . | - . . . | . . . .. + - .. . | ... . . | .. . . . | . .- . | . . . . . | .. . .. | .- . . | - . . . | .. . . . | .. .. . + . - . . | . .. . . | .- . . | . .. . . | . .. . . | .-. . | . - . . | .. . . . | . ... . | . . . . . + - . .. | ..- . | - . . . | . . . .. | .- . . | .. . - | - . - | -. . . | .. . . . | .. . . . + . -. . | ... . . | . .. - | . ... . | . .. . . | ... . . | . . . . . | . . . . . | . . . . . | + . . + . . . . . | .. - . | .... . | .- . . | . . . . . | -- . | .. . . . | - . .. | .. . . . | .. . . . + . . . . . | . - . . | .. . . . | . .. . . | .. - . | . .. . . | .. . . . | -. . . | . .. . . | . . . . . + -. . . | - . . . | . . - . | - . . . | .. . . . | . . . .. | . .. . . | ... - | . . . . . | . . . . . + . . . . . | - . . . | . . . . . | - . . . | . . . . . | - . . . | .. . . . | - . . . | .. . . . | . .. . . + . . . . . | .- . . | . .. . . | . . . . . | . .. . . | . . . . . | . . . . . | . .. . . | .. . . . | - . . . + .- . . | .. . . . | .. . . . | - . . . | . . . . . | . . . .. | .. . . . | .. . . . | ... . . | . .. . . + . . . . . | - . . . | . .. . . | .. . . . | - . . . | . . - . | .. . . . | - . . . | . -. . | - . . . + - . .. | . . . . . | . . . . . | .. . . . | .. . . . | . .. . . | . . . .. | .. . . . | . . . . . | . . . - + .. .. . | . - . . | . . . . . | -. . . | . . . . . | . . - . | .. . . . | . .. . . | . . . . . | . . . . . + .. . . . | . . . . . | . . . . . | . . . . . | - . . . | ... . . | .. . . . | . . . .. | . . . . . | . . . . . + . . . . . | - . . . | . .. . . | . . . . . | . . .. . | - . . . | . . . . . | . .. . . | .. . . . | . . . . . + . . - . | . - . . | .. . - | . . . . . | . . . . . | . . . . . | .. . . . | . . . . . | . . . . . | . . . . . + .. . . . | . - . . | . . .. . | - . . . | . . . . . | .. . . . | . . . . . | . . . . . | -. . . | - . . . + .. . . .
Arturo
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
AutoHotkey
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
AWK
# 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
BASIC
BASIC256
#arraybase 1
print "The first 100 G numbers are:"
col = 1
for n = 4 to 202 step 2
print rjust(string(g(n)), 4);
if col mod 10 = 0 then print
col += 1
next n
print : print "G(1000000) = "; g(1000000)
end
function isPrime(v)
if v <= 1 then return False
for i = 2 to int(sqrt(v))
if v mod i = 0 then return False
next i
return True
end function
function g(n)
cont = 0
if n mod 2 = 0 then
for i = 2 to (1/2) * n
if isPrime(i) = 1 and isPrime(n - i) = 1 then cont += 1
next i
end if
return cont
end function
- Output:
Similar to FreeBASIC entry.
FreeBASIC
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
Gambas
Use "isprime.bas"
Public Sub Main()
Print "The first 100 G numbers are:"
Dim n As Integer, col As Integer = 1
For n = 4 To 202 Step 2
Print Format$(Str(g(n)), "####");
If col Mod 10 = 0 Then Print
col += 1
Next
Print "\nG(1.000.000) = "; g(1000000)
End
Function g(n As Integer) As Integer
Dim i As Integer, count As Integer = 0
If n Mod 2 = 0 Then
For i = 2 To n \ 2 '(1/2) * n
If isPrime(i) And isPrime(n - i) Then count += 1
Next
End If
Return count
End Function
- Output:
Same as FreeBASIC entry.
SmileBASIC
(click the image to view full-size)
![](http://static.miraheze.org/rosettacodewiki/thumb/2/27/GoldbachsComet_SmileBASIC_3DS_screenshot.png/300px-GoldbachsComet_SmileBASIC_3DS_screenshot.png)
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
uBasic/4tH
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
Yabasic
import isprime
print "The first 100 G numbers are:"
col = 1
for n = 4 to 202 step 2
print g(n) using ("####");
if mod(col, 10) = 0 print
col = col + 1
next n
print "\nG(1000000) = ", g(1000000)
end
sub g(n)
count = 0
if mod(n, 2) = 0 then
for i = 2 to (1/2) * n
if isPrime(i) and isPrime(n - i) count = count + 1
next i
fi
return count
end sub
- Output:
Same as FreeBASIC entry.
C++
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <vector>
std::vector<bool> primes;
void initialise_primes(const int32_t& limit) {
primes.resize(limit);
for ( int32_t i = 2; i < limit; ++i ) {
primes[i] = true;
}
for ( int32_t n = 2; n < sqrt(limit); ++n ) {
for ( int32_t k = n * n; k < limit; k += n ) {
primes[k] = false;
}
}
}
int32_t goldbach_function(const int32_t& number) {
if ( number <= 2 || number % 2 == 1 ) {
throw std::invalid_argument("Argument must be even and greater than 2: " + std::to_string(number));
}
int32_t result = 0;
for ( int32_t i = 1; i <= number / 2; ++i ) {
if ( primes[i] && primes[number - i] ) {
result++;
}
}
return result;
}
int main() {
initialise_primes(2'000'000);
std::cout << "The first 100 Goldbach numbers:" << std::endl;
for ( int32_t n = 2; n < 102; ++n ) {
std::cout << std::setw(3) << goldbach_function(2 * n) << ( n % 10 == 1 ? "\n" : "" );
}
std::cout << "\n" << "The 1,000,000th Goldbach number = " << goldbach_function(1'000'000) << std::endl;
}
- Output:
The first 100 Goldbach 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 The 1,000,000th Goldbach number = 5402
Delphi
function GetGoldbachCount(N: integer): integer;
{Count number of prime number combinations add up to N }
var I: integer;
begin
Result:=0;
{Look at all number pairs that add up to N}
{And see if they are prime}
for I:=1 to N div 2 do
if IsPrime(I) and IsPrime(N-I) then Inc(Result);
end;
procedure ShowGoldbachComet(Memo: TMemo);
{Show first 100 Goldback numbers}
var I,N,Cnt,C: integer;
var S: string;
begin
Cnt:=0; N:=2; S:='';
while true do
begin
C:=GetGoldbachCount(N);
if C>0 then
begin
Inc(Cnt);
S:=S+Format('%3d',[C]);
if (Cnt mod 10)=0 then S:=S+CRLF;
if Cnt>=100 then break;
end;
Inc(N,2);
end;
Memo.Lines.Add(S);
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 Elapsed Time: 1.512 ms.
EasyLang
func isprim n .
if n mod 2 = 0 and n > 2
return 0
.
i = 3
sq = sqrt n
while i <= sq
if n mod i = 0
return 0
.
i += 2
.
return 1
.
func goldbach n .
for i = 2 to n div 2
if isprim i = 1
cnt += isprim (n - i)
.
.
return cnt
.
numfmt 0 3
for n = 4 step 2 to 202
write goldbach n
if n mod 20 = 2
print ""
.
.
print goldbach 1000000
J
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).
Java
import java.awt.Color;
import java.awt.Graphics;
import java.awt.image.BufferedImage;
import java.io.File;
import java.io.IOException;
import java.util.List;
import javax.imageio.ImageIO;
public final class GoldbachsComet {
public static void main(String[] aArgs) {
initialisePrimes(2_000_000);
System.out.println("The first 100 Goldbach numbers:");
for ( int n = 2; n < 102; n++ ) {
System.out.print(String.format("%3d%s", goldbachFunction(2 * n), ( n % 10 == 1 ? "\n" : "" )));
}
System.out.println();
System.out.println("The 1,000,000th Goldbach number = " + goldbachFunction(1_000_000));
createImage();
}
private static void createImage() {
final int width = 1040;
final int height = 860;
BufferedImage image = new BufferedImage(width, height, BufferedImage.TYPE_INT_RGB);
Graphics graphics = image.getGraphics();
graphics.setColor(Color.WHITE);
graphics.fillRect(0, 0, width, height);
List<Color> colours = List.of( Color.BLUE, Color.GREEN, Color.RED );
for ( int n = 2; n < 2002; n++ ) {
graphics.setColor(colours.get(n % 3));
graphics.fillOval(n / 2, height - 5 * goldbachFunction(2 * n), 10, 10);
}
try {
ImageIO.write(image, "png", new File("GoldbachsCometJava.png"));
} catch (IOException ioe) {
ioe.printStackTrace();
}
}
private static int goldbachFunction(int aNumber) {
if ( aNumber <= 2 || aNumber % 2 == 1 ) {
throw new AssertionError("Argument must be even and greater than 2: " + aNumber);
}
int result = 0;
for ( int i = 1; i <= aNumber / 2; i++ ) {
if ( primes[i] && primes[aNumber - i] ) {
result += 1;
}
}
return result;
}
private static void initialisePrimes(int aLimit) {
primes = new boolean[aLimit];
for ( int i = 2; i < aLimit; i++ ) {
primes[i] = true;
}
for ( int n = 2; n < Math.sqrt(aLimit); n++ ) {
for ( int k = n * n; k < aLimit; k += n ) {
primes[k] = false;
}
}
}
private static boolean[] primes;
}
- Output:
The first 100 Goldbach 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 The 1,000,000th Goldbach number = 5402
jq
Works with jq and gojq, the C and Go implementations of jq.
Preliminaries
def count(s): reduce s as $_ (0; .+1);
def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;
def nwise($n):
def n: if length <= $n then . else .[0:$n] , (.[$n:] | n) end;
n;
def is_prime:
. as $n
| if ($n < 2) then false
elif ($n % 2 == 0) then $n == 2
elif ($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
($n | sqrt) as $rt
| 23
| until( . > $rt or ($n % . == 0); .+2)
| . > $rt
end;
The Tasks
# emit nothing if . is odd
def G:
select(. % 2 == 0)
| count( range(2; (./2)+1) as $i
| select(($i|is_prime) and ((.-$i)|is_prime)) );
def task1:
"The first 100 G numbers:",
([range(4; 203; 2) | G] | nwise(10) | map(lpad(4)) | join(" "));
def task($n):
$n, 4, 22
|"G(\(.)): \(G)";
task1, "", task(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 G(4): 1 G(22): 3
Julia
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 9The value of G(1000000) is 5402
Lua
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 Language
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]
Nim
To display the Golbach’s comet, we use a library providing a Nim interface to “plotly”. The graph is displayed into a browser.
import std/[math, strformat, strutils, sugar]
import chroma, plotly
const
N1 = 100 # For part 1 of task.
N2 = 1_000_000 # For part 2 of task.
N3 = 2000 # For stretch part.
# Erathostenes sieve.
var isPrime: array[1..N1, bool]
for i in 2..N1: isPrime[i] = true
for n in 2..sqrt(N1.toFloat).int:
for k in countup(n * n, N1, n):
isPrime[k] = false
proc g(n: int): int =
## Goldbach function.
assert n > 2 and n mod 2 == 0, "“n” must be even and greater than 2."
for i in 1..(n div 2):
if isPrime[i] and isPrime[n - i]:
inc result
# Part 1.
echo &"First {N1} G numbers:"
var col = 1
for n in 2..N1:
stdout.write align($g( 2 * n), 3)
stdout.write if col mod 10 == 0: '\n' else: ' '
inc col
# Part 2.
echo &"\nG({N2}) = ", g(N2)
# Stretch part.
const Colors = collect(for name in ["red", "blue", "green"]: name.parseHtmlName())
var x, y: seq[float]
var colors: seq[Color]
for n in 2..N3:
x.add n.toFloat
y.add g(2 * n).toFloat
colors.add Colors[n mod 3]
let trace = Trace[float](type: Scatter, mode: Markers, marker: Marker[float](color: colors), xs: x, ys: y)
let layout = Layout(title: "Goldbach’s comet", width: 1200, height: 400)
Plot[float64](layout: layout, traces: @[trace]).show(removeTempFile = true)
- 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
Perl
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
Phix
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
Picat
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
Python
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
Raku
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
RPL
« IF 2 MOD THEN "GOLDB Error: Odd number" DOERR ELSE 0 2 PICK3 2 / CEIL FOR j IF j ISPRIME? THEN IF OVER j - ISPRIME? THEN 1 + END END NEXT NIP END » 'GOLDB' STO « « n GOLDB » 'n' 4 202 2 SEQ OBJ→ DROP { 10 10 } →ARRY 1000000 GOLDB "G(1000000)" →TAG » 'TASK' STO
- Output:
2: [[ 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 ]] 1: G(1000000): 5402
Ruby
following the J comments, but only generating primes up-to half a million
require 'prime'
n = 100
puts "The first #{n} Godbach numbers are: "
sums = Prime.each(n*2 + 2).to_a[1..].repeated_combination(2).map(&:sum)
sums << 4
sums.sort.tally.values[...n].each_slice(10){|slice| puts "%4d"*slice.size % slice}
n = 1000000
puts "\nThe value of G(#{n}) is #{Prime.each(n/2).count{|pr| (n-pr).prime?} }."
- Output:
The first 100 Godbach 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.
Rust
// [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
Wren
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 "./iterate" 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
XPL0
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
- Programming Tasks
- Prime Numbers
- WikipediaSourced
- 11l
- ALGOL 68
- ALGOL 68-primes
- Arturo
- AutoHotkey
- AWK
- BASIC
- BASIC256
- FreeBASIC
- Gambas
- SmileBASIC
- UBasic/4tH
- Yabasic
- C++
- Delphi
- SysUtils,StdCtrls
- EasyLang
- J
- Java
- Jq
- Julia
- Lua
- Mathematica
- Wolfram Language
- Nim
- Nim-plotly
- Chroma
- Perl
- Ntheory
- Phix
- Phix/pGUI
- Phix/online
- Picat
- Python
- Raku
- RPL
- Ruby
- Rust
- Wren
- DOME
- Wren-math
- Wren-iterate
- Wren-fmt
- Wren-plot
- XPL0