2022.11.04 MATH 3600-001 Goldbach -- Sage

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[100000])

Traceback (click to the left of this block for traceback)
...
IndexError: list index out of range

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_6.py", line 10, in <module>
    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 <module>
    
  File "/tmp/tmplXD5N3/___code___.py", line 3, in <module>
    exec compile(u'print(goldbach_representations[_sage_const_100000 ])
  File "", line 1, in <module>
    
IndexError: list index out of range

How many representations do we get for even numbers?

print(goldbach_representations[78964])

Traceback (click to the left of this block for traceback)
...
IndexError: list index out of range

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_7.py", line 10, in <module>
    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 <module>
    
  File "/tmp/tmpptkO4P/___code___.py", line 3, in <module>
    exec compile(u'print(goldbach_representations[_sage_const_78964 ])
  File "", line 1, in <module>
    
IndexError: list index out of range

print(goldbach_representations[85000])

Traceback (click to the left of this block for traceback)
...
IndexError: list index out of range

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_8.py", line 10, in <module>
    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 <module>
    
  File "/tmp/tmpbfBBPE/___code___.py", line 3, in <module>
    exec compile(u'print(goldbach_representations[_sage_const_85000 ])
  File "", line 1, in <module>
    
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[1]> max_seen: max_seen=x[1] champions.append(x[0])

list_plot(champions)

len(champions)

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[0],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

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.

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[0],5) if a==0: g50.append([x[0],x[1]/(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[0],5) if a==0: g50.append([x[0],x[1]/(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[0],3) b=mod(x[0],5) if a==0: if b==0: gdiv35.append([x[0],x[1]/(8/3)]) else: gdiv3.append([x[0],x[1]/2]) else: if b==0: gdiv5.append([x[0],x[1]/(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[0])) #a=x[0] #print( hl(a) ) gbc.append([x[0],x[1]/hl(x[0])])

hl(6)

m.prime_factors(

Traceback (click to the left of this block for traceback)
...
SyntaxError: invalid syntax

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_41.py", line 10, in <module>
    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 <module>
    
  File "/tmp/tmpPTbJVR/___code___.py", line 3
    
                    ^
SyntaxError: invalid syntax

m=3/2 m=m*2 type(m) m=Integer(m) type(m)

<type 'sage.rings.integer.Integer'>

<type 'sage.rings.integer.Integer'>

A=list_plot(gbc)

var('t')

B=plot(1.7*t/(log(t))^2,t, 2,100000,color='yellow')

show(A+B)

var('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)

2022.11.04 MATH 3600-001 Goldbach

315 days ago by calkin