function test_strass()
n = 16;
A = magic(n);
for i = 1:n
    for j = 1:n
        B(i,j) = cos((i-1)*(j-1)*pi/(n^2));
    end
end
C1 = A*B;
C2 = strass(A,B,4);
ndiff = norm(C2-C1)

function C = strass(A,B,nmin)
% A "divide-and-conquer" recursive algorithm for matrix multiplication
% as n grows, scales like n^2.807 compared to n^3 for direct multiplication
% This version works correctly so long as n = nmin*2^q
[n,k] = size(A);
if n <= nmin
    C = A*B;
else
    m = n/2; u = 1:m; v = m+1:n;
    P1 = strass(A(u,u)+A(v,v),B(u,u)+B(v,v),nmin);
    P2 = strass(A(v,u)+A(v,v),B(u,u),nmin);
    P3 = strass(A(u,u),B(u,v)-B(v,v),nmin);
    P4 = strass(A(v,v),B(v,u)-B(u,u),nmin);
    P5 = strass(A(u,u)+A(u,v),B(v,v),nmin);
    P6 = strass(A(v,u)-A(u,u),B(u,u)+B(u,v),nmin);
    P7 = strass(A(u,v)-A(v,v),B(v,u)+B(v,v),nmin);
    C(u,u) = P1 + P4 - P5 + P7;
    C(u,v) = P3 + P5;
    C(v,u) = P2 + P4;
    C(v,v) = P1 + P3 - P2 + P6;
end
