// <copyright file="FindRoots.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
//
// Copyright (c) 2009-2014 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>

using System;
using IStation.Numerics.RootFinding;
using Complex = System.Numerics.Complex;

namespace IStation.Numerics
{
    public static class FindRoots
    {
        /// <summary>Find a solution of the equation f(x)=0.</summary>
        /// <param name="f">The function to find roots from.</param>
        /// <param name="lowerBound">The low value of the range where the root is supposed to be.</param>
        /// <param name="upperBound">The high value of the range where the root is supposed to be.</param>
        /// <param name="accuracy">Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached. Example: 1e-14.</param>
        /// <param name="maxIterations">Maximum number of iterations. Example: 100.</param>
        public static double OfFunction(Func<double, double> f, double lowerBound, double upperBound, double accuracy = 1e-8, int maxIterations = 100)
        {
            if (!ZeroCrossingBracketing.ExpandReduce(f, ref lowerBound, ref upperBound, 1.6, maxIterations, maxIterations*10))
            {
                throw new NonConvergenceException("The algorithm has failed, exceeded the number of iterations allowed or there is no root within the provided bounds.");
            }

            if (Brent.TryFindRoot(f, lowerBound, upperBound, accuracy, maxIterations, out var root))
            {
                return root;
            }

            if (Bisection.TryFindRoot(f, lowerBound, upperBound, accuracy, maxIterations, out root))
            {
                return root;
            }

            throw new NonConvergenceException("The algorithm has failed, exceeded the number of iterations allowed or there is no root within the provided bounds.");
        }

        /// <summary>Find a solution of the equation f(x)=0.</summary>
        /// <param name="f">The function to find roots from.</param>
        /// <param name="df">The first derivative of the function to find roots from.</param>
        /// <param name="lowerBound">The low value of the range where the root is supposed to be.</param>
        /// <param name="upperBound">The high value of the range where the root is supposed to be.</param>
        /// <param name="accuracy">Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached. Example: 1e-14.</param>
        /// <param name="maxIterations">Maximum number of iterations. Example: 100.</param>
        public static double OfFunctionDerivative(Func<double, double> f, Func<double, double> df, double lowerBound, double upperBound, double accuracy = 1e-8, int maxIterations = 100)
        {
            double root;

            if (RobustNewtonRaphson.TryFindRoot(f, df, lowerBound, upperBound, accuracy, maxIterations, 20, out root))
            {
                return root;
            }

            return OfFunction(f, lowerBound, upperBound, accuracy, maxIterations);
        }

        /// <summary>
        /// Find both complex roots of the quadratic equation c + b*x + a*x^2 = 0.
        /// Note the special coefficient order ascending by exponent (consistent with polynomials).
        /// </summary>
        public static Tuple<Complex, Complex> Quadratic(double c, double b, double a)
        {
            if (b == 0d)
            {
                var t = new Complex(-c/a, 0d).SquareRoot();
                return new Tuple<Complex, Complex>(t, -t);
            }

            var q = b > 0d
                ? -0.5*(b + new Complex(b*b - 4*a*c, 0d).SquareRoot())
                : -0.5*(b - new Complex(b*b - 4*a*c, 0d).SquareRoot());

            return new Tuple<Complex, Complex>(q/a, c/q);
        }

        /// <summary>
        /// Find all three complex roots of the cubic equation d + c*x + b*x^2 + a*x^3 = 0.
        /// Note the special coefficient order ascending by exponent (consistent with polynomials).
        /// </summary>
        public static Tuple<Complex, Complex, Complex> Cubic(double d, double c, double b, double a)
        {
            return RootFinding.Cubic.Roots(d, c, b, a);
        }

        /// <summary>
        /// Find all roots of a polynomial by calculating the characteristic polynomial of the companion matrix
        /// </summary>
        /// <param name="coefficients">The coefficients of the polynomial in ascending order, e.g. new double[] {5, 0, 2} = "5 + 0 x^1 + 2 x^2"</param>
        /// <returns>The roots of the polynomial</returns>
        public static Complex[] Polynomial(double[] coefficients)
        {
            return new Polynomial(coefficients).Roots();
        }

        /// <summary>
        /// Find all roots of a polynomial by calculating the characteristic polynomial of the companion matrix
        /// </summary>
        /// <param name="polynomial">The polynomial.</param>
        /// <returns>The roots of the polynomial</returns>
        public static Complex[] Polynomial(Polynomial polynomial)
        {
            return polynomial.Roots();
        }

        /// <summary>
        /// Find all roots of the Chebychev polynomial of the first kind.
        /// </summary>
        /// <param name="degree">The polynomial order and therefore the number of roots.</param>
        /// <param name="intervalBegin">The real domain interval begin where to start sampling.</param>
        /// <param name="intervalEnd">The real domain interval end where to stop sampling.</param>
        /// <returns>Samples in [a,b] at (b+a)/2+(b-1)/2*cos(pi*(2i-1)/(2n))</returns>
        public static double[] ChebychevPolynomialFirstKind(int degree, double intervalBegin = -1d, double intervalEnd = 1d)
        {
            if (degree < 1)
            {
                return new double[0];
            }

            // transform to map to [-1..1] interval
            double location = 0.5*(intervalBegin + intervalEnd);
            double scale = 0.5*(intervalEnd - intervalBegin);

            // evaluate first kind chebychev nodes
            double angleFactor = Constants.Pi/(2*degree);

            var samples = new double[degree];
            for (int i = 0; i < samples.Length; i++)
            {
                samples[i] = location + scale*Math.Cos(((2*i) + 1)*angleFactor);
            }
            return samples;
        }

        /// <summary>
        /// Find all roots of the Chebychev polynomial of the second kind.
        /// </summary>
        /// <param name="degree">The polynomial order and therefore the number of roots.</param>
        /// <param name="intervalBegin">The real domain interval begin where to start sampling.</param>
        /// <param name="intervalEnd">The real domain interval end where to stop sampling.</param>
        /// <returns>Samples in [a,b] at (b+a)/2+(b-1)/2*cos(pi*i/(n-1))</returns>
        public static double[] ChebychevPolynomialSecondKind(int degree, double intervalBegin = -1d, double intervalEnd = 1d)
        {
            if (degree < 1)
            {
                return new double[0];
            }

            // transform to map to [-1..1] interval
            double location = 0.5*(intervalBegin + intervalEnd);
            double scale = 0.5*(intervalEnd - intervalBegin);

            // evaluate second kind chebychev nodes
            double angleFactor = Constants.Pi/(degree + 1);

            var samples = new double[degree];
            for (int i = 0; i < samples.Length; i++)
            {
                samples[i] = location + scale*Math.Cos((i + 1)*angleFactor);
            }
            return samples;
        }
    }
}