with(LinearAlgebra):

va:=proc(t,s,d) local V,i,j;
V:=<seq(t^i,i=0..d)>;
for j from 2 to s do 
  V:=Matrix(d+1,j,[V,map(x->diff(x,t),Column(V,j-1))]);
 od;
Transpose(V);
end:

#Confluent Vandermonde matrix , T=tau; S=s1,s2,..,sN,  Definition 5:

VandHer:=proc(T,S) local Va,v,i,d;
d:=add( i, i in S)-1;
Va:=subs(t=T[1],va(t,S[1],d));
if nops(T)>1 then for i from 2 to nops(T) do
Va:=<<Va,subs(t=T[i],va(t,S[i],d))>>;od;fi;
Va;
end;

#Given i, compute k(i) y r(i). Definition 6:

kiri:=proc(i,S) local k,r,s,ki,ri;
s:=0;ki:=0;
while s<i do ki:=ki+1;s:=s+S[ki];od;
ri:=-1+i-(s-S[ki]);
return [ki,ri];
end;


#confluent divided differences table in order to compute the derivatives -Definition 8 with P[t_k^a,t_j^b]:

#P :: list 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] ]

tabla:=proc(T,S,P,a,b,k,j) local m,i,c,Tb; 
Tb:= Array(1..a+1,1..b+1);
Tb[2,2]:=(P[add(S[i],i=1..j-1)+1]-P[add(S[i],i=1..k-1)+1])/(T[j]-T[k]);
if a=b and b=1 then return Tb;
 else 
if a>1 then for c from 3 to a+1 do Tb[c,1]:=P[add(S[i],i=1..k-1)+c-2+1]/(c-2)!  ;od;fi;   
if b>1 then for c from 3 to b+1 do Tb[1,c]:=P[add(S[i],i=1..j-1)+c-2+1]/(c-2)!  ;od;fi; 
for i from 3 to a+1 do Tb[i,2]:=(Tb[i-1,2]-Tb[i,1])/(T[j]-T[k]);od; 
for i from 3 to b+1 do Tb[2,i]:=(Tb[1,i]-Tb[2,i-1])/(T[j]-T[k]);od; 
for i from 3 to a+1 do
   for c from 3 to b+1 do Tb[i,c]:=(Tb[i-1,c]-Tb[i,c-1])/(T[j]-T[k]);od;od;
fi;
return Tb;
end:

#Generalized barycentric weights-Proposition 2:

bari:=proc(T,S) local K,f,L,l,i,j; 
f:=1/mul((t-T[i])^S[i],i=1..nops(T));
K:=convert(f,parfrac,t,true);
L:=[];
for i from 1 to nops(T) do  
 l:=[]; 
 for j from 1 to S[i] do
  l:=[op(l),coeff(K,1/(t-T[i])^j)];
 od;
 L:=[op(L),l];
od;
return L;
end:



#Lemma 8: derivatives from divided differences

derivada:=proc(T,S,P,j,orden,bc) local bcjsj,m,A,B,i,rho,k,ta,der;  
if orden <S[j] then der:=P[add(S[i],i=1..j-1)+orden + 1];
else
m:=orden-S[j];  
bcjsj:=bc[j][S[j]]; if bcjsj=0 then error "baricentro j,sj-1 nulo" fi;
if S[j]>1 then  
          B:=add(bc[j][rho+1]*derivada(T,S,P,j,rho+m+1,bc)/(rho+m+1)!,rho=0..S[j]-2);
          else B:=0;
fi;
A:=0;
for k from 1 to nops(T) do if k<>j then
 ta:=tabla(T,S,P,S[k],m+1,k,j);   
 A:= add(bc[k][rho+1]*ta[rho+1+1,m+1+1],rho=0..S[k]-1)+A ; 
                          fi;
od; 
der :=-orden!*(A+B)/bcjsj;
fi;
return der;
end:



#Extended confluent Bezout matrix, using divided differences

