For my CS class there was an extra credit problem in this week’s problem set where we had to write a sieve of Eratosthenes to find prime numbers up to an input integer. I had written the following for fun a few weeks before:
function [ primes, num ] = sieve( n ) % generate list of prime numbers less than or equal to n% start with list of all odd numbers from 3 to n primes = 3:2:n; % check numbers up to floor of sqrt(n) for i = 3:floor(sqrt(n)) % delete all numbers multiple of that primes(mod(primes, i) == 0 & i ~= primes) = ; end % add the 2 back on primes = [2 primes]; num = length(primes); end
So this is a pretty standard way of doing it (although it’s way slower than the built-in primes(), which probably uses a much cleverer method that I can’t figure out).
But we haven’t officially been taught iteration in lecture yet, so this problem had a rule of no using for or while loops. It took me a long time of staring at the command window but eventually I came up with this.
function [ primes ] = sieve( limit ) % first, make a list containing all numbers (potential primes) up to limit % don't need to even bother with the even numbers % but do remember to put the 2 back on at the end primes = 3:2:limit; % divide each number in primes by every other number (but not itself) and find where % that result is a whole number evenlyDivide = floor(primes ./ primes') == primes ./ primes' & floor(primes ./ primes') ~= 1; % anywhere there is a true value means that the number is not prime so % get rid of it primes(any(evenlyDivide)) = ; % remembered to add the 2 back on the end because it was taken off at % the beginning primes = [2 primes]; end
This works because right division (./) in MATLAB works between a 1xN row vector and a Mx1 column vector, where the result is a MxN matrix with every value of the row vector divided by every value of the column vector. This was a little surprising, because my professor said multiple times that you can only right divide vectors and matricies if they are the same size or if one of them is 1×1, which is clearly not the case here.
This method wins for number of lines, but doesn’t win when you run it with a large input like a million:
Clearly, it’s not the most optimal sieve known to man. Firstly, it doesn’t even take into consideration the sqrt(n) limit on numbers you have to check the primes list with, it just divides the entire list of potential primes by itself, which ends up making massive arrays.
By fiddling around with the solution sieve function they give, I discovered that it uses ones() somehow, but I can’t figure out an alternative method, and this homework is due tomorrow so I can’t really be bothered to try and figure out a method that doesn’t error with inputs greater than 100,000. Hopefully they won’t test any inputs that large. Someday, though, I would like to try and work out a faster and less resource-heavy method for doing this, because I like these kind of problems. (Project Euler is fun too! although it’s hard not to just take the easy way out and brute force some of the problems…)