interpolation: rational pole-free, linear spline

Signed-off-by: Christoph Ruegg <git@cdrnet.ch>

src/Managed/Interpolation/Algorithms/BarycentricInterpolation.cs

@@ -0,0 +1,226 @@
+// <copyright file="BarycentricInterpolation.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Interpolation.Algorithms
+{
+    using System;
+    using System.Collections.Generic;
+
+    /// <summary>
+    /// Barycentric Interpolation Algorithm.
+    /// </summary>
+    /// <remarks>
+    /// This algorithm neither supports differentiation nor integration.
+    /// </remarks>
+    public class BarycentricInterpolation : IInterpolation
+    {
+        /// <summary>
+        /// Sample Points t.
+        /// </summary>
+        private IList<double> points;
+
+        /// <summary>
+        /// Sample Values x(t).
+        /// </summary>
+        private IList<double> values;
+
+        /// <summary>
+        /// Barycentric Weights w(t).
+        /// </summary>
+        private IList<double> weights;
+
+        /// <summary>
+        /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative).
+        /// </summary>
+        /// <seealso cref="IInterpolation.Differentiate(double)"/>
+        /// <seealso cref="IInterpolation.Differentiate(double, out double, out double)"/>
+        bool IInterpolation.SupportsDifferentiation
+        {
+            get { return false; }
+        }
+
+        /// <summary>
+        /// Gets a value indicating whether the algorithm supports integration (interpolated quadrature).
+        /// </summary>
+        /// <seealso cref="IInterpolation.Integrate"/>
+        bool IInterpolation.SupportsIntegration
+        {
+            get { return false; }
+        }
+
+        /// <summary>
+        /// Initialize the interpolation method with the given sample set.
+        /// </summary>
+        /// <param name="samplePoints">Sample Points t</param>
+        /// <param name="sampleValues">Sample Values x(t)</param>
+        /// <param name="barycentricWeights">Barycentric weights w(t)</param>
+        public
+        void
+        Initialize(
+            IList<double> samplePoints,
+            IList<double> sampleValues,
+            IList<double> barycentricWeights)
+        {
+            if (null == samplePoints)
+            {
+                throw new ArgumentNullException("samplePoints");
+            }
+
+            if (null == sampleValues)
+            {
+                throw new ArgumentNullException("sampleValues");
+            }
+
+            if (null == barycentricWeights)
+            {
+                throw new ArgumentNullException("barycentricWeights");
+            }
+
+            if (samplePoints.Count < 1)
+            {
+                throw new ArgumentOutOfRangeException("samplePoints");
+            }
+
+            if (samplePoints.Count != sampleValues.Count)
+            {
+                throw new ArgumentException(Properties.Resources.ArgumentVectorsSameLengths);
+            }
+
+            if (samplePoints.Count != barycentricWeights.Count)
+            {
+                throw new ArgumentException(Properties.Resources.ArgumentVectorsSameLengths);
+            }
+
+            this.points = samplePoints;
+            this.values = sampleValues;
+            this.weights = barycentricWeights;
+        }
+
+        /// <summary>
+        /// Interpolate at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated value x(t).</returns>
+        public
+        double
+        Interpolate(double t)
+        {
+            // trivial case: only one sample?
+            if (this.points.Count == 1)
+            {
+                return this.values[0];
+            }
+
+            // evaluate closest point and offset from that point
+            int closestPoint = 0;
+            double offset = t - this.points[0];
+            for (int i = 1; i < this.points.Count; i++)
+            {
+                if (Math.Abs(t - this.points[i]) < Math.Abs(offset))
+                {
+                    offset = t - this.points[i];
+                    closestPoint = i;
+                }
+            }
+
+            // trivial case: on a known sample point?
+            // TODO: Number.AlmostZero(offset) instead of ==
+            if (offset == 0.0)
+            {
+                return this.values[closestPoint];
+            }
+
+            if (Math.Abs(offset) > 1e-150)
+            {
+                // no need to guard against overflow, so use fast formula
+                closestPoint = -1;
+                offset = 1.0;
+            }
+
+            double s1 = 0.0;
+            double s2 = 0.0;
+            for (int i = 0; i < this.points.Count; i++)
+            {
+                if (i != closestPoint)
+                {
+                    double v = offset * this.weights[i] / (t - this.points[i]);
+                    s1 = s1 + (v * this.values[i]);
+                    s2 = s2 + v;
+                }
+                else
+                {
+                    double v = this.weights[i];
+                    s1 = s1 + (v * this.values[i]);
+                    s2 = s2 + v;
+                }
+            }
+
+            return s1 / s2;
+        }
+
+        /// <summary>
+        /// Differentiate at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated first derivative at point t.</returns>
+        /// <seealso cref="IInterpolation.SupportsDifferentiation"/>
+        /// <seealso cref="IInterpolation.Differentiate(double, out double, out double)"/>
+        double IInterpolation.Differentiate(double t)
+        {
+            throw new NotSupportedException();
+        }
+
+        /// <summary>
+        /// Differentiate at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <param name="interpolatedValue">Interpolated value x(t)</param>
+        /// <param name="secondDerivative">Interpolated second derivative at point t.</param>
+        /// <returns>Interpolated first derivative at point t.</returns>
+        /// <seealso cref="IInterpolation.SupportsDifferentiation"/>
+        /// <seealso cref="IInterpolation.Differentiate(double)"/>
+        double IInterpolation.Differentiate(
+            double t,
+            out double interpolatedValue,
+            out double secondDerivative)
+        {
+            throw new NotSupportedException();
+        }
+
+        /// <summary>
+        /// Integrate up to point t.
+        /// </summary>
+        /// <param name="t">Right bound of the integration interval [a,t].</param>
+        /// <returns>Interpolated definite integral over the interval [a,t].</returns>
+        /// <seealso cref="IInterpolation.SupportsIntegration"/>
+        double IInterpolation.Integrate(double t)
+        {
+            throw new NotSupportedException();
+        }
+    }
+}

src/Managed/Interpolation/Algorithms/LinearSplineInterpolation.cs

@@ -0,0 +1,177 @@
+// <copyright file="LinearSplineInterpolation.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Interpolation.Algorithms
+{
+    using System;
+    using System.Collections.Generic;
+
+    /// <summary>
+    /// Linear Spline Interpolation Algorithm.
+    /// </summary>
+    /// <remarks>
+    /// This algorithm supports both differentiation and integration.
+    /// </remarks>
+    public class LinearSplineInterpolation : IInterpolation
+    {
+        /// <summary>
+        /// Internal Spline Interpolation
+        /// </summary>
+        private readonly SplineInterpolation spline;
+
+        /// <summary>
+        /// Initializes a new instance of the LinearSplineInterpolation class.
+        /// </summary>
+        public
+        LinearSplineInterpolation()
+        {
+            this.spline = new SplineInterpolation();
+        }
+
+        /// <summary>
+        /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative).
+        /// </summary>
+        /// <seealso cref="Differentiate(double)"/>
+        /// <seealso cref="Differentiate(double, out double, out double)"/>
+        bool IInterpolation.SupportsDifferentiation
+        {
+            get { return true; }
+        }
+
+        /// <summary>
+        /// Gets a value indicating whether the algorithm supports integration (interpolated quadrature).
+        /// </summary>
+        /// <seealso cref="Integrate"/>
+        bool IInterpolation.SupportsIntegration
+        {
+            get { return true; }
+        }
+
+        /// <summary>
+        /// Initialize the interpolation method with the given spline coefficients.
+        /// </summary>
+        /// <param name="samplePoints">Sample Points t</param>
+        /// <param name="sampleValues">Sample Values x(samplePoints)</param>
+        public
+        void
+        Initialize(
+            IList<double> samplePoints,
+            IList<double> sampleValues)
+        {
+            if (null == samplePoints)
+            {
+                throw new ArgumentNullException("samplePoints");
+            }
+
+            if (null == sampleValues)
+            {
+                throw new ArgumentNullException("sampleValues");
+            }
+
+            if (samplePoints.Count < 2)
+            {
+                throw new ArgumentOutOfRangeException("samplePoints");
+            }
+
+            if (samplePoints.Count != sampleValues.Count)
+            {
+                throw new ArgumentException(Properties.Resources.ArgumentVectorsSameLengths);
+            }
+
+            double[] coefficients = new double[4 * (samplePoints.Count - 1)];
+            double[] sortedPoints = new double[samplePoints.Count];
+            samplePoints.CopyTo(sortedPoints, 0);
+            double[] sortedValues = new double[sampleValues.Count];
+            sampleValues.CopyTo(sortedValues, 0);
+
+            // TODO: Sorting.Sort(sortedPoints, sortedValues);
+
+            for (int i = 0, j = 0; i < sortedPoints.Length - 1; i++, j += 4)
+            {
+                coefficients[j] = sortedValues[i];
+                coefficients[j + 1] = (sortedValues[i + 1] - sortedValues[i]) / (sortedPoints[i + 1] - sortedPoints[i]);
+                coefficients[j + 2] = 0;
+                coefficients[j + 3] = 0;
+            }
+
+            this.spline.Initialize(sortedPoints, coefficients);
+        }
+
+        /// <summary>
+        /// Interpolate at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated value x(t).</returns>
+        public
+        double
+        Interpolate(double t)
+        {
+            return this.spline.Interpolate(t);
+        }
+
+        /// <summary>
+        /// Differentiate at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated first derivative at point t.</returns>
+        /// <seealso cref="IInterpolation.SupportsDifferentiation"/>
+        /// <seealso cref="Differentiate(double, out double, out double)"/>
+        public double Differentiate(double t)
+        {
+            return this.spline.Differentiate(t);
+        }
+
+        /// <summary>
+        /// Differentiate at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <param name="interpolatedValue">Interpolated value x(t)</param>
+        /// <param name="secondDerivative">Interpolated second derivative at point t.</param>
+        /// <returns>Interpolated first derivative at point t.</returns>
+        /// <seealso cref="IInterpolation.SupportsDifferentiation"/>
+        /// <seealso cref="Differentiate(double)"/>
+        public double Differentiate(
+            double t,
+            out double interpolatedValue,
+            out double secondDerivative)
+        {
+            return this.spline.Differentiate(t, out interpolatedValue, out secondDerivative);
+        }
+
+        /// <summary>
+        /// Integrate up to point t.
+        /// </summary>
+        /// <param name="t">Right bound of the integration interval [a,t].</param>
+        /// <returns>Interpolated definite integral over the interval [a,t].</returns>
+        /// <seealso cref="IInterpolation.SupportsIntegration"/>
+        public double Integrate(double t)
+        {
+            return this.spline.Integrate(t);
+        }
+    }
+}

