function [g,MG] = rec(n,nup,x)
%%  REC Recurrence relation for orthonormal Gram polynomials
%
%   [G] = REC(N,NUP,X) returns the orthonormal Gram polynomial of degree N
%   and parameter NUP evaluated at vector X.
%
%   [G,MG] = REC(N,NUP,X) also returns the (N+1)-by-size(X) matrix MG with
%   the Gram polynomials of degrees 0,1,...,N evaluated on X. Its first row
%   contains the polynomial of degree 0, the second row contains the
%   polynomial of degree 1 etc. Its first column are the polynomials on
%   X_1, the second column are the polynomials on X_2 etc.
%
%
%   Copyright 2020 by Dimitar K. Dimitrov and Lourenco L. Peixoto
%   All rights reserved.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   REC implemented by Lourenco L. Peixoto, 2020 - see [1-3].
%
%   Basic references:
%   
%       [1] I. Area, D. K. Dimitrov, E. Godoy, and V. Paschoa, "Approximate
%   calculation of sums I: Bounds for the zeros of Gram polynomials", SIAM
%   J. Numer. Anal., 2014.
%
%       [2] I. Area, D. K. Dimitrov, E. Godoy, and V. Paschoa, "Approximate
%   calculation of sums II: Gaussian type quadrature", SIAM J. Numer.
%   Anal., 2016.
%   
%       [3] D. K. Dimitrov and L. L. Peixoto, "An efficient algorithm for
%   the classical least squares approximation", submitted, 2020.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

comp = length(x);

if n==0
    g = ones(1,comp);
    if (nargout == 2)
        MG = g;
    end
elseif n==1
    cte = sqrt(3*nup^2/(nup^2-1));
    g = cte*x;
    if (nargout == 2)
        MG = [ones(1,comp); g];
    end
else % n>=2
    cte = sqrt(3*nup^2/(nup^2-1));
    g1 = ones(1,comp); g = cte*x;

    k=1:n;
    alpha = sqrt(nup^2*(k.^2-.25)./(k.^2.*(nup^2-k.^2)));

    if (nargout == 1)
        for k=2:n
            g2 = g1; g1 = g; % initialise for polynomial
            g = 2*alpha(k)*x.*g1 - alpha(k)/alpha(k-1).*g2; % polynomial
        end
    else % stores all the polynomials of degree 0,...,n:
        MG = ones(n+1,comp); % preallocating
        MG(2,:) = g; % polynomial of degree 1
        for k=2:n
            g2 = g1; g1 = g; % initialise for polynomial
            g = 2*alpha(k)*x.*g1 - alpha(k)/alpha(k-1).*g2; % polynomial of degree k
            MG(k+1,:) = g; % polynomial of degree k
        end
    end
end

end