// Asymptote Lagrange and Hermite interpolation
// Author : Olivier Guibé
// Thanks : Philippe Ivaldi

// We compute with Asymptote the Lagrange or Hermite
// interpolation polynomial.
// To this end we use --of course-- the Newton's Divided Difference
// together with Horner's rule.
// Let x_0,..,x_n (distincts) (the array xpt) and y_0,...,y_n
// (the array ypt). Diffdiv(xpt,ypt) computes the Newton's Divided Difference
// for Lagrange interpolation.
// Moreover with dy_0,...,dy_n (the array dyp) hdiffdiv(xpt,ypt,dyp)
// computes the Newton's Divided Difference for Hermite interpolation.
// 
// fhorner(xpt, coeff) uses the Horner's rule and
// is the polynomial a_0+a_1(x-x_0)+a_2(x-x_0)(x-x_1)+...+a_n(x-x_0)..(x-x_{n-1})
// coeff is the array a_0,a_1,...,a_n


import graph;

typedef real fhorner(real);

struct horner { // deux tableaux de r\'eels
  // xpt sont les x0,..,xn (non n\'ecessairement distincts)
  // dd sont d0,..,dn. Correspond au polyn^ome
  // d0+(x-x0)d1+(x-x0)(x-x1)d2+....+ (x-x0)..(x-x(n-1))dn
  real [] xpt;
  real [] dd;
}
horner operator init() {return new horner;}

fhorner fhorner( horner sh)
{// algorithme d'\'evaluation de p(x) sous la forme
  // p(x)=d0+(x-x0)(d1+(x-x1).....+(d(n-1)+(x-x(n-1))*dn)))
  // en n(*), n(+)
  return new real(real x)
    {
      int n=sh.xpt.length;
      if (n != sh.dd.length) abort("Horner: nombres de paramètres incorrectes.");
      real s;
      s=sh.dd[n-1];
      for (int k=n-2;k>=0; --k)
	{
	  s=sh.dd[k]+(x-sh.xpt[k])*s;
	}
      return s;
    };
}
      
horner diffdiv(real [] xpt, real [] ypt)
{// algorithme des diff\'erences divis\'ees de Newton, 
  // une fois que c'est fait, n^2/2 (*) op\'erations.
  // retourne s de structure horner.
 int n=xpt.length;
 horner s;
    if (n != ypt.length) abort("Diff Div: nombres de paramètres incorrectes.");
    real [] d;
    for (int i=0;i<n; ++i)
      {
	s.dd[i]=ypt[i];
      }
    for (int k=0; k<n-1; ++k)
      {
	for (int i=n-1;i>k; --i)
	  {
	    s.dd[i]=(s.dd[i]-s.dd[i-1])/(xpt[i]-xpt[i-k-1]);
	  }
      }
    s.xpt=xpt;
    return s;
}

horner hdiffdiv(real [] xpt, real [] ypt, real [] dypt)
{  // algorithme des diff\'erences divis\'ees de Newton 
  // pour l'interpolation d'Hermite simple, i.e. on impose
  // p(xi) et p'(xi). La diff\'erence tient juste dans
  // la 2\`eme \'etape. Retourne s de structure horner.
  int n=xpt.length;
  horner s;
  if (n != ypt.length) abort("Hermite Diff Div: nombres de paramètres incorrectes.");
  if (n != dypt.length) abort("Hermite Diff Div: nombres de paramètres incorrectes.");
  real [] d;
   for (int i=0;i<n; ++i)
      {
	s.dd[2*i]=ypt[i];
	s.dd[2*i+1]=dypt[i];
	s.xpt[2*i]=xpt[i];
	s.xpt[2*i+1]=xpt[i];
      }
   for (int i=n-1; i>0; --i)
     {
       s.dd[2*i]=(s.dd[2*i]-s.dd[2*i-2])/(xpt[i]-xpt[i-1]);
     }
    for (int k=1; k<2*n-1; ++k)
      {
	for (int i=2*n-1;i>k; --i)
	  {
	    s.dd[i]=(s.dd[i]-s.dd[i-1])/(s.xpt[i]-s.xpt[i-k-1]);
	  }
      }
    return s;
}


// Test 1 ou l'effet de Runge pour Lagrange
// et la fameuse fonction 1/(x^2+1)

	

unitsize(2cm);
xlimits(-5,5);
ylimits(-1,1,Crop);
real a,b;
real f ( real x) { return (1/(x^2+1));}
real df( real x) { return (-2*x/(x^2+1)^2);}

//real f (real x) {return (sin(x));}
//real df( real x) {return (cos(x));}
a=-5; b=5;
int n=15;
real [] xpt,ypt,dypt;
fhorner p,ph,ph1;
horner ddh,dd;
xpt=a+(b-a)*sequence(n+1)/n;
ypt=map(f,xpt);
dypt=map(df,xpt);
for(int i=0; i<=n;++i)
{
 dot((xpt[i],ypt[i]),5bp+blue);
}
dd=diffdiv(xpt,ypt);
p=fhorner(dd);
draw(graph(p,a,b,n=500),"$x\longmapsto{}L_{15}$");
draw(graph(f,a,b),red,"$x\longmapsto{}\frac{1}{x^2+1}$");

xaxis("$x$",BottomTop,LeftTicks);
yaxis("$y$",LeftRight,RightTicks);
attach(legend(),point(10S),30S);
shipout("runge1");


// Test 2 ou l'effet de Runge pour Hermite
// et toujours la fameuse fonctin 1/(x^2+1)
erase();


real a,b;
real f ( real x) { return (1/(x^2+1));}
real df( real x) { return (-2*x/(x^2+1)^2);}

//real f (real x) {return (sin(x));}
//real df( real x) {return (cos(x));}
a=-5; b=5;
int n=16;
real [] xpt,ypt,dypt;
fhorner p,ph,ph1;
horner ddh,dd;
xpt=a+(b-a)*sequence(n+1)/n;
ypt=map(f,xpt);
dypt=map(df,xpt);
for(int i=0; i<=n;++i)
{
 dot((xpt[i],ypt[i]),5bp+blue);
}
ddh=hdiffdiv(xpt,ypt,dypt);
ph=fhorner(ddh);
draw(graph(ph,a,b,n=1000),"$x\longmapsto{}H_{16}$");
draw(graph(f,a,b),red,"$x\longmapsto{}\frac{1}{x^2+1}$");

unitsize(2cm);
xlimits(-5,5);
ylimits(-1,5,Crop);

xaxis("$x$",BottomTop,LeftTicks);
yaxis("$y$",LeftRight,RightTicks);
attach(legend(),point(10S),30S);

shipout("runge2");


// Test 3, l'effet de Runge n'est pas r\'eserv\'e
// \`a toute les fonctions. Interpolation de Lagrange
// pour une fonction dont toutes les d\'eriv\'ees successives
// sont born\'ees par une constante M (ici M=1).
// On d\'emontre que \c ca converge

erase();


real a,b;

real f (real x) {return (sin(x));}
real df( real x) {return (cos(x));}
a=-5; b=5;
int n=16;
real [] xpt,ypt,dypt;
fhorner p,ph,ph1;
horner ddh,dd;
xpt=a+(b-a)*sequence(n+1)/n;
ypt=map(f,xpt);
dypt=map(df,xpt);
for(int i=0; i<=n;++i)
{
 dot((xpt[i],ypt[i]),5bp+blue);
}
dd=diffdiv(xpt,ypt);
p=fhorner(dd);

draw(graph(p,a,b,n=500),"$x\longmapsto{}L_{16}$");
draw(graph(f,a,b),red,"$x\longmapsto{}\cos(x)$");

xaxis("$x$",BottomTop,LeftTicks);
yaxis("$y$",LeftRight,RightTicks);
attach(legend(),point(10S),30S);
shipout("runge3");


// Test 4, cependant on remarque ici les limites
// du calcul num\'erique avec une pr\'ecision de 10^(-16)
// ph\'enom\`ene ici purement num\'erique
erase();


real a,b;

real f (real x) {return (sin(x));}
real df( real x) {return (cos(x));}
a=-5; b=5;
int n=72;
real [] xpt,ypt,dypt;
fhorner p,ph,ph1;
horner ddh,dd;
xpt=a+(b-a)*sequence(n+1)/n;
ypt=map(f,xpt);
dypt=map(df,xpt);
for(int i=0; i<=n;++i)
{
 dot((xpt[i],ypt[i]),5bp+blue);
}
dd=diffdiv(xpt,ypt);
p=fhorner(dd);


