using System;
using System.Collections.Generic;
using System.Globalization;
using System.Runtime.Serialization;
using System.Linq;
using Complex = System.Numerics.Complex;
using System.Text;
using IStation.Numerics.LinearAlgebra;
using IStation.Numerics.LinearAlgebra.Double;
using IStation.Numerics.LinearRegression;
using IStation.Numerics.LinearAlgebra.Factorization;

#if !NETSTANDARD1_3
using System.Runtime;
#endif

namespace IStation.Numerics
{
    /// <summary>
    /// A single-variable polynomial with real-valued coefficients and non-negative exponents.
    /// </summary>
    [Serializable]
    [DataContract(Namespace = "urn:IStation/Numerics")]
    public class Polynomial : IFormattable, IEquatable<Polynomial>
#if !NETSTANDARD1_3
        , ICloneable
#endif
    {
        /// <summary>
        /// The coefficients of the polynomial in a
        /// </summary>
        [DataMember(Order = 1)]
        public double[] Coefficients { get; private set; }

        /// <summary>
        /// Only needed for the ToString method
        /// </summary>
        [DataMember(Order = 2)]
        public string VariableName = "x";

        /// <summary>
        /// Degree of the polynomial, i.e. the largest monomial exponent. For example, the degree of y=x^2+x^5 is 5, for y=3 it is 0.
        /// The null-polynomial returns degree -1 because the correct degree, negative infinity, cannot be represented by integers.
        /// </summary>
        public int Degree => EvaluateDegree(Coefficients);

        /// <summary>
        /// Create a zero-polynomial with a coefficient array of the given length.
        /// An array of length N can support polynomials of a degree of at most N-1.
        /// </summary>
        /// <param name="n">Length of the coefficient array</param>
        public Polynomial(int n)
        {
            if (n < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(n), "n must be non-negative");
            }

            Coefficients = new double[n];
        }

        /// <summary>
        /// Create a zero-polynomial
        /// </summary>
        public Polynomial()
        {
#if NET40
            Coefficients = new double[0];
#else
            Coefficients = Array.Empty<double>();
#endif
        }

        /// <summary>
        /// Create a constant polynomial.
        /// Example: 3.0 -> "p : x -> 3.0"
        /// </summary>
        /// <param name="coefficient">The coefficient of the "x^0" monomial.</param>
        public Polynomial(double coefficient)
        {
            if (coefficient == 0.0)
            {
#if NET40
                Coefficients = new double[0];
#else
                Coefficients = Array.Empty<double>();
#endif
            }
            else
            {
                Coefficients = new[] { coefficient };
            }
        }

        /// <summary>
        /// Create a polynomial with the provided coefficients (in ascending order, where the index matches the exponent).
        /// Example: {5, 0, 2} -> "p : x -> 5 + 0 x^1 + 2 x^2".
        /// </summary>
        /// <param name="coefficients">Polynomial coefficients as array</param>
        public Polynomial(params double[] coefficients)
        {
            Coefficients = coefficients;
        }

        /// <summary>
        /// Create a polynomial with the provided coefficients (in ascending order, where the index matches the exponent).
        /// Example: {5, 0, 2} -> "p : x -> 5 + 0 x^1 + 2 x^2".
        /// </summary>
        /// <param name="coefficients">Polynomial coefficients as enumerable</param>
        public Polynomial(IEnumerable<double> coefficients) : this(coefficients.ToArray())
        {
        }

        public static Polynomial Zero => new Polynomial();

        /// <summary>
        /// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k
        /// </summary>
        public static Polynomial Fit(double[] x, double[] y, int order, DirectRegressionMethod method = DirectRegressionMethod.QR)
        {
            var coefficients = Numerics.Fit.Polynomial(x, y, order, method);
            return new Polynomial(coefficients);
        }

        static int EvaluateDegree(double[] coefficients)
        {
            for (int i = coefficients.Length - 1; i >= 0; i--)
            {
                if (coefficients[i] != 0.0)
                {
                    return i;
                }
            }

            return -1;
        }

        #region Evaluation

        /// <summary>
        /// Evaluate a polynomial at point x.
        /// Coefficients are ordered ascending by power with power k at index k.
        /// Example: coefficients [3,-1,2] represent y=2x^2-x+3.
        /// </summary>
        /// <param name="z">The location where to evaluate the polynomial at.</param>
        /// <param name="coefficients">The coefficients of the polynomial, coefficient for power k at index k.</param>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="coefficients"/> is a null reference.
        /// </exception>
        public static double Evaluate(double z, params double[] coefficients)
        {

            // 2020-10-07 jbialogrodzki #730 Since this is public API we should probably
            // handle null arguments? It doesn't seem to have been done consistently in this class though.
            if (coefficients == null)
            {
                throw new ArgumentNullException(nameof(coefficients));
            }

            // 2020-10-07 jbialogrodzki #730 Zero polynomials need explicit handling.
            // Without this check, we attempted to peek coefficients at negative indices!
            int n = coefficients.Length;
            if (n == 0)
            {
                return 0;
            }

            double sum = coefficients[n - 1];
            for (int i = n - 2; i >= 0; --i)
            {
                sum *= z;
                sum += coefficients[i];
            }

            return sum;

        }

        /// <summary>
        /// Evaluate a polynomial at point x.
        /// Coefficients are ordered ascending by power with power k at index k.
        /// Example: coefficients [3,-1,2] represent y=2x^2-x+3.
        /// </summary>
        /// <param name="z">The location where to evaluate the polynomial at.</param>
        /// <param name="coefficients">The coefficients of the polynomial, coefficient for power k at index k.</param>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="coefficients"/> is a null reference.
        /// </exception>
        public static Complex Evaluate(Complex z, params double[] coefficients)
        {

            // 2020-10-07 jbialogrodzki #730 Since this is a public API we should probably
            // handle null arguments? It doesn't seem to have been done consistently in this class though.
            if (coefficients == null)
            {
                throw new ArgumentNullException(nameof(coefficients));
            }

            // 2020-10-07 jbialogrodzki #730 Zero polynomials need explicit handling.
            // Without this check, we attempted to peek coefficients at negative indices!
            int n = coefficients.Length;
            if (n == 0)
            {
                return 0;
            }

            Complex sum = coefficients[n - 1];
            for (int i = n - 2; i >= 0; --i)
            {
                sum *= z;
                sum += coefficients[i];
            }

            return sum;

        }

        /// <summary>
        /// Evaluate a polynomial at point x.
        /// Coefficients are ordered ascending by power with power k at index k.
        /// Example: coefficients [3,-1,2] represent y=2x^2-x+3.
        /// </summary>
        /// <param name="z">The location where to evaluate the polynomial at.</param>
        /// <param name="coefficients">The coefficients of the polynomial, coefficient for power k at index k.</param>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="coefficients"/> is a null reference.
        /// </exception>
        public static Complex Evaluate(Complex z, params Complex[] coefficients)
        {

            // 2020-10-07 jbialogrodzki #730 Since this is a public API we should probably
            // handle null arguments? It doesn't seem to have been done consistently in this class though.
            if (coefficients == null)
            {
                throw new ArgumentNullException(nameof(coefficients));
            }

            // 2020-10-07 jbialogrodzki #730 Zero polynomials need explicit handling.
            // Without this check, we attempted to peek coefficients at negative indices!
            int n = coefficients.Length;
            if (n == 0)
            {
                return 0;
            }

            Complex sum = coefficients[n - 1];
            for (int i = n - 2; i >= 0; --i)
            {
                sum *= z;
                sum += coefficients[i];
            }

            return sum;

        }

