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 2019 by Dimitar K. Dimitrov and Lourenco L. Peixoto
%   All rights reserved.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   GRAM 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Check inputs:
if ( m > floor(2.5*sqrt(nup)) )
    warning('gram: M should be less or equal to floor(2.5*sqrt(NUP)).')
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

[x] = transpose(forster_petras(m)); % initial approximations
dt = inf; 
count = 0;
maxit = 50;
denom = zeros(1,m);
while (norm(dt,inf) > eps && count <= maxit)
    [g] = rec(m,nup,x);  % evaluate the polynomial via recurrence formula
    for k = 1:m
        d = 1;  
        for j = 1:m
            if (k~=j)
                d = d*(x(k)-x(j)); % denominator
            end
        end
        denom(k)=d;
    end
    dt = g./(coef(m,nup)*denom);
    x = x - dt;
    count = count + 1;
    % organize the zeros in ascending order. (It is unnecessary if N>>m).
    x = sort(x); 
end

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


 
%% Compute the weights
[~,MG] = rec(m-1,nup,x);
i=1:m;
w = 1./(sum(MG(:,i).^2.));

x = transpose(x);
    
end
 


%% Compute the leading coefficient of the polynomial
function cl = coef(j,nup)
% compute the leading coefficient of the polynomial of degree J  
k=1:j;
alpha = sqrt(nup^2*(k.^2-.25)./(k.^2.*(nup^2-k.^2))); %alpha_0,...,alpha_{J-1}
cl = 2^j*prod(alpha);
end
 
%% Compute the initial approximations for the zeros of Gram polynomial 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;
    else
        s = 0;
    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