Parallellising the factorisation of triangle numbers

Question

The sequence of triangle numbers is generated by adding the natural numbers. So the 7^th 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. The factors contained in the first four triangle numbers are:

1: 1
3: 1, 3
6: 1, 2, 3, 6
10: 1, 2, 5, 10

We see that 6 is the first triangle number to have over four divisors.

In order to find the first triangle number with over 500 divisors, I wrote the following code:

s1=0;
s2=0;
A=zeros(1,2);
for i=1:4000
    s1=i+s1;
    s2=0;
    for n=1:s1
        if rem(s1,n)==0
            s2=s2+1;
        end
    end
    if s2>=500
        A=[s1 s2];
        disp(A)
    end
end

I would like to speed up this code using parfor instead of for on the outer loop. However, I get the error:

Error: The variable s1 is perhaps intended as a reduction variable, but is actually an uninitialized temporary.

See Parallel for Loops in MATLAB, "Temporary Variables Intended as Reduction Variables".

How can I parallellise this code?

Answer 1

If performance is an issue you can use this code:

ii = 1;
while 1
    tn = ii*(ii+1)/2;                  %compute the next triangular number
    fa = factor(tn);                   %prime factorization
    nb = prod(histc(fa,unique(fa))+1); %number of divisor
    if nb > 500
       fprintf('there is %d divisors for %d\n',nb,tn)
       break
    end
    ii = ii+1;
end

And you will found that 76576500 is the answer.

The prime factorization factor is way faster than finding all the divisors with the divisors function.

So we can compute factor(n) that will output:

x1*n1, x2*n2, x3*n3... 

for example factor(60) give:

2 2 3 5

or

2x2,1x3,1x5

• There is 3 ways to use the prime factor 2: 2^0, 2^1 or 2^2.
• There is 2 ways to use the prime factor 3: 3^0, 3^1
• There is 2 ways to use the prime factor 5: 5^0, 5^1

So we know that 60 have 3*2*2 = 12 divisors, because the prime factor can be combined in 12 different ways.

Answer 2

You cannot parallelise this, as it is recursive: s1 = i + s1.

parfor executes your loop in random order, thus it is not known what value s1 will have in iteration n, as its value explicitly depends on its value in iteration n-1.

Instead of doing a for loop here, I'd go for a while loop with the condition s2<500. Then as soon as you go over the desired threshold, your loop will stop, thus you limit the number of necessary iterations to the bare minimum.

---

Note that parfor is not a magic wand. The overhead of launching the parallel pool will probably be larger than the actual computation time of this short loop. Read this question and this one for background information on how parfor works and how to write suitable code for it.