Highly Composite Numbers

July 12, 2016

[ Today’s exercise was written by Zacharias Voulgaris, based on a Numberphile video. Guest authors are always welcome; contact me if you wish to write an exercise. ]

A highly composite number, also called an anti-prime, is a number n for which d(m) < d(n) for all m < n, where d(x) is the divisor function that gives a count of the number of divisors of x; in other words, a highly composite number has more divisors than any smaller number. Thus, a highly composite number, which has many divisors, is the opposite of a prime number, which has only two divisors. The sequence of highly composite numbers, which begins 1, 2, 4, 6, 12, 24, 36, … (A002182), has been studied by Ramanujan and Erdös, among others, and is a continuing object of study by number theorists. A famous highly composite number, known to Plato, is 5040.

Your task is to write a program that returns all highly composite numbers less than a given limit. 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.

Pages: 1 2

7 Responses to “Highly Composite Numbers”

  1. Jamie Hope said

    This version calculates the number of divisors as suggested by the first comment for [a href=”http://oeis.org/A000005″]A000005[/a]: it calculates the prime factorization of the number, and then calculates the product of (exponent+1). Under Kawa Scheme, at least, this ends up being about 35% faster than the brute force method. (The factorization could be further optimized, too.)

    (define (factor n)
      (let loop ((n n) (p 2) (e 0) (f '()))
        (cond ((> p n) (reverse f))
              ((= p n) (reverse (cons (cons p (+ e 1)) f)))
              ((= 0 (remainder n p))
               (loop (/ n p) p (+ e 1) f))
              ((= 0 e)
               (loop n (if (= p 2) 3 (+ p 2)) 0 f))
              (else
               (loop n (if (= p 2) 3 (+ p 2)) 0 (cons (cons p e) f))))))
    
    (define (nod n)
      (let loop ((d 1) (f (factor n)))
        (if (null? f)
            d
            (loop (* d (+ (cdar f) 1)) (cdr f)))))
    
    (define (hcp limit)
      (let loop ((i 1) (l -1) (ns '()))
        (if (> i limit)
            (reverse ns)
            (let ((pi (nod i)))
              (if (> pi l)
                  (loop (+ i 1) pi (cons (list i pi) ns))
                  (loop (+ i 1) l ns))))))
    
    (display (hcp 1000)) (newline)
    
  2. Rutger said
    highly_composite_list = []
    max_divisors = 0
    for i in range(1, 1000):
        divisors = 2 + sum([not i % d for d in range(2, i/2)])
        if divisors > max_divisors:
            max_divisors = divisors
            highly_composite_list.append(i)
    
    print highly_composite_list
    
  3. matthew said

    So: HCNs are a subset of numbers with non-ascending sequences of prime exponents, we can generate all such sequences using a couple of extension operations ([a,b,c]->[a,b,c,1] and [a,b,c] -> [a,b,c+1]) and use a heap to generate the corresponding numbers in order.

    Runtime is limited by the length of the primes array – with primes up to 119, we get 540 HCNs in about 15 mins (the 540th is 3625096965993854307674502488829776160975104640000 with 10899947520 divisors).

    import heapq
    
    primes = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,
              61,67,71,73,79,83,89,97,101,103,107,109,113,119]
    
    # Construct number and divisor count from sequence of
    # prime exponents.
    def mkentry(s):
        n = 1; d = 1
        for i in range(0,len(s)):
            n *= primes[i]**s[i]
            d *= s[i]+1
        return (n,d,s)
    
    # Starting from [], We can generate all non-ascending sequences
    # of positive integers with these two operations:
    
    def extend1(n,d,s):
        # [a,b,c,d] -> [a,b,c,d,1]
        s1 = s[:]
        s1.append(1)
        return mkentry(s1)
    
    def extend2(n,d,s):
        # [a,b,c,d] -> [a,b,c,d+1] (with d < c)
        if (len(s)) == 0: return None
        if (len(s) > 1 and s[-1] == s[-2]): return None
        s1 = s[:]
        s1[-1] += 1
        return mkentry(s1)
    
    # Generate all candidate numbers in order.
    def genbase():
        q = [(1,1,[])]
        while True:
            (n,d,s) = heapq.heappop(q)
            yield(n,d,s)
            heapq.heappush(q,extend1(n,d,s))
            e = extend2(n,d,s)
            if e: heapq.heappush(q,e)
    
    g = genbase(); max = 0
    while True:
        (n,d,s) = next(g)
        if d > max:
            print("%d %d %s"%(n,d,s))
            max = d
    
  4. Paul said

    @matthew. I like your solution. It is easy to speed it up by a factor of 2, by changing extend1 and extend2 a little to calculate the new n, d from the former ones. By the way, 119 is not a prime.

  5. matthew said

    @Paul: Thanks, much appreciated. Good point on 119, can’t think what I was thinking there (I don’t think it affects the result though).

    Reusing the n and d values is a good idea (also we can reuse some of the arrays, helping with memory allocation – and bytearrays are probably better too).

    There are some nice algorithms linked from https://en.wikipedia.org/wiki/Highly_composite_number that use serious maths to get much larger values.

  6. […] different solution to the previous exercise exploits the form of highly composite numbers, which always consists of small primes to large […]

  7. […] seen programs to compute the sequence of highly composite numbers in two previous exercises. Today, we look at a third algorithm, based on the Sieve of […]

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: