function [ L, D, U ] = CauchyDecomp( x, y )
%CAUCHYDECOMP Accurate LDU factoring [1/(x(i)+y(j)]
%   From Jim Demmel's ISSAC talk 2001
%   See also N. J. Higham, Accuracy and Stability
%   in Numerical Analysis, 2nd ed, SIAM, 2002, p. 514
%   which points to Cho, Math Comp 2(104) 819--825 1968.
%   Cho and Higham give explicit O(n) formulas
%   for entries of L, D, and U.
n = length(x);

U = 1.0./(diag(x)*ones(n)+ones(n)*diag(y));
L = eye(n);
D = eye(n);

for k=1:n,
    D(k,k) = U(k,k);
    L(k:n,k) = U(k:n,k)/D(k,k);
    U(k,k:n) = U(k,k:n)/D(k,k);
    U(k+1:n,k) = zeros(n-k,1);
    for i=k+1:n,
        for j=k+1:n,
            U(i,j) = U(i,j)*(x(i)-x(k))*(y(j)-y(k))/(x(k)+y(j))/(x(i)+y(k));
        end
    end
end
% Shockingly accurate.