        /// <summary>
        /// Evaluate a polynomial at point x.
        /// </summary>
        /// <param name="z">The location where to evaluate the polynomial at.</param>
        public double Evaluate(double z)
        {
            return Evaluate(z, Coefficients);
        }

        /// <summary>
        /// Evaluate a polynomial at point x.
        /// </summary>
        /// <param name="z">The location where to evaluate the polynomial at.</param>
        public Complex Evaluate(Complex z)
        {
            return Evaluate(z, Coefficients);
        }

        /// <summary>
        /// Evaluate a polynomial at points z.
        /// </summary>
        /// <param name="z">The locations where to evaluate the polynomial at.</param>
        public IEnumerable<double> Evaluate(IEnumerable<double> z)
        {
            return z.Select(Evaluate);
        }

        /// <summary>
        /// Evaluate a polynomial at points z.
        /// </summary>
        /// <param name="z">The locations where to evaluate the polynomial at.</param>
        public IEnumerable<Complex> Evaluate(IEnumerable<Complex> z)
        {
            return z.Select(Evaluate);
        }

        #endregion

        #region Calculus

        public Polynomial Differentiate()
        {
            int n = Degree;
            if (n < 0)
            {
                return this;
            }

            if (n == 0)
            {
                // Zero
                return Zero;
            }

            var c = new double[n];
            for (int i = 0; i < c.Length; i++)
            {
                c[i] = Coefficients[i + 1] * (i + 1);
            }

            return new Polynomial(c);
        }

        public Polynomial Integrate()
        {
            int n = Degree;
            if (n < 0)
            {
                return this;
            }

            var c = new double[n + 2];
            for (int i = 1; i < c.Length; i++)
            {
                c[i] = Coefficients[i - 1] / i;
            }

            return new Polynomial(c);
        }

        #endregion

        #region Linear Algebra

        /// <summary>
        /// Calculates the complex roots of the Polynomial by eigenvalue decomposition
        /// </summary>
        /// <returns>a vector of complex numbers with the roots</returns>
        public Complex[] Roots()
        {
            switch (Degree)
            {
                case -1: // Zero-polynomial
                case 0: // Non-zero constant: y = a0
                    return new Complex[0];
                case 1: // Linear: y = a0 + a1*x
                    return new[] { new Complex(-Coefficients[0] / Coefficients[1], 0) };
            }

            DenseMatrix A = EigenvalueMatrix();
            Evd<double> eigen = A.Evd(Symmetricity.Asymmetric);
            return eigen.EigenValues.AsArray();
        }

        /// <summary>
        /// Get the eigenvalue matrix A of this polynomial such that eig(A) = roots of this polynomial.
        /// </summary>
        /// <returns>Eigenvalue matrix A</returns>
        /// <note>This matrix is similar to the companion matrix of this polynomial, in such a way, that it's transpose is the columnflip of the companion matrix</note>
        public DenseMatrix EigenvalueMatrix()
        {
            int n = Degree;
            if (n < 2)
            {
                return null;
            }

            // Negate, and normalize (scale such that the polynomial becomes monic)
            double aN = Coefficients[n];
            double[] p = new double[n];
            for (int i = n - 1; i >= 0; i--)
            {
                p[i] = -Coefficients[i] / aN;
            }

            DenseMatrix A0 = DenseMatrix.CreateDiagonal(n - 1, n - 1, 1.0);
            DenseMatrix A = new DenseMatrix(n);

            A.SetSubMatrix(1, 0, A0);
            A.SetRow(0, p.Reverse().ToArray());

            return A;
        }

        #endregion

        #region Arithmetic Operations

        /// <summary>
        /// Addition of two Polynomials (point-wise).
        /// </summary>
        /// <param name="a">Left Polynomial</param>
        /// <param name="b">Right Polynomial</param>
        /// <returns>Resulting Polynomial</returns>
        public static Polynomial Add(Polynomial a, Polynomial b)
        {
            var ac = a.Coefficients;
            var bc = b.Coefficients;

            var degree = Math.Max(a.Degree, b.Degree);
            var result = new double[degree + 1];

            var commonLength = Math.Min(Math.Min(ac.Length, bc.Length), result.Length);
            for (int i = 0; i < commonLength; i++)
            {
                result[i] = ac[i] + bc[i];
            }

            int acLength = Math.Min(ac.Length, result.Length);
            for (int i = commonLength; i < acLength; i++)
            {
                // no need to add since only one of both applies
                result[i] = ac[i];
            }

            int bcLength = Math.Min(bc.Length, result.Length);
            for (int i = commonLength; i < bcLength; i++)
            {
                // no need to add since only one of both applies
                result[i] = bc[i];
            }

            return new Polynomial(result);
        }

        /// <summary>
        /// Addition of a polynomial and a scalar.
        /// </summary>
        public static Polynomial Add(Polynomial a, double b)
        {
            var ac = a.Coefficients;

            var degree = Math.Max(a.Degree, 0);
            var result = new double[degree + 1];

            var commonLength = Math.Min(ac.Length, result.Length);
            for (int i = 0; i < commonLength; i++)
            {
                result[i] = ac[i];
            }

            result[0] += b;

            return new Polynomial(result);
        }

        /// <summary>
        /// Subtraction of two Polynomials (point-wise).
        /// </summary>
        /// <param name="a">Left Polynomial</param>
        /// <param name="b">Right Polynomial</param>
        /// <returns>Resulting Polynomial</returns>
        public static Polynomial Subtract(Polynomial a, Polynomial b)
        {
            var ac = a.Coefficients;
            var bc = b.Coefficients;

            var degree = Math.Max(a.Degree, b.Degree);
            var result = new double[degree + 1];

            var commonLength = Math.Min(Math.Min(ac.Length, bc.Length), result.Length);
            for (int i = 0; i < commonLength; i++)
            {
                result[i] = ac[i] - bc[i];
            }

            int acLength = Math.Min(ac.Length, result.Length);
            for (int i = commonLength; i < acLength; i++)
            {
                // no need to add since only one of both applies
                result[i] = ac[i];
            }

            int bcLength = Math.Min(bc.Length, result.Length);
            for (int i = commonLength; i < bcLength; i++)
            {
                // no need to add since only one of both applies
                result[i] = -bc[i];
            }

            return new Polynomial(result);
        }

        /// <summary>
        /// Addition of a scalar from a polynomial.
        /// </summary>
        public static Polynomial Subtract(Polynomial a, double b)
        {
            return Add(a, -b);
        }

        /// <summary>
        /// Addition of a polynomial from a scalar.
        /// </summary>
        public static Polynomial Subtract(double b, Polynomial a)
        {
            var ac = a.Coefficients;

            var degree = Math.Max(a.Degree, 0);
            var result = new double[degree + 1];

            var commonLength = Math.Min(ac.Length, result.Length);
            for (int i = 0; i < commonLength; i++)
            {
                result[i] = -ac[i];
            }

            result[0] += b;

            return new Polynomial(result);
        }

        /// <summary>
        /// Negation of a polynomial.
        /// </summary>
        public static Polynomial Negate(Polynomial a)
        {
            var ac = a.Coefficients;

            var degree = a.Degree;
            var result = new double[degree + 1];

            for (int i = 0; i < result.Length; i++)
            {
                result[i] = -ac[i];
            }

            return new Polynomial(result);
        }

