Why is my script so slow? Eratosthenes

I was trying to implement Eratosthenes sieve without "cheating" and looking it up online. What i found was that my script works, i.e. returns the correct numbers, BUT it was sooo slow, took like 25-30 minutes to complete for N=2000000.

This other smaller script I found online was much faster -- timed .6 seconds.

Can someone maybe explain what it is that eats up so much time in script compared to the smaller one?

%%fast script 
N = 2000000;
M = 1:N;
p = M(2);
while p^2 < N
    M(p^2:p:N) = 0;
    p = find(M>p,1);
end

primes = M(M>1);

%%my slow script
clear
E = 2000000;
M = (2:E);

i = 1;
p = M(i);
K = [];
while p^2 < E
    K = [];
    if p ~= 0
        for k = p^2-1:p:length(M)
            if mod(M(k),p) == 0
                K = [K,k];
            end
        end
        M(K) = 0;
    end

      i = i + 1;
      if i <= E-1
          p = M(i);
      else
          break;
      end
  end
  M = M(M>0);

[Adam]

doc profile

will tell you immediately where the time is being spent.

[Martin Hellkvist]

WOW thanks!

I did this and found out that -- surprisingly enough -- the crook in the drama was the vector appending

K = [K,k];

The profiling told me that 99.9% of the time was spent here so I changed it out to

K(end+1) = k;

Which was about 6500 times faster...

Today I learned

[Adam]

Yeah, there are numerous cases like this where certain alternatives are very slow.

I'm having to do rather a lot of profiling at the moment and it certainly does open your eyes to aspects of the programming you may never have even considered relevant previously!

[Philip Borghesani]

That script is not optimal ether:

 >> tic; p=primes(2000000);toc
 Elapsed time is 0.062879 seconds.

Answer 1

[Jan]

Even K(end+1) is not efficient. The iterative growing of an array is a severe mistake. See this example:

K = [];
for i = 1:1e6
  K = [K, i];
  % or K(end+1) = i;  % Same problem!
end

In the first iteration a new array of size 1 is allocated and the contents of i is copied to the last value. In the 2nd iteration a new array of size 2 is allocated, the former values of K are copied and i is inserted as last element, etc. At the end, you did not reserve and copy 1e6 elements (8MB, because a double has 8 bytes), but sum(1:1e6)*8 Bytes: 4e12 or 4TB!

Never let an array grow iteratively. Always allocate the array at once instead. In your case it is not needed to use the vecot K at all:

while p^2 < E
  if p ~= 0
      for k = p^2-1:p:length(M)
          if M(k) && (mod(M(k),p) == 0)
              M(k) = 0;
          end
      end
  end
  ...

Set M(k) directly to 0 instead of collecting the vector of indices at first.

[Martin Hellkvist]

Ok, that's very insightful. Even though I have studied computational science this did not occur to me. I will try with pre allocating, and see what happens.

[Jan]

You have studied computational science and the topic of iterative growing of arrays has not been mentioned? Wow. Then I suggest to read this: https://en.wikipedia.org/wiki/Anti-pattern

Answer 2

[ES]

To the uninformed eyes, there seems to be an unnecessary for loop in the slower code that runs E*E times.

[Martin Hellkvist]

see answer to comment^^