# 2022.11.04 MATH 3600-001 Goldbach

## 315 days ago by calkin

The Goldbach Conjecture: Every even number greater than 2 is the sum of 2 primes in at least one way.

Let us create a list of odd primes.

N=100000 prime_list=[] p=3 while p<N: prime_list.append(p) p=next_prime(p)

goldbach_representations=[0 for i in srange(2*N)] for p in prime_list: for q in prime_list: goldbach_representations[p+q]+=1
print(goldbach_representations)
 Traceback (click to the left of this block for traceback) ... IndexError: list index out of range Traceback (most recent call last): File "", line 1, in File "_sage_input_6.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("cHJpbnQoZ29sZGJhY2hfcmVwcmVzZW50YXRpb25zWzEwMDAwMF0p"),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmplXD5N3/___code___.py", line 3, in exec compile(u'print(goldbach_representations[_sage_const_100000 ]) File "", line 1, in IndexError: list index out of range

How many representations do we get for even numbers?

print(goldbach_representations)
 Traceback (click to the left of this block for traceback) ... IndexError: list index out of range Traceback (most recent call last): File "", line 1, in File "_sage_input_7.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("cHJpbnQoZ29sZGJhY2hfcmVwcmVzZW50YXRpb25zWzc4OTY0XSk="),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmpptkO4P/___code___.py", line 3, in exec compile(u'print(goldbach_representations[_sage_const_78964 ]) File "", line 1, in IndexError: list index out of range
print(goldbach_representations)
 Traceback (click to the left of this block for traceback) ... IndexError: list index out of range Traceback (most recent call last): File "", line 1, in File "_sage_input_8.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("cHJpbnQoZ29sZGJhY2hfcmVwcmVzZW50YXRpb25zWzg1MDAwXSk="),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmpbfBBPE/___code___.py", line 3, in exec compile(u'print(goldbach_representations[_sage_const_85000 ]) File "", line 1, in IndexError: list index out of range
list_plot(goldbach_representations) list_plot(goldbach_representations[:N]) Idea: look at the champions: n for which r(n) is bigger than anything seen before.

goldbach_reps_pairs=[] for i in srange(len(goldbach_representations)/2): if goldbach_representations[i]>0: goldbach_reps_pairs.append([i,goldbach_representations[i]])
list_plot(goldbach_reps_pairs) Let's find the champions in this list.

max_seen=0 champions=[] for x in goldbach_reps_pairs: if x> max_seen: max_seen=x champions.append(x)
list_plot(champions) len(champions)
 52 52
print(champions)
 [6, 8, 10, 16, 22, 24, 34, 36, 48, 60, 78, 84, 90, 114, 120, 168, 180, 210, 300, 330, 390, 420, 510, 630, 780, 840, 990, 1050, 1140, 1260, 1470, 1650, 1680, 1890, 2100, 2310, 2730, 3150, 3570, 3990, 4200, 4410, 4620, 5250, 5460, 6090, 6510, 6930, 7980, 8190, 9030, 9240] [6, 8, 10, 16, 22, 24, 34, 36, 48, 60, 78, 84, 90, 114, 120, 168, 180, 210, 300, 330, 390, 420, 510, 630, 780, 840, 990, 1050, 1140, 1260, 1470, 1650, 1680, 1890, 2100, 2310, 2730, 3150, 3570, 3990, 4200, 4410, 4620, 5250, 5460, 6090, 6510, 6930, 7980, 8190, 9030, 9240]
factor(90)
 2 * 3^2 * 5 2 * 3^2 * 5
factor(114)
 2 * 3 * 19 2 * 3 * 19
factor(120)
 2^3 * 3 * 5 2^3 * 3 * 5
factor(168)
 2^3 * 3 * 7 2^3 * 3 * 7

It appears that champions have x coordinate divisible by 3, 5, 7, usually?

g30=[] g31=[] g32=[] for x in goldbach_reps_pairs: a=mod(x,3) if a==0: g30.append(x) if a==1: g31.append(x) if a==2: g32.append(x)
A30=list_plot(g30,color='red') A31=list_plot(g31,color='blue') A32=list_plot(g32,color='yellow') show(A30+A31+A32) It appears that whether x is 1 or 2 mod 3 makes little difference: however, being divisible by 3 makes a huge difference.

let's redo the plot with just 0 or not mod 3

gdiv=[] # points not divisible by 3 or 5 gdiv3=[] # points divisible by 3 not 5 gdiv5=[] # points divisible by 5 not 3 gdiv35=[] # points divisible by both 3 and 5 for x in goldbach_reps_pairs: a=mod(x,3) b=mod(x,5) if a==0: if b==0: gdiv35.append(x) else: gdiv3.append(x) else: if b==0: gdiv5.append(x) else: gdiv.append(x)
B=list_plot(gdiv, color='yellow') B3=list_plot(gdiv3, color='green') B5=list_plot(gdiv5, color='blue') B35=list_plot(gdiv35, color='red') show(B+B3+B5+B35) Let's redo this, but now, halve the value of the y co-ordinate when x is 0 mod 3.

gdiv=[] # points not divisible by 3 or 5 gdiv3=[] # points divisible by 3 not 5 gdiv5=[] # points divisible by 5 not 3 gdiv35=[] # points divisible by both 3 and 5 for x in goldbach_reps_pairs: a=mod(x,3) b=mod(x,5) if a==0: if b==0: gdiv35.append([x,x/3]) else: gdiv3.append([x,x/2]) else: if b==0: gdiv5.append([x,x*2/3]) else: gdiv.append(x)
B=list_plot(gdiv, color='yellow') B3=list_plot(gdiv3, color='green') B5=list_plot(gdiv5, color='blue') B35=list_plot(gdiv35, color='red') show(B+B3+B5+B35) The last experiment was a little more complicated than we should do yet.  Let's instead just focus on divisibility by 5.

g50=[] g5x=[] for x in goldbach_reps_pairs: a=mod(x,5) if a==0: g50.append([x,x/(5/4)]) else: g5x.append(x)
C0=list_plot(g50,color='blue') C1=list_plot(g5x,color='red') show(C0+C1) It looks as if 5/4 is too small a factor.  Let's try 4/3.

g50=[] g5x=[] for x in goldbach_reps_pairs: a=mod(x,5) if a==0: g50.append([x,x/(4/3)]) else: g5x.append(x)
C0=list_plot(g50,color='blue') C1=list_plot(g5x,color='red') show(C0+C1) show(C1+C0) It appears that divisibility by 5 increases the number of representations by 4/3.

Let's redo the plot with 3 and 5 corrected, with factors of 1, 2, 4/3, 8/3 respectively

gdiv=[] # points not divisible by 3 or 5 gdiv3=[] # points divisible by 3 not 5 gdiv5=[] # points divisible by 5 not 3 gdiv35=[] # points divisible by both 3 and 5 for x in goldbach_reps_pairs: a=mod(x,3) b=mod(x,5) if a==0: if b==0: gdiv35.append([x,x/(8/3)]) else: gdiv3.append([x,x/2]) else: if b==0: gdiv5.append([x,x/(4/3)]) else: gdiv.append(x) B=list_plot(gdiv, color='yellow') B3=list_plot(gdiv3, color='green') B5=list_plot(gdiv5, color='blue') B35=list_plot(gdiv35, color='red') show(B+B3+B5+B35) Let us define a correction factor, the Hardy Littlewood correction factor by

$hl(n) = \prod_{p|n} \frac{p-1}{p-2}$

where $p$ ranges over the odd prime factors of n.

m=9998
m.prime_factors(
 [2, 4999] [2, 4999]

It appears that the method .prime_divisors() is what we want to use.  Let's write the definition of hl(n) so that it would work even if n is odd (what we can do is compute it for 2*n, forcing the input to be even, and then always discard the 2).

def hl(n): m=Integer(2*n) pd= m.prime_factors() hlf=1 for p in pd[1:]: hlf=hlf*(p-1)/(p-2) return(hlf)
hl(14)
 6/5 6/5
hl(30)
 8/3 8/3
gbc= [] #goldbach corrected by hl(n) for x in goldbach_reps_pairs: #print(type(x)) #a=x #print( hl(a) ) gbc.append([x,x/hl(x)])
hl(6)
 2 2
m.prime_factors(
 Traceback (click to the left of this block for traceback) ... SyntaxError: invalid syntax Traceback (most recent call last): File "", line 1, in File "_sage_input_41.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("bS5wcmltZV9mYWN0b3JzKA=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmpPTbJVR/___code___.py", line 3 ^ SyntaxError: invalid syntax
m=3/2 m=m*2 type(m) m=Integer(m) type(m)
  
A=list_plot(gbc)
var('t')
 t t
B=plot(1.7*t/(log(t))^2,t, 2,100000,color='yellow')
show(A+B) var('t')
 t t
C=plot(1.7*t/(log(t))^2-2*sqrt(2*t)/log(t),t, 2,100000,color='red') D=plot(1.7*t/(log(t))^2+2*sqrt(2*t)/log(t),t, 2,100000,color='green') show(A+B+C+D) E=plot(t/log(t), t, 2,100000) F=plot(Li(t),t,2,100000,color='red') show(E+F) 