Programming Praxis

The nearly square divisors exercise has generated a considerable amount of interest, and several excellent solutions in the comments. We looked previously at Matthew Arcus’ solution using a knapsack algorithm, with logarithms to reduce the problem from multiplication to addition, thus allowing languages like C to solve the problem using their native data types instead of switching to big integers.

The knapsack solution works like this: To find the nearly square divisor of n, factor n and form the list of divisors ds. Then use the subset sum algorithm of a previous exercise (a variant of knapsack), taking products rather than sums, to find the maximal product of divisors less than the square root of n. There are several ways to solve the subset sum problem. The standard solution uses dynamic programming. Another solution splits the problem space into two parts. Both solutions require exponential time, though the meet-in-the-middle solution has better asymptotic time of O(n2^n/2). We also studied a solution that takes polynomial time to produce an approximate answer to the subset sum problem, although that solution is not helpful to us because it already takes exponential time to calculate all the divisors.

Today we look at another solution buried in the comments, this one from Paul Hofstra. Here is his solution:

def divisors(facs):
    factors = [(1,) + tuple(accumulate(g, mul)) for _, g in groupby(facs)]
    div = [1]
    for g in factors:
        div = [d * f for d in div for f in g]
    return div
 
def nsd5(number):
    """ output: nearly square divisor
        method: split factors in 2 equal parts
                create divisors with small factors and sort (descending)
                create divisors with large factors (<= ulimit) and sort
                loop over large divisors and small divisors and search
                   for highest product <= ulimit
    """
    ulimit = isqrt(number)
    facs = rho_factors(number)
    mid = len(facs) // 2
    descending = reversed(sorted(divisors(facs[:mid])))
    ascending = iter(sorted(divisors(facs[mid:])))
    best = 0
    desc = next(descending)
    while True:
        for asc in ascending:
            prod = asc * desc
            if prod > best:
                if prod <= ulimit:
                    best = prod
                else:
                    break
        else:
            break  # while
        for desc in descending:
            prod = asc * desc
            if prod <= ulimit:
                if prod > best:
                    best = prod
                break
        else:
            break  # while
    return best

With this solution we are back to using big integers, and we are using the meet-in-the-middle variant of the subset sum solution to calculate the answer. The solution finds the factors, splits them into two halves, calculates the divisors of each half, then uses subset sum to find the maximal divisor less than the square root. Note that we don’t have to compute all the divisors, only the divisors of each of the two halves.

Your task is to implement Hofstra’s solution to the nearly square divisor problem and use it to calculate the nearly square divisor of the product of the primes less than 190. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below:

Pages: 1 2