Project Euler #12 in Julia slower than Python?

Question

My code in Julia, almost identical as the Python code (see below), runs in 4.6 s while the Python version runs in 2.4 s. Obviously there is a lot or room for improvement.

function Problem12()
    #=
     The sequence of triangle numbers is generated by adding the natural
    numbers. So the 7th triangle number would be:
    1 + 2 + 3 + 4 + 5 + 6 + 7 = 28.

    The first ten terms would be:
    1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

    Let us list the factors of the first seven triangle numbers:

     1: 1
     3: 1,3
     6: 1,2,3,6
    10: 1,2,5,10
    15: 1,3,5,15
    21: 1,3,7,21
    28: 1,2,4,7,14,28

    We can see that 28 is the first triangle number to have over five divisors.

    What is the value of the first triangle number to have over five hundred
    divisors?
    =#

    function num_divisors(n)
        res = floor(sqrt(n))
        divs = []
        for i in 1:res
            if n%i == 0
                append!(divs,i)
            end
        end
        if res^2 == n
            pop!(divs)
        end
        return 2*length(divs)
    end

    triangle = 0
    for i in Iterators.countfrom(1)
        triangle += i
        if num_divisors(triangle) > 500
            return string(triangle)
        end
    end
end

Python version below:

import itertools
from math import sqrt, floor


# Returns the number of integers in the range [1, n] that divide n.
def num_divisors(n):
    end = floor(sqrt(n))
    divs = []
    for i in range(1, end + 1):
        if n % i == 0:
            divs.append(i)
    if end**2 == n:
        divs.pop()
    return 2*len(divs)


def compute():
    triangle = 0
    for i in itertools.count(1):
        # This is the ith triangle number, i.e. num = 1 + 2 + ... + i =
        # = i*(i+1)/2
        triangle += i
        if num_divisors(triangle) > 500:
            return str(triangle)

Answer 1

The performance discrepancy in your Julia code comes from

res = floor(sqrt(n))  # is a Float64

which is a floating-point number (specifically, Float64) due to sqrt. This causes the loop variable i to be a Float64 as well, leading to a massive slowdown when you test for divisibility via taking the modulus (if n%i == 0) because Euclidean division (i.e., mod/rem and div/fld/cld) is very slow for floats. It is much faster for integers.

Furthermore, this does not always give the mathematically correct answer due to the inherently limited precision of floats. For example,

julia> Int(floor(sqrt(67108865^2 - 1)))
67108865

Instead, you should call isqrt to take the integer square root

res = isqrt(n)  # is an Int

which also gives the exact answer

julia> isqrt(67108865^2 - 1)
67108864

Avoid Int(floor(sqrt(n))) and floor(Int, sqrt(n)) due to the limited precision outlined above.

Answer 2

You're applying a brute force solution by counting all the divisors through trial division this is excessively slow and we can do better.

First, realize that all integers can be written as a product of prime numbers, this is by definition of prime numbers. In particular as a product of prime numbers to integer powers.

I.e: any integer \$n\$ can be written as \$n = \prod p_i^{k_i}\$ where \$p_i\$ is a prime number and \$k_i\$ is an integer.

Note further that a "triangle number" is just the arithmetic sum from 1 to n.

So we're looking for \$n\$ such that \$n(n+1)/2 = \prod p_i^{k_i}\$ under the condition that \$\prod k_i+1 > 500\$. Note here that product of powers plus one computes the number of ways you can pick zero (the +1) to all of the powers of a prime to contribute to a number, which is the number of divisors you can have.

Your program would then be something like this:

1. generate large list of primes using a prime sieve or hard code.
2. for each n, starting at 1 compute the "triangle number", do prime factorisation using the list of primes, this is our speedup we only need to try dividing by known primes. Then count the prime powers, compute their product to get the divisors and check if larger than 500.

This will be faster than OPs code because we only test division with known primes rather than all integers.

You might be able to do some fancy math to determine a better starting point than n=1 in step 2 by using the expressions above but that's left as an academic exercise for the reader.

Answer 3

You need the number of divisors, not divisors. So you don't need to gather all divisors, just remember the count. Increase it by 1 instead of push!/append, decrease instead of pop/