Goldbach function

goldbach_function.py

import math
import numpy as np
def goldbach(n):
    # for n>17 goldbach function returns vector G with G[i]=g(2i) for  2<=i<=n/2
    Primes5mod6 =np.ones((n//6+1), dtype=bool)
    Primes1mod6 =np.ones((n//6+1), dtype=bool)
    imax=math.isqrt(n)//6+1
    for i in range(1,imax+1):
        if Primes5mod6[i]:
            Primes5mod6[6*i*i::6*i-1]=[False]
            Primes1mod6[6*i*i-2*i::6*i-1]=[False]
        if Primes1mod6[i]:
            Primes5mod6[6*i*i::6*i+1]=[False]
            Primes1mod6[6*i*i+2*i::6*i+1]=[False]
    
    G=[0,0,1,1,1,2]
    for ni in range(12,n-3,6):
        k=ni//6
        G.append(np.count_nonzero(Primes5mod6[1:k]&Primes1mod6[k-1:0:-1]))
        G.append(np.count_nonzero(Primes1mod6[1:math.floor(k/2)+1]&Primes1mod6[k-1:math.ceil(k/2)-1:-1])+Primes5mod6[k])
        G.append(np.count_nonzero(Primes5mod6[1:math.floor((k+1)/2)+1]&Primes5mod6[k:math.ceil((k+1)/2)-1:-1])+Primes1mod6[k])
    return G