        /// <summary>
        /// Multiplies a polynomial by a polynomial (convolution)
        /// </summary>
        /// <param name="a">Left polynomial</param>
        /// <param name="b">Right polynomial</param>
        /// <returns>Resulting Polynomial</returns>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="a"/> or <paramref name="b"/> is a null reference.
        /// </exception>
        public static Polynomial Multiply(Polynomial a, Polynomial b)
        {

            // 2020-10-07 jbialogrodzki #730 Since this is a public API we should probably
            // handle null arguments? It doesn't seem to have been done consistently in this class though.
            if (a == null)
            {
                throw new ArgumentNullException(nameof(a));
            }

            if (b == null)
            {
                throw new ArgumentNullException(nameof(b));
            }

            var ad = a.Degree;
            var bd = b.Degree;

            // 2020-10-07 jbialogrodzki #730 Zero polynomials need explicit handling.
            // Without this check, we attempted to create arrays of negative lengths!
            if (ad < 0 || bd < 0)
            {
                return Polynomial.Zero;
            }

            double[] ac = a.Coefficients;
            double[] bc = b.Coefficients;

            var degree = ad + bd;
            double[] result = new double[degree + 1];

            for (int i = 0; i <= ad; i++)
            {
                for (int j = 0; j <= bd; j++)
                {
                    result[i + j] += ac[i] * bc[j];
                }
            }

            return new Polynomial(result);

        }

        /// <summary>
        /// Scales a polynomial by a scalar
        /// </summary>
        /// <param name="a">Polynomial</param>
        /// <param name="k">Scalar value</param>
        /// <returns>Resulting Polynomial</returns>
        public static Polynomial Multiply(Polynomial a, double k)
        {
            var ac = a.Coefficients;

            var result = new double[a.Degree + 1];
            for (int i = 0; i < result.Length; i++)
            {
                result[i] = ac[i] * k;
            }

            return new Polynomial(result);
        }

        /// <summary>
        /// Scales a polynomial by division by a scalar
        /// </summary>
        /// <param name="a">Polynomial</param>
        /// <param name="k">Scalar value</param>
        /// <returns>Resulting Polynomial</returns>
        public static Polynomial Divide(Polynomial a, double k)
        {
            var ac = a.Coefficients;

            var result = new double[a.Degree + 1];
            for (int i = 0; i < result.Length; i++)
            {
                result[i] = ac[i] / k;
            }

            return new Polynomial(result);
        }

        /// <summary>
        /// Euclidean long division of two polynomials, returning the quotient q and remainder r of the two polynomials a and b such that a = q*b + r
        /// </summary>
        /// <param name="a">Left polynomial</param>
        /// <param name="b">Right polynomial</param>
        /// <returns>A tuple holding quotient in first and remainder in second</returns>
        public static Tuple<Polynomial, Polynomial> DivideRemainder(Polynomial a, Polynomial b)
        {
            var bDegree = b.Degree;
            if (bDegree < 0)
            {
                throw new DivideByZeroException("b polynomial ends with zero");
            }

            var aDegree = a.Degree;
            if (aDegree < 0)
            {
                // zero divided by non-zero is zero without remainder
                return Tuple.Create(a, a);
            }

            if (bDegree == 0)
            {
                // division by scalar
                return Tuple.Create(Divide(a, b.Coefficients[0]), Zero);
            }

            if (aDegree < bDegree)
            {
                // denominator degree higher than nominator degree
                // quotient always be 0 and return c1 as remainder
                return Tuple.Create(Zero, a);
            }

            var c1 = a.Coefficients.ToArray();
            var c2 = b.Coefficients.ToArray();

            var scl = c2[bDegree];
            var c22 = new double[bDegree];
            for (int ii = 0; ii < c22.Length; ii++)
            {
                c22[ii] = c2[ii] / scl;
            }

            int i = aDegree - bDegree;
            int j = aDegree;
            while (i >= 0)
            {
                var v = c1[j];
                for (int k = i; k < j; k++)
                {
                    c1[k] -= c22[k - i] * v;
                }

                i--;
                j--;
            }

            var j1 = j + 1;
            var l1 = aDegree - j;

            var quo = new double[l1];
            for (int k = 0; k < l1; k++)
            {
                quo[k] = c1[k + j1] / scl;
            }

            var rem = new double[j1];
            for (int k = 0; k < j1; k++)
            {
                rem[k] = c1[k];
            }

            return Tuple.Create(new Polynomial(quo), new Polynomial(rem));
        }

