find-highly-composite-number.py

import math
from functools import reduce
from collections import namedtuple

def multiply_all(numbers):
    return reduce(int.__mul__, numbers)

# 1. Highly composite number consists only of the consecutive primes
# 2. Larger primes cannot have larger powers, since the number can always be reduced without effecting
#    number of divisors
primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101]
MAX_NUMBER_ALLOWED = multiply_all(primes)
Solution = namedtuple("Solution", "number, num_divisors, prime_powers")

def find_highly_composite_number(N, prime_index=0, max_allowed_prime_power=None):
    assert 0 <= N < MAX_NUMBER_ALLOWED, f"N={N}"  
    
    if N < primes[prime_index]:
        return Solution(number=1, num_divisors=1, prime_powers=[])
    
    prime = primes[prime_index]
    # Floating point precision errors might make the method less reliable for a very large number
    max_prime_power = math.floor(math.log(N) / math.log(prime))

    if max_allowed_prime_power is not None:
        max_prime_power = min(max_prime_power, max_allowed_prime_power)

    best_solution_num_divisors = 1
    best_solution_prime_powers = []
    
    for power in range(max_prime_power, 0, -1):
        remainder = N / prime**power

        partial_solution = find_highly_composite_number(
            remainder,
            prime_index=prime_index+1,
            max_allowed_prime_power=power,
        )
        
        # Number of divisors equals to the product of numbers which are one larger then power 
        # of each prime: (power_p1 + 1)*(power_p2 + 1)*...*(power_pn + 1)
        if (power + 1) * partial_solution.num_divisors >= best_solution_num_divisors:
            best_solution_prime_powers = [power] + partial_solution.prime_powers
            best_solution_num_divisors = (power + 1) * partial_solution.num_divisors

    number = multiply_all(prime**power for prime, power in zip(primes, best_solution_prime_powers))
    return Solution(number, best_solution_num_divisors, best_solution_prime_powers)

print(find_highly_composite_number(100))
print(find_highly_composite_number(1_000))
print(find_highly_composite_number(100_000))
print(find_highly_composite_number(1_000_000))
# print(find_highly_composite_number(2_000_000_000_000_000_000_000_000_000_000_000_000))