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