        #endregion

        #region Arithmetic Pointwise Operations

        /// <summary>
        /// Point-wise division of two Polynomials
        /// </summary>
        /// <param name="a">Left Polynomial</param>
        /// <param name="b">Right Polynomial</param>
        /// <returns>Resulting Polynomial</returns>
        public static Polynomial PointwiseDivide(Polynomial a, Polynomial b)
        {
            var ac = a.Coefficients;
            var bc = b.Coefficients;

            var degree = a.Degree;
            var result = new double[degree + 1];

            var commonLength = Math.Min(Math.Min(ac.Length, bc.Length), result.Length);
            for (int i = 0; i < commonLength; i++)
            {
                result[i] = ac[i] / bc[i];
            }

            for (int i = commonLength; i < result.Length; i++)
            {
                result[i] = ac[i] / 0.0;
            }

            return new Polynomial(result);
        }

        /// <summary>
        /// Point-wise multiplication of two Polynomials
        /// </summary>
        /// <param name="a">Left Polynomial</param>
        /// <param name="b">Right Polynomial</param>
        /// <returns>Resulting Polynomial</returns>
        public static Polynomial PointwiseMultiply(Polynomial a, Polynomial b)
        {
            var ac = a.Coefficients;
            var bc = b.Coefficients;

            var degree = Math.Min(a.Degree, b.Degree);
            var result = new double[degree + 1];
            for (int i = 0; i < result.Length; i++)
            {
                result[i] = ac[i] * bc[i];
            }

            return new Polynomial(result);
        }

        #endregion

        #region Arithmetic Instance Methods (forwarders)

        /// <summary>
        /// Division of two polynomials returning the quotient-with-remainder of the two polynomials given
        /// </summary>
        /// <param name="b">Right polynomial</param>
        /// <returns>A tuple holding quotient in first and remainder in second</returns>
        public Tuple<Polynomial, Polynomial> DivideRemainder(Polynomial b)
        {
            return DivideRemainder(this, b);
        }

        #endregion

        #region Arithmetic Operator Overloads (forwarders)

        /// <summary>
        /// Addition of two Polynomials (piecewise)
        /// </summary>
        /// <param name="a">Left polynomial</param>
        /// <param name="b">Right polynomial</param>
        /// <returns>Resulting Polynomial</returns>
        public static Polynomial operator +(Polynomial a, Polynomial b)
        {
            return Add(a, b);
        }

        /// <summary>
        /// adds a scalar to a polynomial.
        /// </summary>
        /// <param name="a">Polynomial</param>
        /// <param name="k">Scalar value</param>
        /// <returns>Resulting Polynomial</returns>
        public static Polynomial operator +(Polynomial a, double k)
        {
            return Add(a, k);
        }

        /// <summary>
        /// adds a scalar to a polynomial.
        /// </summary>
        /// <param name="k">Scalar value</param>
        /// <param name="a">Polynomial</param>
        /// <returns>Resulting Polynomial</returns>
        public static Polynomial operator +(double k, Polynomial a)
        {
            return Add(a, k);
        }

        /// <summary>
        /// Subtraction of two polynomial.
        /// </summary>
        /// <param name="a">Left polynomial</param>
        /// <param name="b">Right polynomial</param>
        /// <returns>Resulting Polynomial</returns>
        public static Polynomial operator -(Polynomial a, Polynomial b)
        {
            return Subtract(a, b);
        }

        /// <summary>
        /// Subtracts a scalar from a polynomial.
        /// </summary>
        /// <param name="a">Polynomial</param>
        /// <param name="k">Scalar value</param>
        /// <returns>Resulting Polynomial</returns>
        public static Polynomial operator -(Polynomial a, double k)
        {
            return Subtract(a, k);
        }

