nicely-divisible-numbers.sage

## Find numbers n which have a good score sigma(n)/(n*log(log(n))):
## see https://twitter.com/johncarlosbaez/status/1149700802371608576

## Prime factorization of 2520:
bootstrap = {2: 3, 3: 2, 5: 1, 7: 1}

## Compute the score sigma(n)/(n*log(log(n))) from the prime factorization:
def goodness(d):
    n = prod([p^(d[p]) for p in d.keys()])
    sig = prod([(p^(d[p]+1)-1)/(p-1) for p in d.keys()])
    return N(sig/(n*log(log(n))),prec=1000)

## Compute the number n from its prime factorization:
def unfactor(d):
    n = prod([p^(d[p]) for p in d.keys()])
    return n

## Given a prime factorization, return a list of all candidate "mutations",
## obtained by adding 1 to one of the prime factor exponents (or adding one
## extra prime at the end):
def mutations(d):
    tab = []
    for p in d.keys():
        d2 = d.copy()
        d2[p] += 1
        tab.append(d2)
    p = next_prime(max(d.keys()))
    d2 = d.copy()
    d2[p] = 1
    tab.append(d2)
    return tab

## The bag is a dictionary with the number n itself as key.
## Each value is a list with the following items:
## [0] = the score sigma(n)/(n*log(log(n))) (used for sorting),
## [1] = the number n itself (same as key),
## [2] = the prime factorization (a dictionary),
## [3] = a boolean indicating whether mutations of this n have been added yet.
bag = dict([(unfactor(bootstrap), [goodness(bootstrap), unfactor(bootstrap), bootstrap, False])])
while True:
    tab = sorted([tmp for tmp in bag.values() if tmp[3]==False])
    bst = tab[len(tab)-1]
    # bst is the highest-scoring yet-unmutated number
    print "About to explore %d (goodness=%f)"%(bst[1],bst[0])
    mut = bst[2]
    for urf in mutations(mut):
        n = unfactor(urf)
        if not bag.has_key(n):
            g = goodness(urf)
            bag[n] = [g, n, urf, False]
    bag[bst[1]][3] = True