2022.11.04 MATH 3600-002 Goldbach

315 days ago by calkin

Goldbach's conjecture: every even number greater than 4 can be written as the sum of two odd primes.

 

Let's start by creating a list of odd primes up to N.
N=20 p=3 prime_list=[] while p< N: if is_prime(p): prime_list.append(p) p=p+1 
       
print(prime_list) 
        
[3, 5, 7, 11, 13, 17, 19]
[3, 5, 7, 11, 13, 17, 19]

We could do better.  For a start, we could increment p by 2 instead of 1 each time.
N=100 p=3 prime_list=[] while p< N: if is_prime(p): prime_list.append(p) p=p+2 
       
print(prime_list) 
        
[3, 5, 7, 11, 13, 17, 19]
[3, 5, 7, 11, 13, 17, 19]

We can do even better using the next_prime() method or function.
N=100000 p=3 prime_list=[] while p< N: if is_prime(p): prime_list.append(p) p=p.next_prime() 
       
print(len(prime_list)^2) 
        
1507984
1507984
t=cputime() goldbach_representations=[0 for i in srange(2*N+1)] for p in prime_list: for q in prime_list: goldbach_representations[p+q]+=1 print(cputime()-t) 
        
38.84908
38.84908
#print(goldbach_representations) 
       
plot(goldbach_representations) 
       
list_plot(goldbach_representations) 
       
list_plot(goldbach_representations[:N]) 
       

Let's figure out which points on the plot are champions, i.e. are higher than any point on its left.
gb_reps=[] for i in srange(N): if goldbach_representations[i]>0: gb_reps.append([i,goldbach_representations[i]]) 
       
max_seen=0 champions=[] for x in gb_reps: if x[1]>max_seen: max_seen=x[1] champions.append(x) 
       
len(champions) 
        
52
52
print(champions) 
        
[[6, 1], [8, 2], [10, 3], [16, 4], [22, 5], [24, 6], [34, 7], [36, 8],
[48, 10], [60, 12], [78, 14], [84, 16], [90, 18], [114, 20], [120, 24],
[168, 26], [180, 28], [210, 38], [300, 42], [330, 48], [390, 54], [420,
60], [510, 64], [630, 82], [780, 88], [840, 102], [990, 104], [1050,
114], [1140, 116], [1260, 136], [1470, 146], [1650, 152], [1680, 166],
[1890, 182], [2100, 194], [2310, 228], [2730, 256], [3150, 276], [3570,
308], [3990, 326], [4200, 330], [4410, 342], [4620, 380], [5250, 396],
[5460, 436], [6090, 444], [6510, 482], [6930, 536], [7980, 548], [8190,
584], [9030, 606], [9240, 658]]
[[6, 1], [8, 2], [10, 3], [16, 4], [22, 5], [24, 6], [34, 7], [36, 8], [48, 10], [60, 12], [78, 14], [84, 16], [90, 18], [114, 20], [120, 24], [168, 26], [180, 28], [210, 38], [300, 42], [330, 48], [390, 54], [420, 60], [510, 64], [630, 82], [780, 88], [840, 102], [990, 104], [1050, 114], [1140, 116], [1260, 136], [1470, 146], [1650, 152], [1680, 166], [1890, 182], [2100, 194], [2310, 228], [2730, 256], [3150, 276], [3570, 308], [3990, 326], [4200, 330], [4410, 342], [4620, 380], [5250, 396], [5460, 436], [6090, 444], [6510, 482], [6930, 536], [7980, 548], [8190, 584], [9030, 606], [9240, 658]]
list_plot(champions) 
       

It looks like being divisible by small primes 3,5,7 etc has an impact on the number of representations.  How much of an effect does divisibility by 3 have?
gb30=[] gb31=[] gb32=[] for x in gb_reps: a=mod(x[0],3) if a==0: gb30.append(x) if a==1: gb31.append(x) if a==2: gb32.append(x) 
       
A30=list_plot(gb30,color='blue') A31=list_plot(gb31,color='green') A32=list_plot(gb32,color='orange') show(A30+A32+A31) 
       
show(A30+A32+A31) 
       

Let's redo the plots, but divide the height by 2 when n is divisible by 3.
gb30=[] gb31=[] gb32=[] for x in gb_reps: a=mod(x[0],3) if a==0: gb30.append([x[0],x[1]/2]) if a==1: gb31.append(x) if a==2: gb32.append(x) 
       
