Highly Composite Numbers, Revisited

August 16, 2016

We begin with a list of the prime numbers less than 10,000 in ps, and we assume all of the code of the distinct priority queue is available. The algorithm is exactly identical to the algorithm of the recent exercise.

Function (next ndxs) takes an ndxs list and returns a list of possible successor ndxs lists by incrementing each exponent in turn, if it is possible to do so, including the first (implied) zero exponent following the end of the exponent list. This is exactly the same algorithm as the next function of the previous exercise, except that it returns an ndxs instead of just an n:

(define (next ndxs)
  (let loop ((front (list)) (back (cddr ndxs)) (xs (list)))
    (cond ((null? back)
            (map (lambda (x)
                   (cons (do ((x x (cdr x))
                              (ps ps (cdr ps))
                              (p 1 (* p (expt (car ps) (car x)))))
                             ((null? x) p))
                         (cons (apply * (map add1 x))
                               x)))
                 (cons (reverse (cons 1 front)) xs)))
          ((null? front)
            (loop (cons (car back) front) (cdr back)
                  (cons (cons (+ (car back) 1) (cdr back)) xs)))
          ((< (car back) (car front))
            (loop (cons (car back) front) (cdr back)
                  (cons (append (reverse (cons (+ (car back) 1) front))
                                (cdr back))
                        xs)))
          (else (loop (cons (car back) front) (cdr back) xs)))))

Here’s an example:

