Nearly Square Divisors, Revisited
August 23, 2016
In a recent exercise, we discussed the problem of computing the nearly square divisors of a composite number n = a · b, with a ≥ b, such that the difference a − b is as small as possible. We gave two solutions: The first solution enumerated the divisors one-by-one, by trial division counting from 1, until reaching the square root, which is fast if n is small but slow in general. The second solution factored n, computed the divisors, the picked the two divisors in the middle of the list, which is fast in general but slow if n is highly composite and thus has a large number of divisors.
In a comment on that exercise, Matthew Arcus, who is a regular reader and commentor at Programming Praxis, gave a splendid answer to the problem; it relies on factoring n, but is fast even if n has a large number of divisors. His algorithm reduces the multiplication of the factors to addition of their logarithms, which means that the knapsack algorithm can be used to find the greatest sum less than the logarithm of the square root.
Your task is to implement Matthew’s algorithm to find the nearly-square divisors of a number. 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.
Programming Praxis 
How do this work? With logarithm?
You can use the knapsack algorithm with products, as well as with sums.
On a 64-bit machine all numbers are fixnums, so things go even a bit faster:
> (time (nsd-new 224403121196654400))
(time (nsd-new 224403121196654400))
63 ms real time
64 ms cpu time (64 user, 0 system)
no collections
3280 bytes allocated
no minor faults
no major faults
(473753280 473670855)
> (time (nsd 224403121196654400))
(time (nsd 224403121196654400))
148 ms real time
148 ms cpu time (148 user, 0 system)
77 collections accounting for 62 ms real time (56 user, 0 system)
165921824 bytes allocated
no minor faults
no major faults
(473753280 473670855)
And here’s some code (without sum or factors):
(define (multiply-knap target xs) (if (null? xs) 1 (if (< target (car xs)) (multiply-knap target (cdr xs)) (max (* (car xs) (multiply-knap (quotient target (car xs)) (cdr xs))) (multiply-knap target (cdr xs)))))) (define (nsd-new n) (let* ((fs (reverse (factors n))) (b (multiply-knap (integer-sqrt n) fs)) (a (quotient n b))) (list a b)))Here a description for a version of nsd, which is superfast. It solves Euler 266 in Python in 7 seconds! I implemented a C version of an improved “Matthew” algo, which solved it in about 100 minutes (a comparable Python version would have taken 100 hours).
Split the factors of number n in 2 equal parts, calculate all divisors for the 2 parts and sort both sets of divisors. Loop over both sets of divisors (one ascending and the other descending) and make sure the product is less equal isqrt(n). If the product is too high advance the descending, otherwise the ascending. Keep track of the highest prod below isqrt(n).
@Paul .. Could you give sample snippet or pseudo code of your idea? Especially about the splitting of divisors?
@Rudolf-san. Here is Python code. You need routines for integer square root (isqrt) and factoring (rho_factors).
In fact you do not split the divisors, but you split the factors and then calculate all divisors for them. Every divisor for n is a product of two divisors from either part. You are only interested in these products close to isqrt(n).
def divisors(facs): factors = [(1,) + tuple(accumulate(g, mul)) for _, g in groupby(facs)] # e.g.: [(1, 2), (1, 3), (1, 5) .... (1, 181)] 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