diff --git a/src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs b/src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs
index b97f97e0..a31d910b 100644
--- a/src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs
+++ b/src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs
@@ -30,6 +30,8 @@
 
 using System;
 using System.Collections.Generic;
+using System.Linq;
+using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.Interpolation
 {
@@ -43,31 +45,33 @@ namespace MathNet.Numerics.Interpolation
     /// </remarks>
     public class BulirschStoerRationalInterpolation : IInterpolation
     {
-        /// <summary>
-        /// Sample Points t.
-        /// </summary>
-        IList<double> _points;
-
-        /// <summary>
-        /// Spline Values x(t).
-        /// </summary>
-        IList<double> _values;
+        readonly double[] _x;
+        readonly double[] _y;
 
         /// <summary>
         /// Initializes a new instance of the BulirschStoerRationalInterpolation class.
         /// </summary>
-        public BulirschStoerRationalInterpolation()
+        /// <param name="x">Sample Points t</param>
+        /// <param name="y">Sample Values x(t)</param>
+        public BulirschStoerRationalInterpolation(IEnumerable<double> x, IEnumerable<double> y)
         {
-        }
+            var xx = (x as double[]) ?? x.ToArray();
+            var yy = (y as double[]) ?? y.ToArray();
 
-        /// <summary>
-        /// Initializes a new instance of the BulirschStoerRationalInterpolation class.
-        /// </summary>
-        /// <param name="samplePoints">Sample Points t</param>
-        /// <param name="sampleValues">Sample Values x(t)</param>
-        public BulirschStoerRationalInterpolation(IList<double> samplePoints, IList<double> sampleValues)
-        {
-            Initialize(samplePoints, sampleValues);
+            if (xx.Length != yy.Length)
+            {
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+            }
+
+            if (xx.Length < 1)
+            {
+                throw new ArgumentOutOfRangeException("x");
+            }
+
+            Sorting.Sort(xx, yy);
+
+            _x = xx;
+            _y = yy;
         }
 
         /// <summary>
@@ -86,37 +90,6 @@ namespace MathNet.Numerics.Interpolation
             get { return false; }
         }
 
-        /// <summary>
-        /// Initialize the interpolation method with the given sample pairs.
-        /// </summary>
-        /// <param name="samplePoints">Sample Points t</param>
-        /// <param name="sampleValues">Sample Values x(t)</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 < 1)
-            {
-                throw new ArgumentOutOfRangeException("samplePoints");
-            }
-
-            if (samplePoints.Count != sampleValues.Count)
-            {
-                throw new ArgumentException(Properties.Resources.ArgumentVectorsSameLength);
-            }
-
-            _points = samplePoints;
-            _values = sampleValues;
-        }
-
         /// <summary>
         /// Interpolate at point t.
         /// </summary>
@@ -125,20 +98,20 @@ namespace MathNet.Numerics.Interpolation
         public double Interpolate(double t)
         {
             const double tiny = 1.0e-25;
-            int n = _points.Count;
+            int n = _x.Length;
 
             var c = new double[n];
             var d = new double[n];
 
             int nearestIndex = 0;
-            double nearestDistance = Math.Abs(t - _points[0]);
+            double nearestDistance = Math.Abs(t - _x[0]);
 
             for (int i = 0; i < n; i++)
             {
-                double distance = Math.Abs(t - _points[i]);
+                double distance = Math.Abs(t - _x[i]);
                 if (distance.AlmostEqual(0.0))
                 {
-                    return _values[i];
+                    return _y[i];
                 }
 
                 if (distance < nearestDistance)
@@ -147,18 +120,18 @@ namespace MathNet.Numerics.Interpolation
                     nearestDistance = distance;
                 }
 
-                c[i] = _values[i];
-                d[i] = _values[i] + tiny;
+                c[i] = _y[i];
+                d[i] = _y[i] + tiny;
             }
 
-            double x = _values[nearestIndex];
+            double x = _y[nearestIndex];
 
             for (int level = 1; level < n; level++)
             {
                 for (int i = 0; i < n - level; i++)
                 {
-                    double hp = _points[i + level] - t;
-                    double ho = (_points[i] - t)*d[i]/hp;
+                    double hp = _x[i + level] - t;
+                    double ho = (_x[i] - t)*d[i]/hp;
 
                     double den = ho - c[i + 1];
                     if (den.AlmostEqual(0.0))
diff --git a/src/Numerics/Interpolation/NevillePolynomialInterpolation.cs b/src/Numerics/Interpolation/NevillePolynomialInterpolation.cs
index d2ee878e..baf69009 100644
--- a/src/Numerics/Interpolation/NevillePolynomialInterpolation.cs
+++ b/src/Numerics/Interpolation/NevillePolynomialInterpolation.cs
@@ -30,6 +30,7 @@
 
 using System;
 using System.Collections.Generic;
+using System.Linq;
 using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.Interpolation
@@ -49,31 +50,41 @@ namespace MathNet.Numerics.Interpolation
     /// </remarks>
     public class NevillePolynomialInterpolation : IInterpolation
     {
-        /// <summary>
-        /// Sample Points t.
-        /// </summary>
-        IList<double> _points;
-
-        /// <summary>
-        /// Spline Values x(t).
-        /// </summary>
-        IList<double> _values;
+        readonly double[] _x;
+        readonly double[] _y;
 
         /// <summary>
         /// Initializes a new instance of the NevillePolynomialInterpolation class.
         /// </summary>
-        public NevillePolynomialInterpolation()
+        /// <param name="x">Sample Points t. Optimized for arrays.</param>
+        /// <param name="y">Sample Values x(t). Optimized for arrays.</param>
+        public NevillePolynomialInterpolation(IEnumerable<double> x, IEnumerable<double> y)
         {
-        }
+            var xx = (x as double[]) ?? x.ToArray();
+            var yy = (y as double[]) ?? y.ToArray();
 
-        /// <summary>
-        /// Initializes a new instance of the NevillePolynomialInterpolation class.
-        /// </summary>
-        /// <param name="samplePoints">Sample Points t</param>
-        /// <param name="sampleValues">Sample Values x(t)</param>
-        public NevillePolynomialInterpolation(IList<double> samplePoints, IList<double> sampleValues)
-        {
-            Initialize(samplePoints, sampleValues);
+            if (xx.Length != yy.Length)
+            {
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+            }
+
+            if (xx.Length < 1)
+            {
+                throw new ArgumentOutOfRangeException("x");
+            }
+
+            Sorting.Sort(xx, yy);
+
+            for (var i = 1; i < xx.Length; ++i)
+            {
+                if (xx[i] == xx[i - 1])
+                {
+                    throw new ArgumentException(Resources.Interpolation_Initialize_SamplePointsNotUnique, "x");
+                }
+            }
+
+            _x = xx;
+            _y = yy;
         }
 
         /// <summary>
@@ -92,41 +103,6 @@ namespace MathNet.Numerics.Interpolation
             get { return false; }
         }
 
-        /// <summary>
-        /// Initialize the interpolation method with the given sample pairs.
-        /// </summary>
-        /// <param name="samplePoints">Sample Points t</param>
-        /// <param name="sampleValues">Sample Values x(t)</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 < 1)
-            {
-                throw new ArgumentOutOfRangeException("samplePoints");
-            }
-
-            if (samplePoints.Count != sampleValues.Count)
-            {
-                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
-            }
-
-            for (var i = 1; i < samplePoints.Count; ++i)
-                if (samplePoints[i] == samplePoints[i - 1])
-                    throw new ArgumentException(Resources.Interpolation_Initialize_SamplePointsNotUnique, "samplePoints");
-
-            _points = samplePoints;
-            _values = sampleValues;
-        }
-
         /// <summary>
         /// Interpolate at point t.
         /// </summary>
@@ -134,16 +110,16 @@ namespace MathNet.Numerics.Interpolation
         /// <returns>Interpolated value x(t).</returns>
         public double Interpolate(double t)
         {
-            var x = new double[_values.Count];
-            _values.CopyTo(x, 0);
+            var x = new double[_y.Length];
+            _y.CopyTo(x, 0);
 
             for (int level = 1; level < x.Length; level++)
             {
                 for (int i = 0; i < x.Length - level; i++)
                 {
-                    double hp = t - _points[i + level];
-                    double ho = _points[i] - t;
-                    double den = _points[i] - _points[i + level];
+                    double hp = t - _x[i + level];
+                    double ho = _x[i] - t;
+                    double den = _x[i] - _x[i + level];
                     x[i] = ((hp*x[i]) + (ho*x[i + 1]))/den;
                 }
             }
@@ -158,17 +134,17 @@ namespace MathNet.Numerics.Interpolation
         /// <returns>Interpolated first derivative at point t.</returns>
         public double Differentiate(double t)
         {
-            var x = new double[_values.Count];
-            var dx = new double[_values.Count];
-            _values.CopyTo(x, 0);
+            var x = new double[_y.Length];
+            var dx = new double[_y.Length];
+            _y.CopyTo(x, 0);
 
             for (int level = 1; level < x.Length; level++)
             {
                 for (int i = 0; i < x.Length - level; i++)
                 {
-                    double hp = t - _points[i + level];
-                    double ho = _points[i] - t;
-                    double den = _points[i] - _points[i + level];
+                    double hp = t - _x[i + level];
+                    double ho = _x[i] - t;
+                    double den = _x[i] - _x[i + level];
                     dx[i] = ((hp*dx[i]) + x[i] + (ho*dx[i + 1]) - x[i + 1])/den;
                     x[i] = ((hp*x[i]) + (ho*x[i + 1]))/den;
                 }
@@ -184,18 +160,18 @@ namespace MathNet.Numerics.Interpolation
         /// <returns>Interpolated second derivative at point t.</returns>
         public double Differentiate2(double t)
         {
-            var x = new double[_values.Count];
-            var dx = new double[_values.Count];
-            var ddx = new double[_values.Count];
-            _values.CopyTo(x, 0);
+            var x = new double[_y.Length];
+            var dx = new double[_y.Length];
+            var ddx = new double[_y.Length];
+            _y.CopyTo(x, 0);
 
             for (int level = 1; level < x.Length; level++)
             {
                 for (int i = 0; i < x.Length - level; i++)
                 {
-                    double hp = t - _points[i + level];
-                    double ho = _points[i] - t;
-                    double den = _points[i] - _points[i + level];
+                    double hp = t - _x[i + level];
+                    double ho = _x[i] - t;
+                    double den = _x[i] - _x[i + level];
                     ddx[i] = ((hp*ddx[i]) + (ho*ddx[i + 1]) + (2*dx[i]) - (2*dx[i + 1]))/den;
                     dx[i] = ((hp*dx[i]) + x[i] + (ho*dx[i + 1]) - x[i + 1])/den;
                     x[i] = ((hp*x[i]) + (ho*x[i + 1]))/den;