highly_abundant.py

import math

import mpmath
from fractions import Fraction
from sympy import lcm, divisor_sigma, isprime
mpmath.mp.dps = 100


N = 67 # N = 71 can be refuted with a numerator of 2007 and denominator of 2000
SCORE_NUMERATOR = 251  # original fraction had this 2009, but this works better
SCORE_DENOMINATOR = 250  # original fraction had this 2000, but this works better

print(1 + 1/(N * math.log(N)))
print(SCORE_NUMERATOR / SCORE_DENOMINATOR)

PRIME_BOUND = 2000
PREVIOUS_PRIMES = {}
for n in range(2, 1 + PRIME_BOUND):
    p = n
    while not isprime(p):
        p -= 1
    PREVIOUS_PRIMES[n] = p
LAST_PRIME = PREVIOUS_PRIMES[PRIME_BOUND]

# 5566277615616
# 2783138807808


LN = lcm([n for n in range(1, N + 1)])
SIGMA_LN = divisor_sigma(LN)

# old version of score. This works just fine and should probably be used if generalizing,
# but for the final run the lazy exact arithmetic solution was fast enough, and more convincing
def score(n):
    return mpmath.log(divisor_sigma(n)) - SCORE_NUMERATOR / SCORE_DENOMINATOR * mpmath.log(n)

# exact arithmetic version of score, using 1.0045 = 2009 / 2000
def big_score(n):
    return Fraction(divisor_sigma(n) ** SCORE_DENOMINATOR) / Fraction(n ** SCORE_NUMERATOR)


def score_upper_bound(p, k):
    return (k * (1 - SCORE_NUMERATOR / SCORE_DENOMINATOR) + 1) * mpmath.log(p) - mpmath.log(p - 1)

# exact arithmetic version of score upper bound, using 1.0045 = 2009 / 2000
def big_score_upper_bound(p, k):
    return Fraction(p ** (SCORE_DENOMINATOR - (SCORE_NUMERATOR - SCORE_DENOMINATOR)*k)) / (Fraction((p - 1) ** SCORE_DENOMINATOR))

# TODO quick hack to swap out the scores
# big_score = score
# big_score_upper_bound = score_upper_bound


SIZES = {}

HIGHEST_SCORE = 1
for p in range(1, PRIME_BOUND):
    if isprime(p):
        k = 0
        best = big_score(p**k), k
        while (big_score_upper_bound(p, k), k) >= best:
            k += 1
            best = max(best, (big_score(p**k), k))
        # print(p, best)
        HIGHEST_SCORE *= p ** best[1]
print("Found score maximizer:", HIGHEST_SCORE)
# print(score(HIGHEST_SCORE) - score(L67))

for p in range(1, PRIME_BOUND):
    if isprime(p):
        temp = HIGHEST_SCORE
        exponent = 0
        while temp % p == 0:
            temp //= p
            exponent += 1
        lower_bound = exponent
        while lower_bound > 0 and big_score(HIGHEST_SCORE // (p ** (exponent - (lower_bound - 1)))) >= big_score(LN):
            lower_bound -= 1

        upper_bound = exponent
        while big_score(HIGHEST_SCORE * p ** (upper_bound - exponent + 1)) >= big_score(LN):
            upper_bound += 1

        SIZES[p] = lower_bound, upper_bound

print(SIZES)
assert SIZES[LAST_PRIME] == (0, 0)  # I only search primes up to 200, so if generalizing one should change this bound


# The old score(n) function could be incorporated here to throw out solutions much faster
def search(n, sigma, p):
    if n > LN:  # stop searching this branch if the partial product is too large already.
        return
    if p < 2:
        if sigma > SIGMA_LN:
            assert divisor_sigma(n) == sigma
            print("="*30)
            print("NEW:", n, divisor_sigma(n))
            print("OLD:", LN, SIGMA_LN)
        return
    p = PREVIOUS_PRIMES[p]

    lb, ub = SIZES[p]
    assert sigma % ((p ** (lb + 1) - 1) // (p - 1) ) == 0
    base_sigma = sigma // ((p ** (lb + 1) - 1) // (p - 1) )

    # n *= p ** lb
    # sigma *= (p ** (lb + 1) - 1) // (p - 1)  # update sigma and n as if our number is divisible by p**lb
    for index in range(lb, ub + 1):
        # n is currtly divisible by p**index
        search(n, sigma, p - 1)
        n *= p
        if n > LN: # another early exit
            return
        sigma = p * sigma + base_sigma  # update sigma and n as if our number is divisible by p**(index + 1)

total = 1
for (l, u) in SIZES.values():
    total *= (u - l + 1)


print(f"Naive search checks {total} possibilities")
search_base = 1
for p, (l, u) in SIZES.items():
    search_base *= p ** l


search(search_base, divisor_sigma(search_base), LAST_PRIME + 1)
print("Done!")