gfadd.m

Proakis《contemporarycommunication systems using matlab》matlab源代码

function c = gfadd(a, b, p, len)
%GFADD  Add two GF(P) polynomials or two GF(P^M) elements.
%       C = GFADD(A, B) adds two GF(2) polynomials A and B. The resulting
%       GF(2) polynomial C will keep the larger length between A and B.
%
%       C = GFADD(A, B, P) adds two polynomials A and B in GF(P) when P is
%       a scalar prime number.
%       When P is a matrix that contains the tuple of all elements in 
%       GF(Q^M), this function takes A and B as indices (power number of the
%       exponential form) of GF(Q^M) elements, the output C is 
%       alpha^C = alpha^A + alpha^B in GF(Q^M). The computation is 
%       element-by-element computation. You can generate the tuple of all
%       elements in GF(Q^M) by P = GFTUPLE([-1:Q^M-2]', M, Q).  In power
%       representation form, [-Inf, 0, 1, 2, ...] represents 
%       [0, 1, alpha, alpha^2, ...] in GF(p^m).
%
%       C = GFADD(A, B, P, LEN) adds two GF(P) polynomials A and B. The
%       resulting GF(P) polynomial C keeps the given length LEN. When LEN is a
%       negative number, the length of C equals degree(C) + 1. P must a scalar
%       integer in using this format.
%
%       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, GFTUPLE, GFPLUS.

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

% 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(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(p(a + 2, :) + p(b + 2, :), prim);
            cal_len = len_a;
        elseif  (len_a == 1)
            tmp = rem(ones(len_b,1)*p(a + 2, :) + p(b + 2, :), prim);
            cal_len = len_b;
        elseif  (len_b == 1)
            tmp = rem(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(p(a(1:cal_len) + 2, :) + p(b(1:cal_len) + 2, :), prim);
        end;
        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 GFADD---