%
% CONDVL  Frobenius condition number and gradient
%    of a Vandermonde-like matrix involving orthonormal
%    polynomials relative to a measure defined by
%    the array ab of recurrence coefficients (for
%    the monic orthogonal polynomials)
%
function [c,g]=condVl(n,x,ab)
%function c=condVl(n,x,ab)
s1=0;
for nu=1:n
  p=ab(1,2); sp=1/p;
  p0=0; p1=1;
  for k=1:n-1
    p=ab(k+1,2)*p;
    pm1=p0; p0=p1;
    p1=(x(nu)-ab(k,1))*p0-ab(k,2)*pm1;
    sp=sp+p1*p1/p;
  end
  s1=s1+sp;
end
xw=gauss(n,ab);
s2=0;
for nu=1:n
  sg=0;
  for k=1:n
    mu=(1:n)';
    mmn=find(mu-nu);
    p=prod((xw(k,1)-x(mmn))./(x(nu)-x(mmn)));
    sg=sg+xw(k,2)*p*p;
  end
  s2=s2+sg;
end
c=sqrt(s1*s2);
%return
rho=sqrt(s2/s1);
g=zeros(n,1);
for i=1:n
  p=ab(1,2); sp=0;
  p0=0; p1=1; dp0=0; dp1=0;
  for k=1:n-1
    p=ab(k+1,2)*p;
    pm1=p0; p0=p1;
    dpm1=dp0; dp0=dp1;
    p1=(x(i)-ab(k,1))*p0-ab(k,2)*pm1;
    dp1=p0+(x(i)-ab(k,1))*dp0-ab(k,2)*dpm1;
    sp=sp+p1*dp1/p;
  end
  s=0;
  for nu=1:n
    sg=0;
    for k=1:n
      mu=(1:n)';
      mmn=find(mu-nu);
      p1=prod((xw(k,1)-x(mmn))./(x(nu)-x(mmn)));
      if nu~=i
        mu=(1:n)';
        mmni=find(mu-nu&mu-i);
        p2=prod((xw(k,1)-x(mmni))./(x(nu)-x(mmni)));
        p2=(xw(k,1)-x(nu))*p2/((x(nu)-x(i))^2);
      else
        l=(1:n)';
        lmn=find(l-nu);
        ss=sum(1./(x(nu)-x(lmn)));
        p2=-p1*ss;
      end
      sg=sg+xw(k,2)*p1*p2;
    end
    s=s+sg;
  end
  g(i)=rho*sp+s/rho;
end
gnorm=norm(g,2);
