% MCCHEBYSHEV Multiple-component discretized Chebyshev algorithm.
%
%    This is an alternative to the routine chebyshev. It performs
%    a sequence of discretizations of the given weight function
%    (or measure), each discretization being used to compute
%    discretized modified moments, which in turn are passed
%    as input to the modified Chebyshev algorithmi in order
%    to produce approximations to the desired recurrence
%    coefficients. The fineness of the discretization is
%    characterized by a discretization parameter M. The support 
%    of the continuous part of the weight function is decomposed 
%    into a given number mc of subintervals (some or all of which 
%    may be identical). The routine then applies to each 
%    subinterval an M-point quadrature rule to discretize the
%    weight function on that subinterval. The discrete part of the
%    weight function (if there is any) is added on to the
%    discretized continuous weight function. The sequence of 
%    discretizations, if chosen judiciously, leads to convergence
%    of the recurrence coefficients for the discretized measures 
%    to those of the given measure. If convergence to within a
%    prescribed accuracy eps0 occurs before M reaches its maximum
%    allowed value Mmax, then the value of M that yields 
%    convergence is output as Mcap, and so is the number of 
%    iterations, kount. If there is no convergence, the routine 
%    displays the message "Mcap exceeds Mmax in mccheb" prior to
%    exiting.
%
%    The details of the discretization are to be specified prior 
%    to calling the procedure. They are embodied in the following 
%    global parameters:
%
%    mc     = the number of component intervals
%    mp     = the number of points in the discrete part of the 
%             measure (mp=0 if there is none)
%    iq     = a parameter to be set equal to 1, if the user
%             provides his or her own quadrature routine, and 
%             different from 1 otherwise
%    idelta = a parameter whose default value is 1, but is 
%             preferably set equal to 2, if iq=1 and the user 
%             provides Gauss-type quadrature routines
% 
%    The component intervals have to be specified (in the order 
%    left to right) by a global mcx2 array AB=[[a1 b1];[a2 b2];
%    ...;[amc bmc]], where for infinite extreme intervals a1=-Inf
%    resp. bmc=Inf. The discrete spectrum (if mp>0) is similarly
%    specified by a global mpx2 array DM=[[x1 y1];[x2 y2];...;
%    [xmp ymp]] containing the abscissae and jumps.
%
%    If the user provides his or her own quadrature routine 
%    "quadown", the routine mcdis must be called with the input 
%    parameter "quad" replaced by "@quadown", otherwise with 
%    "quad" replaced by "@quadgp", a general-purpose routine
%    provided in the package. The quadrature routine must have 
%    the form
%
%                             function xw=quad(M,i)
%
%    where M is the number of nodes and i identifies the interval
%    to which the routine is to be applied.
%
%    The routine mcchebyshev also applies to measures given
%    originally in multi-component form.
%
function [ab,Mcap,kount]=mcchebyshev(N,eps0,quad,Mmax)
global mc mp iq idelta DM uv AB 
global om2 abm
if N<1, error('Input variable N out of range'), end
iMcap=1; kount=-1;
ab(:,2)=zeros(N,1); b=ones(N,1);
Mcap=floor((2*N-1)/idelta);
while any(abs(ab(:,2)-b)>eps0*abs(ab(:,2)))
  b=ab(:,2);
  kount=kount+1;
  if kount>1, iMcap=2^(floor(kount/5))*N; end
  Mcap=Mcap+iMcap;
  if Mcap>Mmax, error('Mcap exceeds Mmax in mccheb'), end
  mtMcap=mc*Mcap;
  if iq~=1, uv=fejer(Mcap); end
  for i=1:mc
    im1tn=(i-1)*Mcap;
    xw=feval(quad,Mcap,i);
    xwm(im1tn+1:im1tn+Mcap,:)=xw;
  end
  if mp~=0, xwm(mtMcap+1:mtMcap+mp,:)=DM; end
  for k=1:2*N
    p1=zeros(mtMcap+mp,1);
    p=ones(mtMcap+mp,1);
    if k>1
      for l=1:k-1
        pm1=p1;
        p1=p;
        p=(xwm(:,1)-abm(l,1)).*p1-abm(l,2)*pm1;
      end
    end
    mom(k)=xwm(:,2)'*p;
  end
  ab=chebyshev(N,mom,abm);
end