koetter_interp.m

在matlab中

function Qout = koetter_interp(L, alphas, M, k)

% cword = koetter_interp(L, alphas, M, k)
%
% Interpolates a bivariate polynomial with zeros at
% (alphas(i), j) of multiplicity M(j,i).  M is a multiplicity
% matrix of size q x n.
%
% The finite field over which the interpolation is done is GF(q), 
% where q is assumed to be a power of 2, e.g. q = 2^m
% The primitive polynomial is whatever MATLAB sets the default to be
% for that Galois field.
% 
% (q-1,k) are the parameters of the Reed-Solomon code.


q = size(M,1);
m = log2(q);     % field is GF(2^m)
n = size(M,2);

% initialize G matrix (xdegree, ydegree, L+1)
% so that G(j,1,j) = y^j
G = gf(zeros(L+1,L+1,L+1), m);
for ind1=0:L
   G(1, ind1+1, ind1+1) = 1;
end

% allocate discrepancy vector
Delta = gf(zeros(1,L+1),m);


% constraints are of the form Q_{r,s}(alpha_i, beta_i) = 0
% for 0 <= r+s <= m_i
% we have to impose them in (m_ij - 1, 1) lex order
[Is, Js, Ms] = find(M);
for ind1=1:length(Is)
  i = Is(ind1);
  j = Js(ind1);
  mij = Ms(ind1);
  alpha = alphas(j);
  beta = gf(i-1,m);

  % loop over the nonzero multiplicities and constraints
  for cons=1:mij
     for r=0:(mij-1)
        for s=0:(mij-1-r)

           % calculate the discrepancies for this constraint
           for ind2=0:L
              gj = reshape(G(:,:,ind2+1),size(G,1),size(G,2));
              Delta(ind2+1) = hasse_deriv2(gj,alpha,beta,r,s,m);
           end

           % find nonzero discrepancies
           Dfix = find(Delta~=0);

           % check if there are any discrepancies
           if (length(Dfix) ~= 0)
              % find the smallest j <= L such that there is a
              % discrepancy, call that jstar and the associated
              % polynomial gstar and discrepancy Dstar
 
              jstar = minwdegree(G,Dfix,1,k-1);
              gstar = reshape(G(:,:,jstar),size(G,1),size(G,2));
              Dstar = Delta(jstar);

              % loop over the guys with discrepancies
              for ind3=Dfix;
                % for the non-jstar entries
                if (ind3~=jstar)
                  gj = reshape(G(:,:,ind3),size(G,1),size(G,2));
                  tmp = Dstar*gj - Delta(ind3)*gstar;
                  G(:,:,ind3) = tmp;
                % for the jstar entry
                elseif (ind3==jstar)
                  xf = [zeros(1,size(gstar,2)); gstar];
                  nf = [gstar; zeros(1,size(gstar,2))];

                  gstar = nf;
                  
                  % replace that polynomial in jstar
                  gj = Dstar*xf - ...
		      hasse_deriv2(xf,alpha,beta,r,s,m)*nf;

                  Gnew = gf(zeros(size(G,1)+1,size(G,2),size(G,3)),m);
                  Gnew(1:size(G,1), 1:size(G,2), 1:size(G,3)) = G;
                  G = Gnew;

                  G(:,:,ind3) = gj;
                end
              end
           end
	end
     end
  end
end
       
qind = minwdegree(G,(0:L)+1,1,k-1);

Qout = reshape(G(:,:,qind),size(G,1),size(G,2));

[Is, Js] = find(Qout ~= 0);

Qout = Qout(1:max(Is), 1:max(Js));