bezH:=proc(T,S,P,Q) local B,d,i,j,alpha,beta,ri,rj,ki,kj,nu,bc;
d:=add(i, i in S)-1;
B:=Matrix(d+1,d+1,shape=symmetric);
bc:=bari(T,S);
for i from 1 to d+1  do
 for j from i to d+1 do  
   ri:=kiri(i,S)[2];rj:=kiri(j,S)[2]; ki:=kiri(i,S)[1];kj:=kiri(j,S)[1];
   if ki<> kj then
     B[i,j]:=add(add((-1)^(ri+alpha)*binomial(ri, alpha)*binomial(rj, beta)*(ri+rj-alpha-beta)!*(P[add(S[r],r=1..ki-1)+alpha+1]*Q[add(S[r],r=1..kj-1)+beta+1]-P[add(S[r],r=1..kj-1)+beta+1]*Q[add(S[r],r=1..ki-1)+alpha+1])/((T[ki]-T[kj])^(ri+rj-alpha-beta+1)) ,beta=0..rj),alpha=0..ri); 
  else  
      B[i,j]:=0; 
     for nu from  ri+1 to d do 
      for beta from 0 to ri+rj-nu+1  do 
     B[i,j]:= (-1)^(nu-ri-1)*rj!*(derivada(T,S,P,ki,nu+beta,bc )*derivada(T,S,Q,ki,ri+rj-nu-beta+1,bc )-derivada(T,S,P,ki,ri+rj-nu-beta+1,bc)*derivada(T,S,Q,ki,nu+beta,bc) )/(nu*(beta!)*(ri+rj-nu-beta+1)!*(nu-ri-1)!)  +B[i,j]; 
     od;
      od;
   fi;
 od;
od;
return B;
end:

#The Confluent Bezout matrix , using divided differences

bezHred:=proc(T,S,P,Q) local B,d,i,j,alpha,beta,ri,rj,ki,kj,nu,bc;
d:=add(i, i in S)-1;
B:=Matrix(d,d,shape=symmetric);
bc:=bari(T,S);
for i from 1 to d  do
 for j from i to d do  
   ri:=kiri(i,S)[2];rj:=kiri(j,S)[2]; ki:=kiri(i,S)[1];kj:=kiri(j,S)[1];
   if ki<> kj then
     B[i,j]:=add(add((-1)^(ri+alpha)*binomial(ri, alpha)*binomial(rj, beta)*(ri+rj-alpha-beta)!*(P[add(S[r],r=1..ki-1)+alpha+1]*Q[add(S[r],r=1..kj-1)+beta+1]-P[add(S[r],r=1..kj-1)+beta+1]*Q[add(S[r],r=1..ki-1)+alpha+1])/((T[ki]-T[kj])^(ri+rj-alpha-beta+1)) ,beta=0..rj),alpha=0..ri); 
  else  
      B[i,j]:=0; 
     for nu from  ri+1 to d do 
      for beta from 0 to ri+rj-nu+1  do 
        
     B[i,j]:= (-1)^(nu-ri-1)*rj!*(derivada(T,S,P,ki,nu+beta,bc )*derivada(T,S,Q,ki,ri+rj-nu-beta+1,bc )-derivada(T,S,P,ki,ri+rj-nu-beta+1,bc)*derivada(T,S,Q,ki,nu+beta,bc) )/(nu*(beta!)*(ri+rj-nu-beta+1)!*(nu-ri-1)!)  +B[i,j]; 
     od;
      od;
   fi;
 od;
od;
return B;
end:

#####################################
################# derivatives computed from the differentiation matrix##############
#
#The user must download the code "BHIP.mpl" from http://www.nfillion.com/index.php/resources/ginm-code-repository
# 
#
#P :: 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]] ]



read"BHIP.mpl.txt":

derivadas_matrix:=proc(T,S,P) local DM,L,N,Pf,Pt,Pn,Pv,i,j,k,r,om,sm,vi,vn;#Pt: lista en Taylor form
N:=nops(T);
Pt:=[]; 
for i from 1 to N do
 L:=[]:
 for j from 1 to nops(P[i]) do
  L:=[op(L),P[i][j]/(j-1)!];
 od;
 Pt:=[op(Pt),L];
od;
Pf:=P;  
DM:= BHIP( Pt, T,t,Dmat=true)[3];
sm:=max(S);
om:=2*max(S)-1; 
vi:=convert(Pt,Vector);#vector inicial
Pn:=Pt;
for i from 1 to sm do
 vn:=DM.vi; 
 for j from 1 to N do
   Pf[j]:=[op(Pf[j]), (S[j]+1-2)!*vn[add(S[k],k=1..j)]]; 
   Pn[j]:=subsop(1=NULL,Pn[j]  );
   Pn[j]:=[op(Pn[j]), vn[add(S[k],k=1..j)]]; 
 od; 
 vi:=vn;
