mrdivide.m

几个关于多小波的程序

function [Q,method] = mrdivide(P1,P2)

% MRDIVIDE -- right division of matrix polynomials
%
%        Q = P1 / P2
%        Q = mrdivide(P1,P2)
%
% This function is not meant to be called by the user. It is called by
% Matlab if an expression of the form P1/P2 is encountered, where at
% least one of the two operands is a matrix polynomial. P2 must be
% square, either of a size compatible with P1, or a scalar.
%
% This routine knows how to handle 3 special cases at the moment,
% otherwise it gives up:
%
% 1. P2 has a non-singular first or last coefficient, in which case
%    we try long division; if it fails, we know that division is
%    not possible.
%
% 2. P2 is a projection factor, which has a known inverse
%
% 3. P2 is upper or lower triangular; this will succeed if the diagonal 
%    consists of nonzero monomials, but it may succeed even in other
%    cases.
%
% The METHOD output parameter is strictly for testing purposes. It 
% indicates which of the special cases was applied.

% Copyright (c) 2004 by Fritz Keinert ([e-mail]),
% Dept. of Mathematics, Iowa State University, Ames, IA 50011.
% This software may be freely used and distributed for non-commercial
% purposes, provided this copyright statement is preserved, and
% appropriate credit for its use is given.
%
% Last update: Feb 20, 2004

% make sure P1, P2 are both matrix polynomials

if ~isa(P1,'mpoly')
    P1 = mpoly(P1);
end
if ~isa(P2,'mpoly')
    P2 = mpoly(P2);
end
P1 = trim(P1);
P2 = trim(P2);
P1max = get(P1,'max');

% check that P2 is square and compatible with P1

[m1,m2] = size(P1);
[n1,n2] = size(P2);
if (n1 ~= n2 | (n1 ~= 1 & n1 ~= m2))
    error('incompatible sizes in division');
end

% Case 1: check if either the first or last coefficient of P2 is
% non-singular. If yes, attempt long division. 
%
% If the first coefficient is nonsingular, we divide left to right. 
% If the last coefficient is nonsingular, but not the first, we would
% have to do the division right to left. To save on programming
% effort, we instead reverse the order of the coefficients before and
% after the division.
%
% If both the first and last coefficients of P2 are nonsingular 
% (which includes the cases P2 = scalar and P2 = monomial), then 
% degree of Q = degree of P1 - degree of P2. 
%
% Otherwise, the degree of Q could be larger. The worst case I can
% think of that could be useful is when P2 is unit lower or upper
% triangular, in which case the degree of P2inv could be as large as
% (n-1)*(degree of P2), n = size of P2. If the division does not
% terminate at that point, we give up.

% Algorithm:
%
% 1. Pad P1 with P1extra zero matrices, where P1extra=0 if the first
% and last coefficient of P2 are both non-singular, and 
% P1extra=n*(length(P2)-1) otherwise. 
%
% The earliest possibility for the division to terminate
% successfully is when we reach the end of the original P1, which is
% after 
%
%        kmin = length(P1) - length(P2) + 1
%
% division steps. The last possibility we consider is when we reach
% the end of (P1 + P1extra zeros), which happens after
%
%        kmax = kmin + P1extra
%
% steps.
%
% We start with (kmin-1) division steps without checking the
% remainder, then we start checking.

method = 1;

first = ~issingular(P2.coef(:,:,1));
last  = ~issingular(P2.coef(:,:,end));

if (last & ~first)
    P1 = reverse(P1);
    P2 = reverse(P2);
end
    
if (~first & size(P2.coef,3) == 1)
% we are trying to divide by a single singular matrix
    error('right division failed');
end

if (first | last)
    P1length = size(P1.coef,3);
    P2length = size(P2.coef,3);
    if (first & last)
	if (P1length < P2length)
	    error('right division failed');
	end
	P1extra = 0;
    else
	P1extra = n1 * (P2length - 1);
    end
    Qlengthmin = P1length - P2length + 1;
    Qlength = Qlengthmin + P1extra;
% pad P1 with P1extra zeros
    P1 = trim(P1,[P1.min,P1max+P1extra]);
