%
% Copyright (C) May/June 2005 by Henrik Pagenkopf
%
function P = product(A, B, zero, rec_depth) % approximates A* B via interval sectioning with "regula-falsi" for square matrices.
%
% - cond: A and B are square matrices and have same sizes > 0.
%
sa = size(A);
sb = size(B);
P = zeros(sa(1)); % assign output arguments in case of recursion.
if sa(1) >= 1 && sa(1)== sa(2) && sb(1)== sb(2) && sa(1)== sb(1) % => sa(2)== sb(2)
s = sa(1);
[rA,cA,dA]= max_coord_dist(A);
[rB,cB,dB]= max_coord_dist(B);
if dA <= zero | dB <= zero
P = zeros(s);
else
sum = 0;
for j = 1 :s
for k = 1 :s
sum = sum + A(j,k)+ B(j,k);
end;
end; % for
avg = sum/(2* s);
C = 0.5 *(A + B); % "average" matrix
[rC,cC,dC]= max_coord_dist(C);
S = avg * C; % =approx.=(C ^2) =approx.=(A* B) if A =approx.= B
S = correct_PR(C, C, S, rC, cC, rC, cC); % correcture on C^2
if rec_depth > 0 % make use of the distributive laws in the ring of square matrices of constant size s.
Delta = avg * eye(s) - C;
[rD,cD,dD]= max_coord_dist(Delta);
if dD <= min(1, abs(zero)^1.5) % exemplary constant for exponent!
% disp('Delta :=(avg* I)- C is almost zero.');
else % -- recursion: P = S + (C - avg* I)* C + (A - C)* C + A *(B - C);
S = S + product(C - avg * eye(s), C, zero, rec_depth - 1);
S = S + intersect_PR(real(A), real(B), zero, rec_depth); % -> the real part
S = S - intersect_PR(imag(A), imag(B), zero, rec_depth);
S = S + i* intersect_PR(real(A), imag(B), zero, rec_depth); % -> i*(the imag part)
S = S + i* intersect_PR(imag(A), real(B), zero, rec_depth);
P = correct_PR(A, B, S, rA, cA, rB, cB);
end; % if(dD)
end; % if(rec_depth)
end; % if(dA)
else % error
end
% return.
