3
\$\begingroup\$

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)
\$\endgroup\$
4
  • \$\begingroup\$ What are compiler options for julia? Are you using interactive mode? \$\endgroup\$ Jul 22 '21 at 8:16
  • \$\begingroup\$ Not sure about that. I'm just running vanilla Julia in VSCode \$\endgroup\$ Jul 22 '21 at 8:32
  • 1
    \$\begingroup\$ Check it. Also how do you measure time? \$\endgroup\$ Jul 22 '21 at 9:40
  • \$\begingroup\$ The big performance killer is due to res being a floating-point number. Replace res = floor(sqrt(n)) by res = isqrt(n) and you should get an order of magnitude improvement in speed. I'll write up an answer that includes this. \$\endgroup\$
    – Vincent Yu
    Jul 23 '21 at 1:34
5
\$\begingroup\$

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.

\$\endgroup\$
2
  • \$\begingroup\$ The biggest performance issue is not isqrt vs sqrt, it's using trial division instead of a better factorisation method. \$\endgroup\$
    – Emily L.
    Jul 24 '21 at 12:42
  • 1
    \$\begingroup\$ @EmilyL. Thanks, I've edited my answer to clarify what it answers. Note that the OP's question is why their translated Julia code is slower than their Python code. The answer as explained is that the inner loop in the Julia code is taking floating-point modulus (via libc's fmod) whereas the Python code is taking integer modulus, ultimately because floor(sqrt(n)) is a float in Julia. The solution is to use isqrt(n). \$\endgroup\$
    – Vincent Yu
    Jul 24 '21 at 13:56
1
\$\begingroup\$

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.

\$\endgroup\$
0
\$\begingroup\$

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/

\$\endgroup\$
1
  • \$\begingroup\$ Thx for your answer! Although your suggestions do help in terms of allocations, they don't help in terms of performance :( \$\endgroup\$ Jul 22 '21 at 8:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.