Programming Praxis

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.

Pages: 1 2