m = 4;
n = 3;
p = 5;
k = 500000;
A = rand(m,n);
B = rand(n*k, p);

% Loop solution
tic
result1 = zeros(m*k, p);
for i = 1 : k
    result1(i*m-(m-1):i*m,:) = A*B(i*n-(n-1):i*n,:);
end
fprintf('loop time: %.4f\n', toc);

tic;
[m,n]=size(A);
[r,c]=size(B);
BB = permute(reshape(B.', p, n, []), [2 1 3]);
%[a,b,c]=size(BB);
BBB = reshape(BB,[n k*p]);
A_times_B_slices=A*BBB;
C = permute(reshape(A_times_B_slices,m,p,[]),[2 1 3]);
result2 = reshape(C, p, m*k).';
fprintf('matrix manipulation time: %.4f\n', toc);

tic;
result3 = block_product(A, B);
toc;

%result2 = C(:,:).'; %<<<-------------------------***had to change this line***

equality = isequal(result1,result3)