src/Managed/Interpolation/Algorithms/RationalPoleFreeInterpolation.cs

@@ -0,0 +1,248 @@
+// <copyright file="RationalPoleFreeInterpolation.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Interpolation.Algorithms
+{
+    using System;
+    using System.Collections.Generic;
+
+    /// <summary>
+    /// Barycentric Rational Interpolation without poles, using Floater and Hormann's Algorithm.
+    /// </summary>
+    /// <remarks>
+    /// This algorithm neither supports differentiation nor integration.
+    /// </remarks>
+    public class RationalPoleFreeInterpolation : IInterpolation
+    {
+        /// <summary>
+        /// Internal Barycentric Interpolation
+        /// </summary>
+        private readonly BarycentricInterpolation barycentric;
+
+        /// <summary>
+        /// Initializes a new instance of the RationalPoleFreeInterpolation class.
+        /// </summary>
+        public
+        RationalPoleFreeInterpolation()
+        {
+            this.barycentric = new BarycentricInterpolation();
+        }
+
+        /// <summary>
+        /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative).
+        /// </summary>
+        /// <seealso cref="IInterpolation.Differentiate(double)"/>
+        /// <seealso cref="IInterpolation.Differentiate(double, out double, out double)"/>
+        bool IInterpolation.SupportsDifferentiation
+        {
+            get { return false; }
+        }
+
+        /// <summary>
+        /// Gets a value indicating whether the algorithm supports integration (interpolated quadrature).
+        /// </summary>
+        /// <seealso cref="IInterpolation.Integrate"/>
+        bool IInterpolation.SupportsIntegration
+        {
+            get { return false; }
+        }
+
+        /// <summary>
+        /// Initialize the interpolation method with the given sample set.
+        /// </summary>
+        /// <remarks>
+        /// The interpolation scheme order will be set to 3.
+        /// </remarks>
+        /// <param name="samplePoints">Sample Points t</param>
+        /// <param name="sampleValues">Sample Values x(t)</param>
+        public
+        void
+        Initialize(
+            IList<double> samplePoints,
+            IList<double> sampleValues)
+        {
+            this.Initialize(samplePoints, sampleValues, Math.Min(3, samplePoints.Count - 1));
+        }
+
+        /// <summary>
+        /// Initialize the interpolation method with the given sample set.
+        /// </summary>
+        /// <param name="samplePoints">Sample Points t</param>
+        /// <param name="sampleValues">Sample Values x(t)</param>
+        /// <param name="order">
+        /// Order of the interpolation scheme, 0 &lt;= order &lt;= N.
+        /// In most cases a value between 3 and 8 gives good results.
+        /// </param>
+        public
+        void
+        Initialize(
+            IList<double> samplePoints,
+            IList<double> sampleValues,
+            int order)
+        {
+            if (null == samplePoints)
+            {
+                throw new ArgumentNullException("samplePoints");
+            }
+
+            if (null == sampleValues)
+            {
+                throw new ArgumentNullException("sampleValues");
+            }
+
+            if (samplePoints.Count < 1)
+            {
+                throw new ArgumentOutOfRangeException("samplePoints");
+            }
+
+            if (samplePoints.Count != sampleValues.Count)
+            {
+                throw new ArgumentException(Properties.Resources.ArgumentVectorsSameLengths);
+            }
+
+            if (0 > order || samplePoints.Count <= order)
+            {
+                throw new ArgumentOutOfRangeException("order");
+            }
+
+            double[] sortedWeights = new double[sampleValues.Count];
+            double[] sortedPoints = new double[samplePoints.Count];
+            samplePoints.CopyTo(sortedPoints, 0);
+
+            // order: odd -> negative, even -> positive
+            double sign = ((order & 0x1) == 0x1) ? -1.0 : 1.0;
+
+            // init permutation vector
+            int[] perm = new int[sortedWeights.Length];
+            for (int i = 0; i < perm.Length; i++)
+            {
+                perm[i] = i;
+            }
+
+            // sort and update permutation vector
+            for (int i = 0; i < perm.Length - 1; i++)
+            {
+                for (int j = i + 1; j < perm.Length; j++)
+                {
+                    if (sortedPoints[j] < sortedPoints[i])
+                    {
+                        double s = sortedPoints[i];
+                        sortedPoints[i] = sortedPoints[j];
+                        sortedPoints[j] = s;
+                        int k = perm[i];
+                        perm[i] = perm[j];
+                        perm[j] = k;
+                    }
+                }
+            }
+
+            // compute barycentric weights
+            for (int k = 0; k < sortedWeights.Length; k++)
+            {
+                double s = 0;
+                for (int i = Math.Max(k - order, 0); i <= Math.Min(k, sortedWeights.Length - 1 - order); i++)
+                {
+                    double v = 1;
+                    for (int j = i; j <= i + order; j++)
+                    {
+                        if (j != k)
+                        {
+                            v = v / Math.Abs(sortedPoints[k] - sortedPoints[j]);
+                        }
+                    }
+
+                    s = s + v;
+                }
+
+                sortedWeights[k] = sign * s;
+                sign = -sign;
+            }
+
+            // reorder back to original order, based on the permutation vector.
+            double[] weights = new double[sortedWeights.Length];
+            for (int i = 0; i < weights.Length; i++)
+            {
+                weights[perm[i]] = sortedWeights[i];
+            }
+
+            this.barycentric.Initialize(samplePoints, sampleValues, weights);
+        }
+
+        /// <summary>
+        /// Interpolate at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated value x(t).</returns>
+        public
+        double
+        Interpolate(double t)
+        {
+            return this.barycentric.Interpolate(t);
+        }
+
+        /// <summary>
+        /// Differentiate at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated first derivative at point t.</returns>
+        /// <seealso cref="IInterpolation.SupportsDifferentiation"/>
+        /// <seealso cref="IInterpolation.Differentiate(double, out double, out double)"/>
+        double IInterpolation.Differentiate(double t)
+        {
+            throw new NotSupportedException();
+        }
+
+        /// <summary>
+        /// Differentiate at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <param name="interpolatedValue">Interpolated value x(t)</param>
+        /// <param name="secondDerivative">Interpolated second derivative at point t.</param>
+        /// <returns>Interpolated first derivative at point t.</returns>
+        /// <seealso cref="IInterpolation.SupportsDifferentiation"/>
+        /// <seealso cref="IInterpolation.Differentiate(double)"/>
+        double IInterpolation.Differentiate(
+            double t,
+            out double interpolatedValue,
+            out double secondDerivative)
+        {
+            throw new NotSupportedException();
+        }
+
+        /// <summary>
+        /// Integrate up to point t.
+        /// </summary>
+        /// <param name="t">Right bound of the integration interval [a,t].</param>
+        /// <returns>Interpolated definite integral over the interval [a,t].</returns>
+        /// <seealso cref="IInterpolation.SupportsIntegration"/>
+        double IInterpolation.Integrate(double t)
+        {
+            throw new NotSupportedException();
+        }
+    }
+}

