clear;
n = 5; t_its = 30;
%First we test if x^t [A + B] x = x^T A x + x^T B x for all x, A, B.
max_e = 0;
for i=1:t_its
  x = rand(n,1);
  A = rand(n,n);
  B = randn(n,n);
  max_e = max(max_e,abs(x'*A*x + x'*B*x - x'*(A+B)*x));
end%for <- with for for Octave
fprintf('Maximum error in MatSum-test = %g. Strike a key to continue...\n',max_e)
pause
max_e = 0;
for i=1:t_its
  x = rand(n,1);
  v = rand(n,1);
  A = v*v';
  max_e = max(max_e,abs((x'*v)^2 - x'*A*x));
end%for
fprintf('Maximum error in Rewrite-test = %g. Strike a key to continue...\n',max_e)
pause
%Third test concerns the definition and positive definiteness (positivity
%of the smallest eigenvalue) of C.
N = 13;
max_e = 0; min_l = 1e14;
for i = 1:t_its
  x = rand(N,n); c = cov(x); v = x - repmat_own(mean(x),N,1);
  cv = v'*v/(N-1);
  max_e = max(max_e,norm(c - cv));
  [u,lambda] = eig(c);
  min_l = min(min_l,lambda(1,1)); %Eigenvalues always in ascending order
end%for
fprintf('Cov-test: max. difference in definition = %g, smallest eigenvalue = %g.\n\n',max_e,min_l)

