<commit_before>﻿// <copyright file="Operations.Algebraics.cs" company="Shkyrockett" >
//    Copyright © 2005 - 2019 Shkyrockett. All rights reserved.
// </copyright>
// <author id="shkyrockett">Shkyrockett</author>
// <license>
//    Licensed under the MIT License. See LICENSE file in the project root for full license information.
// </license>
// <summary></summary>
// <remarks></remarks>

// <copyright company="kevlindev" >
//     Many of the Roots methods were adapted from Kevin Lindsey's site http://www.kevlindev.com/gui/math/intersection/.
//     Copyright © 2000 - 2003 <NAME>. All rights reserved.
// </copyright>
// <author id="thelonious"><NAME></author>
// <license>
//     Licensed under the BSD-3-Clause https://github.com/thelonious/kld-intersections/blob/development/LICENSE
// </license>

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Runtime.CompilerServices;
using static System.Math;
using static Engine.Mathematics;

namespace Engine
{
    /// <summary>
    /// Algebraic Operations class.
    /// </summary>
    public static partial class Operations
    {
        /// <summary>
        /// The inverse sqrt.
        /// </summary>
        /// <param name="number">The number.</param>
        /// <returns>The <see cref="double"/>.</returns>
        [DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static double InverseSqrt(double number) => 1d / Sqrt(number);

        /// <summary>
        /// Returns the specified root a specified number.
        /// </summary>
        /// <param name="x">A double-precision floating-point number to find the specified root of.</param>
        /// <param name="y">A double-precision floating-point number that specifies a root.</param>
        /// <returns>The y root of the number x.</returns>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static double Root(double x, double y) => (x < 0d && Math.Abs((y % 2d) - 1d) < double.Epsilon) ? -Pow(-x, 1d / y) : Pow(x, 1d / y);

        /// <summary>
        /// Cube root equivalent of the sqrt function. (note that there are actually
        /// three roots: one real, two complex, and we don't care about the latter):
        /// </summary>
        /// <param name="value"></param>
        /// <returns></returns>
        /// <acknowledgment>
        /// http://stackoverflow.com/questions/26823024/cubic-bezier-reverse-getpoint-equation-float-for-vector-vector-for-float?answertab=active#tab-top
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static double Crt(double value) => value < 0 ? -Pow(-value, OneThird) : Pow(value, OneThird);

        /// <summary>
        /// The inverse cube root.
        /// </summary>
        /// <param name="number">The number.</param>
        /// <returns>The <see cref="double"/>.</returns>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static double InverseCrt(double number) => 1d / Crt(number);

        /// <summary>
        /// Calculates the real order or degree of the polynomial.
        /// </summary>
        /// <param name="coefficients">The coefficients.</param>
        /// <param name="epsilon">The <paramref name="epsilon"/> or minimal value to represent a change.</param>
        /// <returns>Returns a <see cref="PolynomialDegree"/> value representing the order of degree of the polynomial.</returns>
        /// <remarks><para>Primarily used to locate where to trim off any leading zero coefficients of the internal coefficients array.</para></remarks>
        /// <acknowledgment>
        /// A hodge-podge helper method based on Simplify from of: http://www.kevlindev.com/
        /// as well as Trim and RealOrder from: https://github.com/superlloyd/Poly
        /// </acknowledgment>
        [DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static PolynomialDegree DegreeRealOrder(IList<double> coefficients, double epsilon = double.Epsilon)
        {
            var pos = 1;
            var count = coefficients?.Count;

            // Monomial can be a zero constant, skip them and check the rest.
            if (count > 1)
            {
                // Count the number of leading zeros. Because the coefficients array is reversed, start at the end because there should generally be fewer leading zeros than other coefficients.
                for (var i = count.Value - 1; i >= 1 /* Monomials can be 0. */; i--)
                {
                    // Check if coefficient is a leading zero.
                    if (Math.Abs(coefficients[i]) <= epsilon)
                    {
                        pos++;
                    }
                    else
                    {
                        // Break early if a non zero value was found. This indicates the end of any leading zeros.
                        break;
                    }
                }
            }

            // If coefficients is empty return constant, otherwise return the calculated order of degree of the polynomial.
            return (PolynomialDegree)(count - pos ?? 0);
        }

        /// <summary>
        /// Align points to a line.
        /// </summary>
        /// <param name="points">The points to align.</param>
        /// <param name="x1">The x1.</param>
        /// <param name="y1">The y1.</param>
        /// <param name="x2">The x2.</param>
        /// <param name="y2">The y2.</param>
        /// <returns>The <see cref="List{T}"/>.</returns>
        /// <acknowledgment>
        /// https://pomax.github.io/bezierinfo/#aligning
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<Point2D> AlignPoints(IList<Point2D> points, double x1, double y1, double x2, double y2)
        {
            //var angle = -Atan2(y2 - y1, x2 - x1);
            //var sinA = Sin(angle);
            //var cosA = Cos(angle);

            // Atan2, Sin and Cos are kind of slow. In theory this should be faster.
            var dx = x2 - x1;
            var dy = y2 - y1;
            var det = (dx * dx) + (dy * dy);

            // I believe det should only be 0 if the line is a point. Not sure what the correct value shoud be for overlapping points.
            var sinA = det == 0 ? 0 : -dy / det;
            var cosA = det == 0 ? 1 : -dx / det;

            var results = new List<Point2D>();
            if (!(points is null))
            {
                foreach (var point in points)
                {
                    results.Add(new Point2D(
                        ((point.X - x1) * cosA) - ((point.Y - y1) * sinA),
                        ((point.X - x1) * sinA) + ((point.Y - y1) * cosA))
                        );
                }
            }

            return results;
        }

        #region D Root Finding
        /// <summary>
        /// The d roots.
        /// </summary>
        /// <param name="coefficients">The coefficients.</param>
        /// <param name="epsilon">The <paramref name="epsilon"/> or minimal value to represent a change.</param>
        /// <returns>The <see cref="Array"/>.</returns>
        /// <acknowledgment>
        /// http://pomax.github.io/bezierinfo/
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> DRoots(IList<double> coefficients, double epsilon = double.Epsilon)
            // ToDo: What are DRoots?
            => coefficients?.Count switch
            {
                2 => LinearDRoots(coefficients[0], coefficients[1], epsilon),
                3 => QuadraticDRoots(coefficients[0], coefficients[1], coefficients[2], epsilon),
                _ => Array.Empty<double>(),
            };

        /// <summary>
        /// The linear d roots.
        /// </summary>
        /// <param name="a">The a.</param>
        /// <param name="b">The b.</param>
        /// <param name="epsilon">The <paramref name="epsilon"/> or minimal value to represent a change.</param>
        /// <returns>The <see cref="IList{T}"/>.</returns>
        /// <acknowledgment>
        /// http://pomax.github.io/bezierinfo/
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> LinearDRoots(double a, double b, double epsilon = double.Epsilon)
        {
            // ToDo: What are DRoots?
            _ = epsilon;
            return a != b ? (new double[] { a / (a - b) }) : Array.Empty<double>();
        }

        /// <summary>
        /// The quadratic d roots.
        /// </summary>
        /// <param name="a">The a.</param>
        /// <param name="b">The b.</param>
        /// <param name="c">The c.</param>
        /// <param name="epsilon">The <paramref name="epsilon"/> or minimal value to represent a change.</param>
        /// <returns>The <see cref="IList{T}"/>.</returns>
        /// <acknowledgment>
        /// http://pomax.github.io/bezierinfo/
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> QuadraticDRoots(double a, double b, double c, double epsilon = double.Epsilon)
        {
            // ToDo: What are DRoots?
            _ = epsilon;
            var det = a - (2d * b) + c;
            if (det != 0)
            {
                // Negative square root discriminant. Missing the 4?
                var sqrtd = -Sqrt((b * b) - (a * c));
                var m2 = b - a;
                var v1 = -(m2 + sqrtd) / det;
                var v2 = -(m2 - sqrtd) / det;
                return new double[] { v1, v2 };
            }
            else if (b != c && det == 0)
            {
                return new double[] { ((2d * b) - c) / (2d * (b - c)) };
            }

            return Array.Empty<double>();
        }
        #endregion D Root Finding

        #region Root Finding
        /// <summary>
        /// Find the Roots of up to Quintic degree <see cref="Polynomial"/>s.
        /// </summary>
        /// <param name="coefficients">The coefficients.</param>
        /// <param name="epsilon">The <paramref name="epsilon"/> or minimal value to represent a change.</param>
        /// <returns>The <see cref="Array"/>.</returns>
        /// <acknowledgment>
        /// http://www.kevlindev.com/geometry/2D/intersections/
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> Roots(IList<double> coefficients, double epsilon = double.Epsilon)
            => DegreeRealOrder(coefficients) switch
            {
                PolynomialDegree.Constant => (coefficients is null) ? Array.Empty<double>() : new double[] { coefficients[0] },
                PolynomialDegree.Linear => LinearRoots(coefficients[1], coefficients[0], epsilon),
                PolynomialDegree.Quadratic => QuadraticRoots(coefficients[2], coefficients[1], coefficients[0], epsilon),
                PolynomialDegree.Cubic => CubicRoots(coefficients[3], coefficients[2], coefficients[1], coefficients[0], epsilon),
                PolynomialDegree.Quartic => QuarticRoots(coefficients[4], coefficients[3], coefficients[2], coefficients[1], coefficients[0], epsilon),
                PolynomialDegree.Quintic => QuinticRoots(coefficients[5], coefficients[4], coefficients[3], coefficients[2], coefficients[1], coefficients[0], epsilon),
                PolynomialDegree.Sextic => Array.Empty<double>(), //SexticRoots(epsilon),
                PolynomialDegree.Septic => Array.Empty<double>(), //SepticRoots(epsilon),
                PolynomialDegree.Octic => Array.Empty<double>(), //OcticRoots(epsilon),
                _ => Array.Empty<double>(), // Should try Newton's method and/or bisection
            };

        /// <summary>
        /// The linear roots.
        /// </summary>
        /// <param name="a">The a.</param>
        /// <param name="b">The b.</param>
        /// <param name="epsilon">The <paramref name="epsilon"/> or minimal value to represent a change.</param>
        /// <returns>The <see cref="List{T}"/>.</returns>
        /// <acknowledgment>
        /// http://www.kevlindev.com/geometry/2D/intersections/
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> LinearRoots(double a, double b, double epsilon = double.Epsilon)
        {
            var result = new HashSet<double>();
            if (Math.Abs(a) < epsilon)
            {
                if (!(Math.Abs(b) < epsilon))
                {
                    result.Add(b);
                }
                else
                {
                    return result.ToList();
                }
            }
            result.Add(-b / a);

            return result.ToList();
        }

        /// <summary>
        /// The quadratic roots.
        /// </summary>
        /// <param name="a">The a.</param>
        /// <param name="b">The b.</param>
        /// <param name="c">The c.</param>
        /// <param name = "epsilon"> The minimal value to represent a change.</param>
        /// <returns>The <see cref="List{T}"/>.</returns>
        /// <acknowledgment>
        /// http://www.kevlindev.com/geometry/2D/intersections/
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> QuadraticRoots(double a, double b, double c, double epsilon = double.Epsilon)
        {
            // Is the coefficient of the highest term zero?
            if (Math.Abs(a) < epsilon)
            {
                // If the highest term coefficient is 0, then it is a lower degree polynomial.
                return LinearRoots(b, c, epsilon);
            }

            var b_ = b / a;
            var c_ = c / a;

            // Polynomial discriminant
            var discriminant = (b_ * b_) - (4d * c_);

            // ToDo: May need to switch from a hash set to a list for scan-beams.
            var results = new HashSet<double>();

            if (Math.Abs(discriminant) <= epsilon)
            {
                discriminant = 0;
            }

            if (discriminant > 0)
            {
                // Complex or duplicate roots
                var e = Sqrt(discriminant);
                results.Add(OneHalf * (-b_ + e));
                results.Add(OneHalf * (-b_ - e));
            }
            else if (discriminant == 0)
            {
                // really two roots with same value, but we only return one
                results.Add(OneHalf * -b_);
            }

            return results.ToList();
        }

        /// <summary>
        /// Cubic Roots
        /// </summary>
        /// <param name="a">t^3</param>
        /// <param name="b">t^2</param>
        /// <param name="c">t</param>
        /// <param name="d">1</param>
        /// <param name="epsilon">The <paramref name="epsilon"/> or minimal value to represent a change.</param>
        /// <returns></returns>
        /// <acknowledgment>
        /// based on http://abecedarical.com/javascript/script_exact_cubic.html
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> CubicRoots(double a, double b, double c, double d, double epsilon = double.Epsilon)
        {
            // Is the coefficient of the highest term zero?
            if (Math.Abs(a) < epsilon)
            {
                // If the highest term coefficient is 0, then it is a lower degree polynomial.
                return QuadraticRoots(b, c, d, epsilon);
            }

            var A = b / a;
            var B = c / a;
            var C = d / a;

            var Q = ((3d * B) - (A * A)) / 9d;
            var R = (-(2d * A * A * A) + (9d * A * B) - (27d * C)) / 54d;

            var offset = A * OneThird;

            // Polynomial discriminant
            var discriminant = (R * R) + (Q * Q * Q);

            // ToDo: May need to switch from a hash set to a list for scan-beams.
            var results = new HashSet<double>();

            if (Math.Abs(discriminant) <= epsilon)
            {
                discriminant = 0d;
            }

            if (discriminant == 0d)
            {
                var t = Sign(R) * Pow(Math.Abs(R), OneThird);

                // Real root.
                results.Add(-offset + (t + t));

                // Real part of complex root.
                results.Add(-offset - ((t + t) * OneHalf));
            }
            if (discriminant > 0)
            {
                var s = Sign(R + Sqrt(discriminant)) * Pow(Math.Abs(R + Sqrt(discriminant)), OneThird);
                var t = Sign(R - Sqrt(discriminant)) * Pow(Math.Abs(R - Sqrt(discriminant)), OneThird);

                // Real root.
                results.Add(-offset + (s + t));

                // Complex part of root pair.
                var Im = Math.Abs(Sqrt(3d) * (s - t) * OneHalf);
                if (Im == 0d)
                {
                    // Real part of complex root.
                    results.Add(-offset - ((s + t) * OneHalf));
                }
            }
            else if (discriminant < 0)
            {
                // Distinct real roots.
                var th = Acos(R / Sqrt(-Q * Q * Q));

                results.Add((2d * Sqrt(-Q) * Cos(th * OneThird)) - offset);
                results.Add((2d * Sqrt(-Q) * Cos((th + Tau) * OneThird)) - offset);
                results.Add((2d * Sqrt(-Q) * Cos((th + (4d * PI)) * OneThird)) - offset);
            }

            return results.ToList();
        }

        /// <summary>
        /// The quartic roots.
        /// </summary>
        /// <param name="a">The a.</param>
        /// <param name="b">The b.</param>
        /// <param name="c">The c.</param>
        /// <param name="d">The d.</param>
        /// <param name="e">The e.</param>
        /// <param name = "epsilon"> The minimal value to represent a change.</param>
        /// <returns>The <see cref="List{T}"/>.</returns>
        /// <remarks>
        /// <para>ToDo: Translate code found at: https://web.archive.org/web/20150504111126/http://abecedarical.com/javascript/script_quintic.html
        /// and http://jwezorek.com/2015/01/my-code-for-doing-two-things-that-sooner-or-later-you-will-want-to-do-with-bezier-curves/:
        /// This method computes complex and real roots for any quintic polynomial. Then returns the real roots.
        /// It applies the Lin-Bairstow algorithm which iteratively solves for the
        /// roots starting from random guesses for a solution.
        /// The calculator is designed to solve for the roots of a quintic polynomial
        /// with the form: x⁵ + ax⁴ + bx³ + cx² + dx + e = 0
        /// ⁰¹²³⁴⁵⁶⁷⁸⁹</para>
        /// </remarks>
        /// <acknowledgment>
        /// http://www.kevlindev.com/geometry/2D/intersections/
        /// https://web.archive.org/web/20150504111126/http://abecedarical.com/javascript/script_quintic.html
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> QuarticRoots(double a, double b, double c, double d, double e, double epsilon = double.Epsilon)
        {
            // Is the coefficient of the highest term zero?
            if (Math.Abs(a) < epsilon)
            {
                // If the highest term coefficient is 0, then it is a lower degree polynomial.
                return CubicRoots(b, c, d, e, epsilon);
            }

            // ToDo: Translate code found at: https://web.archive.org/web/20150504111126/http://abecedarical.com/javascript/script_quintic.html
            // and http://jwezorek.com/2015/01/my-code-for-doing-two-things-that-sooner-or-later-you-will-want-to-do-with-bezier-curves/

            var A = b / a;
            var B = c / a;
            var C = d / a;
            var D = e / a;

            var resolveRoots = CubicRoots(
                1d,
                -B,
                (A * C) - (4d * D),
                (-A * A * D) + (4d * B * D) - (C * C),
                epsilon);
            var y = resolveRoots[0];
            var discriminant = (A * A * OneQuarter) - B + y;

            // ToDo: May need to switch from a hash set to a list for scan-beams.
            var results = new HashSet<double>();

            if (Math.Abs(discriminant) <= epsilon)
            {
                discriminant = 0d;
            }

            if (discriminant > 0d)
            {
                var ee = Sqrt(discriminant);
                var t1 = (3d * A * A * OneQuarter) - (ee * ee) - (2d * B);
                var t2 = ((4d * A * B) - (8d * C) - (A * A * A)) / (4d * ee);
                var plus = t1 + t2;
                var minus = t1 - t2;
                if (Math.Abs(plus) <= epsilon)
                {
                    plus = 0d;
                }

                if (Math.Abs(minus) <= epsilon)
                {
                    minus = 0d;
                }

                if (plus >= 0d)
                {
                    var f = Sqrt(plus);
                    results.Add((-A * OneQuarter) + ((ee + f) * OneHalf));
                    results.Add((-A * OneQuarter) + ((ee - f) * OneHalf));
                }
                if (minus >= 0d)
                {
                    var f = Sqrt(minus);
                    results.Add((-A * OneQuarter) + ((f - ee) * OneHalf));
                    results.Add((-A * OneQuarter) - ((f + ee) * OneHalf));
                }
            }
            else if (discriminant < 0d)
            {
            }
            else
            {
                var t2 = (y * y) - (4d * D);
                if (t2 >= -epsilon)
                {
                    if (t2 < 0)
                    {
                        t2 = 0d;
                    }

                    t2 = 2d * Sqrt(t2);
                    var t1 = (3d * A * A * OneQuarter) - (2d * B);
                    if (t1 + t2 >= epsilon)
                    {
                        var d0 = Sqrt(t1 + t2);
                        results.Add((-A * OneQuarter) + (d0 * OneHalf));
                        results.Add((-A * OneQuarter) - (d0 * OneHalf));
                    }
                    if (t1 - t2 >= epsilon)
                    {
                        var d1 = Sqrt(t1 - t2);
                        results.Add((-A * OneQuarter) + (d1 * OneHalf));
                        results.Add((-A * OneQuarter) - (d1 * OneHalf));
                    }
                }
            }

            return results.ToList();
        }

        /// <summary>
        /// The quintic roots.
        /// </summary>
        /// <param name="a">The a.</param>
        /// <param name="b">The b.</param>
        /// <param name="c">The c.</param>
        /// <param name="d">The d.</param>
        /// <param name="e">The e.</param>
        /// <param name="f">The f.</param>
        /// <param name = "epsilon"> The minimal value to represent a change.</param>
        /// <returns>The <see cref="List{T}"/>.</returns>
        /// <acknowledgment>
        /// This is a Copy and paste port of the method found at:
        /// https://web.archive.org/web/20150504111126/http://abecedarical.com/javascript/script_quintic.html
        /// http://www.convertalot.com/quintic_root_calculator.html
        /// There has been little attempt to fix it up and get it working correctly.
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> QuinticRoots(double a, double b, double c, double d, double e, double f, double epsilon = double.Epsilon)
        {
            // Is the coefficient of the highest term zero?
            if (Math.Abs(a) < epsilon)
            {
                // If the highest term coefficient is 0, then it is a lower degree polynomial.
                return QuarticRoots(b, c, d, e, f, epsilon);
            }

            var coeff = new List<double> { a, b, c, d, e, f };

            // Order
            var n = 4; // 5;
            var n1 = 5; // 6;
            var n2 = 6; // 7;

            var a_ = new List<double> { 0d, 0d, 0d, 0d, 0d, 0d };
            var b_ = new List<double> { 0d, 0d, 0d, 0d, 0d, 0d };
            //var c_ = new List<double> { 0d, 0d, 0d, 0d, 0d, 0d };
            var d_ = new List<double> { 0d, 0d, 0d, 0d, 0d, 0d };
            var real = new List<double> { 0d, 0d, 0d, 0d, 0d, 0d };
            var imag = new List<double> { 0d, 0d, 0d, 0d, 0d, 0d };

            //  Copy into working array
            for (var i = 0; i <= n; i++)
            {
                a_[a_.Count - 1 - i] = coeff[i];
            }

            // Initialize root counter
            var count = 0;

            // Start the main Lin-Bairstow iteration loop
            do
            {
                // Initialize the counter and guesses for the coefficients of quadratic factor: p(x) = x^2 + alfa1*x + beta1
                // ToDo: The random alphas make this method non-deterministic. Need a better guess method.
                var alfa1 = Mathematics.Random(OneHalf, 1d);
                var beta1 = Mathematics.Random(OneHalf, 1d);
                var limit = 1000;

                double delta1;
                do
                {
                    b_[0] = 0d;
                    d_[0] = 0d;
                    b_[1] = 1d;
                    d_[1] = 1d;

                    for (int i = 2, j = 1, k = 0; i < a_.Count; i++)
                    {
                        b_[i] = a_[i] - (alfa1 * b_[j]) - (beta1 * b_[k]);
                        d_[i] = b_[i] - (alfa1 * d_[j]) - (beta1 * d_[k]);
                        j += 1;
                        k += 1;
                    }

                    double alfa2;

                    double beta2;
                    {
                        var j = n - 1;
                        var k = n - 2;
                        delta1 = (d_[j] * d_[j]) - ((d_[n] - b_[n]) * d_[k]);
                        alfa2 = ((b_[n] * d_[j]) - (b_[n1] * d_[k])) / delta1;
                        beta2 = ((b_[n1] * d_[j]) - ((d_[n] - b_[n]) * b_[n])) / delta1;
                        alfa1 += alfa2;
                        beta1 += beta2;
                    }

                    if (--limit < 0)
                    {
                        // Cannot solve
                        return new List<double>();
                    }

                    if (Math.Abs(alfa2) < epsilon && Math.Abs(beta2) < epsilon)
                    {
                        break;
                    }
                }
                while (true);

                delta1 = (alfa1 * alfa1) - (4d * beta1);

                double delta2;
                // Imaginary roots
                if (delta1 < 0)
                {
                    delta2 = Sqrt(Math.Abs(delta1)) * OneHalf;
                    var delta3 = -alfa1 * OneHalf;

                    real[count] = delta3;
                    imag[count] = delta2;

                    real[count + 1] = delta3;
                    // Sign is inverted on display
                    imag[count + 1] = delta2;
                }
                else
                {
                    // Roots are real
                    delta2 = Sqrt(delta1);

                    real[count] = (delta2 - alfa1) * OneHalf;
                    imag[count] = 0;

                    real[count + 1] = (delta2 + alfa1) * -OneHalf;
                    imag[count + 1] = 0;
                }

                // Update root counter
                count += 2;

                // Reduce polynomial order
                n -= 2;
                n1 -= 2;
                n2 -= 2;

                // For n >= 2 calculate coefficients of
                // The new polynomial
                if (n >= 2)
                {
                    for (var i = 1; i <= n1; i++)
                    {
                        a_[i] = b_[i];
                    }
                }

                if (n < 2)
                {
                    break;
                }
            }
            while (true);

            if (n == 1)
            {
                // Obtain last single real root
                real[count] = -b_[2];
                imag[count] = 0;
            }

            return real;
        }
        #endregion Root Finding

        /// <summary>
        /// Newton's (Newton-Raphson) method for finding Real roots on univariate function. <br/>
        /// When using bounds, algorithm falls back to secant if newton goes out of range.
        /// Bisection is fall-back for secant when determined secant is not efficient enough.
        /// </summary>
        /// <param name="x0">Initial root guess</param>
        /// <param name="f">Function which root we are trying to find</param>
        /// <param name="df">Derivative of function f</param>
        /// <param name="maxIterations">Maximum number of algorithm iterations</param>
        /// <param name="min">Left bound value</param>
        /// <param name="max">Right bound value</param>
        /// <returns>root</returns>
        /// <remarks>
        /// <para>https://github.com/thelonious/kld-polynomial
        /// http://en.wikipedia.org/wiki/Newton%27s_method
        /// http://en.wikipedia.org/wiki/Secant_method
        /// http://en.wikipedia.org/wiki/Bisection_method</para>
        /// </remarks>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static double NewtonSecantBisection(double x0, Func<double, double> f, Func<double, double> df, int maxIterations, double? min = null, double? max = null)
        {
            var prev_dfx = 0d;
            var prev_x_ef_correction = 0d;
            var y_atmin = 0d;
            var y_atmax = 0d;
            var x = x0;
            const int ACCURACY = 14;
            var min_correction_factor = Pow(10, -ACCURACY);
            var isBounded = min != null && max != null;
            if (isBounded)
            {
                if (min > max)
                {
                    throw new Exception("newton root finding: min must be greater than max");
                }

                y_atmin = f(min.Value);
                y_atmax = f(max.Value);
                if (Sign(y_atmin) == Sign(y_atmax))
                {
                    throw new Exception("newton root finding: y values of bounds must be of opposite sign");
                }
            }

            double x_correction;
            bool isEnoughCorrection()<commit_msg>Fixing Ellipse Intersection Edge case<commit_after>﻿// <copyright file="Operations.Algebraics.cs" company="Shkyrockett" >
//    Copyright © 2005 - 2019 Shkyrockett. All rights reserved.
// </copyright>
// <author id="shkyrockett">Shkyrockett</author>
// <license>
//    Licensed under the MIT License. See LICENSE file in the project root for full license information.
// </license>
// <summary></summary>
// <remarks></remarks>

// <copyright company="kevlindev" >
//     Many of the Roots methods were adapted from Kevin Lindsey's site http://www.kevlindev.com/gui/math/intersection/.
//     Copyright © 2000 - 2003 <NAME>. All rights reserved.
// </copyright>
// <author id="thelonious"><NAME></author>
// <license>
//     Licensed under the BSD-3-Clause https://github.com/thelonious/kld-intersections/blob/development/LICENSE
// </license>

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Runtime.CompilerServices;
using static System.Math;
using static Engine.Mathematics;

namespace Engine
{
    /// <summary>
    /// Algebraic Operations class.
    /// </summary>
    public static partial class Operations
    {
        /// <summary>
        /// The inverse sqrt.
        /// </summary>
        /// <param name="number">The number.</param>
        /// <returns>The <see cref="double"/>.</returns>
        [DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static double InverseSqrt(double number) => 1d / Sqrt(number);

        /// <summary>
        /// Returns the specified root a specified number.
        /// </summary>
        /// <param name="x">A double-precision floating-point number to find the specified root of.</param>
        /// <param name="y">A double-precision floating-point number that specifies a root.</param>
        /// <returns>The y root of the number x.</returns>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static double Root(double x, double y) => (x < 0d && Math.Abs((y % 2d) - 1d) < double.Epsilon) ? -Pow(-x, 1d / y) : Pow(x, 1d / y);

        /// <summary>
        /// Cube root equivalent of the sqrt function. (note that there are actually
        /// three roots: one real, two complex, and we don't care about the latter):
        /// </summary>
        /// <param name="value"></param>
        /// <returns></returns>
        /// <acknowledgment>
        /// http://stackoverflow.com/questions/26823024/cubic-bezier-reverse-getpoint-equation-float-for-vector-vector-for-float?answertab=active#tab-top
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static double Crt(double value) => value < 0 ? -Pow(-value, OneThird) : Pow(value, OneThird);

        /// <summary>
        /// The inverse cube root.
        /// </summary>
        /// <param name="number">The number.</param>
        /// <returns>The <see cref="double"/>.</returns>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static double InverseCrt(double number) => 1d / Crt(number);

        /// <summary>
        /// Calculates the real order or degree of the polynomial.
        /// </summary>
        /// <param name="coefficients">The coefficients.</param>
        /// <param name="epsilon">The <paramref name="epsilon"/> or minimal value to represent a change.</param>
        /// <returns>Returns a <see cref="PolynomialDegree"/> value representing the order of degree of the polynomial.</returns>
        /// <remarks><para>Primarily used to locate where to trim off any leading zero coefficients of the internal coefficients array.</para></remarks>
        /// <acknowledgment>
        /// A hodge-podge helper method based on Simplify from of: http://www.kevlindev.com/
        /// as well as Trim and RealOrder from: https://github.com/superlloyd/Poly
        /// </acknowledgment>
        [DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static PolynomialDegree DegreeRealOrder(IList<double> coefficients, double epsilon = double.Epsilon)
        {
            var pos = 1;
            var count = coefficients?.Count;

            // Monomial can be a zero constant, skip them and check the rest.
            if (count > 1)
            {
                // Count the number of leading zeros. Because the coefficients array is reversed, start at the end because there should generally be fewer leading zeros than other coefficients.
                for (var i = count.Value - 1; i >= 1 /* Monomials can be 0. */; i--)
                {
                    // Check if coefficient is a leading zero.
                    if (Math.Abs(coefficients[i]) <= epsilon)
                    {
                        pos++;
                    }
                    else
                    {
                        // Break early if a non zero value was found. This indicates the end of any leading zeros.
                        break;
                    }
                }
            }

            // If coefficients is empty return constant, otherwise return the calculated order of degree of the polynomial.
            return (PolynomialDegree)(count - pos ?? 0);
        }

        /// <summary>
        /// Align points to a line.
        /// </summary>
        /// <param name="points">The points to align.</param>
        /// <param name="x1">The x1.</param>
        /// <param name="y1">The y1.</param>
        /// <param name="x2">The x2.</param>
        /// <param name="y2">The y2.</param>
        /// <returns>The <see cref="List{T}"/>.</returns>
        /// <acknowledgment>
        /// https://pomax.github.io/bezierinfo/#aligning
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<Point2D> AlignPoints(IList<Point2D> points, double x1, double y1, double x2, double y2)
        {
            //var angle = -Atan2(y2 - y1, x2 - x1);
            //var sinA = Sin(angle);
            //var cosA = Cos(angle);

            // Atan2, Sin and Cos are kind of slow. In theory this should be faster.
            var dx = x2 - x1;
            var dy = y2 - y1;
            var det = (dx * dx) + (dy * dy);

            // I believe det should only be 0 if the line is a point. Not sure what the correct value shoud be for overlapping points.
            var sinA = det == 0 ? 0 : -dy / det;
            var cosA = det == 0 ? 1 : -dx / det;

            var results = new List<Point2D>();
            if (!(points is null))
            {
                foreach (var point in points)
                {
                    results.Add(new Point2D(
                        ((point.X - x1) * cosA) - ((point.Y - y1) * sinA),
                        ((point.X - x1) * sinA) + ((point.Y - y1) * cosA))
                        );
                }
            }

            return results;
        }

        #region D Root Finding
        /// <summary>
        /// The d roots.
        /// </summary>
        /// <param name="coefficients">The coefficients.</param>
        /// <param name="epsilon">The <paramref name="epsilon"/> or minimal value to represent a change.</param>
        /// <returns>The <see cref="Array"/>.</returns>
        /// <acknowledgment>
        /// http://pomax.github.io/bezierinfo/
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> DRoots(IList<double> coefficients, double epsilon = double.Epsilon)
            // ToDo: What are DRoots?
            => coefficients?.Count switch
            {
                2 => LinearDRoots(coefficients[0], coefficients[1], epsilon),
                3 => QuadraticDRoots(coefficients[0], coefficients[1], coefficients[2], epsilon),
                _ => Array.Empty<double>(),
            };

        /// <summary>
        /// The linear d roots.
        /// </summary>
        /// <param name="a">The a.</param>
        /// <param name="b">The b.</param>
        /// <param name="epsilon">The <paramref name="epsilon"/> or minimal value to represent a change.</param>
        /// <returns>The <see cref="IList{T}"/>.</returns>
        /// <acknowledgment>
        /// http://pomax.github.io/bezierinfo/
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> LinearDRoots(double a, double b, double epsilon = double.Epsilon)
        {
            // ToDo: What are DRoots?
            _ = epsilon;
            return a != b ? (new double[] { a / (a - b) }) : Array.Empty<double>();
        }

        /// <summary>
        /// The quadratic d roots.
        /// </summary>
        /// <param name="a">The a.</param>
        /// <param name="b">The b.</param>
        /// <param name="c">The c.</param>
        /// <param name="epsilon">The <paramref name="epsilon"/> or minimal value to represent a change.</param>
        /// <returns>The <see cref="IList{T}"/>.</returns>
        /// <acknowledgment>
        /// http://pomax.github.io/bezierinfo/
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> QuadraticDRoots(double a, double b, double c, double epsilon = double.Epsilon)
        {
            // ToDo: What are DRoots?
            _ = epsilon;
            var det = a - (2d * b) + c;
            if (det != 0)
            {
                // Negative square root discriminant. Missing the 4?
                var sqrtd = -Sqrt((b * b) - (a * c));
                var m2 = b - a;
                var v1 = -(m2 + sqrtd) / det;
                var v2 = -(m2 - sqrtd) / det;
                return new double[] { v1, v2 };
            }
            else if (b != c && det == 0)
            {
                return new double[] { ((2d * b) - c) / (2d * (b - c)) };
            }

            return Array.Empty<double>();
        }
        #endregion D Root Finding

        #region Root Finding
        /// <summary>
        /// Find the Roots of up to Quintic degree <see cref="Polynomial"/>s.
        /// </summary>
        /// <param name="coefficients">The coefficients.</param>
        /// <param name="epsilon">The <paramref name="epsilon"/> or minimal value to represent a change.</param>
        /// <returns>The <see cref="Array"/>.</returns>
        /// <acknowledgment>
        /// http://www.kevlindev.com/geometry/2D/intersections/
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> Roots(IList<double> coefficients, double epsilon = double.Epsilon)
            => DegreeRealOrder(coefficients) switch
            {
                PolynomialDegree.Constant => (coefficients is null) ? Array.Empty<double>() : new double[] { coefficients[0] },
                PolynomialDegree.Linear => LinearRoots(coefficients[1], coefficients[0], epsilon),
                PolynomialDegree.Quadratic => QuadraticRoots(coefficients[2], coefficients[1], coefficients[0], epsilon),
                PolynomialDegree.Cubic => CubicRoots(coefficients[3], coefficients[2], coefficients[1], coefficients[0], epsilon),
                PolynomialDegree.Quartic => QuarticRoots(coefficients[4], coefficients[3], coefficients[2], coefficients[1], coefficients[0], epsilon),
                PolynomialDegree.Quintic => QuinticRoots(coefficients[5], coefficients[4], coefficients[3], coefficients[2], coefficients[1], coefficients[0], epsilon),
                PolynomialDegree.Sextic => Array.Empty<double>(), //SexticRoots(epsilon),
                PolynomialDegree.Septic => Array.Empty<double>(), //SepticRoots(epsilon),
                PolynomialDegree.Octic => Array.Empty<double>(), //OcticRoots(epsilon),
                _ => Array.Empty<double>(), // Should try Newton's method and/or bisection
            };

        /// <summary>
        /// The linear roots.
        /// </summary>
        /// <param name="a">The a.</param>
        /// <param name="b">The b.</param>
        /// <param name="epsilon">The <paramref name="epsilon"/> or minimal value to represent a change.</param>
        /// <returns>The <see cref="List{T}"/>.</returns>
        /// <acknowledgment>
        /// http://www.kevlindev.com/geometry/2D/intersections/
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> LinearRoots(double a, double b, double epsilon = double.Epsilon)
        {
            var result = new HashSet<double>();
            if (Math.Abs(a) < epsilon)
            {
                if (!(Math.Abs(b) < epsilon))
                {
                    result.Add(b);
                }
                else
                {
                    return result.ToList();
                }
            }
            result.Add(-b / a);

            return result.ToList();
        }

        /// <summary>
        /// The quadratic roots.
        /// </summary>
        /// <param name="a">The a.</param>
        /// <param name="b">The b.</param>
        /// <param name="c">The c.</param>
        /// <param name = "epsilon"> The minimal value to represent a change.</param>
        /// <returns>The <see cref="List{T}"/>.</returns>
        /// <acknowledgment>
        /// http://www.kevlindev.com/geometry/2D/intersections/
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> QuadraticRoots(double a, double b, double c, double epsilon = double.Epsilon)
        {
            // Is the coefficient of the highest term zero?
            if (Math.Abs(a) < epsilon)
            {
                // If the highest term coefficient is 0, then it is a lower degree polynomial.
                return LinearRoots(b, c, epsilon);
            }

            var ba = b / a;
            var ca = c / a;

            // Polynomial discriminant
            var discriminant = (ba * ba) - (4d * ca);

            // ToDo: May need to switch from a hash set to a list for scan-beams.
            var results = new HashSet<double>();

            if (Math.Abs(discriminant) <= epsilon)
            {
                discriminant = 0;
            }

            if (discriminant > 0)
            {
                // Real Roots.
                var e = Sqrt(discriminant);
                results.Add(OneHalf * (-ba + e));
                results.Add(OneHalf * (-ba - e));
            }
            else if (discriminant == 0)
            {
                // Technically two roots with same value, but we only need the one.
                results.Add(OneHalf * -ba);
            }
            else
            {
                // Imaginary roots? Technically incorrect, but it seems to fix the problem where ellipses intersections were broken... 
                var e = -Sqrt(-discriminant);
                results.Add(OneHalf * (-ba + e));
                results.Add(OneHalf * (-ba - e));
            }

            return results.ToList();
        }

        /// <summary>
        /// Cubic Roots
        /// </summary>
        /// <param name="a">t^3</param>
        /// <param name="b">t^2</param>
        /// <param name="c">t</param>
        /// <param name="d">1</param>
        /// <param name="epsilon">The <paramref name="epsilon"/> or minimal value to represent a change.</param>
        /// <returns></returns>
        /// <acknowledgment>
        /// based on http://abecedarical.com/javascript/script_exact_cubic.html
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> CubicRoots(double a, double b, double c, double d, double epsilon = double.Epsilon)
        {
            // Is the coefficient of the highest term zero?
            if (Math.Abs(a) < epsilon)
            {
                // If the highest term coefficient is 0, then it is a lower degree polynomial.
                return QuadraticRoots(b, c, d, epsilon);
            }

            var ba = b / a;
            var ca = c / a;
            var da = d / a;

            var q = ((3d * ca) - (ba * ba)) / 9d;
            var r = (-(2d * ba * ba * ba) + (9d * ba * ca) - (27d * da)) / 54d;

            var offset = ba * OneThird;

            // Polynomial discriminant
            var discriminant = (r * r) + (q * q * q);

            // ToDo: May need to switch from a hash set to a list for scan-beams.
            var results = new HashSet<double>();

            if (Math.Abs(discriminant) <= epsilon)
            {
                discriminant = 0d;
            }

            if (discriminant == 0d)
            {
                var t = Sign(r) * Pow(Math.Abs(r), OneThird);

                // Real root.
                results.Add(-offset + (t + t));

                // Real part of complex root.
                results.Add(-offset - ((t + t) * OneHalf));
            }
            if (discriminant > 0)
            {
                var s = Sign(r + Sqrt(discriminant)) * Pow(Math.Abs(r + Sqrt(discriminant)), OneThird);
                var t = Sign(r - Sqrt(discriminant)) * Pow(Math.Abs(r - Sqrt(discriminant)), OneThird);

                // Real root.
                results.Add(-offset + (s + t));

                // Complex part of root pair.
                var Im = Math.Abs(Sqrt(3d) * (s - t) * OneHalf);
                if (Im == 0d)
                {
                    // Real part of complex root.
                    results.Add(-offset - ((s + t) * OneHalf));
                }
            }
            else if (discriminant < 0)
            {
                // Distinct real roots.
                var th = Acos(r / Sqrt(-q * q * q));

                results.Add((2d * Sqrt(-q) * Cos(th * OneThird)) - offset);
                results.Add((2d * Sqrt(-q) * Cos((th + Tau) * OneThird)) - offset);
                results.Add((2d * Sqrt(-q) * Cos((th + (4d * PI)) * OneThird)) - offset);
            }

            return results.ToList();
        }

        /// <summary>
        /// The quartic roots.
        /// </summary>
        /// <param name="a">The a.</param>
        /// <param name="b">The b.</param>
        /// <param name="c">The c.</param>
        /// <param name="d">The d.</param>
        /// <param name="e">The e.</param>
        /// <param name = "epsilon"> The minimal value to represent a change.</param>
        /// <returns>The <see cref="List{T}"/>.</returns>
        /// <remarks>
        /// <para>ToDo: Translate code found at: https://web.archive.org/web/20150504111126/http://abecedarical.com/javascript/script_quintic.html
        /// and http://jwezorek.com/2015/01/my-code-for-doing-two-things-that-sooner-or-later-you-will-want-to-do-with-bezier-curves/:
        /// This method computes complex and real roots for any quintic polynomial. Then returns the real roots.
        /// It applies the Lin-Bairstow algorithm which iteratively solves for the
        /// roots starting from random guesses for a solution.
        /// The calculator is designed to solve for the roots of a quintic polynomial
        /// with the form: x⁵ + ax⁴ + bx³ + cx² + dx + e = 0
        /// ⁰¹²³⁴⁵⁶⁷⁸⁹</para>
        /// </remarks>
        /// <acknowledgment>
        /// http://www.kevlindev.com/geometry/2D/intersections/
        /// https://web.archive.org/web/20150504111126/http://abecedarical.com/javascript/script_quintic.html
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> QuarticRoots(double a, double b, double c, double d, double e, double epsilon = double.Epsilon)
        {
            // Is the coefficient of the highest term zero?
            if (Math.Abs(a) < epsilon)
            {
                // If the highest term coefficient is 0, then it is a lower degree polynomial.
                return CubicRoots(b, c, d, e, epsilon);
            }

            // ToDo: Translate code found at: https://web.archive.org/web/20150504111126/http://abecedarical.com/javascript/script_quintic.html
            // and http://jwezorek.com/2015/01/my-code-for-doing-two-things-that-sooner-or-later-you-will-want-to-do-with-bezier-curves/

            var ba = b / a;
            var ca = c / a;
            var da = d / a;
            var ea = e / a;

            var resolveRoots = CubicRoots(
                1d,
                -ca,
                (ba * da) - (4d * ea),
                (-ba * ba * ea) + (4d * ca * ea) - (da * da),
                epsilon);
            var y = resolveRoots[0];
            var discriminant = (ba * ba * OneQuarter) - ca + y;

            // ToDo: May need to switch from a hash set to a list for scan-beams.
            var results = new HashSet<double>();

            if (Math.Abs(discriminant) <= epsilon)
            {
                discriminant = 0d;
            }

            if (discriminant > 0d)
            {
                // Real Roots.
                var ee = Sqrt(discriminant);
                var t1 = (3d * ba * ba * OneQuarter) - (ee * ee) - (2d * ca);
                var t2 = ((4d * ba * ca) - (8d * da) - (ba * ba * ba)) / (4d * ee);
                var plus = t1 + t2;
                var minus = t1 - t2;
                if (Math.Abs(plus) <= epsilon)
                {
                    plus = 0d;
                }

                if (Math.Abs(minus) <= epsilon)
                {
                    minus = 0d;
                }

                if (plus >= 0d)
                {
                    var f = Sqrt(plus);
                    results.Add((-ba * OneQuarter) + ((ee + f) * OneHalf));
                    results.Add((-ba * OneQuarter) + ((ee - f) * OneHalf));
                }
                if (minus >= 0d)
                {
                    var f = Sqrt(minus);
                    results.Add((-ba * OneQuarter) + ((f - ee) * OneHalf));
                    results.Add((-ba * OneQuarter) - ((f + ee) * OneHalf));
                }
            }
            else if (discriminant == 0d)
            {
                // Origin Roots.
                var t2 = (y * y) - (4d * ea);
                if (t2 >= -epsilon)
                {
                    if (t2 < 0)
                    {
                        t2 = 0d;
                    }

                    t2 = 2d * Sqrt(t2);
                    var t1 = (3d * ba * ba * OneQuarter) - (2d * ca);
                    if (t1 + t2 >= epsilon)
                    {
                        var d0 = Sqrt(t1 + t2);
                        results.Add((-ba * OneQuarter) + (d0 * OneHalf));
                        results.Add((-ba * OneQuarter) - (d0 * OneHalf));
                    }
                    if (t1 - t2 >= epsilon)
                    {
                        var d1 = Sqrt(t1 - t2);
                        results.Add((-ba * OneQuarter) + (d1 * OneHalf));
                        results.Add((-ba * OneQuarter) - (d1 * OneHalf));
                    }
                }
            }
            else if (discriminant < 0d)
            {
                // Imaginary roots?
            }

            return results.ToList();
        }

        /// <summary>
        /// The quintic roots.
        /// </summary>
        /// <param name="a">The a.</param>
        /// <param name="b">The b.</param>
        /// <param name="c">The c.</param>
        /// <param name="d">The d.</param>
        /// <param name="e">The e.</param>
        /// <param name="f">The f.</param>
        /// <param name = "epsilon"> The minimal value to represent a change.</param>
        /// <returns>The <see cref="List{T}"/>.</returns>
        /// <acknowledgment>
        /// This is a Copy and paste port of the method found at:
        /// https://web.archive.org/web/20150504111126/http://abecedarical.com/javascript/script_quintic.html
        /// http://www.convertalot.com/quintic_root_calculator.html
        /// There has been little attempt to fix it up and get it working correctly.
        /// </acknowledgment>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static IList<double> QuinticRoots(double a, double b, double c, double d, double e, double f, double epsilon = double.Epsilon)
        {
            // Is the coefficient of the highest term zero?
            if (Math.Abs(a) < epsilon)
            {
                // If the highest term coefficient is 0, then it is a lower degree polynomial.
                return QuarticRoots(b, c, d, e, f, epsilon);
            }

            var coeff = new List<double> { a, b, c, d, e, f };

            // Order
            var n = 4; // 5;
            var n1 = 5; // 6;
            var n2 = 6; // 7;

            var a_ = new List<double> { 0d, 0d, 0d, 0d, 0d, 0d };
            var b_ = new List<double> { 0d, 0d, 0d, 0d, 0d, 0d };
            //var c_ = new List<double> { 0d, 0d, 0d, 0d, 0d, 0d };
            var d_ = new List<double> { 0d, 0d, 0d, 0d, 0d, 0d };
            var real = new List<double> { 0d, 0d, 0d, 0d, 0d, 0d };
            var imag = new List<double> { 0d, 0d, 0d, 0d, 0d, 0d };

            //  Copy into working array
            for (var i = 0; i <= n; i++)
            {
                a_[a_.Count - 1 - i] = coeff[i];
            }

            // Initialize root counter
            var count = 0;

            // Start the main Lin-Bairstow iteration loop
            do
            {
                // Initialize the counter and guesses for the coefficients of quadratic factor: p(x) = x^2 + alfa1*x + beta1
                // ToDo: The random alphas make this method non-deterministic. Need a better guess method.
                var alfa1 = Mathematics.Random(OneHalf, 1d);
                var beta1 = Mathematics.Random(OneHalf, 1d);
                var limit = 1000;

                double delta1;
                do
                {
                    b_[0] = 0d;
                    d_[0] = 0d;
                    b_[1] = 1d;
                    d_[1] = 1d;

                    for (int i = 2, j = 1, k = 0; i < a_.Count; i++)
                    {
                        b_[i] = a_[i] - (alfa1 * b_[j]) - (beta1 * b_[k]);
                        d_[i] = b_[i] - (alfa1 * d_[j]) - (beta1 * d_[k]);
                        j += 1;
                        k += 1;
                    }

                    double alfa2;

                    double beta2;
                    {
                        var j = n - 1;
                        var k = n - 2;
                        delta1 = (d_[j] * d_[j]) - ((d_[n] - b_[n]) * d_[k]);
                        alfa2 = ((b_[n] * d_[j]) - (b_[n1] * d_[k])) / delta1;
                        beta2 = ((b_[n1] * d_[j]) - ((d_[n] - b_[n]) * b_[n])) / delta1;
                        alfa1 += alfa2;
                        beta1 += beta2;
                    }

                    if (--limit < 0)
                    {
                        // Cannot solve
                        return new List<double>();
                    }

                    if (Math.Abs(alfa2) < epsilon && Math.Abs(beta2) < epsilon)
                    {
                        break;
                    }
                }
                while (true);

                delta1 = (alfa1 * alfa1) - (4d * beta1);

                double delta2;
                // Imaginary roots
                if (delta1 < 0)
                {
                    delta2 = Sqrt(Math.Abs(delta1)) * OneHalf;
                    var delta3 = -alfa1 * OneHalf;

                    real[count] = delta3;
                    imag[count] = delta2;

                    real[count + 1] = delta3;
                    // Sign is inverted on display
                    imag[count + 1] = delta2;
                }
                else
                {
                    // Roots are real
                    delta2 = Sqrt(delta1);

                    real[count] = (delta2 - alfa1) * OneHalf;
                    imag[count] = 0;

                    real[count + 1] = (delta2 + alfa1) * -OneHalf;
                    imag[count + 1] = 0;
                }

                // Update root counter
                count += 2;

                // Reduce polynomial order
                n -= 2;
                n1 -= 2;
                n2 -= 2;

                // For n >= 2 calculate coefficients of
                // The new polynomial
                if (n >= 2)
                {
                    for (var i = 1; i <= n1; i++)
                    {
                        a_[i] = b_[i];
                    }
                }

                if (n < 2)
                {
                    break;
                }
            }
            while (true);

            if (n == 1)
            {
                // Obtain last single real root
                real[count] = -b_[2];
                imag[count] = 0;
            }

            return real;
        }
        #endregion Root Finding

        /// <summary>
        /// Newton's (Newton-Raphson) method for finding Real roots on univariate function. <br/>
        /// When using bounds, algorithm falls back to secant if newton goes out of range.
        /// Bisection is fall-back for secant when determined secant is not efficient enough.
        /// </summary>
        /// <param name="x0">Initial root guess</param>
        /// <param name="f">Function which root we are trying to find</param>
        /// <param name="df">Derivative of function f</param>
        /// <param name="maxIterations">Maximum number of algorithm iterations</param>
        /// <param name="min">Left bound value</param>
        /// <param name="max">Right bound value</param>
        /// <returns>root</returns>
        /// <remarks>
        /// <para>https://github.com/thelonious/kld-polynomial
        /// http://en.wikipedia.org/wiki/Newton%27s_method
        /// http://en.wikipedia.org/wiki/Secant_method
        /// http://en.wikipedia.org/wiki/Bisection_method</para>
        /// </remarks>
        //[DebuggerStepThrough]
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static double NewtonSecantBisection(double x0, Func<double, double> f, Func<double, double> df, int maxIterations, double? min = null, double? max = null)
        {
            if (f is null)
            {
                throw new ArgumentNullException(nameof(f));
            }

            var prev_dfx = 0d;
            var prev_x_ef_correction = 0d;
            var y_atmin = 0d;
            var y_atmax = 0d;
            var x = x0;
            const int ACCURACY = 14;
            var min_correction_factor = Pow(10, -ACCURACY);
            var isBounded = min != null && max != null;
            if (isBounded)
            {
                if (min > max)
                {
                    throw new Exception("newton root finding: min must be greater than max");
                }

                y_atmin = f(min.Value);
                y_atmax = f(max.Value);
                if (Sign(y_atmin) == Sign(y_atmax))
                {
                    throw new Exception("newton root finding: y values of bounds must be of opposite sign");
                }
            }

            double x_correction;
            bool isEnoughCorrection()
{'fraction_non_alphanumeric': 0.13414672944769662, 'fraction_numerical': 0.01321901247136676, 'mean_word_length': 1.3735180850260515, 'pattern_counts': {'":': 0, '<': 627, '<?xml version=': 0, '>': 639, 'https://': 16, 'lorem ipsum': 0, 'www.': 14, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7387313', 'n_tokens_mistral': 19008, 'n_tokens_neox': 17374, 'n_words': 4778}