src/Managed/Interpolation/Algorithms/SplineInterpolation.cs

@@ -0,0 +1,249 @@
+// <copyright file="SplineInterpolation.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Interpolation.Algorithms
+{
+    using System;
+    using System.Collections.Generic;
+
+    /// <summary>
+    /// Third-Degree Spline Interpolation Algorithm.
+    /// </summary>
+    /// <remarks>
+    /// This algorithm supports both differentiation and integration.
+    /// </remarks>
+    public class SplineInterpolation : IInterpolation
+    {
+        /// <summary>
+        /// Sample Points t.
+        /// </summary>
+        private IList<double> points;
+
+        /// <summary>
+        /// Spline Coefficients c(t).
+        /// </summary>
+        private IList<double> coefficients;
+
+        /// <summary>
+        /// Number of samples.
+        /// </summary>
+        private int sampleCount;
+
+        /// <summary>
+        /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative).
+        /// </summary>
+        /// <seealso cref="Differentiate(double)"/>
+        /// <seealso cref="Differentiate(double, out double, out double)"/>
+        bool IInterpolation.SupportsDifferentiation
+        {
+            get { return true; }
+        }
+
+        /// <summary>
+        /// Gets a value indicating whether the algorithm supports integration (interpolated quadrature).
+        /// </summary>
+        /// <seealso cref="Integrate"/>
+        bool IInterpolation.SupportsIntegration
+        {
+            get { return true; }
+        }
+
+        /// <summary>
+        /// Initialize the interpolation method with the given spline coefficients.
+        /// </summary>
+        /// <param name="samplePoints">Sample Points t (length: N)</param>
+        /// <param name="splineCoefficients">Spline Coefficients (length: 4*(N-1))</param>
+        public
+        void
+        Initialize(
+            IList<double> samplePoints,
+            IList<double> splineCoefficients)
+        {
+            if (null == samplePoints)
+            {
+                throw new ArgumentNullException("samplePoints");
+            }
+
+            if (null == splineCoefficients)
+            {
+                throw new ArgumentNullException("splineCoefficients");
+            }
+
+            if (samplePoints.Count < 1)
+            {
+                throw new ArgumentOutOfRangeException("samplePoints");
+            }
+
+            if (splineCoefficients.Count != 4 * (samplePoints.Count - 1))
+            {
+                throw new ArgumentOutOfRangeException("splineCoefficients");
+            }
+
+            this.points = samplePoints;
+            this.coefficients = splineCoefficients;
+            this.sampleCount = samplePoints.Count;
+        }
+
+        /// <summary>
+        /// Interpolate at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated value x(t).</returns>
+        public
+        double
+        Interpolate(double t)
+        {
+            // Binary search in the [ t[0], ..., t[n-2] ] (t[n-1] is not included)
+            int low = 0;
+            int high = this.sampleCount - 1;
+            while (low != high - 1)
+            {
+                int middle = (low + high) / 2;
+                if (this.points[middle] > t)
+                {
+                    high = middle;
+                }
+                else
+                {
+                    low = middle;
+                }
+            }
+
+            // Interpolation
+            t = t - this.points[low];
+            int k = low << 2;
+            return this.coefficients[k] + (t * (this.coefficients[k + 1] + (t * (this.coefficients[k + 2] + (t * this.coefficients[k + 3])))));
+        }
+
+        /// <summary>
+        /// Differentiate at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated first derivative at point t.</returns>
+        /// <seealso cref="IInterpolation.SupportsDifferentiation"/>
+        /// <seealso cref="Differentiate(double, out double, out double)"/>
+        public double Differentiate(double t)
+        {
+            // Binary search in the [ t[0], ..., t[n-2] ] (t[n-1] is not included)
+            int low = 0;
+            int high = this.sampleCount - 1;
+            while (low != high - 1)
+            {
+                int middle = (low + high) / 2;
+                if (this.points[middle] > t)
+                {
+                    high = middle;
+                }
+                else
+                {
+                    low = middle;
+                }
+            }
+
+            // Differentiation
+            t = t - this.points[low];
+            int k = low << 2;
+            return this.coefficients[k + 1] + (2 * t * this.coefficients[k + 2]) + (3 * t * t * this.coefficients[k + 3]);
+        }
+
+        /// <summary>
+        /// Differentiate at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <param name="interpolatedValue">Interpolated value x(t)</param>
+        /// <param name="secondDerivative">Interpolated second derivative at point t.</param>
+        /// <returns>Interpolated first derivative at point t.</returns>
+        /// <seealso cref="IInterpolation.SupportsDifferentiation"/>
+        /// <seealso cref="Differentiate(double)"/>
+        public double Differentiate(
+            double t,
+            out double interpolatedValue,
+            out double secondDerivative)
+        {
+            // Binary search in the [ t[0], ..., t[n-2] ] (t[n-1] is not included)
+            int low = 0;
+            int high = this.sampleCount - 1;
+            while (low != high - 1)
+            {
+                int middle = (low + high) / 2;
+                if (this.points[middle] > t)
+                {
+                    high = middle;
+                }
+                else
+                {
+                    low = middle;
+                }
+            }
+
+            // Differentiation
+            t = t - this.points[low];
+            int k = low << 2;
+            interpolatedValue = this.coefficients[k] + (t * (this.coefficients[k + 1] + (t * (this.coefficients[k + 2] + (t * this.coefficients[k + 3])))));
+            secondDerivative = (2 * this.coefficients[k + 2]) + (6 * t * this.coefficients[k + 3]);
+            return this.coefficients[k + 1] + (2 * t * this.coefficients[k + 2]) + (3 * t * t * this.coefficients[k + 3]);
+        }
+
+        /// <summary>
+        /// Integrate up to point t.
+        /// </summary>
+        /// <param name="t">Right bound of the integration interval [a,t].</param>
+        /// <returns>Interpolated definite integral over the interval [a,t].</returns>
+        /// <seealso cref="IInterpolation.SupportsIntegration"/>
+        public double Integrate(double t)
+        {
+            // Binary search in the [ t[0], ..., t[n-2] ] (t[n-1] is not included)
+            int low = 0;
+            int high = this.sampleCount - 1;
+            while (low != high - 1)
+            {
+                int middle = (low + high) / 2;
+                if (this.points[middle] > t)
+                {
+                    high = middle;
+                }
+                else
+                {
+                    low = middle;
+                }
+            }
+
+            // Integration
+            double result = 0;
+            for (int i = 0, j = 0; i < low; i++, j += 4)
+            {
+                double w = this.points[i + 1] - this.points[i];
+                result += w * (this.coefficients[j] + ((w * (this.coefficients[j + 1] * 0.5)) + (w * ((this.coefficients[j + 2] / 3) + (w * this.coefficients[j + 3] * 0.25)))));
+            }
+
+            t = t - this.points[low];
+            int k = low << 2;
+            return result + (t * (this.coefficients[k] + ((t * (this.coefficients[k + 1] * 0.5)) + (t * (this.coefficients[k + 2] / 3)) + (t * this.coefficients[k + 3] * 0.25))));
+        }
+    }
+}

src/Managed/Interpolation/Interpolation.cs

@@ -30,6 +30,7 @@ namespace MathNet.Numerics.Interpolation
 {
     using System;
     using System.Collections.Generic;
+    using Algorithms;
 
     /// <summary>
     /// Interpolation Factory.
@@ -37,7 +38,7 @@ namespace MathNet.Numerics.Interpolation
     public static class Interpolation
     {
         /// <summary>
-        /// Create a rational pole-free interpolation based on arbitrary points. This is the default interpolation scheme.
+        /// Creates an interpolation based on arbitrary points.
         /// </summary>
         /// <param name="points">The sample points t. Supports both lists and arrays.</param>
         /// <param name="values">The sample point values x(t). Supports both lists and arrays.</param>
@@ -50,7 +51,7 @@ namespace MathNet.Numerics.Interpolation
             IList<double> points,
             IList<double> values)
         {
-            throw new NotImplementedException();
+            return CreateRationalPoleFree(points, values);
         }
 
         /// <summary>
@@ -67,7 +68,9 @@ namespace MathNet.Numerics.Interpolation
             IList<double> points,
             IList<double> values)
         {
-            throw new NotImplementedException();
+            LinearSplineInterpolation method = new LinearSplineInterpolation();
+            method.Initialize(points, values);
+            return method;
         }
 
         /// <summary>
@@ -84,7 +87,9 @@ namespace MathNet.Numerics.Interpolation
             IList<double> points,
             IList<double> values)
         {
-            throw new NotImplementedException();
+            RationalPoleFreeInterpolation method = new RationalPoleFreeInterpolation();
+            method.Initialize(points, values);
+            return method;
         }
     }
 }

src/Managed/Managed.csproj

@@ -54,6 +54,10 @@
     <Compile Include="Distributions\IContinuousDistribution.cs" />
     <Compile Include="Distributions\IDiscreteDistribution.cs" />
     <Compile Include="Distributions\IDistribution.cs" />
+    <Compile Include="Interpolation\Algorithms\BarycentricInterpolation.cs" />
+    <Compile Include="Interpolation\Algorithms\LinearSplineInterpolation.cs" />
+    <Compile Include="Interpolation\Algorithms\RationalPoleFreeInterpolation.cs" />
+    <Compile Include="Interpolation\Algorithms\SplineInterpolation.cs" />
     <Compile Include="Interpolation\IInterpolation.cs" />
     <Compile Include="Interpolation\Interpolation.cs" />
     <Compile Include="Precision.cs" />

src/Native/Native.csproj

@@ -71,6 +71,18 @@
     <Compile Include="..\Managed\Distributions\IDistribution.cs">
       <Link>Distributions\IDistribution.cs</Link>
     </Compile>
+    <Compile Include="..\Managed\Interpolation\Algorithms\BarycentricInterpolation.cs">
+      <Link>Interpolation\Algorithms\BarycentricInterpolation.cs</Link>
+    </Compile>
+    <Compile Include="..\Managed\Interpolation\Algorithms\LinearSplineInterpolation.cs">
+      <Link>Interpolation\Algorithms\LinearSplineInterpolation.cs</Link>
+    </Compile>
+    <Compile Include="..\Managed\Interpolation\Algorithms\RationalPoleFreeInterpolation.cs">
+      <Link>Interpolation\Algorithms\RationalPoleFreeInterpolation.cs</Link>
+    </Compile>
+    <Compile Include="..\Managed\Interpolation\Algorithms\SplineInterpolation.cs">
+      <Link>Interpolation\Algorithms\SplineInterpolation.cs</Link>
+    </Compile>
     <Compile Include="..\Managed\Interpolation\IInterpolation.cs">
       <Link>Interpolation\IInterpolation.cs</Link>
     </Compile>