diff --git a/src/Numerics.Tests/FractionalCalculusTests/RiemannLiouvilleTests.cs b/src/Numerics.Tests/FractionalCalculusTests/RiemannLiouvilleTests.cs
new file mode 100644
index 00000000..ca0593d1
--- /dev/null
+++ b/src/Numerics.Tests/FractionalCalculusTests/RiemannLiouvilleTests.cs
@@ -0,0 +1,75 @@
+using NUnit.Framework;
+using System;
+
+namespace MathNet.Numerics.Tests.FractionalCalculus
+{
+    [TestFixture, Category("DifferIntegration")]
+    public class RiemannLiouvilleTests
+    {
+        [TestCase(3.0, 3.0, 0.0)]                               // D(3) x^2 = 0
+        [TestCase(2.5, 3.0, 0.651470015870559895448512830451)]  // D(5/2) x^2 = Gamma(3)/Gamma(1/2)*x^(-1/2)
+        [TestCase(2.0, 3.0, 2.0)]                               // D(2) x^2 = 2
+        [TestCase(1.5, 3.0, 3.90882009522335937269107698271)]   // D(3/2) x^2 = Gamma(3)/Gamma(3/2)*x^(1/2)        
+        [TestCase(1.0, 3.0, 6.0)]                               // D(1) x^2 = 2 x
+        [TestCase(0.5, 3.0, 7.81764019044671874538215396541)]   // D(1/2) x^2 = Gamma(3)/Gamma(5/2)*x^(3/2)        
+        [TestCase(0.0, 3.0, 9.0)]                               // D(0) x^2 = x^2
+        [TestCase(-0.5, 3.0, 9.38116822853606249445858475850)]  // D(-1/2) x^2 = Gamma(3)/Gamma(7/2)*x^(5/2)        
+        [TestCase(-1.0, 3.0, 9.0)]                              // D(-1) x^2 = x^3/3
+        [TestCase(-1.5, 3.0, 8.04100133874519642382164407871)]  // D(-3/2) x^2 = Gamma(3)/Gamma(9/2)*x^(7/2)        
+        [TestCase(-2.0, 3.0, 6.75)]                             // D(-2) x^2 = x^4/12
+        [TestCase(-2.5, 3.0, 5.36066755916346428254776271914)]  // D(-5/2) x^2 = Gamma(3)/Gamma(11/2)*x^(9/2)        
+        [TestCase(-3.0, 3.0, 4.05)]                             // D(-3) x^2 = x^5/60
+        public void TestPolynomial(double n, double x, double expected)
+        {
+            Func<double, double> f = (t) => Math.Pow(t, 2);
+
+            Assert.LessOrEqual(
+                (expected - DifferIntegrate.DoubleExponential(f, x, n, x0: 0.0)) / expected,
+                1E-7,
+                "D({0}) x^2 (Double-Exponetial) where x = {1}", n, x);
+
+            Assert.LessOrEqual(
+                (expected - DifferIntegrate.GaussLegendre(f, x, n, x0: 0.0)) / expected,
+                1E-6,
+                "D({0}) x^2 (Gauss-Legendre) where x = {1}", n, x);
+
+            Assert.LessOrEqual(
+                (expected - DifferIntegrate.GaussKronrod(f, x, n, x0: 0.0)) / expected,
+                1E-4,
+                "D({0}) x^2 (Gauss-Kronrod) where x = {1}", n, x);
+        }
+
+        [TestCase(3.0, -1.5, 0.278538406900792139127369841989)]     // D(3) exp(π x) = exp(π x) π^3,
+        [TestCase(2.5, -1.5, 0.157148467791413402880359987678)]     // D(5/2) exp(π x) = exp(π x) π^(2.5),        
+        [TestCase(2.0, -1.5, 0.0886615285984055199686857617477)]    // D(2) exp(π x) = exp(π x) π^2,
+        [TestCase(1.5, -1.5, 0.0500219108966418944762345595265)]    // D(3/2) exp(π x) = exp(π x) π^(1.5),        
+        [TestCase(1.0, -1.5, 0.0282218410770393627236690193887)]    // D(1) exp(π x) = exp(π x) π,
+        [TestCase(0.5, -1.5, 0.0159224687642057998088504370906)]    // D(1/2) exp(π x) = exp(π x) π^(0.5),        
+        [TestCase(0.0, -1.5, 0.00898329102112942788966495190793)]   // D(0) exp(π x) = exp(π x)
+        [TestCase(-0.5, 3.5, 33631.1952282529497668583487591)]      // D(-1/2) exp(π x) = exp(π x)/π^(0.5)        
+        [TestCase(-1.0, 3.5, 18974.3700300413201713375621666)]      // D(-1) exp(π x) = exp(π x)/π
+        [TestCase(-1.5, 3.5, 10705.1419253300403750707800464)]      // D(-3/2) exp(π x) = exp(π x)/π^(1.5)        
+        [TestCase(-2.0, 3.5, 6039.72956467158140885534431539)]      // D(-2) exp(π x) = exp(π x)/π^2
+        [TestCase(-2.5, 3.5, 3407.55250783313088752769495213)]      // D(-5/2) exp(π x) = exp(π x)/π^(2.5) 
+        [TestCase(-3.0, 3.5, 1922.50563031148665828996231149)]      // D(-3) exp(π x) = exp(π x)/π^3
+        public void TestExponential(double n, double x, double expected)
+        {
+            Func<double, double> f = (t) => Math.Exp(Math.PI * t);
+
+            Assert.LessOrEqual(
+               (expected - DifferIntegrate.DoubleExponential(f, x, n, x0: double.NegativeInfinity)) / expected,
+               1E-9, 
+               "D({0}) exp(π x) (Double-Exponential) where x = {1}", n, x);
+
+            Assert.LessOrEqual(
+               (expected - DifferIntegrate.GaussLegendre(f, x, n, x0: double.NegativeInfinity)) / expected,
+               1E-6,
+               "D({0}) exp(π x) (Gauss-Legendre) where x = {1}", n, x);
+
+            Assert.LessOrEqual(
+               (expected - DifferIntegrate.GaussKronrod(f, x, n, x0: double.NegativeInfinity)) / expected,
+               1E-7,
+               "D({0}) exp(π x) (Gauss-Kronrod) where x = {1}", n, x);
+        }
+    }
+}
diff --git a/src/Numerics/DifferIntegrate.cs b/src/Numerics/DifferIntegrate.cs
new file mode 100644
index 00000000..158b2414
--- /dev/null
+++ b/src/Numerics/DifferIntegrate.cs
@@ -0,0 +1,134 @@
+using System;
+
+namespace MathNet.Numerics
+{
+    public static class DifferIntegrate
+    {
+        /// <summary>
+        /// Evaluates the Riemann-Liouville fractional derivative that uses the double exponential integration.
+        /// </summary>
+        /// <remarks>
+        /// <para>order = 1.0 : normal derivative</para>
+        /// <para>order = 0.5 : semi-derivative</para>
+        /// <para>order = -0.5 : semi-integral</para>
+        /// <para>order = -1.0 : normal integral</para>
+        /// </remarks>
+        /// <param name="f">The analytic smooth function to differintegrate.</param>
+        /// <param name="x">The evaluation point.</param>
+        /// <param name="order">The order of fractional derivative.</param>
+        /// <param name="x0">The reference point of integration.</param>
+        /// <param name="targetAbsoluteError">The expected relative accuracy of the Double-Exponential integration.</param>
+        /// <returns>Approximation of the differintegral of order n at x.</returns>
+        public static double DoubleExponential(Func<double, double> f, double x, double order, double x0 = 0, double targetAbsoluteError = 1E-10)
+        {
+            // The Riemann–Liouville fractional derivative of f(x) of order n is defined as
+            //    \,_{x_0}{\mathbb{D}}^n_xf(x) = \frac{1}{\Gamma(m-n)} \frac{d^m}{dx^m} \int_{x_0}^{x} (x-t)^{m-n-1} f(t) dt
+            //    where m is the smallest interger greater than n.
+            // see https://en.wikipedia.org/wiki/Differintegral
+
+            if (Math.Abs(order) < double.Epsilon)
+            {
+                return f(x);
+            }
+            else if (order > 0 && Math.Abs(order - (int)order) < double.Epsilon)
+            {
+                return Differentiate.Derivative(f, x, (int)order);
+            }
+            else
+            {
+                int m = (int)Math.Ceiling(order) + 1;
+                if (m < 1) m = 1;
+                double r = m - order - 1;
+                Func<double, double> g = (v) => Integrate.DoubleExponential((t) => Math.Pow(v - t, r) * f(t), x0, v, targetAbsoluteError: targetAbsoluteError);
+                double numerator = Differentiate.Derivative(g, x, m);
+                double denominator = SpecialFunctions.Gamma(m - order);
+                return numerator / denominator;
+            }
+        }
+
+        /// <summary>
+        /// Evaluates the Riemann-Liouville fractional derivative that uses the Gauss-Legendre integration.
+        /// </summary>
+        /// <remarks>
+        /// <para>order = 1.0 : normal derivative</para>
+        /// <para>order = 0.5 : semi-derivative</para>
+        /// <para>order = -0.5 : semi-integral</para>
+        /// <para>order = -1.0 : normal integral</para>
+        /// </remarks>
+        /// <param name="f">The analytic smooth function to differintegrate.</param>
+        /// <param name="x">The evaluation point.</param>
+        /// <param name="order">The order of fractional derivative.</param>
+        /// <param name="x0">The reference point of integration.</param>
+        /// <param name="gaussLegendrePoints">The number of Gauss-Legendre points.</param>
+        /// <returns>Approximation of the differintegral of order n at x.</returns>
+        public static double GaussLegendre(Func<double, double> f, double x, double order, double x0 = 0, int gaussLegendrePoints = 128)
+        {
+            // The Riemann–Liouville fractional derivative of f(x) of order n is defined as
+            //    \,_{x_0}{\mathbb{D}}^n_xf(x) = \frac{1}{\Gamma(m-n)} \frac{d^m}{dx^m} \int_{x_0}^{x} (x-t)^{m-n-1} f(t) dt
+            //    where m is the smallest interger greater than n.
+            // see https://en.wikipedia.org/wiki/Differintegral
+
+            if (Math.Abs(order) < double.Epsilon)
+            {
+                return f(x);
+            }
+            else if (order > 0 && Math.Abs(order - (int)order) < double.Epsilon) 
+            {
+                return Differentiate.Derivative(f, x, (int)order);
+            }
+            else
+            {
+                int m = (int)Math.Ceiling(order) + 1;
+                if (m < 1) m = 1;
+                double r = m - order - 1;
+                Func<double, double> g = (v) => Integrate.GaussLegendre((t) => Math.Pow(v - t, r) * f(t), x0, v, order: gaussLegendrePoints);
+                double numerator = Differentiate.Derivative(g, x, m);
+                double denominator = SpecialFunctions.Gamma(m - order);
+                return numerator / denominator;
+            }
+        }
+
+        /// <summary>
+        /// Evaluates the Riemann-Liouville fractional derivative that uses the Gauss-Kronrod integration.
+        /// </summary>
+        /// <remarks>
+        /// <para>order = 1.0 : normal derivative</para>
+        /// <para>order = 0.5 : semi-derivative</para>
+        /// <para>order = -0.5 : semi-integral</para>
+        /// <para>order = -1.0 : normal integral</para>
+        /// </remarks>
+        /// <param name="f">The analytic smooth function to differintegrate.</param>
+        /// <param name="x">The evaluation point.</param>
+        /// <param name="order">The order of fractional derivative.</param>
+        /// <param name="x0">The reference point of integration.</param>
+        /// <param name="targetRelativeError">The expected relative accuracy of the Gauss-Kronrod integration.</param>
+        /// <param name="gaussKronrodPoints">The number of Gauss-Kronrod points. Pre-computed for 15, 21, 31, 41, 51 and 61 points.</param>
+        /// <returns>Approximation of the differintegral of order n at x.</returns>
+        public static double GaussKronrod(Func<double, double> f, double x, double order, double x0 = 0, double targetRelativeError = 1E-10, int gaussKronrodPoints = 15)
+        {
+            // The Riemann–Liouville fractional derivative of f(x) of order n is defined as
+            //    \,_{x_0}{\mathbb{D}}^n_xf(x) = \frac{1}{\Gamma(m-n)} \frac{d^m}{dx^m} \int_{x_0}^{x} (x-t)^{m-n-1} f(t) dt
+            //    where m is the smallest interger greater than n.
+            // see https://en.wikipedia.org/wiki/Differintegral
+
+            if (Math.Abs(order) < double.Epsilon)
+            {
+                return f(x);
+            }
+            else if (order > 0 && Math.Abs(order - (int)order) < double.Epsilon) 
+            {
+                return Differentiate.Derivative(f, x, (int)order);
+            }
+            else
+            {
+                int m = (int)Math.Ceiling(order) + 1;
+                if (m < 1) m = 1;
+                double r = m - order - 1;
+                Func<double, double> g = (v) => Integrate.GaussKronrod((t) => Math.Pow(v - t, r) * f(t), x0, v, targetRelativeError: targetRelativeError, order: gaussKronrodPoints);
+                double numerator = Differentiate.Derivative(g, x, m);
+                double denominator = SpecialFunctions.Gamma(m - order);
+                return numerator / denominator;
+            }
+        }
+    }
+}