Diferenças divididas e polinómio de Newton

> #restart;
with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> #-------------------------------------------

> DDivididas:=proc(i0,i1,x::vector,f::vector,d::posint)
local s;
Digits:=d;
if i0=i1 then
s:=evalf(f[i0])
else s:=evalf((DDivididas(i0+1,i1,x,f,d)-DDivididas(i0,i1-1,x,f,d))/(x[i1]-x[i0]));
fi;
RETURN(s);
end:

> #-------------------------------------------

> PoliNewton:=proc(X::vector,fX::vector,d::posint)
local i,s,p,n;
Digits:=d;
n:=vectdim(X);
s:=0;
p:=1;
for i from 1 to n do
s:=s+DDivididas(1,i,X,fX,d)*p;
p:=p*evalf(x-X[i]);
od;
RETURN(simplify(s));
end:

> #-------------------------------------------

> PoliNewtonPt:=proc(pt::numeric,X::vector,fX::vector,d::posint)
local polinomio, res;
Digits:=d;
polinomio:=PoliNewton(X,fX,d);
res:=subs(x=pt,polinomio);
RETURN(evalf(res));
end: