1 files changed, 140 insertions, 0 deletions
diff --git a/Master/texmf/asymptote/interpolate.asy b/Master/texmf/asymptote/interpolate.asy
new file mode 100644
index 00000000000..a6b820bc7d5
--- /dev/null
+++ b/Master/texmf/asymptote/interpolate.asy
@@ -0,0 +1,140 @@
+// Lagrange and Hermite interpolation in Asymptote
+// Author: Olivier Guibé
+// Acknowledgements: Philippe Ivaldi
+
+// diffdiv(x,y) computes Newton's Divided Difference for
+// Lagrange interpolation with distinct values {x_0,..,x_n} in the array x
+// and values y_0,...,y_n in the array y,
+
+// hdiffdiv(x,y,dyp) computes Newton's Divided Difference for
+// Hermite interpolation where dyp={dy_0,...,dy_n}. 
+// 
+// fhorner(x,coeff) uses Horner's rule to compute 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}),
+// where coeff={a_0,a_1,...,a_n}.
+
+// fspline does standard cubic spline interpolation of a function f
+// on the interval [a,b].
+// The points a=x_1 < x_2 < .. < x_n=b form the array x;
+// the points y_1=f(x_1),....,y_n=f(x_n) form the array y
+// We use the Hermite form for the spline.
+
+// The syntax is:
+// s=fspline(x,y);                            default not_a_knot condition
+// s=fspline(x,y,natural);                    natural spline
+// s=fspline(x,y,periodic);                   periodic spline
+// s=fspline(x,y,clamped(1,1));               clamped spline
+// s=fspline(x,y,monotonic);                  piecewise monotonic spline
+
+// Here s is a real function that is constant on (-infinity,a] and [b,infinity).
+
+private import math;
+import splinetype;
+
+typedef real fhorner(real);
+
+struct horner {
+  // x={x0,..,xn}(not necessarily distinct)
+  // a={a0,..,an} corresponds to the polyonmial
+  // a_0+a_1(x-x_0)+a_2(x-x_0)(x-x_1)+...+a_n(x-x_0)..(x-x_{n-1}),
+  real[] x;
+  real[] a;
+}
+
+// Evaluate p(x)=d0+(x-x0)(d1+(x-x1)+...+(d(n-1)+(x-x(n-1))*dn)))
+// via Horner's rule: n-1 multiplications, 2n-2 additions.
+fhorner fhorner(horner sh)
+{
+  int n=sh.x.length;
+  checklengths(n,sh.a.length);
+  return new real(real x) {
+    real s=sh.a[n-1];
+    for(int k=n-2; k >= 0; --k)
+      s=sh.a[k]+(x-sh.x[k])*s;
+    return s;
+  };
+}
+      
+// Newton's Divided Difference method: n(n-1)/2 divisions, n(n-1) additions.
+horner diffdiv(real[] x, real[] y)
+{
+  int n=x.length;
+  horner s;
+  checklengths(n,y.length);
+  for(int i=0; i < n; ++i)
+    s.a[i]=y[i];
+  for(int k=0; k < n-1; ++k) {
+    for(int i=n-1; i > k; --i) {
+      s.a[i]=(s.a[i]-s.a[i-1])/(x[i]-x[i-k-1]);
+    }
+  }
+  s.x=x;
+  return s;
+}
+
+// Newton's Divided Difference for simple Hermite interpolation,
+// where one specifies both p(x_i) and p'(x_i).
+horner hdiffdiv(real[] x, real[] y, real[] dy)
+{
+  int n=x.length;
+  horner s;
+  checklengths(n,y.length);
+  checklengths(n,dy.length);
+  for(int i=0; i < n; ++i) {
+    s.a[2*i]=y[i];
+    s.a[2*i+1]=dy[i];
+    s.x[2*i]=x[i];
+    s.x[2*i+1]=x[i];
+  }
+
+  for(int i=n-1; i > 0; --i)
+    s.a[2*i]=(s.a[2*i]-s.a[2*i-2])/(x[i]-x[i-1]);
+
+  int stop=2*n-1;
+  for(int k=1; k < stop; ++k) {
+    for(int i=stop; i > k; --i) {
+      s.a[i]=(s.a[i]-s.a[i-1])/(s.x[i]-s.x[i-k-1]);
+    }
+  }
+  return s;
+}
+
+typedef real realfunction(real);
+
+// piecewise Hermite interpolation:
+// return the piecewise polynomial p(x), where on [x_i,x_i+1], deg(p) <= 3,
+// p(x_i)=y_i, p(x_{i+1})=y_i+1, p'(x_i)=dy_i, and p'(x_{i+1})=dy_i+1.
+// Outside [x_1,x_n] the returned function is constant: y_1 on (infinity,x_1]
+// and y_n on [x_n,infinity).
+realfunction pwhermite(real[] x, real[] y, real[] dy)
+{
+  int n=x.length;
+  checklengths(n,y.length);
+  checklengths(n,dy.length);
+  if(n < 2) abort(morepoints);
+  if(!increasing(x,strict=true)) abort("array x is not strictly increasing");
+  return new real(real t) {
+    int i=search(x,t);
+    if(i == n-1) {
+      i=n-2;
+      t=x[n-1];
+    } else if(i == -1) {
+      i=0;
+      t=x[0];
+    }
+    real h=x[i+1]-x[i];
+    real delta=(y[i+1]-y[i])/h;
+    real e=(3*delta-2*dy[i]-dy[i+1])/h;
+    real f=(dy[i]-2*delta+dy[i+1])/h^2;
+    real s=t-x[i];
+    return y[i]+s*(dy[i]+s*(e+s*f));
+  };
+}
+
+realfunction fspline(real[] x, real[] y, splinetype splinetype=notaknot)
+{
+  real[] dy=splinetype(x,y);
+  return new real(real t) {
+    return pwhermite(x,y,dy)(t);
+  }; 
+}