1 files changed, 108 insertions, 0 deletions
diff --git a/src/interpolation.jl b/src/interpolation.jl
new file mode 100644
index 0000000..fa3bcdb
--- /dev/null
+++ b/src/interpolation.jl
@@ -0,0 +1,108 @@
+## Copyright 2021 Ian Jauslin
+## 
+## Licensed under the Apache License, Version 2.0 (the "License");
+## you may not use this file except in compliance with the License.
+## You may obtain a copy of the License at
+## 
+##     http://www.apache.org/licenses/LICENSE-2.0
+## 
+## Unless required by applicable law or agreed to in writing, software
+## distributed under the License is distributed on an "AS IS" BASIS,
+## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+## See the License for the specific language governing permissions and
+## limitations under the License.
+
+# linear interpolation: given vectors x,y, compute a linear interpolation for y(x0)
+# assume x is ordered
+@everywhere function linear_interpolation(x0,x,y)
+  # if x0 is beyond all x's, then return the corresponding boundary value.
+  if x0>x[length(x)]
+    return y[length(y)]
+  elseif x0<x[1]
+    return y[1]
+  end
+
+  # find bracketing interval
+  i=bracket(x0,x,1,length(x))
+
+  # interpolate
+  return y[i]+(y[i+1]-y[i])*(x0-x[i])/(x[i+1]-x[i])
+end
+@everywhere function bracket(x0,x,a,b)
+  i=floor(Int64,(a+b)/2)
+  if x0<x[i]
+    return bracket(x0,x,a,i)
+  elseif x0>x[i+1]
+    return bracket(x0,x,i+1,b)
+  else
+    return i
+  end
+end
+
+
+# polynomial interpolation of a family of points
+@everywhere function poly_interpolation(x,y)
+  # init for recursion
+  rec=Array{Polynomial{Float64}}(undef,length(x))
+  for i in 1:length(x)
+    rec[i]=Polynomial([1.])
+  end
+
+  # compute \prod (x-x_i)
+  poly_interpolation_rec(rec,x,1,length(x))
+
+  # sum terms together
+  out=0.
+  for i in 1:length(y)
+    out+=rec[i]/rec[i](x[i])*y[i]
+  end
+
+  return out
+end
+# recursive helper function
+@everywhere function poly_interpolation_rec(out,x,a,b)
+  if a==b
+    return
+  end
+  # mid point
+  mid=floor(Int64,(a+b)/2)
+  # multiply left and right of mid
+  right=Polynomial([1.])
+  for i in mid+1:b
+    right*=Polynomial([-x[i],1.])
+  end
+  left=Polynomial([1.])
+  for i in a:mid
+    left*=Polynomial([-x[i],1.])
+  end
+  
+  # multiply into left and right
+  for i in a:mid
+    out[i]*=right
+  end
+  for i in mid+1:b
+    out[i]*=left
+  end
+
+  # recurse
+  poly_interpolation_rec(out,x,a,mid)
+  poly_interpolation_rec(out,x,mid+1,b)
+
+  return
+end
+## the following does the same, but has complexity N^2, the function above has N*log(N)
+#@everywhere function poly_interpolation(x,y)
+#  out=Polynomial([0.])
+#  for i in 1:length(x)
+#    prod=Polynomial([1.])
+#    for j in 1:length(x)
+#      if j!=i
+#	prod*=Polynomial([-x[j]/(x[i]-x[j]),1/(x[i]-x[j])])
+#      end
+#    end
+#    out+=prod*y[i]
+#  end
+#
+#  return out
+#end
+