function [w,D,cap] = rationalwts( sigma, tau, sin, C, Taylor )
%
% [w,D] = rationalwts( sigma, tau, s, <optional capacity C>, <optional Taylor=true> )
% Generalized  Barycentric weights : rational case
% (c)  Robert M. Corless, 2007
% and (optionally) Differentiation matrix
% in (optionally) Taylor form
%
% sigma is the Hermite data for the denominator
%
% on distinct nodes tau with integer confluencies s
% size(w) = [n,smax]
% [n,1] = size(tau(:))
% smax = max(s)

[n,dummy] = size(tau(:));

if max(size(sin))==1,
    s = sin*ones(n,1);
else
    s = sin;
end;
smax = max(s);

%  The following computes a skew-symmetric matrix of differences
%  between nodes. We need both sums and products of these.
taudiff = tau(:)*ones(1,n) - ones(n,1)*(tau(:).');

%  The following heuristic scaling seems to work well both for real
%  intervals and for the unit circle. Other sets need testing.  The purpose
%  of the scaling is to prevent underflow or overflow.
if nargin>=4 & C>0,
    capacity = C;
elseif n==1,
    capacity = 1;
else
    % max geometric mean of differences
    capacity = max(abs(prod(eye(n)+taudiff)).^(1/(n-1)));
end;
if nargout>2,
    cap = capacity;
end;
taudiff = taudiff/capacity;

taudiff = taudiff + eye(n);
iz = find( taudiff==0 );
if size(iz) ~= [0,0],
    error('Interpolation nodes must be explicitly distinct. Use confluency vector s to indicate repetition.');
end;

u  = zeros(n,smax);
nu = zeros(n,smax);

nu(:,1) = ones(n,1);

taudiffrecip = 1.0./taudiff;
taupow = ones(n,n);

for m=1:smax-1,
    taupow = taupow.*taudiffrecip;
    u(:,m) = sum( diag(s)*(taupow-eye(n)) ).';
    % Cauchy convolution for series recurrence relation.
    for i=1:n,
        nu(i,m+1) = u(i,1:m)*(nu(i,m:-1:1).') /m;
    end;
end;

% raise each col to s(j) power
for j=1:n,
    taudiffrecip(:,j) = taudiffrecip(:,j).^s(j);
end;

w = diag( prod( taudiffrecip, 2 ) )*nu;
% now reverse the order!
for i=1:n
    w(i,:)=[w(i,(s(i):-1:1)),zeros(1,smax-s(i))];
end;

ik = 1+cumsum(s(:)');
ik = [1,ik(1:n-1)]; % [1,s(1)+1,s(1)+s(2)+1,...,] points at y-vals
gam = zeros(n,smax);
for i=1:n,
    for j=1:s(i),
        gam(i,j) = dot( w(i,j:s(i)), sigma(ik(i):(ik(i)+s(i)-j)) );
    end;
    if abs(gam(i,s(i)))<=1.0e-6*norm(sigma,inf),
        warning('Potentially unattainable point tau(%d) = %12.4f',i,tau(i))
    end;
end;
w = gam;

if nargout > 1,
    %compute differentiation matrix in theta form if requested, H form
    %otherwise
    H = zeros(n,n,smax);
    for L=1:n,
        for i=[1:L-1,L+1:n],
            td = (tau(L)-tau(i))/capacity;
            for j=1:s(i),
                H(L,i,j) = 0;
                for k=s(i)-j:-1:0,
                    H(L,i,j) = w(i,j+k) + H(L,i,j)/td;
                end;
                if nargin>=5 & Taylor,
                    H(L,i,j) = s(L)*H(L,i,j)/w(L,s(L))/td;
                else
                    H(L,i,j) = factorial(s(L))*H(L,i,j)/w(L,s(L))/td/factorial(j-1);
                end;
            end;
        end;
        for mu=1:s(L),
            H(L,L,mu) = 0;
            for i=[1:L-1,L+1:n],
                td = (tau(i)-tau(L))/capacity;
                for j=1:min(mu,s(i)),
                    if nargin>=5 & Taylor,
                        H(L,L,mu) = H(L,L,mu) - td^(mu-j)*H(L,i,j)*prod(mu-1:-1:mu-j+1)/factorial(j-1);
                    else
                        H(L,L,mu) = H(L,L,mu) - td^(mu-j)*H(L,i,j)/factorial(mu-j);
                    end;
                end;
            end;
        end;
    end;
    
    dp1 = sum(s);
    D = zeros(dp1,dp1);
    irow = 0;
    for i=1:n,
        for j=2:s(i),
            irow = irow + 1;
            if nargin>=5 & Taylor,
                inc = j-1;
            else
                inc = 1;
            end;
            D(irow,irow+1) = inc;
        end;
        irow = irow + 1;
    end;
    
    irow = 0;
    for i=1:n,
        irow = irow + s(i);
        icol = 0;
        for j=1:n,
            for k=1:s(j),
                icol = icol + 1;
                D(irow,icol) = H(i,j,k);
            end;
        end;
    end;
end;