How do you implement the divisor function in code?

Question

Overall Problem: Project Euler 12 - What is the value of the first triangle number to have over five hundred divisors?

Focus of problem: The divisor function

Language: Python

Description: The function I used is brute and the time it take for the program to find a number with more divisors than x increases almost exponentially with each 10 or 20 numbers highers. I need to get to 500 or more divisors. I've identified that the divisor function is what is hogging down the program. The research I did lead me to divisor functions and specifically the divisor function which is supposed to be a function that will count all the divisors of any integer. Every page I've looked at seems to be directed toward mathematics majors and I only have high-school maths. Although I did come across some page that mentioned allot about primes and the Sieve of Atkins but I could not make the connection between primes and finding all the divisors of any integer nor find anything on the net about it.

Main Question: Could someone explain how to code the divisor function or even provide a sample? Maths concepts make more sense to me when I look at them with code. So much appreciated.

brute force divisor function:

def countdiv(a):
    count = 0
    for i in range(1,(a/2)+1):
        if a % i == 0:
            count += 1
    return count + 1    # +1 to account for number itself as a divisor

Answer 1

If you need a bruteforce function to calculate Number of Divisors (also known as tau(n))

Here's what it looks like

def tau(n):
        sqroot,t = int(n**0.5),0
        for factor in range(1,sqroot+1):
                if n % factor == 0:
                        t += 2 # both factor and N/factor
        if sqroot*sqroot == n: t = t - 1 # if sqroot is a factor then we counted it twice, so subtract 1
        return t

The second method involves a decomposing n into its prime factors (and its exponents).

tau(n) = (e1+1)(e2+1)....(em+1) where n = p1^e1 * p2^e2 .... pm^em and p1,p2..pm are primes

More info here

The third method and much more simpler to understand is simply using a Sieve to calculate tau.

def sieve(N):
        t = [0]*(N+1)
        for factor in range(1,N+1):
                for multiple in range(factor,N+1,factor):
                        t[multiple]+=1
        return t[1:]

Here's it in action at ideone

Answer 2

I agree with the two other answers submitted here in that you will only need to search up to the square root of the number. I have one thing to add to this however. The solutions offered will get you the correct answer in a reasonable amount of time. But when the problems start getting tougher, you will need an even more powerful function.

Take a look at Euler's Totient function. Though it only indirectly applies here, it is incredibly useful in later problems. Another related concept is that of Prime Factorization.

A quick way to improve your algorithm is to find the prime factorization of the number. In the Wikipedia article, they use 36 as an example, whose prime factorization is 2^2 * 3^2. Therefore, knowing this, you can use combinatorics to find the number of factors of 36. With this, you will not actually be computing each factor, plus you'd only have to check divisors 2 and 3 before you're complete.

Answer 3

When searching for divisors of n you never have to search beyond the square root of the number n. Whenever you find a divisor that's less than sqrt(n) there is exactly one matching divisor which is greater than the root, so you can increment your count by 2 (if you find divisor d of n then n/d will be the counterpart).

Watch out for square numbers, though. :) The root will be a divisor that doesn't count twice, of course.

Answer 4

If you're going to solve the Project Euler problems you need some functions that deal with prime numbers and integer factorization. Here is my modest library, which provides primes(n), is_prime(n) and factors(n); the focus is on simplicity, clarity and brevity at the expense of speed, though these functions should be sufficient for Project Euler:

def primes(n):
    """
        list of primes not exceeding n in ascending
        order; assumes n is an integer greater than
        1; uses Sieve of Eratosthenes
    """
    m = (n-1) // 2
    b = [True] * m
    i, p, ps = 0, 3, [2]
    while p*p < n:
        if b[i]:
            ps.append(p)
            j = 2*i*i + 6*i + 3
            while j < m:
                b[j] = False
                j = j + 2*i + 3
        i += 1; p += 2
    while i < m:
        if b[i]:
            ps.append(p)
        i += 1; p += 2
    return ps

def is_prime(n):
    """
        False if n is provably composite, else
        True if n is probably prime; assumes n
        is an integer greater than 1; uses
        Miller-Rabin test on prime bases < 100
    """
    ps = [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]
    def is_spsp(n, a):
        d, s = n-1, 0
        while d%2 == 0:
            d /= 2; s += 1
        if pow(a,d,n) == 1:
            return True
        for r in xrange(s):
            if pow(a, d*pow(2,r), n) == n-1:
                return True
        return False
    if n in ps: return True
    for p in ps:
        if not is_spsp(n,p):
            return False
    return True

def factors(n):
    """
        list of prime factors of n in ascending
        order; assumes n is an integer, may be
        positive, zero or negative; uses Pollard's
        rho algorithm with Floyd's cycle finder
    """
    def gcd(a,b):
        while b: a, b = b, a%b
        return abs(a)
    def facts(n,c,fs):
        f = lambda(x): (x*x+c) % n
        if is_prime(n): return fs+[n]
        t, h, d = 2, 2, 1
        while d == 1:
            t = f(t); h = f(f(h))
            d = gcd(t-h, n)
        if d == n:
            return facts(n, c+1, fs)
        if is_prime(d):
            return facts(n//d, c+1, fs+[d])
        return facts(n, c+1, fs)
    if -1 <= n <= 1: return [n]
    if n < -1: return [-1] + factors(-n)
    fs = []
    while n%2 == 0:
        n = n//2; fs = fs+[2]
    if n == 1: return fs
    return sorted(facts(n,1,fs))

Once you know how to factor a number, it is easy to count the number of divisors. Consider 76576500 = 2^2 * 3^2 * 5^3 * 7^1 * 11^1 * 13^1 * 17^1. Ignore the bases and look at the exponents, which are 2, 2, 3, 1, 1, 1, and 1. Add 1 to each exponent, giving 3, 3, 4, 2, 2, 2, and 2. Now multiply that list to get the number of divisors of the original number 76576500: 3 * 3 * 4 * 2 * 2 * 2 * 2 = 576. Here's the function:

def numdiv(n):
    fs = factors(n)
    f = fs.pop(0); d = 1; x = 2
    while fs:
        if f == fs[0]:
            x += 1
        else:
            d *= x; x = 2
        f = fs.pop(0)
    return d * x

You can see these functions at work at http://codepad.org/4j8qp60u, and learn more about how they work at my blog. I'll leave it to you to work out the solution to Problem 12.