function n = sieve(k) decom = zeros(k^4, 2) for a = 1:k-1 for b = 1:k-a for c = 1:k-1 for d = 1:k-c A = [a,b;c,d] for e = 1:min( (k-b)/a, (k-d)/c ) for f = 1:min( (k-b)/a, (k-d)/c ) for g = 1:min ( (k-a*e)/b, (k-c*e)/d ) for h = 1:min ( (k-a*f)/b, (k-c*f)/d ) B = [e,f;g,h] C = A*B // Note that C is a matrix such that max(C) <= k r = (C(1,1)-1)*k^3 + (C(1,2)-1)*k^2 + (C(2,1)-1)*k + C(2,2)-1; decom(r,:) = [1, max(C)] // We mark that matrix as decomposable, with its maximal value end end end end end end end end n = zeros(k,1) for i = 1:size(decom, 1) // We look for the decomposable matrices computed if decom(i,1) == 1 then n(decom(i,2):k) = n(decom(i,2):k) + 1 end end n = [1:k]'.^4 - n // n(i) is the number of undecomposable matrix with max() <= i endfunction