gfsub.m

数字通信第四版原书的例程

function c = gfsub(a, b, p, len)
%GFSUB  Subtract computation of two GF(P) polynomials or GF(P^M) elements.
%       C = GFSUB(A, B) subtracts GF(2) polynomial B from A. The resulting
%       GF(2) polynomial C keeps the larger length between A and B. In this
%       case, the result is the same as using GFADD.
%
%       C = GFSUB(A, B, P) subtracts polynomial B from A in GF(P) for a
%       scalar prime number P. When P is a matrix that lists all elements of
%       a GF(Q^M) field, this function takes A and B as GF(Q^M) exponential
%       form index. The output C is alpha^C = alpha^A - alpha^B in GF(Q^M).
%       The output C is a column vector. A list of all element in GF(Q^M)
%       can be generated by P = GFTUPLE([-1:Q^M-2]',M,Q), 
%
%       C = BFSUB(A, B, P, LEN) subtracts Galois field polynomials B from A.
%       The resulting GF(P) polynomial C keeps the given length LEN. When 
%       LEN is a negative number, the length of C will be degree(C) + 1.
%
%       In polynomial computation, A, B, and C are in ascending order, i.e.,
%       A = [a_0, a_1, a_2,..., a_(n-1), a_n] represents
%       A(X) = a_0 + a_1 X + a_2 X^2 +...+ a_(n-1) X^(n-1) + a_n X^n
%       a_i must be a element in GF(P).
%
%       In power representation form, [-Inf, 0, 1, 2, ...] represents 
%       [0, 1, alpha, alpha^2, ...] in GF(p^m).
%
%       See also GFMUL, GFDIV, GFADD, GFTUPLE

%        Wes Wang 6/8/94, 7/18/94, 10/7/95.
%        Copyright (c) 1995-96 by The MathWorks, Inc.
%        $Revision: 1.1 $  $Date: 1996/04/01 18:00:07 $

% input validation
a=a(:)'; b=b(:)';

if nargin < 2
       error('Not enough input for GFADD')
elseif nargin < 3
    p = 2;
end

[n_p, m_p] = size(p);

if n_p <= 1
    %the case of polynomial computation
    if ~isempty([find(a~=floor(a)), find(a<0), find(a>=p)])
     error('The polynomial coeficients must be in GF(P)')
    elseif ~isempty([find(b~=floor(b)), find(b<0), find(b>=p)])
     error('The polynomial coeficients must be in GF(P)')
    end;

    % match the length of vectors.
    na = length(a);
    nb = length(b);
    if na > nb
      b = [b zeros(1, na-nb)];
    elseif nb > na
      a = [a zeros(1, nb-na)];
    end;
    c = rem(p + a - b, p);

    % truncate the resul as requested.
    if (nargin > 3)
        nc = max(na, nb);
        if isempty(len)
             c = gftrunc(c);
        elseif len < 0
           c = gftrunc(c);
        elseif len <= nc
            c = c(1:len);
        else
            c = [c zeros(1, len-nc)];
        end;
    end;
else
    % GF(Q^M) field calculation.
    % The first row of P is -Inf, second is 1 and so on.

    % handle vector case
    indx = find(a > n_p - 2);
    if ~isempty(indx)
        a(indx) = rem(a(indx), n_p - 1);
    end;
    indx = find(b > n_p - 2);
    if ~isempty(indx)
        b(indx) = rem(b(indx), n_p - 1);
    end;
    if (a < 0)
        c = b;
    elseif (b < 0)
        c = a;
    else
        prim = max(max(p)) + 1;
        indx = find(~(a >= 0));
        a(indx) = - ones(1, length(indx));
        indx = find(~(b >= 0));
        b(indx) = - ones(1, length(indx));
        len_a = length(a);
        len_b = length(b);
        if (len_a == len_b) 
            tmp = rem(prim + p(a + 2, :) - p(b + 2, :), prim);
            cal_len = len_a;
        elseif  (len_a == 1)
            tmp = rem(prim + ones(len_b,1) * p(a + 2, :) - p(b + 2, :), prim);
            cal_len = len_b;
        elseif  (len_b == 1)
            tmp = rem(prim + p(a + 2, :) - ones(len_a, 1)*p(b + 2, :), prim);
            cal_len = len_a;
        else
            cal_len = min(len_a, len_b);
            tmp = rem(prim + p(a(1:cal_len) + 2, :) - p(b(1:cal_len) + 2, :), prim);
        end;

        % find the match
        for i = 1: cal_len
            indx = find(sum(rem(p + prim - ones(n_p, 1) * tmp(i,:), prim)') == 0);
            if length(indx) == 1
                c(i) = indx - 2;
            else
                if isempty(indx)
                    error('The list of Galois field is not a complete one');
                else
                    error('The list of Galois field has repeat elements');
                end;
            end;
        end;
    end;
    indx = find(c < 0);
    if ~isempty(indx)
        c(indx) = indx - Inf;
    end;
    c = c(:)';
    if cal_len < len_a
        c = [c a(cal_len + 1 : len_a)];
    elseif cal_len < len_b
        c = [c b(cal_len + 1 : len_b)];
    end;
end;
%--end of GFSUB---