
function rule = cub_interval_fejer_gegenbauer(n,lambda,ftype)

%--------------------------------------------------------------------------
% Object
%--------------------------------------------------------------------------
% Fejer I, FejerII, Clenshaw-Curtis quadrature rule on (-1,1) with 
% cardinality "n".
%
% The rule is w.r.t. Gegenbauer weight function.
%
%                  w(x)=(1-x^2)^(lambda-0.5),    lambda > -0.5.
%
% Here "n" is the cardinality of the interpolatory rule described by
% the variable "ftype".
%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
% n: cardinality of the cubature rule
% ftype: 1: Fejer I, 2: Fejer II, 3: Clenshaw-Curtis.
%--------------------------------------------------------------------------
% Output:
%--------------------------------------------------------------------------
% rule: n x 2 matrix of nodes "rule(:,1)" and weights "rule(:,2)".
%--------------------------------------------------------------------------
% Routines:
%--------------------------------------------------------------------------
% 1. fejerI
% 2. fejerII
% 3. clenshaw-curtis
% 4. moms_cheb1type_gegenbauer
% 
% as well as all their subroutines. All the necessary files are attached to
% this function.
%--------------------------------------------------------------------------
% Reference:
%--------------------------------------------------------------------------
% J. Waldvogel
% "Fast Construction of Fejer and Clenshaw-Curtis Quadrature rules".
% BIT, 2003, Vol. 43, No. 1, pp. 001–018
%
% Fast Construction of Fejer and Clenshaw-Curtis rules for general weight 
% functions,
% Computers Mathematics with Applications, Volume 65, Issue 4, 
% February 2013, Pages 682-693.
%--------------------------------------------------------------------------
% License:
%--------------------------------------------------------------------------
% Copyright (C) 2026- Alvise Sommariva.
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
%
% Authors:
%
% Alvise Sommariva <alvise@math.unipd.it>
%
% Date: January 11, 2026
%--------------------------------------------------------------------------

if nargin < 1, n=10; end
if nargin < 2, lambda=[]; end
if isempty(lambda), lambda=0.5; end
if nargin < 3, ftype=3; end

% .............................. Nodes/weights  ...........................

moms=moms_cheb1type_gegenbauer(n-1,lambda);

switch ftype

    case 1  % Fejer1
        [x,w]=fejerI(n,moms);

    case 2 % Fejer2
        [x,w]=fejerII(n,moms);

    case 3 % Clenshaw-Curtis
        [x,w]=clenshaw_curtis(n,moms);

end

rule=[x w];









%--------------------------------------------------------------------------
% Subroutines.
%--------------------------------------------------------------------------

function [x,w]=fejerI(n,moms)

% OBJECT: FEJER RULE. TYPE I.

% INPUTS:
% n: NUMBER OF POINTS. TESTED WITH n <= 14 MILLION NODES.
% moms: WEIGHT FUNCTION MODIFIED MOMENTS (w.r.t. CHEB. POLY. TYPE I).

% OUTPUTS:
% x: QUADRATURE NODES IN THE INTERVAL (-1,1). COLUMN VECTOR.
% w: QUADRATURE WEIGHTS. COLUMN VECTOR.

% REFERENCE:
% ALVISE SOMMARIVA
% "Fast Construction of Fej?er and Clenshaw-Curtis rules for general weight
% functions".

% AUTHOR:
% ALVISE SOMMARIVA, JULY 23, 2012.

if nargin < 2
    moms=moms_cheb1type_legendre(n-1);
end

moms(1)=sqrt(1/2)*moms(1);
x=cos((2*(1:n)'-1)*pi/(2*n));
w=sqrt(2/n)*idct(moms);









function [x,w]=fejerII(n,moms)

% OBJECT: FEJER RULE. TYPE II.

% INPUTS:
% n: NUMBER OF POINTS. TESTED WITH n <= 14 MILLION NODES.
% moms: WEIGHT FUNCTION MODIFIED MOMENTS (w.r.t. CHEB. POLY. TYPE I).

% OUTPUTS:
% x: QUADRATURE NODES IN THE INTERVAL (-1,1). COLUMN VECTOR.
% w: QUADRATURE WEIGHTS. COLUMN VECTOR.

% REFERENCE:
% ALVISE SOMMARIVA
% "Fast Construction of Fej?er and Clenshaw-Curtis rules for general weight
% functions".

% AUTHOR:
% ALVISE SOMMARIVA, JULY 23, 2012.

if nargin < 2
    moms=moms_cheb1type_legendre(n-1);
end

momsII=compute_moments_IIw(n-1,moms);
theta=(1:n)*pi/(n+1); theta=theta'; x=cos(theta);
w=(2*sin(theta)/(n+1)).*dst(momsII);









function [x,w]=clenshaw_curtis(n,moms)

% OBJECT: CLENSHAW-CURTIS RULE.

% INPUTS:
% n: NUMBER OF POINTS. TESTED WITH n <= 14 MILLION NODES.
% moms: WEIGHT FUNCTION MODIFIED MOMENTS (w.r.t. CHEB. POLY. TYPE I).

% OUTPUTS:
% x: QUADRATURE NODES IN THE INTERVAL [-1,1]. COLUMN VECTOR.
% w: QUADRATURE WEIGHTS. COLUMN VECTOR.

% REFERENCE:
% ALVISE SOMMARIVA
% "Fast Construction of Fej?er and Clenshaw-Curtis rules for general weight
% functions".

% AUTHOR:
% ALVISE SOMMARIVA, JULY 23, 2012.

if nargin < 2
    moms=moms_cheb1type_legendre(n-1);
end

momsIIcc=compute_moments_IIwcc(n-3,moms);
theta=(1:n-2)'*pi/(n-1); xx=cos(theta);
w=((2*sin(theta)/(n-1)).*dst(momsIIcc))./(1-xx.^2);
w1=(2*sum(moms)-moms(1)-moms(end))/(2*(n-1));
wn=moms(1)-w1-sum(w);
x=[1;xx;-1]; w=[w1;w;wn];









function [momsII,momsI]=compute_moments_IIwcc(n,momsI)

% COMPUTATION OF MOMENTS W.R.T. CHEBYSHEV SECOND TYPE FROM CHEBYSHEV FIRST
% TYPE.

momsI=momsI(1:n+3);
momsII=0.5*(momsI(1:end-2)-momsI(3:end));









function momsII=compute_moments_IIw(n,momsI)

% COMPUTATION OF MOMENTS W.R.T. CHEBYSHEV SECOND TYPE FROM CHEBYSHEV FIRST
% TYPE.

momsI=momsI(1:n+1);

momsIe=momsI(1:2:end); momsIo=momsI(2:2:end);
momsII=zeros(size(momsI));
momsII(1:2:end)=2*cumsum(momsIe)-momsI(1);
momsII(2:2:end)=2*cumsum(momsIo);









function moms=moms_cheb1type_legendre(n)

moms=zeros(n+1,0);
m=0:2:n; m=m';
moms_even_degree=2./(1-m.^2);
moms(m+1,1)=moms_even_degree;









function moms=moms_cheb1type_gegenbauer(n,lambda)

moms=zeros(n+1,1);

% CHOOSE LAMBDA > -0.5. WEIGHT: (1-x^2)^(LAMBDA-0.5).
moms(1)=gamma(lambda+0.5)*sqrt(pi)/gamma(lambda+1);

nn=floor(n/2);
term1=gamma(lambda+0.5)*sqrt(pi)/gamma(lambda+1);
j=1:nn; j=j';
v=(j-1-lambda)./(j+lambda);
G=cumprod(v);
moms(3:2:end)=term1*G;