% do division
    for i = 1:Qlength
	Qcoef(:,:,i) = P1.coef(:,:,i) / P2.coef(:,:,1);
	for j = 2:P2length
	    P1.coef(:,:,i+j-1) = P1.coef(:,:,i+j-1) - Qcoef(:,:,i)*P2.coef(:,:,j);
	end
	if (i >= Qlengthmin) % check remainder
	    remainder = mpoly(P1.coef(:,:,i+1:i+P2length-1));
	    if iszero(remainder) % success!
		Q = mpoly(Qcoef,P1.min-P2.min);
		Q = trim(Q);
		[type,m,r] = match_type(P1,P2);
		Q = set(Q,'type',type,'m',m,'r',r);
		if (last & ~first)
		    Q = reverse(Q);
		end
		return
	    end
	end
    end
% If we get to this point, long division failed.
    error('right division failed');
end
    
% case 2: check if P2 is a generalized projection factor
% That means we try to factor P2 as
%
%                              s   k
%        P2 = P2  [ A + (I-A) z ] z
%               0
%       2
% with A  = A. If this succeeds, then
%
%          -1               -s    -1  -k
%        P2   = [A + (I-A) z  ] P2   z
%                                 0

method = 2;

P20 = moment(P2);
if ~issingular(P20)
    F = P2;
    for i = 1:size(F.coef,3)
	F.coef(:,:,i) = F.coef(:,:,i) / P20;
    end
    A = F.coef(:,:,1);
    if iszero(A^2 - A)
	if iszero(F.coef(:,:,end) - eye(m1) + A)
	    k = F.min;
	    s = F.min + size(F.coef,3) - 1 - k;
	    middle = mpoly(F.coef(:,:,2:end-1));
	    if iszero(middle) % success!
		P2inv = F;
		P2inv.coef(:,:,1) = eye(m1) - A;
		P2inv.coef(:,:,end) = A;
		P2inv.min = - k - s;
		for i = 1:size(P2inv.coef,3)
		    P2inv.coef(:,:,i) = P20 \ P2inv.coef(:,:,i);
		end
		Q = trim(P1 * P2inv);
		[type,m,r] = match_type(P1,P2);
		Q = set(Q,'type',type,'m',m,'r',r);
		return
	    end
	end
    end
end

% case 3: check if the matrix is upper or lower triangular.
% If it is, we attempt forward or back substitution.

% We don't check for division by zero, since that would lead
% to an error message, anyway.

method = 3;

if iszero(triu(P2,1)) % matrix is lower triangular
    Q = mpoly(zeros(n1,n2),0);
    if (isa(P1.coef,'sym') | isa(P2.coef,'sym'))
	Q = sym(Q);
    end
% Q(:,n2) = P1(:,n2) / P2(n2,n2);
    Q = subsasgn(Q,substruct('()',{':',n2}),...
		 subsref(P1,substruct('()',{':',n2})) / ...
	subsref(P2,substruct('()',{n2,n2}))); 
    for j = n2-1:-1:1
	% Q(:,j) = (P1(:,j) - Q(:,j+1:n2)*P2(j+1:n2,j)) / P2(j,j);
	Q = subsasgn(Q,substruct('()',{':',j}),...
	    (subsref(P1,substruct('()',{':',j})) - ...
	    subsref(Q,substruct('()',{':',j+1:n2})) * ...
	    subsref(P2,substruct('()',{j+1:n2,j}))) / ...
	    subsref(P2,substruct('()',{j,j})));
    end
    Q = trim(Q);
    [type,m,r] = match_type(P1,P2);
    Q = set(Q,'type',type,'m',m,'r',r);
 return
elseif iszero(tril(P2,-1)) % matrix is upper triangular
    Q = mpoly(zeros(n1,n2),0);
    % Q(:,1) = P1(:,1) / P2(1,1);
    Q = subsasgn(Q,substruct('()',{':',1}),...
	subsref(P1,substruct('()',{':',1})) / ...
	subsref(P2,substruct('()',{1,1}))); 
    for j = 2:n2
	% Q(:,j) = (P1(:,j) - Q(:,1:j-1)*P2(1:j-1,j)) / P2(j,j);
	Q = subsasgn(Q,substruct('()',{':',j}),...
	    (subsref(P1,substruct('()',{':',j})) - ...
	    subsref(Q,substruct('()',{':',1:j-1})) * ...
	    subsref(P2,substruct('()',{1:j-1,j}))) / ...
	    subsref(P2,substruct('()',{j,j})));
    end
    Q = trim(Q);
    [type,m,r] = match_type(P1,P2);
    Q = set(Q,'type',type,'m',m,'r',r);
    return
end

% if we get to this point, division has failed

error('right division failed');