format long e
n=64; 
q=ones(1,n); e=65536*ones(1,n-1);
disp('Differential QD-Algorithm')
for p=1:3
  d=q(1);
  for k=1:n-1
    q(k)=d+e(k);
    e(k)=e(k)*(q(k+1)/q(k));
    d=d*(q(k+1)/q(k));
  end
  q(n)=d;
  p, min(q)
end

disp('Progressive QD-Algorithm')
q=ones(1,n); e=65536*ones(1,n-1);
for p=1:3
  e(n)=0;
  q(1)=q(1)+e(1);
  for k=2:n,
    e(k-1)=(e(k-1)/q(k-1))*q(k);
    q(k)=(q(k)-e(k-1)) + e(k);
  end
  p, min(q)
end

disp('Matlab eig')
q=ones(1,n); e=65536*ones(1,n-1);
R=diag(q)+diag(sqrt(e),1); B=R'*R; 
min(eig(B))