        /// <summary>
        /// Subtracts a polynomial from a scalar.
        /// </summary>
        /// <param name="k">Scalar value</param>
        /// <param name="a">Polynomial</param>
        /// <returns>Resulting Polynomial</returns>
        public static Polynomial operator -(double k, Polynomial a)
        {
            return Subtract(k, a);
        }

        /// <summary>
        /// Negates a polynomial.
        /// </summary>
        /// <param name="a">Polynomial</param>
        /// <returns>Resulting Polynomial</returns>
        public static Polynomial operator -(Polynomial a)
        {
            return Negate(a);
        }

        /// <summary>
        /// Multiplies a polynomial by a polynomial (convolution).
        /// </summary>
        /// <param name="a">Left polynomial</param>
        /// <param name="b">Right polynomial</param>
        /// <returns>resulting Polynomial</returns>
        public static Polynomial operator *(Polynomial a, Polynomial b)
        {
            return Multiply(a, b);
        }

        /// <summary>
        /// Multiplies a polynomial by a scalar.
        /// </summary>
        /// <param name="a">Polynomial</param>
        /// <param name="k">Scalar value</param>
        /// <returns>Resulting Polynomial</returns>
        public static Polynomial operator *(Polynomial a, double k)
        {
            return Multiply(a, k);
        }

        /// <summary>
        /// Multiplies a polynomial by a scalar.
        /// </summary>
        /// <param name="k">Scalar value</param>
        /// <param name="a">Polynomial</param>
        /// <returns>Resulting Polynomial</returns>
        public static Polynomial operator *(double k, Polynomial a)
        {
            return Multiply(a, k);
        }

        /// <summary>
        /// Divides a polynomial by scalar value.
        /// </summary>
        /// <param name="a">Polynomial</param>
        /// <param name="k">Scalar value</param>
        /// <returns>Resulting Polynomial</returns>
        public static Polynomial operator /(Polynomial a, double k)
        {
            return Divide(a, k);
        }

        #endregion

        #region ToString

        /// <summary>
        /// Format the polynomial in ascending order, e.g. "4.3 + 2.0x^2 - x^3".
        /// </summary>
        public override string ToString()
        {
            return ToString("G", CultureInfo.CurrentCulture);
        }

        /// <summary>
        /// Format the polynomial in descending order, e.g. "x^3 + 2.0x^2 - 4.3".
        /// </summary>
        public string ToStringDescending()
        {
            return ToStringDescending("G", CultureInfo.CurrentCulture);
        }

        /// <summary>
        /// Format the polynomial in ascending order, e.g. "4.3 + 2.0x^2 - x^3".
        /// </summary>
        public string ToString(string format)
        {
            return ToString(format, CultureInfo.CurrentCulture);
        }

        /// <summary>
        /// Format the polynomial in descending order, e.g. "x^3 + 2.0x^2 - 4.3".
        /// </summary>
        public string ToStringDescending(string format)
        {
            return ToStringDescending(format, CultureInfo.CurrentCulture);
        }

        /// <summary>
        /// Format the polynomial in ascending order, e.g. "4.3 + 2.0x^2 - x^3".
        /// </summary>
        public string ToString(IFormatProvider formatProvider)
        {
            return ToString("G", formatProvider);
        }

        /// <summary>
        /// Format the polynomial in descending order, e.g. "x^3 + 2.0x^2 - 4.3".
        /// </summary>
        public string ToStringDescending(IFormatProvider formatProvider)
        {
            return ToStringDescending("G", formatProvider);
        }

