Highly Composite Numbers

July 12, 2016

We give an answer in Julia:

function nod(x::Int64) # number of divisors
	n = 2 # every number is divisible by 1 and by itself
	
	for z = 2:div(x,2)
		if x % z == 0
			n += 1
		end
	end
	
	return n
end

function hcp(N::Int64 = 1000) # Highly Composite Numbers
	y = [1] # list of highly composite numbers
	z = [1] # list of corresponding number of divisors

	for x = 2:N
		n = nod(x)
		
		if n > z[end]
			push!(y, x)
			push!(z, n)
		end
	end

	return hcat(y, z)
end

The calculation of the number of divisors is brute force, checking each z from 1 to n − 1 for divisibility. Function hcp iterates from 2 to the limit and stores each number whose divisor count exceeds the current maximum.

Since ideone does not provide Julia, a corresponding Scheme program is available at http://ideone.com/3EjR09.

Posted by programmingpraxis

Filed in Exercises

8 Comments »

8 Responses to “Highly Composite Numbers”

Jamie Hope said

July 12, 2016 at 7:07 PM

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)

Rutger said

July 13, 2016 at 12:18 PM

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

matthew said

July 14, 2016 at 6:23 AM

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

Paul said
July 15, 2016 at 7:41 AM
@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.
matthew said
July 15, 2016 at 8:21 AM
@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.
Highly Composite Numbers, Using Priority Queues | Programming Praxis said
July 15, 2016 at 9:01 AM
[…] different solution to the previous exercise exploits the form of highly composite numbers, which always consists of small primes to large […]
Highly Composite Numbers, A Sieving Approach | Programming Praxis said
July 19, 2016 at 9:01 AM
[…] 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 […]
Highly Abundant Numbers | Programming Praxis said
December 20, 2016 at 9:01 AM
[…] 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 […]

Programming Praxis