A30=list_plot(gb30,color='blue') A31=list_plot(gb31,color='green') A32=list_plot(gb32,color='orange') show(A30+A31+A32) 
       
show(A32+A31+A30) 
       

Let's try something similar for divisibility by 5
 
       
gb50=[] gb51=[] gb52=[] gb53=[] gb54=[] for x in gb_reps: a=mod(x[0],5) if a==0: gb50.append([x[0],x[1]*3/4]) if a==1: gb51.append(x) if a==2: gb52.append(x) if a==3: gb53.append(x) if a==4: gb54.append(x) 
       
A50=list_plot(gb50,color='blue') A51=list_plot(gb51,color='green') A52=list_plot(gb52,color='orange') A53=list_plot(gb53,color='red') A54=list_plot(gb54,color='yellow') show(A50+A51+A52+A53+A54) 
       

It appears experimentally that divisibility by 3 increases the number of representations by a factor of 2, and divisibility by 5 increases by a factor of 4/3.

Let's combine them.
gb=[] #not divisible by 3 or 5 gb3=[] #divisible by 3, not 5 gb5=[] #divisible by 5, not 3 gb35=[] #divisible by 3 and 5 for x in gb_reps: a=mod(x[0],3) b=mod(x[0],5) if a==0: if b==0: gb35.append(x) else: gb3.append(x) else: if b==0: gb5.append(x) else: gb.append(x) 
       
C=list_plot(gb,color='red') C3=list_plot(gb3,color='yellow') C5=list_plot(gb5,color='green') C35=list_plot(gb35,color='blue') show(C+C3+C5+C35) 
       

Let's take a wild guess that the correction for divisibilty by $p$ turns out to be $(p-1)/(p-2)$, and that the corrections act roughly independently:

define $hl(n)$ by 

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

over odd primes dividing $n$.
def hl(n): m=2*n pd=m.prime_divisors() hlf=1 for p in pd[1:]: hlf=hlf*(p-1)/(p-2) return(hlf) 
       
hl(105) 
        
16/5
16/5
n.prime_divisors( 
        
<type 'sage.rings.integer.Integer'>
<type 'sage.rings.integer.Integer'>
gbc = [] # Goldbach corrected for x in gb_reps: gbc.append([x[0], x[1]/hl(x[0])]) 
       
A=list_plot(gb_reps,color='red') B=list_plot(gbc, color='blue') show(A+B) 
       
A=list_plot(gbc) 
       
var('t') a=.75 B=plot(a*f(t),t,2,100000,color='red') 
       
show(A+B) 
       

Let's throw in some upper and lower curves on either side, to see if they look about the right distance away from the curve.
B=plot(a*f(t),t,2,100000,color='red') C=plot(a*t/log(t)^2-2*sqrt(a*t)/log(t),t,2,100000,color='red') D=plot(a*t/log(t)^2+2*sqrt(a*t)/log(t),t,2,100000,color='red') show(A+B+C+D) 
       

Let's try Li(x) instead of x/log(x).
C=plot(a*t/log(t)^2-2*sqrt(a*t)/log(t),t,2,100000,color='red') D=plot(a*t/log(t)^2+2*sqrt(a*t)/log(t),t,2,100000,color='red') show(A+B+C+D) 
       
def f(t): return(Li(t)*2*(Li(2*t)-Li(t))/t ) 
       
C=plot(a*f(t)-2*sqrt(f(t)),t,2,100000,color='red') D=plot(a*f(t)+2*sqrt(f(t)),t,2,100000,color='red') show(A+B+C+D) 
       
f(10) 
        
Traceback (click to the left of this block for traceback)
...
TypeError: 'sage.rings.integer.Integer' object is not callable
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_48.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("ZigxMCk="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpoCHWmL/___code___.py", line 3, in <module>
    exec compile(u'f(_sage_const_10 )
  File "", line 1, in <module>
    
  File "/tmp/tmpKFp1ox/___code___.py", line 4, in f
    return(Li(t)*_sage_const_2 (Li(_sage_const_2 *t)-Li(t))/t   )
TypeError: 'sage.rings.integer.Integer' object is not callable