od; 
return Pf;
end:

#extended confluent Bezout matrix

bezH_DM:=proc(T,S,P,Q) local B,d,i,j,alpha,beta,ri,rj,ki,kj,nu,PP,QQ,derivadasP,derivadasQ;
d:=add(i, i in S)-1;
B:=Matrix(d+1,d+1,shape=symmetric);
derivadasP:=derivadas_matrix(T,S,P);
derivadasQ:=derivadas_matrix(T,S,Q);
PP:=convert(convert(P,Vector),list);
QQ:=convert(convert(Q,Vector),list); 
for i from 1 to d+1  do
 for j from i to d+1 do 
   ri:=kiri(i,S)[2];rj:=kiri(j,S)[2]; ki:=kiri(i,S)[1];kj:=kiri(j,S)[1];
   if ki<> kj then
     B[i,j]:=add(add((-1)^(ri+alpha)*binomial(ri, alpha)*binomial(rj, beta)*(ri+rj-alpha-beta)!*(PP[add(S[r],r=1..ki-1)+alpha+1]*QQ[add(S[r],r=1..kj-1)+beta+1]-PP[add(S[r],r=1..kj-1)+beta+1]*QQ[add(S[r],r=1..ki-1)+alpha+1])/((T[ki]-T[kj])^(ri+rj-alpha-beta+1)) ,beta=0..rj),alpha=0..ri); 
  else   
      B[i,j]:=0; 
     for nu from  ri+1 to d do 
      for beta from 0 to ri+rj-nu+1  do 
         B[i,j]:= (-1)^(nu-ri-1)*rj!*(derivadasP[ki][nu+beta+1]*derivadasQ[ki][ri+rj-nu-beta+1+1]-derivadasP[ki][ri+rj-nu-beta+1+1]*derivadasQ[ki][nu+beta+1])/(nu*(beta!)*(ri+rj-nu-beta+1)!*(nu-ri-1)!)  +B[i,j];
     od;
      od;
   fi;
 od;
od;
return B;
end:


#Confluent Bezout matrix


bezH_DM_red:=proc(T,S,P,Q) local B,d,i,j,alpha,beta,ri,rj,ki,kj,nu,PP,QQ,derivadasP,derivadasQ;
d:=add(i, i in S)-1;
B:=Matrix(d,d,shape=symmetric);
derivadasP:=derivadas_matrix(T,S,P);
derivadasQ:=derivadas_matrix(T,S,Q);
PP:=convert(convert(P,Vector),list);
QQ:=convert(convert(Q,Vector),list); 
for i from 1 to d  do
 for j from i to d do 
   ri:=kiri(i,S)[2];rj:=kiri(j,S)[2]; ki:=kiri(i,S)[1];kj:=kiri(j,S)[1];
   if ki<> kj then
     B[i,j]:=add(add((-1)^(ri+alpha)*binomial(ri, alpha)*binomial(rj, beta)*(ri+rj-alpha-beta)!*(PP[add(S[r],r=1..ki-1)+alpha+1]*QQ[add(S[r],r=1..kj-1)+beta+1]-PP[add(S[r],r=1..kj-1)+beta+1]*QQ[add(S[r],r=1..ki-1)+alpha+1])/((T[ki]-T[kj])^(ri+rj-alpha-beta+1)) ,beta=0..rj),alpha=0..ri); 
  else   
      B[i,j]:=0; 
     for nu from  ri+1 to d do 
      for beta from 0 to ri+rj-nu+1  do 
B[i,j]:= (-1)^(nu-ri-1)*rj!*(derivadasP[ki][nu+beta+1]*derivadasQ[ki][ri+rj-nu-beta+1+1]-derivadasP[ki][ri+rj-nu-beta+1+1]*derivadasQ[ki][nu+beta+1])/(nu*(beta!)*(ri+rj-nu-beta+1)!*(nu-ri-1)!)  +B[i,j];
     od;
      od;
   fi;
 od;
od;
return B;
end: