#
# BHIP: Barycentric Hermite Interpolation Program
#
# (c) Robert M. Corless, December 2007
#
# Compute the barycentric form of the unique Hermite interpolant 
# of the polynomial given by values and derivative values of 
# p(t) at the nodes tau.
#
# CALLING SEQUENCES
#
#   ( p, gam ) := BHIP( flist, tau, t );
#   ( p, gam ) := BHIP( ftayl, tau, t, 'Taylor' = true, 'Denominator' = q );
#
# Input:
#    flist :: list of lists of derivative values
#             [ [f[1,0],f[1,1],...,f[1,s[1]-1]], 
#               [f[2,0],f[2,1],...,f[2,s[2]-1]],
#                ... 
#               [f[n,0],f[n,1],...,f[n,s[n]-1]] ]  
#    or, equivalently, when 'Taylor' is true,
#    ftayl :: list of lists of Taylor Coefficients, 
#             same as above but divided by factorials
#             [ [f[1,0]/0!,f[1,1]/1!,...,f[1,s[1]-1]/(s[1]-1)!],
#               [f[2,0]/0!,f[2,1]/1!,...,f[2,s[2]-1]/(s[2]-1)!],
#                ...
#               [f[n,0]/0!,f[n,1]/1!,...,f[n,s[n]-1]/(s[n]-1)!] ] 
#      
#    As a convenience, singleton lists may be entered without []    
#    If plist[i] = [] is empty, then nothing
#    is known about f at t = tau[i], and 
#    tau[i] will not figure in the output.
#    Sometimes this is useful.  Not all flist[i] may
#    be empty.  
#           
#    tau :: list of distinct complex nodes
#           [tau[1],tau[2],...,tau[n]]
#
#    t :: name of the variable to use
#
#    Taylor :: boolean, 'true' if input is in Taylor form
#
#    Denominator :: polynomial in t, default 1 if absent
#                :: list of lists of values and derivatives 
#                   or Taylor coefficients if Taylor is true
#
#    N.B.   f[i,j] = (D@@j)(f)(tau[i])  0 <= j <= s[i]-1
#          = flist[i][j+1]
#
# Output:
#         p :: The unique (rational) Hermite interpolant, in barycentric
#              form.  Namely, 
#              w(t)*add(add(gam[i,j]/(t-tau[i])^(1+j)*add(p[i,j]/j!*(t-tau[i])^k)
#              This can be converted to distributed form by calling
#              distrib(p);
#
#         gam :: The Array(1..n, 0..smax-1) of coefficients
#
#         Here w(t) = mul( ( t-tau[i])^s[i], i = 1..n )
#              smax = max( op(s) )
#         and  s = list of sizes of each list in plist,
#              so s[i] is the number of pieces of 
#              information at node tau[i].
#
# Processing: local Laurent series.
#             This approach is different to that of
# Reference: C. Schneider & W.Werner, "Hermite
#            Interpolation: The Barycentric Approach",
#            Computing 46, 1991, pp 35-51.
#
BHIP := proc( pin::list, tau::list, t::name, 
              {Taylor::truefalse:=false}, {Conditioning::truefalse:=false},
              {Denominator::{algebraic,list}:=1} )
        local denr, dens, gam, dgam, dr, ghat,h, i, j, k, n, numr, nums,  
              p, P, q, r, rs, rt, s, smax, sq;
                  
        n := nops(tau);
        if nops(pin) <> n then
           error "Mismatched size of node list and data list"
        end if;
        
        if nops(convert(tau,set)) < n then
           error "Nodes must be distinct, with confluency explicitly specified."
        end if;
        
        p := map(t -> `if`(t::list,t,[t]),pin); # singletons ok
        
        s := map(nops,p); # confluency
        smax := max(op(s));
        if smax = 0 then 
          error "At least one piece of data is necessary."
        end if;
        
        p := `if`( Taylor, p, [seq([seq(p[i][j]/(j-1)!,j=1..s[i])],i=1..n)] );
         
        gam := Array( 1..n, 0..smax-1 ); #default 0
        if Conditioning then
           dgam := Array( 1..n, 0..smax-1, 1..n ); #default 0
        end if;
              
        # The following works for n>=1
        for i to n do 
            if s[i] > 0 then # ignore empty lists
               h[i] := mul( (t-tau[j])^s[j], j = 1..i-1 )*
                       mul( (t-tau[j])^s[j], j=i+1..n );
               r[i] := series( 1/h[i], t=tau[i], s[i] );
               for j to s[i] do
                  gam[i,s[i]-j] := coeff( r[i], t-tau[i], j-1) ; #op( 2*j-1, r[i] );
               end do;
               if Conditioning then
                  # We could compose a series for 1/(t-tau[k]) with 
                  # what we know, but using the kernel function "series"
                  # is likely faster.
                  for k to i-1 do
                     dr[i,k] := series( s[k]/h[i]/(t-tau[k]), t=tau[i], s[i] );
                     for j to s[i] do
                       dgam[i, s[i]-j, k] := coeff( dr[i,k], t-tau[i], j-1 );
                     end do;
                  end do;
                  # We could reuse earlier series, and do one O(n^2)
                  # computation to get gam[i,-1], but it's simpler to
                  # use series (and likely faster because series is in the kernel)
                  dr[i,i] := series( 1/h[i], t=tau[i], s[i]+1 );
                  # We implicitly divide this by t-tau[i], and take
                  # coefficients one higher.
                  for j to s[i] do
                     dgam[i,s[i]-j,i] := j*coeff( dr[i,i], t-tau[i], j );
                  end do;
                  for k from i+1 to n do
                     dr[i,k] := series( s[k]/h[i]/(t-tau[k]), t=tau[i], s[i] );
                     for j to s[i] do
                       dgam[i, s[i]-j, k] := coeff( dr[i,k], t-tau[i], j-1 );
                     end do;
                  end do;
               end if;
            end if;
        end do;
        
        if not (Denominator::algebraic and Denominator=1) then
           # adjust gam by folding in q
           if Denominator::list then
              if nops(Denominator)<>n then
                 error "Denominator list (q) has the wrong length."
              end if;
              q :=`if`( Taylor, q, [seq([seq(q[i][j]/(j-1)!,j=1..s[i])],i=1..n)] );
           else
              ghat := Array( 1..n );
              for i to n do
                 sq := series(Denominator,t=tau[i],s[i]);
                 ghat[i] :=[seq(coeff(sq,t-tau[i],j),j=0..s[i]-1)];
              end do;
              q := [seq(ghat[i],i=1..n)];
           end if;
           ghat := Array( 1..n, 0..smax-1 );
           for i to n do
              for j from 0 to s[i]-1 do
                 ghat[i,j] := add( gam[i,j+k]*q[i][k+1], k=0..s[i]-j-1 );
              end do;
           end do;
           gam := ghat;
        end if;
        
        P := mul( (t-tau[i])^s[i],i=1..n)*
             add(add(gam[i,j]/(t-tau[i])^(1+j)*
                     add(p[i][1+k]*(t-tau[i])^k, k=0..j),
                     j=0..s[i]-1),
                 i=1..n );
        return P, gam, `if`(Conditioning,dgam,NULL) ;
