What is the best algorithm for finding Goldbug Numbers?

Question

A Goldbug Number of order k is an even number 2n for which there exists some k order subset of the prime non-divisors of n $2 < p_1 < p_2 < p_3 < \cdots < p_k < n$ such that $(2n-p_1)(2n-p_2)(2n-p_3)\cdots (2n-p_k)$ only has $p_1,p_2,p_3,...,p_k$ as factors.

More information can be found at the related OEIS entry. https://oeis.org/A306746

Answer 1

I am posting the code from the forum written by UAU. This is the best implementation to date and identifies 6 Goldbug numbers under 100k. Keep in mind for speed this is searching for the maximal set and there can be subsets of the maximal set which satisfy the Goldbug property. I have run the code and verified it produces 6 Goldbugs under 100k (128,1718,1862,1928,2200,6142). For more details see the original post:

https://www.mersenneforum.org/showpost.php?p=501714&postcount=23

Thanks UAU!

#!/usr/bin/python3
import sys

def sieve(n):
    n //= 2
    r = [[] for _ in range(n)]
    for x in range(1, n):
        if r[x]:
            continue
        r[x] = None
        y = 3*x+1
        while y < n:
            r[y].append(2*x+1)
            y += x+x+1
    return r

def check(n, r):
    s = {}
    for i in range(3, n//2, 2):
        if r[i>>1] is not None or not n % i:
            continue
        for f in r[(n-i)>>1] or [n-i]:
            s.setdefault(f, []).append(i)
    bad = [x for x in s if x >= n//2]
    bad_seen = set(bad)
    ok = set(i for i in range(3, n // 2, 2) if r[i>>1] is None and n % i)
    while bad:
        x = s.get(bad.pop(), ())
        for p in x:
            if p not in bad_seen:
                bad.append(p)
                bad_seen.add(p)
                ok.discard(p)
    return ok

def search(lim):
    r = sieve(lim)
    for i in range(2, lim, 2):
        x = check(i, r)
        if x:
             print(i, sorted(x))

search(int(sys.argv[1]))

Answer 2

I am posting sage code (rather than pure py code), see also

www.sagemath.org

which is testing if a given number is a Goldburg number, corrected version.

(In the first version primes till $N=2n$ were included, i definitively need glasses. Sorry for the first version.)

Python is maybe too weak to start programming number theory from the scratch in it. This is the reason why using sage, which is python plus batteries. Here is my code, second try, some Goldburg numbers were found...

def deliverGoldburgSolution(N, k):
    """The number N shoud be even.
    This routine tests if there are prime numbers p1, p2, ... , pk
    2 < p1 < p2 < ... < pk < n (*** bug fix ***) 
    so that for any pj in the list 
    N - pj
    involves as factors only the primes p1, p2, ... , pk themselves again.
    This is a lazy quick implementation.
    """

    n = ZZ(N/2)
    D = Set(n.divisors())

    # *** bug fix below w.r.t. my first post *** using only primes < n
    # (the older code went till N)
    allowedPrimes = Set([ p for p in primes(3, n) ])

    for p in allowedPrimes:
        P = [p, ]    # init start, we soon extend
        stillSearch = True

        while stillSearch:
            stillSearch = False
            P_New = [ prime_number
                      for q in P
                      for prime_number in (N-q).prime_divisors() ]

            P_New = list(Set(P_New + P))

            if ( len(P_New) > k
                 or D.intersection(Set(P_New))
                 or not Set(P_New).issubset(allowedPrimes)
                 ):
                continue    # with the next p

            # else
            if len(P) == len(P_New):
                if len(P) == k:
                    # print "Solution: {}".format(P)
                    return P

            else:
                P = P_New[:]
                stillSearch = True

for N in range(2, 10000, 2):
    for k in (2,3):
        sol = deliverGoldburgSolution(N, k)
        if sol:
            print( "GOLDBURG NUMBER: N={} k={} SOLUTION FOR THE PRIMES {}"
                   .format(N, k, sol) )
    if N % 500 == 0:
        print "\t... N={} so far".format(N)

Results:

        ... N=500 so far
        ... N=1000 so far
        ... N=1500 so far
        ... N=2000 so far
GOLDBURG NUMBER: N=2200 k=2 SOLUTION FOR THE PRIMES [3, 13]
        ... N=2500 so far
        ... N=3000 so far
        ... N=3500 so far
        ... N=4000 so far
        ... N=4500 so far
        ... N=5000 so far
        ... N=5500 so far
        ... N=6000 so far
        ... N=6500 so far
        ... N=7000 so far
        ... N=7500 so far
        ... N=8000 so far
        ... N=8500 so far
        ... N=9000 so far
        ... N=9500 so far

So there is only one $k$-Goldburg number, $k=2,3$ (allowed) up to $9999$, which is $2200$ for $k=2$ and the primes $3,13$.

The code is no longer building "clusters", so that in the "peculiar case" a number is both a $k_1$- and a $k_2$-Goldburg number for different sets of primes (so that it is then also a $(k_1+k_2)$-Goldburg number) the code will not detect it.