function [x,w] = gram(m,nup)
%%  GRAM    Nodes and weights of Gaussian type quadrature for sums.
%   
%   [X] = GRAM(M, NUP) returns the zeros X of the orthonormal Gram
%   polynomial of degree M with parameter NUP.
%
%   [X, W] = GRAM(M, NUP) returns both X and a row vector W of
%   weights for the Gaussian type quadrature.
%     
%   Calls REC.
%
%   Copyright 2021 by Dimitar K. Dimitrov and Lourenco L. Peixoto
%   All rights reserved.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   GRAM implemented by Lourenco L. Peixoto, 2021 - 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", SIAM J. Sci. Comput., 2020.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Check inputs:
if ( m > floor(2.5*sqrt(nup)) )
    warning('gram: M = %i. M must be less than floor(2.5*sqrt(NUP)).', m)
end

% Deal with trivial cases:
if (m<=3)
    if (m<1)
        error('gram:inputs:m', 'First input should be a positive integer.');
    end
    if (m==1)
        x = 0;
        w = 1;
        return
    elseif (m==2) 
        x = sqrt((nup^2-1)/(3*nup^2));
        x = [-x x]';
        w = [.5 .5];
        return
    else % m=3
        x = sqrt((3*nup^2-7)/(5*nup^2));
        x = [-x 0 x]';
        w = 4/3*(nup^2-4)/(3*nup^2-7);
        w = [.5*(1-w) w .5*(1-w)];
        return
    end
end
    
 
 
%% Compute the zeros with WDDK method

[x1] = transpose(forster_petras(m)); % initial approximations

m_odd = (mod(m,2));
if m_odd % m is odd number
    x1(floor(m/2)+1) = 0;
end

dt = inf; 
count = 0;
maxit = 50;
denom = zeros(1,m);

% alpha_0,...,alpha_{M-1}:
k=1:m;
alpha = sqrt(nup^2*(k.^2-.25)./(k.^2.*(nup^2-k.^2))); 

% compute the leading coefficient of the polynomial of degree M:
coef_lead = 2^m*prod(alpha);

% initialise:
x = x1; 
lim = 1-1/nup; % interval limit

while (norm(dt,inf) > eps && count <= maxit)
    g = rec(m,nup,x,alpha);  % evaluate the polynomial via recurrence formula
    
    % WDDK:
    for k = 1:m
        d = coef_lead;
        for j = 1:m
            if (k~=j)
                d = d*(x(k)-x(j)); % denominator
                if d==0 % underflow:
                    error('Underflow in WDDK.')
                end
            end
        end
        denom(k) = d;
    end
    dt = g./denom;
    x = x - dt;
    
    if (x(1) < -lim || x(end) > lim)
        error('Nodes are out of range (-1+1/NUP, 1-1/NUP). The value M = %i is too large.', m)
    end
    
    if (m_odd) % m is an odd number
        x(floor(m/2)+1) = 0;
    end
    
    count = count + 1; 
end

if (count == maxit+1 && norm(dt,inf) > eps)
     warning('gram: %i iterations in WDDK method. ', maxit)
end



%% Compute the weights

% only non negative nodes:
z = x(floor(m/2)+1:m);

cte = sqrt(3*nup^2/(nup^2-1));
g1 = ones(1,length(z)); g = cte*z;

w = g1.^2 + g.^2;
for k=2:m-1
    % initialise for polynomial:
    g2 = g1; g1 = g; 
    % polynomial of degree k:
    g = 2*alpha(k)*z.*g1 - alpha(k)/alpha(k-1).*g2; 
    w = w + g.^2;
end
w = 1./w;

% Reflect for symmetric values:
if mod(m,2) 
    w = [flip(w), w(2:end)];
else
    w = [flip(w), w];
end

x = transpose(x);

end







%% Compute initial approximations for the zeros of Gram pol. of degree M
function [x] = forster_petras(m)
% X = FORSTER_PETRAS(M) returns the approximations in Forster and Petras
% (1993, Theorem 1) for the zeros of Legendre polynomial.
if m==1
    x=0;
    s=1;
else
    if ( mod(m,2) )
        s = 1; % odd m
    else
        s = 0; % even m
    end
    k = ((m+s)/2:-1:1).';
    
    theta = (k-.25)/(m+.5)*pi;
    x0 = cos(theta);
    cos2 = x0.*x0;
    % Sharp initial guess (Forster and Petras, 1993):
    x = cos( theta + .25/(2*(m+.5)^2)*(1-(6+.25*(9-2*cos2))./(12*(m+.5)^2*...
    (1-cos2))).*cot(theta) );
end
x = [-x(end:-1:1+s) ; x]; % reflect negative zeros

end