draw(graph(p,a,b,n=500),"$x\longmapsto{}L_{72}$");
draw(graph(f,a,b),red,"$x\longmapsto{}\cos(x)$");
ylimits(-1,5,Crop);
xaxis("$x$",BottomTop,LeftTicks);
yaxis("$y$",LeftRight,RightTicks);
attach(legend(),point(10S),30S);


shipout("runge4");

// Test 5, avec les points de Tchebychev
// ce n'est que du bonheur, ou presque
// les points sont compliqu\'es \`a calculer
erase();	

unitsize(2cm);
xlimits(-5,5);
ylimits(-1,2,Crop);
real a,b;
real f ( real x) { return (1/(x^2+1));}
real df( real x) { return (-2*x/(x^2+1)^2);}

//real f (real x) {return (sin(x));}
//real df( real x) {return (cos(x));}
a=-5; b=5;
int n=16;
real [] xpt,ypt,dypt;
fhorner p,ph,ph1;
horner ddh,dd;
for(int i=0; i<=n;++i)
{
  xpt[i]=(a+b)/2+(b-a)/2*cos((2*i+1)/(2*n+2)*pi);
}
ypt=map(f,xpt);
dypt=map(df,xpt);
for(int i=0; i<=n;++i)
{
 dot((xpt[i],ypt[i]),5bp+blue);
}
dd=diffdiv(xpt,ypt);
p=fhorner(dd);

draw(graph(p,a,b,n=500),"$x\longmapsto{}L_{16}$");
draw(graph(f,a,b),red,"$x\longmapsto{}\frac{1}{x^2+1}$");

xaxis("$x$",BottomTop,LeftTicks);
yaxis("$y$",LeftRight,RightTicks);
attach(legend(),point(10S),30S);

shipout("runge5");


// Test 6, avec les points de Tchebychev
// ce n'est que du bonheur, ou presque
// les points sont compliqu\'es \`a calculer

erase();	
unitsize(2cm);
xlimits(-5,5);
ylimits(-1,2,Crop);
real a,b;
real f ( real x) { return (1/(x^2+1));}
real df( real x) { return (-2*x/(x^2+1)^2);}

//real f (real x) {return (sin(x));}
//real df( real x) {return (cos(x));}
a=-5; b=5;
int n=26;
real [] xpt,ypt,dypt;
fhorner p,ph,ph1;
horner ddh,dd;
for(int i=0; i<=n;++i)
{
  xpt[i]=(a+b)/2+(b-a)/2*cos((2*i+1)/(2*n+2)*pi);
}
ypt=map(f,xpt);
dypt=map(df,xpt);
for(int i=0; i<=n;++i)
{
 dot((xpt[i],ypt[i]),5bp+blue);
}
dd=diffdiv(xpt,ypt);
p=fhorner(dd);
draw(graph(p,a,b,n=500),"$x\longmapsto{}L_{26}$");
draw(graph(f,a,b),red,"$x\longmapsto{}\frac{1}{x^2+1}$");

xaxis("$x$",BottomTop,LeftTicks);
yaxis("$y$",LeftRight,RightTicks);
attach(legend(),point(10S),30S);


shipout("runge6");


// Test 7, avec les points de Tchebychev
// ce n'est que du bonheur, ou presque
// les points sont compliqu\'es \`a calculer

erase();	
unitsize(2cm);
xlimits(-2,2);
ylimits(-0.5,2,Crop);
real a,b;
real f ( real x) { return (sqrt(abs(x-1)));}


//real f (real x) {return (sin(x));}
//real df( real x) {return (cos(x));}
a=-2; b=2;
int n=30;
real [] xpt,ypt,dypt;
fhorner p,ph,ph1;
horner ddh,dd;
for(int i=0; i<=n;++i)
{
  xpt[i]=(a+b)/2+(b-a)/2*cos((2*i+1)/(2*n+2)*pi);
}
ypt=map(f,xpt);
dypt=map(df,xpt);
for(int i=0; i<=n;++i)
{
 dot((xpt[i],ypt[i]),5bp+blue);
}
dd=diffdiv(xpt,ypt);
p=fhorner(dd);
draw(graph(p,a,b,n=500),"$x\longmapsto{}L_{30}$");
draw(graph(f,a,b),red,"$x\longmapsto{}\sqrt{|x-1|}$");

xaxis("$x$",BottomTop,LeftTicks);
yaxis("$y$",LeftRight,RightTicks);
attach(legend(),point(10S),30S);


shipout("runge7");