% SCHEBYSHEV_SOB Symbolic counterpart of chebyshev_sob.m
%
function [B,normsq]=schebyshev_sob(dig,n,mom,abm)
syms sig mu x y sgn
digits(dig);
if n<=0, error('n out of range'), end
if n>size(mom,2)/2, n=size(mom,2)/2; end
if nargin<3, abm=vpa(zeros(2*n-1,2)); end
if n>(size(abm,1)+1)/2, n=(size(abm,1)+1)/2; end
smom=vpa(mom,dig); sabm=vpa(abm,dig);
B=vpa(zeros(n)); sig=vpa(zeros(n,2*n));
% 
% Initialization
%
sig(1,:)=smom(1,:); B(1,1)=sabm(1,1)+sig(1,2)/sig(1,1);
normsq(1)=sig(1,1);
if n==1, return, end
mu(1,1,1,:)=smom(2,:);
mu(1,2,1,1)=vpa(0); mu(2,1,1,1)=vpa(0); mu(2,1,1,2)=vpa(0); 
mu(1,2,1,2)=smom(2,1);
tau=vpa(zeros(3,2*n-1)); tau(3,1)=vpa(1);
for k=2:2*n-2
  for l=1:k-1
    tau(1,l)=tau(2,l); tau(2,l)=tau(3,l);
  end
  s=vpa(0);
  tau(3,1)=tau(2,2)*sabm(2,2)+tau(2,1)*(sabm(1,1)-sabm(k,1))...
    -tau(1,1)*sabm(k,2);
  s=s+tau(3,1)*smom(2,1);
  for j=2:k
    tau(3,j)=tau(2,j-1)+tau(2,j+1)*sabm(j+1,2)+tau(2,j)*...
      (sabm(j,1)-sabm(k,1))-tau(1,j)*sabm(k,2);
    if j==k, tau(3,j)=tau(3,j)+1; end
    s=s+tau(3,j)*smom(2,j);
  end
  mu(1,2,1,k+1)=s; mu(2,1,1,k+1)=vpa(0);
end
%
% Continuation
%
for k=2:n
  for l=k:2*n-k+1
    x=B(1:k-1,k-1); y=sig(k-1:-1:1,l);
    s=sum(x.*y);
    sig(k,l)=sig(k-1,l+1)+sabm(l,1)*sig(k-1,l)...
      +sabm(l,2)*sig(k-1,l-1)+mu(1,2,k-1,l)...
      -mu(2,1,k-1,l)-s;
    if(k==l), normsq(k)=sig(k,l); end
  end
  for j=1:k
    B(k-j+1,k)=(sig(j,k+1)+sabm(k,1)*sig(j,k))/sig(j,j);
    if j<k
      B(k-j+1,k)=B(k-j+1,k)+sabm(k,2)*sig(j,k-1)/sig(j,j);
      s=vpa(0);
      for l=j:k-1
        s=s+B(l-j+1,l)*sig(l,k)/sig(l,l);
      end
      B(k-j+1,k)=B(k-j+1,k)-s;
    end
    if j>1
      B(k-j+1,k)=B(k-j+1,k)-sig(j-1,k)/sig(j-1,j-1);
    end
  end
  if k<n
    for l=k:2*n-k
      for i2=1:2
        for i1=1:2
          if i1*i2<4
            s=0;
            for j=1:k-1
              s=s+B(j,k-1)*mu(i1,i2,k-j,l);
            end
            mu(i1,i2,k,l)=mu(i1,i2,k-1,l+1)+sabm(l,1)*...
              mu(i1,i2,k-1,l)+sabm(l,2)*mu(i1,i2,k-1,l-1)-s;
            if i1~=i2
              if i2==1
                sgn=vpa(1);
              else
                sgn=vpa(-1);
              end
              mu(i1,i2,k,l)=mu(i1,i2,k,l)+sgn*mu(1,1,k-1,l);
            end
          end
        end
      end
    end
  end
end

