function D=diff_div(X,Y)


%	diff_div : calcul des différences divisées pour des points de supports quelconques
%
% *********************************************************
%
%	D=diff_div(X,Y)
%     les différences divisées sont calculées par l'algorithme 
%     pyramidal donné en TD.
%
%	variables d'entrées : 
%     * X : contient les valeurs x_i, pour 0 <=i<=n (dans l'ordre croissant)
%     * Y : contient les valeurs f^(r)(x_i), pour 0 <=i<=n, de telle sorte que :
%             si x_k est multiple avec x_{k-1}~=x_k=x_{k+1}=...=x_{k+r}~=x_{k+r+1} alors
%                y_k=f(x_k), y_{k+1}=f'(x_k), ..., y_{k+r}=f^(r)(x_k).
%
%	variables de sortie
%     * D : contient les différences divisées (généralisée à support quelconque)	 
%           f[x_0], f[x_0,x_1], ...., f[x_0,x_1,...,x_n]
%
%
%
% ************ Fonctions auxiliaires utilisées ************
%
%   val_derivee_multiple     
%
% *********************************************************


% Contrôles d'entrée
%
% nombre d'arguments
if nargin~=2
   error('nombre d''arguments de la fonction incorrect');
end
% taille des deux variables 
n=length(X)-1;
np=length(Y)-1;
if n~=np
   error('les deux tableaux n''ont pas la même taille');
end
% vérification de l'ordre croissant des x_i 
% (si X non formel)
if isnumeric(X)
   if max(abs(sort(X)-X))~=0
      error('support non ordonné');
   end
end

% Corps d'algorithme
i=0;
while (i<=n)
   [alpha,valfder]=val_derivee_multiple(i,X,Y);
   D(i+1:i+alpha)=valfder(1)*ones(1,alpha);
   i=i+alpha;
end
factoi=1;
for i=1:n;
   factoi=factoi*i;
   for j=n:-1:i;
       dx=X(j+1)-X(j-i+1);
       if (dx==0)
          [alpha,valfder]=val_derivee_multiple(j-i,X,Y);
          D(j+1)=(valfder(i+1))/factoi;
       else
          D(j+1)=(D(j+1)-D(j))/(dx);
       end
    end
end
