Nearly Square Divisors, Meet In The Middle
September 6, 2016
Before we look at code, let’s look at an example so we are sure to have the algorithm clearly fixed in our minds; we calculate the nearly square divisor of 360 = 2 * 2 * 2 * 3 * 3 * 5. There are six factors, so we split them into halves and calculate their divisors: the low half are the factors 2, 2 and 2, with divisors 1, 2, 4 and 8, and the high half are the factors 3, 3, and 5, with divisors 45, 15, 9, 5, 3 and 1 in descending order. The square root of 360 is approximately 18.973666, so the upper limit on the nearly square divisor is 18.
Now we iterate the meet-in-the-middle algorithm. We repeatedly multiply the head of the low divisors list by the head of the high divisors list. If the product exceeds the limit, we discard the head of the high divisors list. Otherwise we discard the head of the low divisors list, after first checking if the product is a new maximum. The table shows an index, the current maximum, the low divisors, and the high divisors in reverse order:
0 0 (1 2 4 8) (45 15 9 5 3 1) 1 0 (1 2 4 8) (15 9 5 3 1) 2 15 (2 4 8) (15 9 5 3 1) 3 15 (2 4 8) (9 5 3 1) 4 18 (4 8) (9 5 3 1) 5 18 (4 8) (5 3 1) 6 18 (4 8) (3 1) 7 18 (8) (3 1) 8 18 (8) (1) 9 18 () (1)
Row 0 shows the initial state, with a current maximum of 0 and the two lists of divisors. The first calculation takes 1 × 45 = 45, which exceeds the limit, so Row 1 shows 45 removed from the list of high divisors. The second calculation takes 1 × 15 = 15, which doesn’t exceed the limit, so Row 2 shows a new best of 15 and 1 removed from the list of low divisors. The third calculation takes 2 × 15 = 30, which exceeds the limit, so 15 is removed from the list of high divisors. The fourth calculation takes 2 × 9 = 18, which doesn’t exceed the limit, so we save the new best and remove 2 from the list of low divisors. And so on, stopping after the list of low divisors is exhausted, The final solution is 18. Note that 360 has 24 divisors, but only 9 are used (divisor 1 is used twice); the divisors 6, 10, 12, 18, 20, 24, 30, 36, 40, 60, 72, 90, 120, 180 and 360 never appear.
Now we can look at our version of Hofstra’s solution:
(define (nsd n) (let* ((limit (isqrt n)) (fs (factors n)) (mid (quotient (length fs) 2)) (los (divisors (apply * (take mid fs)))) (his (reverse (divisors (apply * (drop mid fs)))))) (let loop ((best 0) (los los) (his his)) (if (or (null? los) (null? his)) best (let ((t (* (car los) (car his)))) (cond ((< limit t) (loop best los (cdr his))) ((< best t) (loop t (cdr los) his)) (else (loop best (cdr los) his)))))))) > (time (nsd (apply * (primes 190)))) (time (nsd (apply * ...))) 67 collections 2.066603685s elapsed cpu time, including 0.950848951s collecting 2.072700979s elapsed real time, including 0.954510719s collecting 626128608 bytes allocated, including 601787408 bytes reclaimed 2323218950659046189161096883702440585
This is reasonably fast, even when n has a large number of factors (which is a relatively rare occurrence). First, n must be factored, which may be the hard part of the algorithm. Then the two lists of divisors are created. The iteration removes one divisor at each step, and performs at most two comparisons. Solving Project Euler 266 takes about 2 seconds on my machine. Sloane’s sequence A060798 finds the nearly square divisors of the primorials:
(define (primorial-nsd n) (let ((p 1)) (do ((ps (primes n) (cdr ps)) (i 1 (+ i 1))) ((null? ps)) (set! p (* p (car ps))) (display i) (display ": ") (display (nsd p)) (newline)))) > (primorial-nsd 190) 1: 1 2: 2 3: 5 4: 14 5: 42 6: 165 7: 714 8: 3094 9: 14858 10: 79534 11: 447051 12: 2714690 13: 17395070 14: 114371070 15: 783152070 16: 5708587335 17: 43848093003 18: 342444658094 19: 2803119896185 20: 23619540863730 21: 201813981102615 22: 1793779293633437 23: 16342050964565645 24: 154170926013430326 25: 1518409177581024365 26: 15259825120722837478 27: 154870329715038713659 28: 1601991088154989174258 29: 16725278985898957454695 30: 177792163538134124432895 31: 2003615963401499430234017 32: 22932432897807928504246930 33: 268417244825173788264575798 34: 3164592642338712034986945585 35: 38628776202949578678974208410 36: 474678348768025365932004176547 37: 5947702664563379291963314981629 38: 75935184293205358249943157743631 39: 981298843695756278498820457991905 40: 12906971131408359840474658782214165 41: 172683504643732793885778158650669178 42: 2323218950659046189161096883702440585
You can run the program at http://ideone.com/V0KZDB, where you will also see the auxiliary functions that are used.
I’d come across this way of solving the problem but have to say I’m impressed with how much faster than branch-and-bound it is.
How about a simulated annealing solution?