> (next '(293318625600 5040 6 4 2 2 1 1 1 1))
 ((6746328388800 10080 6 4 2 2 1 1 1 1 1)
  (3226504881600 7560 6 4 2 2 2 1 1 1)
  (1466593128000 6720 6 4 3 2 1 1 1 1)
  (879955876800 6048 6 5 2 2 1 1 1 1)
  (586637251200 5760 7 4 2 2 1 1 1 1))

Now the algorithm to compute highly composite numbers is also exactly the same as the algorithm of the previous exercise: Initialize a distinct priority queue with the ndxs (1 1), the number 1 with 1 divisor, and also a counter ctr to track the index of the highly composite numbers and a variable record that tracks the current record for the number of divisors of the highly composite number. Repeatedly pop the first item in the distinct priority queue; if it creates a new record, update the counter and record and print the item. Then, whether or not a new record has been found, compute the next possible highly composite numbers and, for any that are greater than the current record, add them to the distinct priority queue. This solves two problems with the previous exercise: first, items that could not possibly be records were added to the queue, and second, multiple copies of the same item could be added to the queue. Here’s the code:

(define (hcn)
  (let ((dpq (make-dpq (lambda (xs ys) (< (car xs) (car ys)))))
        (ctr 0) (record 0))
    (let loop1 ((dpq (dpq-insert (list 1 1) dpq)))
      (let* ((candidate (dpq-first dpq)) (dpq (dpq-rest dpq)))
        (when (< record (cadr candidate))
          (set! ctr (+ ctr 1)) (set! record (cadr candidate))
          (display ctr) (display " ") (display candidate) (newline))
        (let loop2 ((xs (next candidate)) (dpq dpq))
          (cond ((null? xs) (loop1 dpq))
                ((< record (cadar xs))
                  (loop2 (cdr xs) (dpq-insert (car xs) dpq)))
                (else (loop2 (cdr xs) dpq))))))))

When you run that, highly composite numbers flash past on the screen faster than you can see them; here are the first hundred (we started counting from 1, as some people do; other people start counting from 0, so depending on where you look, your index might be off by 1):

> (hcn)
1 (1 1)
2 (2 2 1)
3 (4 3 2)
4 (6 4 1 1)
5 (12 6 2 1)
6 (24 8 3 1)
7 (36 9 2 2)
8 (48 10 4 1)
9 (60 12 2 1 1)
10 (120 16 3 1 1)
11 (180 18 2 2 1)
12 (240 20 4 1 1)
13 (360 24 3 2 1)
14 (720 30 4 2 1)
15 (840 32 3 1 1 1)
16 (1260 36 2 2 1 1)
17 (1680 40 4 1 1 1)
18 (2520 48 3 2 1 1)
19 (5040 60 4 2 1 1)
20 (7560 64 3 3 1 1)
21 (10080 72 5 2 1 1)
22 (15120 80 4 3 1 1)
23 (20160 84 6 2 1 1)
24 (25200 90 4 2 2 1)
25 (27720 96 3 2 1 1 1)
26 (45360 100 4 4 1 1)
27 (50400 108 5 2 2 1)
28 (55440 120 4 2 1 1 1)
29 (83160 128 3 3 1 1 1)
30 (110880 144 5 2 1 1 1)
31 (166320 160 4 3 1 1 1)
32 (221760 168 6 2 1 1 1)
33 (277200 180 4 2 2 1 1)
34 (332640 192 5 3 1 1 1)
35 (498960 200 4 4 1 1 1)
36 (554400 216 5 2 2 1 1)
37 (665280 224 6 3 1 1 1)
38 (720720 240 4 2 1 1 1 1)
39 (1081080 256 3 3 1 1 1 1)
40 (1441440 288 5 2 1 1 1 1)
41 (2162160 320 4 3 1 1 1 1)
42 (2882880 336 6 2 1 1 1 1)
43 (3603600 360 4 2 2 1 1 1)
44 (4324320 384 5 3 1 1 1 1)
45 (6486480 400 4 4 1 1 1 1)
46 (7207200 432 5 2 2 1 1 1)
47 (8648640 448 6 3 1 1 1 1)
48 (10810800 480 4 3 2 1 1 1)
49 (14414400 504 6 2 2 1 1 1)
50 (17297280 512 7 3 1 1 1 1)
51 (21621600 576 5 3 2 1 1 1)
52 (32432400 600 4 4 2 1 1 1)
53 (36756720 640 4 3 1 1 1 1 1)
54 (43243200 672 6 3 2 1 1 1)
55 (61261200 720 4 2 2 1 1 1 1)
56 (73513440 768 5 3 1 1 1 1 1)
57 (110270160 800 4 4 1 1 1 1 1)
58 (122522400 864 5 2 2 1 1 1 1)
59 (147026880 896 6 3 1 1 1 1 1)
60 (183783600 960 4 3 2 1 1 1 1)
61 (245044800 1008 6 2 2 1 1 1 1)
62 (294053760 1024 7 3 1 1 1 1 1)
63 (367567200 1152 5 3 2 1 1 1 1)
64 (551350800 1200 4 4 2 1 1 1 1)
65 (698377680 1280 4 3 1 1 1 1 1 1)
66 (735134400 1344 6 3 2 1 1 1 1)
67 (1102701600 1440 5 4 2 1 1 1 1)
68 (1396755360 1536 5 3 1 1 1 1 1 1)
69 (2095133040 1600 4 4 1 1 1 1 1 1)
70 (2205403200 1680 6 4 2 1 1 1 1)
71 (2327925600 1728 5 2 2 1 1 1 1 1)
72 (2793510720 1792 6 3 1 1 1 1 1 1)
73 (3491888400 1920 4 3 2 1 1 1 1 1)
74 (4655851200 2016 6 2 2 1 1 1 1 1)
75 (5587021440 2048 7 3 1 1 1 1 1 1)
76 (6983776800 2304 5 3 2 1 1 1 1 1)
77 (10475665200 2400 4 4 2 1 1 1 1 1)
78 (13967553600 2688 6 3 2 1 1 1 1 1)
79 (20951330400 2880 5 4 2 1 1 1 1 1)
80 (27935107200 3072 7 3 2 1 1 1 1 1)
81 (41902660800 3360 6 4 2 1 1 1 1 1)
82 (48886437600 3456 5 3 2 2 1 1 1 1)
83 (64250746560 3584 6 3 1 1 1 1 1 1 1)
84 (73329656400 3600 4 4 2 2 1 1 1 1)
85 (80313433200 3840 4 3 2 1 1 1 1 1 1)
86 (97772875200 4032 6 3 2 2 1 1 1 1)
87 (128501493120 4096 7 3 1 1 1 1 1 1 1)
88 (146659312800 4320 5 4 2 2 1 1 1 1)
89 (160626866400 4608 5 3 2 1 1 1 1 1 1)
90 (240940299600 4800 4 4 2 1 1 1 1 1 1)
91 (293318625600 5040 6 4 2 2 1 1 1 1)
92 (321253732800 5376 6 3 2 1 1 1 1 1 1)
93 (481880599200 5760 5 4 2 1 1 1 1 1 1)
94 (642507465600 6144 7 3 2 1 1 1 1 1 1)
95 (963761198400 6720 6 4 2 1 1 1 1 1 1)
96 (1124388064800 6912 5 3 2 2 1 1 1 1 1)
97 (1606268664000 7168 6 3 3 1 1 1 1 1 1)
98 (1686582097200 7200 4 4 2 2 1 1 1 1 1)
99 (1927522396800 7680 7 4 2 1 1 1 1 1 1)
100 (2248776129600 8064 6 3 2 2 1 1 1 1 1)

That closes the book on highly composite numbers, at least as far as I am concerned; I am satisfied that we finally have a good solution. I am fascinated that a clever algorithm and a proper data structure can accomplish so much in such little time. The distinct priority queue is critical to the success of this algorithm. When we used exactly the same algorithm in a previous exercise, but using a normal priority queue instead of the distinct priority queue, we quickly reached a breaking point due to lack of memory. Now we are able to compute highly composite numbers to a very high limit, which is a very satisfying result.

You can run the program at http://ideone.com/s1vq5T.

Pages: 1 2

10 Responses to “Highly Composite Numbers, Revisited”

  1. Paul said

    I implemented the proposed scheme in Python. It works indeed superfast. I am not sure if I made a mistake, but I notice that I am missing some hcn’s. The first is at n=815. The correct hcn is 10 6 4 2^5 1^25 (notation of , but the answer I get is 10 5 4 3 2^2 1^28 (which nr 816 on the “official” list). Could you please, run your code and check, if I made a mistake?

  2. Paul said

    I implemented the proposed scheme in Python. It works indeed superfast. I am not sure if I made a mistake, but I notice that I am missing some hcn’s. The first is at n=815. The correct hcn is 10 6 4 2^5 1^25 (notation of , but the answer I get is 10 5 4 3 2^2 1^28 (which nr 816 on the “official” list). Could you please, run your code and check, if I made a mistake?
    (I hope this format better)

  3. Paul said

    I implemented the proposed scheme in Python. It works indeed superfast. I am not sure if I made a mistake, but I notice that I am missing some hcn’s. The first is at n=815. The correct hcn is 10 6 4 2^5 1^25 (notation of , but the answer I get is 10 5 4 3 2^2 1^28 (which is nr 816 on the “official” list). Could you please, run your code and check, if I made a mistake?
    (I hope this format better and sorry for creating this mess)

  4. Paul said

    I think, I found my mistake. I was only adding successors for every hcn found. It has to be done for every number tested, I guess. This reduces speed considerably and much more memory is needed. It is surprising that the incorrect way of adding successors gives correct hcn’s up to number 815.

  5. Paul said

    I found a hack, that works: only add successors if the ratio between the number of divisors and the maximum number of divisors is higher than a cut-off (I tried 0.95). Then you can calculate the first 10000 hcn’s correctly in 668 secs, which is still great for such a simple scheme. If I do not use the hack and add successors only when finding an hcn the time is 21 secs, but the hcn’s are not correct. I will post code in Python later.

  6. Paul said

    I posted Python3 code on ideone. The results were checked upto 30000 hcns. Running times are 13 secs for 5000, 75 secs for 10000 and 1250 secs for 30000.The cut-off value is set to 0.99. That means that every number that has a number of divisors, that is at least 0.99 times the maximum number of divisors found sofar will be used to add successors.

  7. Paul said

    I did a few optimizations. The code runs new incredibly fast. It calculates the first 100,000 hcn’s in about 5 minutes. A dramatic speedup compared to the 15 minutes for 540 hcn’s that @matthew reported in his first post.

  8. Paul said

    I checked to results with the data of Achim Flammenkamp up to 770,000, which is about as high as his list gets. Still correct. Amazing.

  9. Ajit said

    Is there a C++ implementation?

  10. […] studied highly composite numbers in a series of several previous exercises, and had a lot of fun. In today’s exercise we look at a related concept: highly abundant […]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: