using System;

namespace IStation.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;
            }
        }
    }
}