Nearly Square Divisors, Revisited
August 23, 2016
The knap function takes a target, which is the logarithm of the square root of n, and a list xs of the logarithms of the factors of n, with multiplicity, in descending order, and returns the logarithm of the largest divisor of n less than the logarithm of the square root. The function is recursive, with the target decreasing and logarithms of the factors being removed from the list. The base case occurs when the list of logarithms of the factors of n is empty, when the function returns 0. If the largest logarithm is larger than the target, the recursion continues with the remaining logarithms. Otherwise, the algorithm returns the maximum of two values: the sum of the largest logarithm and the reduced knapsack problem with the largest logarithm subtracted from the target and removed from the list, or the recursion with the remaining logarithms. The code is simpler than the explanation:
(define (knap target xs)
(if (null? xs) 0
(if (< target (car xs))
(knap target (cdr xs))
(max (+ (car xs) (knap (- target (car xs)) (cdr xs)))
(knap target (cdr xs))))))
Then the nearly square divisors are computed by factoring n, arranging the logarithms in descending order, calculating the target as half the sum of the logarithms (which is the logarithm of the square root of n), calling the knap function, taking the anti-logarithm, and finally dividing to find the other divisor:
(define (nsd n)
(let* ((fs (map log (reverse (factors n))))
(b (inexact->exact (round (exp
(knap (/ (sum fs) 2) fs)))))
(a (/ n b)))
(list a b)))
Here’s the example from the previous exercise:
> (time (nsd 224403121196654400))
(time (nsd 224403121196654400))
8 collections
0.121862989s elapsed cpu time, including 0.000410511s collecting
0.124334069s elapsed real time, including 0.000441672s collecting
65259936 bytes allocated, including 67360576 bytes reclaimed
(473753280 473670855)
This sadly runs into a problem; the logarithms lose precision, and if n is too large, the program fails to find the correct solution. It will take something like the unlimited-precision floating-point functions in gmp to make it work for very large numbers.
You can run the program at http://ideone.com/n5Wcja, where you will also see the sum and factors functions.
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