end proc:

#Condition number
BBHIP := proc( pin::list, tau::list, t::name, {Taylor::boolean:=false} )
   local B, gam, n, i, j, k, p, s;
   n := nops(tau);
   p := map(t -> `if`(t::list,t,[t]),pin); # singletons ok
        
   s := map(nops,p); # confluency

   ( p, gam ) := BHIP( pin, tau, t, ':-Taylor'=Taylor );
   
   return abs(mul((t-tau[i])^s[i],i=1..n))*add( add( 
          abs( add( gam[i,j+k]/k!/(t-tau[i])^(j+1),j=0..s[i]-1-k) ),
          k=0..s[i]-1),
          i=1..n )
end proc:

# Generalized companion matrix pencil
CMP := proc( pin::list, tau::list, t::name, {Taylor::boolean:=false} )
   local C0, C1, gam, n, i, j, k, p, P, s, d;
   n := nops(tau);
   p := map(t -> `if`(t::list,t,[t]),pin); # singletons ok
   s := map(nops,p); # confluency
   p := `if`(Taylor, p, [seq([seq(p[i][j]/(j-1)!,j=1..s[i])],i=1..n)]);
        
   d := -1 + add( s[i], i=1..n );

   # P is irrelevant here, bar shape
   ( P, gam ) := BHIP( p, tau, t, ':-Taylor'=true );
   
   C1 := Matrix( d+2, d+2, (i,j)->`if`(i<>j,0,`if`(i=d+2,0,1)) );
   C0 := Matrix( d+2, d+2, 0 );
   k := 0;
   for i to n do
      for j to s[i] do
         k := k+1;
         if j < s[i] then
            C0[k+1,k] := 1;
         end if;
         C0[k,k]   :=  tau[i];
         C0[k,d+2] :=  p[i][j];
         C0[d+2,k] := -gam[i,j-1];
      end do;
   end do;
   
   return P, C0, C1
   
end proc:

# distribution utility
distrib := proc( p )
  local w;
  if op(0,p) <> `*` then
     return p
  elif op(-1,p)::`+` then
     w := p/op(-1,p);
     return add( w*t, t in [op(op(-1,p))] )
  else
     return p
  end;
end: