#!/usr/bin/octave
disp("Testing Divided Differences for sin(x) sin(4x) example [-5,5]  f(x)=1/(1+x^2)");


 
 
function [a,errkn]=RoznDzielInterLagr(x,y,n=length(x))
#Function computes divided differences of f with y_k=f(x_k)
#Input 
#x=[x_1...x_n] n interpolation knots 
#y=[y_1...y_n] values of f at  interpolation knots 
#Output
#a=[a_1 ... a_n] divided differences of f at x_k
#errkn Max of absolute value of Error L_{n-1}f - f at  knots (should be zero up to rounding error)

if(length(y)!=length(x))
 printf("Wrong input. No of knots differs from no of values at knots\n");
 return;  
endif

a=y;
for k=1:n-1,
%Slower loop version
% for j=n:-1:k+1,
%    h=x(j)-x(j-k);
%    if(abs(h)>1e-10)
%      a(j)=(a(j)-a(j-1))/h;
%    endif
% endfor #j  
#vector version
  h=x(k+1:n)-x(k:n-1);
  if(norm(h,'inf')>1e-8)
    a(k+1:n)=(a(k+1:n)-a(k:n-1))./h;
  endif
endfor #k
err=HornerNewton(a,x,x)-y;
errkn=norm(err,'inf');
if(errkn>1e-8)
 printf("Warning. max(|y_k-p(x_k)})>1e-8\n"); 
endif

endfunction


function y=HornerNewton(a,x,t)
#Schemat Hornera dla wiel w bazie newtona a(1)+a(2)(x-x1)+a(3)(x-x1)(x-x2) itd
#a wspolczynniki w bazie 
#x - wezly (co najmniej n-1)
#t - wartosc w ktorej chcemy policzyc wartosc

n=length(a);
if(length(x) < n-1) 
printf("za malo wezlow");
end

y=a(n);
for k=n-1:-1:1,
y=y.*(t.-x(k)) + a(k);
endfor

endfunction



########################TESTING sin#####################################
a=-2*pi;b=pi;
rf=@(x) sin(x);
s=linspace(a,b,200); #equidistant mesh (200 points)
rfs=rf(s);
N=2;
for k=1:3,
 N*=2
#equidistant knots 
 x=linspace(a,b,N+1);
 y=rf(x);
 c=RoznDzielInterLagr(x,y);
 li=HornerNewton(c,x,s);
 err=norm(rfs-li,'inf')
#Chebyshev knots 
 xc=0.5*(a+b)+ 0.5*(b-a)*cos(pi*(0.5+(0:N))/(N+1));
 yc=rf(xc);
 cc=RoznDzielInterLagr(xc,yc);
 lic=HornerNewton(cc,xc,s);
 errc=norm(rfs-lic,'inf')
 plot(s,rfs,";sin f;",s,li,[";LI ", num2str(N+1), " equidist knots;"],s,lic,[";LI ", num2str(N+1), " Chebys knots;"]);
 pause(1);

endfor

rf=@(x) sin(4*x);
s=linspace(a,b,200); #equidistant mesh (200 points)
rfs=rf(s);
N=2;
for k=1:4,
 N*=2
#equidistant knots 
 x=linspace(a,b,N+1);
 y=rf(x);
 c=RoznDzielInterLagr(x,y);
 li=HornerNewton(c,x,s);
 err=norm(rfs-li,'inf')
#Chebyshev knots 
 xc=0.5*(a+b)+ 0.5*(b-a)*cos(pi*(0.5+(0:N))/(N+1));
 yc=rf(xc);
 cc=RoznDzielInterLagr(xc,yc);
 lic=HornerNewton(cc,xc,s);
 errc=norm(rfs-lic,'inf')
 plot(s,rfs,";sin(4x) f;",s,li,[";LI ", num2str(N+1), " equidist knots;"],s,lic,[";LI ", num2str(N+1), " Chebys knots;"]);
 pause(1);

endfor



#now polyfit for same N

disp("Testing octave interpolation function polyfit( vs our divided differences algorithm");
N
cc=polyfit(xc,yc,N+1);
norm(rfs-polyval(cc,s),'inf')
disp("and with fewer  knots?");
N=24
 xc=0.5*(a+b)+ 0.5*(b-a)*cos(pi*(0.5+(0:N))/(N+1));
 yc=rf(xc);
 cc=RoznDzielInterLagr(xc,yc);
 lic=HornerNewton(cc,xc,s);
 errc=norm(rfs-lic,'inf')
ac=polyfit(xc,yc,N+1);
norm(rfs-polyval(ac,s),'inf')
#w wezlach...

norm(yc-polyval(ac,xc),'inf')
printf("why polyfit does not work? ...
Look at the warning after calling polyfit..... polyfit calls backslash \n");
#and solver for linear systems fails - matrix is numerically singular - i.e. condition of the matrix too large
