% STURAN Symbolic Gauss-Turan quadrature formula.
%
function [sxw,it,condmax]=sturan(dig,n,s,eps0,sab0,hom)
digits(dig);
if hom==1, s0=1; else s0=s; end
sxw=vpa(zeros(n,2*s+2));
sn=(s+1)*n; sn2=(s+1)*n^2;
if size(sab0,1)<sn, error('array sab0 too short'), end
sxw0=sgauss(dig,sn,sab0(1:sn,:));
f=vpa(zeros(2*n-1,1)); fjac=vpa(zeros(2*n-1));
A=vpa(zeros(n+1,sn2)); B=vpa(zeros(n+1,sn2)); 
D=vpa(zeros(n,sn2)); P=vpa(zeros(n,sn2)); 
Q=vpa(zeros(n,sn2)); PI=vpa(zeros(n+1,sn)); 
PI(1,:)=vpa(ones(1,sn));
sab=sab0(1:n,:);
for sp=s0:s
  rho=[sab(:,1);sab(2:n,2)]; rhop=vpa(zeros(2*n-1,1));
%
% Newton iteration
%
  it=0; condmax=0;
  while any(subs(abs(rho-rhop))>eps0*subs(abs(rho)+1))
    rhop=rho;
    it=it+1;
    if it>50
      fprintf('condmax=%11.4e\n',condmax)
      error('Newton iteration does not converge') 
    end
    sabp=[rhop(1:n) [0;rhop(n+1:2*n-1)]];
    PI(2,:)=sxw0(:,1)'-sabp(1,1)*PI(1,:);
    for k=2:n
      PI(k+1,:)=(sxw0(:,1)'-sabp(k,1)).*PI(k,:)-sabp(k,2)*PI(k-1,:);
    end
%
% Computation of f
%
    for nu=1:n
      f(2*nu-1)=((sabp(nu,1)-sxw0(:,1)').*PI(nu,:).^2.*PI(n+1,:).^( ...
        2*sp))*sxw0(:,2);
      if nu<n
        f(2*nu)=((sabp(nu+1,2)*PI(nu,:).^2-PI(nu+1,:).^2).*PI( ...
          n+1,:).^(2*sp))*sxw0(:,2);
      end
    end  
%
% Computation of auxiliary matrices
%
    for mu=1:n
      m=mu:n:sn2;
      A(mu+1,m)=-PI(mu,:);
      for nu=mu+1:n
        A(nu+1,m)=(sxw0(:,1)'-sabp(nu,1)).*A(nu,m)-sabp(nu,2).*A( ...
          nu-1,m);
      end
      if mu>1
        B(mu+1,m)=-PI(mu-1,:);
        for nu=mu+1:n
          B(nu+1,m)=(sxw0(:,1)'-sabp(nu,1)).*B(nu,m)-sabp(nu,2)*B( ...
            nu-1,m);
        end
      end
    end
    for nu=1:n
      for mu=1:n
        m=mu:n:sn2;
        if nu==mu, D(nu,m)=vpa(1); end
        P(nu,m)=PI(nu,:).*(A(nu,m).*PI(n+1,:)+vpa(sp)*A(n+1,m).*PI( ...
          nu,:));
        Q(nu,m)=PI(nu,:).*(B(nu,m).*PI(n+1,:)+vpa(sp)*B(n+1,m).*PI( ...
          nu,:));
      end
    end
%
% Computation of the Jacobian of f
%
    for nu=1:n
      for mu=1:n
        m=mu:n:sn2;
        fjac(2*nu-1,mu)=vpa(2)*(PI(n+1,:).^(2*sp-1).*((sabp(nu,1) ...
          -sxw0(:,1)').*P(nu,m)+vpa(1/2)*D(nu,m).*PI(nu,:).^2.*PI( ...
          n+1,:))*sxw0(:,2));
        if mu>1
          fjac(2*nu-1,mu+n-1)=vpa(2)*(PI(n+1,:).^(2*sp-1).*(sabp(nu,1) ...
            -sxw0(:,1)').*Q(nu,m))*sxw0(:,2);
        end
        if nu<n
          fjac(2*nu,mu)=vpa(2)*(PI(n+1,:).^(2*sp-1).*(sabp(nu+1,2) ...
            *P(nu,m)-P(nu+1,m)))*sxw0(:,2);
          if mu>1
            fjac(2*nu,mu+n-1)=vpa(2)*(PI(n+1,:).^(2*sp-1).*((sabp(nu+1,2) ...
              *Q(nu,m)-Q(nu+1,m))+vpa(1/2)*D(nu+1,mu).*PI(nu,:).^2.*PI( ...
                n+1,:))*sxw0(:,2));
          end
        end
      end
    end
%
% Newton correction
%
    if cond(subs(fjac))>condmax, condmax=cond(subs(fjac)); end
    d=fjac\f; rho=rhop-d;
  end
  sab=[rho(1:n) [0;rho(n+1:2*n-1)]];
end
%it 
sxwg=sgauss(dig,n,sab);
sxw(:,1)=sxwg(:,1);
%
% Computation of the weights
%
Ahat=vpa(zeros(2*s+1));
dtau=vpa(zeros(n,1)); u=vpa(zeros(1,2*s));
for nu=1:n
  prod=vpa(ones(sn,1));
  for mu=1:n
    if mu~=nu
      prod=(sxw0(:,1)-sxw(mu,1)).*prod/(sxw(nu,1)-sxw(mu,1));
    end
  end
  prod=prod.^(2*s+1);
  for sig=0:2*s
    muhat(sig+1)=((sxw0(:,1)'-sxw(nu,1)).^sig.*prod')*sxw0(:,2);
  end
  dtau=sxw(:,1)-sxw(nu,1);
  for l=1:2*s
    if n==1
      u(l)=vpa(0);
    else
      u(l)=sum(vpa(1)./dtau(find(subs(dtau))).^l);
    end
  end
  for k=1:2*s+1
    Ahat(k,k)=vpa(1);
  end
  for j=1:2*s
    for k=1:2*s+1-j
      Ahat(k,k+j)=-vpa((2*s+1))*(u(1:j)*Ahat(1:j,j))/vpa(j);
    end
  end
  sxw(nu,2:2*s+2)=(Ahat\muhat')';
  fac=vpa(1);
  for sig=0:2*s
    sxw(nu,sig+2)=sxw(nu,sig+2)/fac;
    fac=vpa((sig+1))*fac;
  end
end   