        /// <summary>
        /// Format the polynomial in ascending order, e.g. "4.3 + 2.0x^2 - x^3".
        /// </summary>
        public string ToString(string format, IFormatProvider formatProvider)
        {
            if (Degree < 0)
            {
                return "0";
            }

            var sb = new StringBuilder();
            bool first = true;
            for (int i = 0; i < Coefficients.Length; i++)
            {
                double c = Coefficients[i];
                if (c == 0.0)
                {
                    continue;
                }

                if (first)
                {
                    sb.Append(c.ToString(format, formatProvider));
                    if (i > 0)
                    {
                        sb.Append(VariableName);
                    }

                    if (i > 1)
                    {
                        sb.Append("^");
                        sb.Append(i);
                    }

                    first = false;
                }
                else
                {
                    if (c < 0.0)
                    {
                        sb.Append(" - ");
                        sb.Append((-c).ToString(format, formatProvider));
                    }
                    else
                    {
                        sb.Append(" + ");
                        sb.Append(c.ToString(format, formatProvider));
                    }

                    if (i > 0)
                    {
                        sb.Append(VariableName);
                    }

                    if (i > 1)
                    {
                        sb.Append("^");
                        sb.Append(i);
                    }
                }
            }

            return sb.ToString();
        }

        /// <summary>
        /// Format the polynomial in descending order, e.g. "x^3 + 2.0x^2 - 4.3".
        /// </summary>
        public string ToStringDescending(string format, IFormatProvider formatProvider)
        {
            if (Degree < 0)
            {
                return "0";
            }

            var sb = new StringBuilder();
            bool first = true;
            for (int i = Coefficients.Length - 1; i >= 0; i--)
            {
                double c = Coefficients[i];
                if (c == 0.0)
                {
                    continue;
                }

                if (first)
                {
                    sb.Append(c.ToString(format, formatProvider));
                    if (i > 0)
                    {
                        sb.Append(VariableName);
                    }

                    if (i > 1)
                    {
                        sb.Append("^");
                        sb.Append(i);
                    }

                    first = false;
                }
                else
                {
                    if (c < 0.0)
                    {
                        sb.Append(" - ");
                        sb.Append((-c).ToString(format, formatProvider));
                    }
                    else
                    {
                        sb.Append(" + ");
                        sb.Append(c.ToString(format, formatProvider));
                    }

                    if (i > 0)
                    {
                        sb.Append(VariableName);
                    }

                    if (i > 1)
                    {
                        sb.Append("^");
                        sb.Append(i);
                    }
                }
            }

            return sb.ToString();
        }

        #endregion

        #region Equality

        public bool Equals(Polynomial other)
        {
            if (ReferenceEquals(null, other)) return false;
            if (ReferenceEquals(this, other)) return true;

            int n = Degree;
            if (n != other.Degree)
            {
                return false;
            }

            for (var i = 0; i <= n; i++)
            {
                if (!Coefficients[i].Equals(other.Coefficients[i]))
                {
                    return false;
                }
            }

            return true;
        }

        public override bool Equals(object obj)
        {
            if (ReferenceEquals(null, obj)) return false;
            if (ReferenceEquals(this, obj)) return true;
            if (obj.GetType() != typeof(Polynomial)) return false;
            return Equals((Polynomial)obj);
        }

        public override int GetHashCode()
        {
            var hashNum = Math.Min(Degree + 1, 25);
            int hash = 17;
            unchecked
            {
                for (var i = 0; i < hashNum; i++)
                {
                    hash = hash * 31 + Coefficients[i].GetHashCode();
                }
            }

            return hash;
        }

        #endregion

        #region Clone

        public Polynomial Clone()
        {
            int degree = EvaluateDegree(Coefficients);
            var coefficients = new double[degree + 1];
            Array.Copy(Coefficients, coefficients, coefficients.Length);
            return new Polynomial(coefficients);
        }

#if !NETSTANDARD1_3
        /// <summary>
        /// Creates a new object that is a copy of the current instance.
        /// </summary>
        /// <returns>
        /// A new object that is a copy of this instance.
        /// </returns>
        object ICloneable.Clone()
        {
            return Clone();
        }
#endif

        #endregion
    }
}