function [fek,leb] = symmfek(a,n)

% by Marco Vianello, University of Padova
% October 1 2013

% computes n+1 approximate Fekete points for polynomial interpolation 
% of degree n on the union of two symmetric disjoint intervals, 
% [-1,-a] U [a,1], 0<a<1
% and estimates their Lebesgue constant (that appears to be much lower 
% for odd than for even n)

% extraction mesh and corresponding rectangular Vandermonde-like matrix 
m=ceil(n^2/(1-a));
x=linspace(a,1,m);
y=-x;
pts=[y,x]; 
% recurrence coeffs of the monic Chebyshev polynomials
abm(:,1)=zeros(2*n+1,1);
abm(:,2)=0.25*ones(2*n+1,1); abm(1,2)=pi; abm(2,2)=0.5;
% modified moments of the monic Chebyshev polynomials
mom=zeros(1,2*n+2);
j=(0:2:2*n+1);
mom(1:2:end)=(j./(j.^2-1)).*(cos((j+1).*acos(1))-cos((j+1).*acos(a)));
mom(1:2:end)=mom(1:2:end)+(1./(j-1)).*(a*cos(j.*acos(a))-cos(j.*acos(1)));
mom(1:2:end)=2*mom(1:2:end);
mom(3:2:end)=mom(3:2:end)./2.^(j(2:end)-1);
% recurrence coefficients of orthogonal polynomials on [-1,-a] U [a,1]
ab = chebyshev(n+1,mom,abm);
% computing the Vandermonde matrix at pts in the orthonormal basis
B=eval_op(n,ab,pts');
xw=gauss(n+1,ab);
D=diag(sqrt(xw(:,2)));
W=eval_op(n,ab,xw(:,1));
DD=diag(1./sqrt(diag((D*W)'*(D*W))));
B=B*DD;

% preconditioning the polynomial basis
[Q,R]=qr(B,0);
B=B/R;
[Q,R1]=qr(B,0);
B=B/R1;
% [new basis] = [old basis]*inv(R)*inv(R1)

% computing the approximate Fekete points
mm=ones(n+1,1); 
w=B'\mm;
ind=find(abs(w)>0);
fek=pts(ind);

% roughly estimating their Lebesgue constant
V=eval_op(n,ab,fek');
V=V*DD/R/R1;
leb=norm(V'\B',1);



function V=eval_op(deg,ab,x)

%--------------------------------------------------------------------------
% OBJECT.
%----------
%
% deg: degree.
% ab: orthogonal polynomials three terms recurrence.
% x: points in which evaluate the polynomial of degree <= deg.
%    P_{n+1}(x)=(x-a_{n})*P_n(x)-b_{n}*P_{n-1}(x)
%
% V: Vandermonde matrix evaluated at x of degree <= deg.

L=size(x,1);

aa=ab(:,1);
bb=ab(:,2);

V=[zeros(L,1) ones(L,1)];


for ii=1:deg
   Vloc=(x-aa(ii)).*V(:,ii+1)-bb(ii)*V(:,ii);
   V=[V Vloc];
end

V=V(:,2:end);



% CHEBYSHEV Modified Chebyshev algorithm.
%
%    Given a weight function w encoded by its first 2n modified
%    moments, stored in the (row) vector mom, relative to monic 
%    polynomials defined by the (2n-1)x2 array abm of their
%    recurrence coefficients, [ab,normsq]=CHEBYSHEV(n,mom,abm)
%    generates the array ab of the first n recurrence coefficients
%    of the orthogonal polynomials for the weight function w, and 
%    the vector normsq of their squared norms. The n alpha-
%    coefficients are stored in the first column, the n beta-
%    coefficients in the second column, of the nx2 array ab. The
%    call [ab,normsq]=CHEBYSHEV(n,mom) does the same, but using the 
%    classical Chebyshev algorithm. If n is larger than the sizes
%    of mom and abm warrant, then n is reduced accordingly.
%
function [ab,normsq]=chebyshev(N,mom,abm)
if N<=0, error('N out of range'), end
if N>size(mom,2)/2, N=size(mom,2)/2; end
if nargin<3, abm=zeros(2*N-1,2); end
if N>(size(abm,1)+1)/2, N=(size(abm,1)+1)/2; end
ab(1,1)=abm(1,1)+mom(2)/mom(1); ab(1,2)=mom(1);
if N==1, normsq(1)=mom(1); return, end
sig(1,1:2*N)=0; sig(2,:)=mom(1:2*N);
for n=3:N+1
  for m=n-1:2*N-n+2
    sig(n,m)=sig(n-1,m+1)-(ab(n-2,1)-abm(m,1))*sig(n-1,m) ...
      -ab(n-2,2)*sig(n-2,m)+abm(m,2)*sig(n-1,m-1);
  end
  ab(n-1,1)=abm(n-1,1)+sig(n,n)/sig(n,n-1)-sig(n-1,n-1)/ ...
    sig(n-1,n-2);
  ab(n-1,2)=sig(n,n-1)/sig(n-1,n-2);
end
for n=1:N, normsq(n)=sig(n+1,n); end; normsq=normsq'; 


% GAUSS Gauss quadrature rule.
%
%    Given a weight function w encoded by the nx2 array ab of the 
%    first n recurrence coefficients for the associated orthogonal
%    polynomials, the first column of ab containing the n alpha-
%    coefficients and the second column the n beta-coefficients, 
%    the call xw=GAUSS(n,ab) generates the nodes and weights xw of
%    the n-point Gauss quadrature rule for the weight function w.
%    The nodes, in increasing order, are stored in the first 
%    column, the n corresponding weights in the second column, of
%    the nx2 array xw.
%
function xw=gauss(N,ab)
N0=size(ab,1); if N0<N, error('input array ab too short'), end
J=zeros(N);
for n=1:N, J(n,n)=ab(n,1); end
for n=2:N
  J(n,n-1)=sqrt(ab(n,2));
  J(n-1,n)=J(n,n-1);
end
[V,D]=eig(J);
[D,I]=sort(diag(D));
V=V(:,I);
xw=[D ab(1,2)*V(1,:)'.^2];


