Skalowanie obrazu z wykorzystaniem filtrów.

Skalowanie w dwóch wymiarach możemy zrealizować jako superpozycję skalowania w poszczególnych wymiarach osobno. Np. najpierw w poziomie, a później w pionie. Samo skalowanie to splot odpowiedzi impulsowej filtra z pikselami poszczególnych linii i kolumn. Splot po stronie wymiaru przestrzennego oznacza mnożenie po stronie częstotliwości.

Przy pomniejszaniu każdy piksel wyjściowy rekontrujemy z pewnego zbioru pikseli wejściowych. Jeśli promień filtra to 2, a współczynnik pomniejszania to 3 to bierzemy pod uwagę 12 pikseli źródła.

Przy powiększaniu bierzemy 4 piksele.

Na te 12 i 4 piksele rozciągamy odpowiednio filtr.

Ponieważ dla każdego wiersza i każdej kolumny współczynniki filtra są takie same warto je wstępnie policzyć.

Pojawia się też pytanie: co zrobić dla skrajnych pikseli, dla których obszar do filtrowania wychodzi poza bitmapę. Tutaj mamy generalnie 3 rozwiązania: kopiować skrajny piksel, dokonać lustrzanego odbicia lub zawinąć. To ostanie przydaje się kiedy zamierzamy teksturować przedmiot wielokrotnie tą samą teksturą. Pierwsze nadaje zbyt dużego znaczenia skrajnym pikselowi. Ja z reguły dla samplowania bitmap wykorzystuję metodę drugą.

Szerokość efektywna filtra (nasze 12 i 4) może być liczbą zmiennoprzecinkową. Sam środek naszego filtra także. Bardzo ważne jest tutaj takie dobranie jego szerokości całkowitej, by wartości filtra wewnątrz obszaru były różne od zera, a poza równe mu.

Na samym końcu warto jeszcze raz przejrzeć wagi i wyeliminować z nich te o wadze 0, które mimo wszystko mogą się pojawić. Poza tym możemy także znormalizować nasz filtr tak by pole powierzchni pod nim była równa jeden (suma wag filtra), dzięki temu unikniemy małych zmian w jasności.

Samego filtrowania powinniśmy dokonywać w liniowej przestrzeni kolorów.

Kod metody wyliczającej wagi:

public class PixelWeight : IComparable<PixelWeight>
{
    public double Weight;
    public int Pos;

    public override string ToString()
    {
        return String.Format("pos: {0}; weight: {1}", Pos, Weight);
    }

    public int CompareTo(PixelWeight a_other)
    {
        return Pos - a_other.Pos;
    }
}

public static void PrecalculateWeights(int a_src_width, int a_src_height,
    int a_dest_width, int a_dest_height, 
    out List<List<PixelWeight>> a_horz_weights,
    out List<List<PixelWeight>> a_vert_weights, Filter a_filter, 
    bool a_optimize = true)
{
    double scalex = 1.0 * a_dest_width / a_src_width;
    double scaley = 1.0 * a_dest_height / a_src_height;

    double src_ray_x = a_filter.Ray / scalex;
    double src_ray_y = a_filter.Ray / scaley;

    a_horz_weights = new List<List<PixelWeight>>(a_dest_width);

    if (scalex < 1)
    {
        OverlayCorrector overlay_corrector = OverlayCorrector.Create(
            OverlayMethod.Mirror,
            a_src_width, a_src_height, (int)(src_ray_x + 2));

        for (int x = 0; x < a_dest_width; x++)
        {
            double src_center_x = x / scalex;

            int start_xx = (int)Math.Floor(src_center_x - src_ray_x);
            int end_xx = (int)Math.Ceiling(src_center_x + src_ray_x);

            if (a_filter.FilterType == FilterType.NearestNeighbour)
                src_center_x = src_center_x.Round();

            List<PixelWeight> weights = new List<PixelWeight>();

            for (int xx = start_xx; xx <= end_xx; xx++)
            {
                weights.Add(new PixelWeight()
                {
                    Pos = overlay_corrector.CorrectX(xx),
                    Weight = a_filter.Evaluate((src_center_x - xx) * scalex) * scalex
                });
            }

            Debug.Assert(a_filter.Evaluate(
                (src_center_x - (start_xx - 1)) * scalex).IsAlmostEquals(0));
            Debug.Assert(a_filter.Evaluate(
                (src_center_x - (end_xx + 1)) * scalex).IsAlmostEquals(0));

            a_horz_weights.Add(weights);
        }
    }
    else
    {
        OverlayCorrector overlay_corrector = OverlayCorrector.Create(
            OverlayMethod.Mirror,
            a_src_width, a_src_height, (int)(a_filter.Ray + 2));

        for (int x = 0; x < a_dest_width; x++)
        {
            double src_center_x = x / scalex;

            int start_xx = (int)Math.Floor(src_center_x - a_filter.Ray);
            int end_xx = (int)Math.Ceiling(src_center_x + a_filter.Ray);

            List<PixelWeight> weights = new List<PixelWeight>();

            if (a_filter.FilterType == FilterType.NearestNeighbour)
                src_center_x = src_center_x.Round();

            for (int xx = start_xx; xx <= end_xx; xx++)
            {
                weights.Add(new PixelWeight()
                {
                    Pos = overlay_corrector.CorrectX(xx),
                    Weight = a_filter.Evaluate(src_center_x - xx)
                });
            }

            Debug.Assert(a_filter.Evaluate(
                src_center_x - (start_xx - 1)).IsAlmostEquals(0));
            Debug.Assert(a_filter.Evaluate(
                src_center_x - (end_xx + 1)).IsAlmostEquals(0));

            a_horz_weights.Add(weights);
        }
    }

    a_vert_weights = new List<List<PixelWeight>>(a_dest_height);

    if (scaley < 1)
    {
        OverlayCorrector overlay_corrector = OverlayCorrector.Create(
            OverlayMethod.Mirror,
            a_src_width, a_src_height, (int)(src_ray_y + 2));

        for (int y = 0; y < a_dest_height; y++)
        {
            double src_center_y = y / scaley;

            int start_yy = (int)Math.Floor(src_center_y - src_ray_y);
            int end_yy = (int)Math.Ceiling(src_center_y + src_ray_y);

            List<PixelWeight> weights = new List<PixelWeight>();

            if (a_filter.FilterType == FilterType.NearestNeighbour)
                src_center_y = src_center_y.Round();

            for (int yy = start_yy; yy <= end_yy; yy++)
            {
                weights.Add(new PixelWeight()
                {
                    Pos = overlay_corrector.CorrectY(yy),
                    Weight = a_filter.Evaluate((src_center_y - yy) * scaley) * scaley
                });
            }

            Debug.Assert(a_filter.Evaluate(
                (src_center_y - (start_yy - 1)) * scaley).IsAlmostEquals(0));
            Debug.Assert(a_filter.Evaluate(
                (src_center_y - (end_yy + 1)) * scaley).IsAlmostEquals(0));

            a_vert_weights.Add(weights);
        }
    }
    else
    {
        OverlayCorrector overlay_corrector = OverlayCorrector.Create(
            OverlayMethod.Mirror,
            a_src_width, a_src_height, (int)(a_filter.Ray + 2));

        for (int y = 0; y < a_dest_height; y++)
        {
            double src_center_y = y / scaley;

            int start_yy = (int)Math.Floor(src_center_y - a_filter.Ray);
            int end_yy = (int)Math.Ceiling(src_center_y + a_filter.Ray);

            List<PixelWeight> weights = new List<PixelWeight>();

            if (a_filter.FilterType == FilterType.NearestNeighbour)
                src_center_y = src_center_y.Round();

            for (int yy = start_yy; yy <= end_yy; yy++)
            {
                weights.Add(new PixelWeight()
                {
                    Pos = overlay_corrector.CorrectY(yy),
                    Weight = a_filter.Evaluate(src_center_y - yy)
                });
            }

            Debug.Assert(a_filter.Evaluate(
                src_center_y - (start_yy - 1)).IsAlmostEquals(0));
            Debug.Assert(a_filter.Evaluate(
                src_center_y - (end_yy + 1)).IsAlmostEquals(0));

            a_vert_weights.Add(weights);
        }
    }

    if (a_optimize)
    {
        foreach (var ws in a_vert_weights.Concat(a_horz_weights))
        {
            double sum = ws.Sum(w => w.Weight);

            foreach (var w in ws)
                w.Weight /= sum;

            for (int i = ws.Count - 1; i >= 0; i--)
            {
                if (ws[i].Weight.IsAlmostEquals(0))
                    ws.RemoveAt(i);
            }

            ws.Sort();
        }
    }
}

Kod samej metody reskalującej:

public void Resize(ColorArrayFloat a_dest,
    FilterType a_filter)
{
    List<List<PixelWeight>> horz_weights, vert_weights;
    ResizerResampler.PrecalculateWeights(Width, Height, a_dest.Width, a_dest.Height, 
        out horz_weights, out vert_weights, Filter.Create(a_filter));

    ColorArrayFloat temp_bd = new ColorArrayFloat(a_dest.Width, Height);

    for (int x = 0; x < temp_bd.Width; x++)
    {
        List<PixelWeight> weights = horz_weights[x];

        for (int y = 0; y < temp_bd.Height; y++)
        {
            ColorFloat color = new ColorFloat();

            foreach (PixelWeight pw in weights)
                color += GetColor(pw.Pos, y) * pw.Weight;

            temp_bd.SetColor(x, y, color.AboveZero);
        }
    }

    for (int y = 0; y < a_dest.Height; y++)
    {
        List<PixelWeight> weights = vert_weights[y];

        for (int x = 0; x < a_dest.Width; x++)
        {
            ColorFloat color = new ColorFloat();

            foreach (PixelWeight pw in weights)
                color += temp_bd.GetColor(x, pw.Pos) * pw.Weight;

            a_dest.SetColor(x, y, color.AboveZero);
        }
    }
}

Przykład klasy filtra:

public class LanczosFilter : Filter
{
    public double Tau = 3;

    public override double Ray
    {
        get
        {
            return Tau;
        }
    }

    public override double Evaluate(double a_value)
    {
        double v = Math.Abs(a_value);

        if (v > Tau)
            return 0;

        if (v < Constants.DOUBLE_PRECISION)
            return 1;

        v = v * MathExtensions.PI;
        double p3 = v / Tau;
        return Math.Sin(v) * Math.Sin(p3) / (v * p3);
    }
}

Root Finder - podsumowanie

Najgorzej sprawuje się metoda Secant, która potrafi zrenderować niepoprawnie torus. Pozostałe metody renderują go bezbłędnie. Miejsca w których metody nie polegające na pochodnej zawodzą to renderowanie prosto ze wzoru brył typu superelipsoida, supertoroid, superegg. Wydaje się, że o wiele lepiej powinny się tam sprawdzić metody polegające na pochodnej.

Porównanie czasów renderowania typowej sceny, w której funkcja dana jest wielomianem (intersekcja z torusem):

Bisection: 7.9
Modified Regula Falsi: 7.4
Secant: 7.7
Newton: 7.4
Halley: 7.7

Widzimy więc, że najszybsza metoda Modified Regula Falsi jest zarazem tą najmniej problematyczną - zawsze znajduje rozwiązanie. Po zaimplementowaniu bryły Superquadric mogłem trochę dokładniej przetestować algorytmy znajdowania miejsc zerowych. Okazało się, że Newton i Halley i prawdopodobnie wszystkie oparte na analizie pochodnych nie działają. Raz, że nie powinny, gdyż bryła ta nie ma powierzchni ciągłej. Dwa, że Root Finder-y tego typu bardzo łatwo potrafią wyskoczyć z lokalnego minimum i zacząć aproksymować do sąsiedniego miejsca zerowego. Na razie je zostawiłem, ale generalnie nie nadają się one w tej chwili do niczego. Próbowałem także zaimplementować Root Finder, który możemy nazwać Brute Forcem. Idziemy co krok rzędu 1e-4 i jeśli namierzymy zmianę znaku to szukamy w takim małym przedziale miejsca zerowego. Działa, ale jest niezwykle wolny, niepraktyczny.

Wszystkie metody znajdowania miejsc zerowych wykorzystują bezwzględny test na przerwanie iteracji. Test względny wymagał by znajomości wielkości obiektu. W moim przypadku większość renderowanych brył jest ograniczona w małym przedziale, typowo (-1,1), dzięki czemu RootFinder-y nie wymagają podawania względnych warunków na zakończenie testów.

Secant Root Finder

Kolejna metoda na znajdowanie pierwiastków równania. Do rozpoczęcia wymagane są dwa punkty a i b, które powinny być dostatecznie blisko miejsca zerowego. Nie jest wymagane by miejsce zerowe znajdowało się wewnątrz (a,b), ale wtedy wzrasta ryzyko, że metoda może nie znaleźć miejsca zerowego. Podobnie dzieje się jeśli w zakresie (a,b) istnieje ekstremum nie będące miejscem zerowym. Tak więc możemy powiedzieć, że zakres (a,b) jest mniejszy tym lepiej.

Przybliżanie następuje rekurencyjnie. W dowolnym momencie iteracji mamy dwa punkty $y_0=f(x_0)$ i $y_1=f(y_1)$. Wyznaczamy równanie linii przez nie przechodzące:

$y=\displaystyle\frac{f(x_1)-f(x_0)}{x_1-x_0}(x-x_1) + f(x_1)$

Wyznaczym punkt przecięcia się naszej funkcji z linią:

$0=\displaystyle\frac{f(x_1)-f(x_0)}{x_1-x_0}(x-x_1) + f(x_1)$

$x_2=\displaystyle x_1-f(x_1)\frac{x_1-x_0}{f(x_1)-f(x_0)}$

Bierzemy nowy przedział $(x_1, x_2)$ i powtarzamy iterację. Ponieważ metoda nie zawsze znajduje rozwiązanie powinniśmy wyznaczyć górną ilość iteracji. Normalnie metoda powinna zakończyć działanie po znalezieniu takie x dla, którego f(x) jest w granicach precyzji.

Kod:

public override double FindRoot(double a_a, double a_b, Func<double, double> a_func)
{
    Debug.Assert(a_a < a_b);

    double ya = a_func(a_a);
    double yb = a_func(a_b);

    if (Math.Abs(ya) < Constants.ROOT_FINDER_PRECISION)
        return a_a;

    if (Math.Abs(yb) < Constants.ROOT_FINDER_PRECISION)
        return a_b;

    if (Math.Sign(ya) * Math.Sign(yb) > 0.0)
        return Double.PositiveInfinity;

    double x = 0;

    for (int i = 0; i < Constants.ROOT_FINDER_MAXIMUM_ITERATIONS; i++)
    {
        x = a_a - ya * (a_a - a_b) / (ya - yb);

        Debug.Assert(x.IsNumber());

        double y = a_func(x);

        if (Math.Abs(x) > Constants.ROOT_FINDER_PRECISION)
        {
            if (Math.Abs(y / x) < Constants.ROOT_FINDER_PRECISION)
                return x;
        }
        else if (Math.Abs(y) < Constants.ROOT_FINDER_PRECISION)
            return x;

        yb = ya;
        a_b = a_a;

        ya = y;
        a_a = x;

        if (ya == yb)
            return Double.PositiveInfinity;
    }

    return x;
}

Halley Root Finder

Kolejna metoda znajdowania miejsc zerowych. W założeniach bardzo podobna do metody Newtona. Od funkcji wymaga się ciągłości pierwszej i drugiej pochodnej. Musimy znać wzory na te pochodne, by je ciągle przeliczać. Jest to metoda iteracyjna. Krok iteracji ma postać:

$x_1 = x_0 - \displaystyle \frac {2 f(x_0) f'(x_0)} {2 {[f'(x_0)]}^2 - f(x_0) f''(x_0)}$

Uzasadnienie jest dla mnie dość pokrętne. Weźmy funkcję:

$g(x) = \displaystyle\frac {f(x)} {\sqrt{|f'(x)|}}$

Miejsca zerowe g(x) to miejsca zerowe f(x). Miejsca zerowe f(x) o ile nie są miejscami zerowymi $f'(x)$ są miejscami zerowymi g(x). Funkcja g(x) ma wartości nieoznaczone w ekstremach. Funkcja ta została zbudowana w taki sposób nieprzypadkowo. Jeśli podstawimy ją do wzoru na metodę Newtona to otrzymamy:

$x_1=x_0-\displaystyle\frac{g(x_0)}{g'(x_0)}=\frac{f(x_0)}{\sqrt{\left|f'(x_0)\right|}}\left(\frac{\sqrt{\left|f'(x_0)\right|}}{f(x_0)}\right)'$
$x_1=\displaystyle\frac{f(x_0)}{\sqrt{\left|f'(x_0)\right|}} \frac{2 f'(x) \sqrt{|f'(x)|}} {2 {[f'(x)]}^2 - f(x) f''(x)}$

Podobnie jak w metodzie Newtona staramy się poszukiwania tą metodą przeprowadzić dwa razy, za każdym razem za punkt startowy obierając inny koniec zakresu (a,b). I podobnie jak w metodzie Newtona nie ma gwarancji na znalezienie miejsca zerowego.

public class CentralDifferenceMethod
{
    private double m_h;
    private Func<double, double> m_func;

    public CentralDifferenceMethod()
        : base()
    {
    }
        
    public CentralDifferenceMethod(Func<double, double> a_func, 
        double a_h = Constants.DOUBLE_PRECISION * 2)
    {
        m_h = a_h;
        m_func = a_func;
    }

    public double Derive1(double a_x)
    {
        return (m_func(a_x + m_h) - m_func(a_x - m_h)) / (2 * m_h);
    }

    public double Derive2(double a_x)
    {
        return (m_func(a_x + m_h) - 2 * m_func(a_x) + 
            m_func(a_x - m_h)) / (m_h * m_h);
    }

    public static  Func<double, double> Derive1(
        Func<double, double> a_func, double a_h)
    {
        CentralDifferenceMethod cdm = 
            new CentralDifferenceMethod(a_func, a_h);

        return (x) =>
        {
            return cdm.Derive1(x);
        };
    }

    public static Func<double, double> Derive2(
        Func<double, double> a_func, double a_h)
    {
        CentralDifferenceMethod cdm = 
            new CentralDifferenceMethod(a_func, a_h);

        return (x) =>
        {
            return cdm.Derive2(x);
        };
    }

    public static Func<double, double> Derive2(
        Func<double, double> a_func, 
        Func<double, double> a_d1, double a_h)
    {
        CentralDifferenceMethod cdm = 
            new CentralDifferenceMethod(a_d1, a_h);

        return (x) =>
        {
            return cdm.Derive1(x);
        };
    }
}

public abstract class RootFinder
{
    public int MaxIterations;
    public double AbsoluteError;

    public RootFinder()
    {
        MaxIterations = Constants.ROOT_FINDER_MAXIMUM_ITERATIONS;
        AbsoluteError = Constants.ROOT_FINDER_ABSOLUTE_ERROR;
    }

    public static RootFinder Create(RootFindingMethod a_method)
    {
        switch (a_method)
        {
            case RootFindingMethod.Bisection: 
                return new BisectionRootFinder();
            case RootFindingMethod.ModifedRegulaFalsi: 
                return new ModifiedRegulaFalsiRootFinder();
            case RootFindingMethod.RegulaFalsi: 
                return new RegulaFalsiRootFinder();
            case RootFindingMethod.Secant: 
                return new SecantRootFinder();
            case RootFindingMethod.Newton: 
                return new NewtonRootFinder();
            case RootFindingMethod.Halley: 
                return new HalleyRootFinder();
            default: throw new NotImplementedException();
        }
    }

    public abstract double FindRoot(double a_a, double a_b, 
        Polynomial a_poly);
    public abstract double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func);
    public abstract double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func, Func<double, double> a_d1);
    public abstract double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func, Func<double, double> a_d1, 
        Func<double, double> a_d2);
}

public abstract class DerivativeRootFinder : RootFinder
{
    protected struct RootState
    {
        public double Root;
        public bool Finded;

        public RootState(double a_x, bool a_finded = true)
        {
            Debug.Assert(a_x.IsNumber());

            Root = a_x;
            Finded = a_finded;
        }
    }

    protected double FindRootFromBothSide(double a_a, double a_b, 
        Func<double, double> a_func, Func<double, 
        RootState> a_root_finder)
    {
        Debug.Assert(a_a < a_b);
        Debug.Assert(a_a.IsNumber());
        Debug.Assert(a_b.IsNumber());

        double ya = a_func(a_a);
        double yb = a_func(a_b);

        if (Math.Abs(ya) < AbsoluteError)
            return a_a;
        if (Math.Abs(yb) < AbsoluteError)
            return a_b;

        if (Math.Sign(ya) * Math.Sign(yb) > 0)
            return Double.PositiveInfinity;

        RootState x1 = a_root_finder(a_a);
        if (x1.Finded && (x1.Root > a_a) && (x1.Root < a_b))
            return x1.Root;

        RootState x2 = a_root_finder(a_b);
        if (x2.Finded && (x2.Root > a_a) && (x2.Root < a_b))
            return x2.Root;

        return Double.PositiveInfinity;
    }
}

public class HalleyRootFinder : DerivativeRootFinder
{
    public override double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func)
    {
        return FindRoot(a_a, a_b, a_func, 
            CentralDifferenceMethod.Derive1(a_func, AbsoluteError),
            CentralDifferenceMethod.Derive2(a_func, AbsoluteError));
    }

    public override double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func, Func<double, double> a_d1)
    {
        return FindRoot(a_a, a_b, a_func, a_d1, 
            CentralDifferenceMethod.Derive2(a_func, a_d1, AbsoluteError));
    }

    public override double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func, Func<double, double> a_d1, 
        Func<double, double> a_d2)
    {
        return FindRootFromBothSide(a_a, a_b, a_func, 
            (s) => FindRoot(s, a_func, a_d1, a_d2));
    }

    private RootState FindRoot(double a_x, Func<double, double> a_func,
        Func<double, double> a_d1, Func<double, double> a_d2)
    {
        double prev_y = a_func(a_x);

        for (int i = 0; i < MaxIterations; i++)
        {
            double dy = a_d1(a_x);
            double ddy = a_d2(a_x);

            double den = 2 * dy * dy - prev_y * ddy;
            if (den.IsAlmostEquals(0))
                return new RootState(a_x);

            a_x = a_x - 2 * prev_y * dy / den;

            double y = a_func(a_x);
            prev_y = y;

            if (Math.Abs(a_x) > AbsoluteError)
            {
                if (Math.Abs(y / a_x) < AbsoluteError)
                    return new RootState(a_x);
            }
            else if (Math.Abs(y) < AbsoluteError)
                return new RootState(a_x);
        }

        return new RootState(a_x, false);
    }

    public override double FindRoot(double a_a, double a_b, 
        Polynomial a_poly)
    {
        Polynomial d1 = a_poly.Differentiate();
        Polynomial d2 = d1.Differentiate();
        return FindRoot(a_a, a_b, (x) => a_poly.Evaluate(x),
            (x) => d1.Evaluate(x), (x) => d2.Evaluate(x));
    }
}

Newton Root Finder

Nie wymaga zakresu do poszukiwania miejsca zerowego, tylko punktu startowego. Wymaga za to obliczania wartości pochodnej, czyli także ciągłości funkcji. Metodę tą ciężko zastosować jeśli nie możemy podać funkcji na pochodną. Jest to także metoda iteracyjna. Krok iteracji ma postać:

$x_{1} = x_0 - \frac{f(x_0)}{f'(x_0)}$

Metoda ta prawie na pewno nie poradzi sobie z ekstremami funkcji - zacznie wokół nich oscylować. Może się zdarzyć, że zostanie znaleziony pierwiastek z poza zakresu, gdyż w jego kierunku iteracja zacznie podążać wzdłuż krzywej.

Poniższy kod stara się poszukuje miejsca zerowego dwa razy, obierając jako punkt startu a i b zakresu poszukiwań. Tym samym zwiększamy trochę szansę znalezienia pierwiastków. Ale ciągle nie radzi sobie z ekstremami.

Kod:

public class CentralDifferenceMethod
{
    private double m_h;
    private Func<double, double> m_func;

    public CentralDifferenceMethod()
        : base()
    {
    }
        
    public CentralDifferenceMethod(Func<double, double> a_func, 
        double a_h = Constants.DOUBLE_PRECISION * 2)
    {
        m_h = a_h;
        m_func = a_func;
    }

    public double Derive1(double a_x)
    {
        return (m_func(a_x + m_h) - m_func(a_x - m_h)) / (2 * m_h);
    }

    public double Derive2(double a_x)
    {
        return (m_func(a_x + m_h) - 2 * m_func(a_x) + 
            m_func(a_x - m_h)) / (m_h * m_h);
    }

    public static  Func<double, double> Derive1(
        Func<double, double> a_func, double a_h)
    {
        CentralDifferenceMethod cdm = 
            new CentralDifferenceMethod(a_func, a_h);

        return (x) =>
        {
            return cdm.Derive1(x);
        };
    }

    public static Func<double, double> Derive2(
        Func<double, double> a_func, double a_h)
    {
        CentralDifferenceMethod cdm = 
            new CentralDifferenceMethod(a_func, a_h);

        return (x) =>
        {
            return cdm.Derive2(x);
        };
    }

    public static Func<double, double> Derive2(
        Func<double, double> a_func, 
        Func<double, double> a_d1, double a_h)
    {
        CentralDifferenceMethod cdm = 
            new CentralDifferenceMethod(a_d1, a_h);

        return (x) =>
        {
            return cdm.Derive1(x);
        };
    }
}

public abstract class RootFinder
{
    public int MaxIterations;
    public double AbsoluteError;

    public RootFinder()
    {
        MaxIterations = Constants.ROOT_FINDER_MAXIMUM_ITERATIONS;
        AbsoluteError = Constants.ROOT_FINDER_ABSOLUTE_ERROR;
    }

    public static RootFinder Create(RootFindingMethod a_method)
    {
        switch (a_method)
        {
            case RootFindingMethod.Bisection: 
                return new BisectionRootFinder();
            case RootFindingMethod.ModifedRegulaFalsi: 
                return new ModifiedRegulaFalsiRootFinder();
            case RootFindingMethod.RegulaFalsi: 
                return new RegulaFalsiRootFinder();
            case RootFindingMethod.Secant: 
                return new SecantRootFinder();
            case RootFindingMethod.Newton: 
                return new NewtonRootFinder();
            case RootFindingMethod.Halley: 
                return new HalleyRootFinder();
            default: throw new NotImplementedException();
        }
    }

    public abstract double FindRoot(double a_a, double a_b, 
        Polynomial a_poly);
    public abstract double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func);
    public abstract double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func, Func<double, double> a_d1);
    public abstract double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func, Func<double, double> a_d1, 
        Func<double, double> a_d2);
}

public abstract class DerivativeRootFinder : RootFinder
{
    protected struct RootState
    {
        public double Root;
        public bool Finded;

        public RootState(double a_x, bool a_finded = true)
        {
            Debug.Assert(a_x.IsNumber());

            Root = a_x;
            Finded = a_finded;
        }
    }

    protected double FindRootFromBothSide(double a_a, double a_b, 
        Func<double, double> a_func, Func<double, 
        RootState> a_root_finder)
    {
        Debug.Assert(a_a < a_b);
        Debug.Assert(a_a.IsNumber());
        Debug.Assert(a_b.IsNumber());

        double ya = a_func(a_a);
        double yb = a_func(a_b);

        if (Math.Abs(ya) < AbsoluteError)
            return a_a;
        if (Math.Abs(yb) < AbsoluteError)
            return a_b;

        if (Math.Sign(ya) * Math.Sign(yb) > 0)
            return Double.PositiveInfinity;

        RootState x1 = a_root_finder(a_a);
        if (x1.Finded && (x1.Root > a_a) && (x1.Root < a_b))
            return x1.Root;

        RootState x2 = a_root_finder(a_b);
        if (x2.Finded && (x2.Root > a_a) && (x2.Root < a_b))
            return x2.Root;

        return Double.PositiveInfinity;
    }
}

public class NewtonRootFinder : DerivativeRootFinder
{
    public override double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func)
    {
        return FindRoot(a_a, a_b, a_func, 
            CentralDifferenceMethod.Derive1(a_func, AbsoluteError));
    }

    public override double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func, Func<double, double> a_d1, 
        Func<double, double> a_d2)
    {
        return FindRoot(a_a, a_b, a_func, a_d1);
    }

    public override double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func, Func<double, double> a_d1)
    {
        return FindRootFromBothSide(a_a, a_b, a_func, 
            (s) => FindRoot(s, a_func, a_d1));
    }

    private RootState FindRoot(double a_x, Func<double, double> a_func,
        Func<double, double> a_d1)
    {
        double prev_y = a_func(a_x);

        for (int i = 0; i < MaxIterations; i++)
        {
            double dy = a_d1(a_x);
            if (dy.IsAlmostEquals(0))
                return new RootState(a_x);

            a_x = a_x - prev_y / dy;
            double y = a_func(a_x);
            prev_y = y;

            if (Math.Abs(a_x) > AbsoluteError)
            {
                if (Math.Abs(y / a_x) < AbsoluteError)
                    return new RootState(a_x);
            }
            else if (Math.Abs(y) < AbsoluteError)
                return new RootState(a_x);
        }

        return new RootState(a_x, false);
    }

    public override double FindRoot(double a_a, double a_b, 
        Polynomial a_poly)
    {
        Polynomial d1 = a_poly.Differentiate();
        return FindRoot(a_a, a_b, (x) => a_poly.Evaluate(x),
            (x) => d1.Evaluate(x));
    }
}

Root Finder Modifed Regula Falsi

Metoda szybsza od dwóch pierwszej Bisection. Zawsze gwarantująca rozwiązanie w przeciwieństwie do metody Secant.

Wymaga przedziału startowego (a,b) takiego, że znaki f(a) i f(b) są różne, co oznacza, że w przedziale (a,b) istnieje miejsce zerowe. By tak się stało funkcja powinna być ciągła.

Jest to metoda rekurencyjna, w której stopniowo przybliżamy miejsce zerowe. W przedziale (a,b) mogą istnieć ekstrema nie będące miejsca miejscami zerowymi albo parzyste miejsca zerowe. Metoda ta nie zawsze gwarantuje nam znalezienie rozwiązania.

W dowolnym kroku iteracji mamy przedział $(x_0, x_1)$, podobnie jak w metodzie Secant wyznaczamy linię przebiegającą przez punkty $(x_0,f(x_0))$ i $(x_1,f(x_1))$. Znajdujemy jej miejsce zerowe:

$x_2=x_1-f(x_1)\frac{x_1-x_0}{f(x_1)-f(x_0)}=\frac{f(y_1)x_0-f(y_0)x_1}{f(x_1)-f(x_0)}$

Teraz wybieramy nowy zakres tak by znaki funkcji na końcach zakresu miały przeciwny znak. Cały proces powtarzamy.

Kod:

public override double FindRoot(
    double a_a, double a_b, Func a_func)
{
    Debug.Assert(a_a < a_b);

    double ya = a_func(a_a);
    double yb = a_func(a_b);

    if (Math.Abs(ya) < Constants.ROOT_FINDER_PRECISION)
        return a_a;

    if (Math.Abs(yb) < Constants.ROOT_FINDER_PRECISION)
        return a_b;

    if (Math.Sign(ya) * Math.Sign(yb) > 0.0)
        return Double.PositiveInfinity;

    double x = 0;

    for (int i = 0; i < Constants.ROOT_FINDER_MAXIMUM_ITERATIONS; i++)
    {
        x = (yb * a_a - ya * a_b) / (yb - ya);
        double y = a_func(x);

        if (Math.Abs(x) > Constants.ROOT_FINDER_PRECISION)
        {
            if (Math.Abs(y / x) < Constants.ROOT_FINDER_PRECISION)
                return x;
        }
        else if (Math.Abs(y) < Constants.ROOT_FINDER_PRECISION)
            return x;

        if ((Math.Sign(ya) * Math.Sign(y)) < 0)
        {
            a_b = x;
            yb = y;
        }
        else
        {
            a_a = x;
            ya = y;
        }
    }

    return x;
}

Metoda ta nie zawsze gwarantuje nam zbieżność. Istnieje jej modyfikacja która zawsze gwarantuje nam zbieżność. To dzielenie przez 2 ma uzasadnienie teorii.

public override double FindRoot(
    double a_a, double a_b, Func<double, double> a_func)
{
    Debug.Assert(a_a < a_b);

    double ya = a_func(a_a);
    double yb = a_func(a_b);

    if (Math.Sign(ya) * Math.Sign(yb) > 0.0)
        return Double.PositiveInfinity;

    if (Math.Abs(ya) < Constants.ROOT_FINDER_PRECISION)
        return a_a;

    if (Math.Abs(yb) < Constants.ROOT_FINDER_PRECISION)
        return a_b;

    double prev_y = ya;
    double x = 0;

    for (int i = 0; i < Constants.ROOT_FINDER_MAXIMUM_ITERATIONS; i++)
    {
        x = (yb * a_a - ya * a_b) / (yb - ya);
        double y = a_func(x);

        if (Math.Abs(x) > Constants.ROOT_FINDER_PRECISION)
        {
            if (Math.Abs(y / x) < Constants.ROOT_FINDER_PRECISION)
                return x;
        }
        else if (Math.Abs(y) < Constants.ROOT_FINDER_PRECISION)
            return x;

        if ((Math.Sign(ya) * Math.Sign(y)) < 0)
        {
            a_b = x;
            yb = y;
            if ((prev_y * y) > 0)
                ya /= 2;
        }
        else
        {
            a_a = x;
            ya = y;
            if ((prev_y * y) > 0)
                yb /= 2;
        }

        prev_y = y;
    }

    return x;
}

Szybka transformata Fouriera Radix-2

FFT czyli dyskretna transformata Fouriera w wersji szybkiej.

Dyskretną transformatę Fouriera definiujemy jako:

$X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i 2 \pi \frac{k}{N} n}$

Definicja odwrotnej dyskretnej transformaty jest bardzo podobna:

$x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k \cdot e^{+i 2 \pi \frac{k}{N} n}$

W praktyce obie transformaty korzystają z tych samych procedur z różnie zadanymi parametrami. My będziemy tylko zajmować się FFT.

Wersja naiwna wprost z definicji:

public static Complex[] fft(Complex[] a_ar, int a_length = -1)
{
    if (a_length == -1)
        a_length = a_ar.Length;

    Complex[] result = new Complex[a_length];

    for (int k = 0; k < result.Length; k++)
    {
        for (int n = 0; n < a_ar.Length; n++)
        {
            double phase = -2 * MathExtensions.PI * k * n / a_ar.Length;
            result[k] += a_ar[n] * new Complex(Math.Cos(phase), Math.Sin(phase));
        }
    }

    return result;
}

Czas wykonania: 55s. Złożoność obliczeniowa to $O(n^2)$. Test to FFT w dwóch wymiarach.

Pierwszą rzeczą jaką możemy wykonać to zrezygnować z wyliczania za każdym razem wartości new Complex(Math.Cos(phase), Math.Sin(phase)). Zauważmy, że:

$coeff_{n+1}=e^{-i 2 \pi \frac{k}{N} (n+1)}=e^{-i 2 \pi \frac{k}{N} n}e^{-i 2 \pi \frac{k}{N} }=coeff_n \cdot e^{-i 2 \pi \frac{k}{N} }$

Implementacja w kodzie:

public static Complex[] fft(Complex[] a_ar, int a_length = -1)
{
    if (a_length == -1)
        a_length = a_ar.Length;

    Complex[] result = new Complex[a_length];

    for (int k = 0; k < result.Length; k++)
    {
        Complex step = Complex.FromPolar(-2 * MathExtensions.PI * k / a_ar.Length);
        Complex coeff = Complex.ONE;

        for (int n = 0; n < a_ar.Length; n++)
        {
            result[k] += a_ar[n] * coeff;
            coeff *= step;
        }
    }

    return result;
}

Czas wykonania: 16s. Całkiem niezłe przyspieszenie.

Czas na FFT w dwóch wymiarach. Jeśli zbiór danych potraktujemy jako tablicę dwuwymiarową, to FFT wykonujemy na rzędach, a później kolumnach, lub odwrotnie. Kod:

public static Complex[,] fft2(Complex[,] a)
{
    int w = a.GetLength(0);
    int h = a.GetLength(1);

    Complex[,] result = new Complex[w, h];

    for (int y = 0; y < h; y++)
    {
        Complex[] row = new Complex[w];

        for (int x = 0; x < w; x++)
            row[x] = a[x, y];

        Complex[] frow = fft(row);

        for (int x = 0; x < w; x++)
            result[x, y] = frow[x];
    }

    for (int x = 0; x < w; x++)
    {
        Complex[] col = new Complex[h];

        for (int y = 0; y < h; y++)
            col[y] = result[x, y];

        Complex[] fcol = fft(col);

        for (int y = 0; y < h; y++)
            result[x, y] = fcol[y];
    }

    return result;
}

Wynik zwracany przez fft(), fft2() warto przesunąć, tak by element o f=0 znajdował się w środku, a nie na początku tablicy.

public static T[,] fftshift<T>(T[,] a_ar)
{
    int deltax = (a_ar.GetLength(0) / 2.0).Floor();
    int deltay = (a_ar.GetLength(1) / 2.0).Floor();

    T[,] result = new T[a_ar.GetLength(0), a_ar.GetLength(1)];

    for (int x = 0; x < a_ar.GetLength(0); x++)
    {
        int cx = OverlapOverlayCorrector.Correct_XY(x + deltax, a_ar.GetLength(0));

        for (int y = 0; y < a_ar.GetLength(1); y++)
        {
            int cy = OverlapOverlayCorrector.Correct_XY(y + deltay, 
               a_ar.GetLength(1));
            result[cx, cy] = a_ar[x, y];
        }
    }

    return result;
}

Warto zwrócić uwagę, że dla parzystej długości tablicy, podwójne wywołanie fftshift() daje dane początkowe. Tak nie jest w przypadku tablicy o długości nieparzystej. Dla odwracalności definiujemy:

public static T[,] ifftshift<T>(T[,] a_ar)
{
    int deltax = (a_ar.GetLength(0) / 2.0).Ceiling();
    int deltay = (a_ar.GetLength(1) / 2.0).Ceiling();
    return generic_shift(a_ar, deltax, deltay);
}

Kompletny kod funkcji na przesuwanie:

public static T[] fftshift<T>(T[] a_ar)
{
    int delta = (a_ar.Length / 2.0).Floor();
    return generic_shift(a_ar, delta);
}

public static T[] ifftshift<T>(T[] a_ar)
{
    int delta = (a_ar.Length / 2.0).Ceiling();
    return generic_shift(a_ar, delta);
}

private static T[] generic_shift<T>(T[] a_ar, int delta)
{
    T[] result = new T[a_ar.Length];

    for (int i = 0; i < a_ar.Length; i++)
    {
        int ci = OverlapOverlayCorrector.Correct_XY(i + delta, a_ar.Length);
        result[ci] = a_ar[i];
    }

    return result;
}

public static T[,] fftshift<T>(T[,] a_ar)
{
    int deltax = (a_ar.GetLength(0) / 2.0).Floor();
    int deltay = (a_ar.GetLength(1) / 2.0).Floor();
    return generic_shift(a_ar, deltax, deltay);
}

public static T[,] ifftshift<T>(T[,] a_ar)
{
    int deltax = (a_ar.GetLength(0) / 2.0).Ceiling();
    int deltay = (a_ar.GetLength(1) / 2.0).Ceiling();
    return generic_shift(a_ar, deltax, deltay);
}

private static T[,] generic_shift<T>(T[,] a_ar, int deltax, int deltay)
{
    T[,] result = new T[a_ar.GetLength(0), a_ar.GetLength(1)];

    for (int x = 0; x < a_ar.GetLength(0); x++)
    {
        int cx = OverlapOverlayCorrector.Correct_XY(x + deltax, a_ar.GetLength(0));

        for (int y = 0; y < a_ar.GetLength(1); y++)
        {
            int cy = OverlapOverlayCorrector.Correct_XY(y + deltay, 
                a_ar.GetLength(1));
            result[cx, cy] = a_ar[x, y];
        }
    }

    return result;
}

Podana wcześniej wersja fft2() jest dla liczb zespolonych, tymczasem nasza tablica wejściowa to zbiór liczb rzeczywistych. Tak więc stwórzmy osobną wersję na to:

public static Complex[,] fft2(Complex[,] a)
{
    int w = a.GetLength(0);
    int h = a.GetLength(1);

    Complex[,] result = new Complex[w, h];

    for (int y = 0; y < h; y++)
    {
        Complex[] row = new Complex[w];

        for (int x = 0; x < w; x++)
            row[x] = a[x, y];

        Complex[] frow = fft(row);

        for (int x = 0; x < w; x++)
            result[x, y] = frow[x];
    }

    for (int x = 0; x < w; x++)
    {
        Complex[] col = new Complex[h];

        for (int y = 0; y < h; y++)
            col[y] = result[x, y];

        Complex[] fcol = fft(col);

        for (int y = 0; y < h; y++)
            result[x, y] = fcol[y];
    }

    return result;
}

Naiwna implementacja FFT dla liczb rzeczywistych:

public static Complex[] fft(double[] a_ar, int a_length = -1)
{
    Complex[] ar = new Complex[a_ar.Length];
    for (int i = 0; i < a_ar.Length; i++)
        ar[i] = new Complex(a_ar[i]);
    return fft(ar, a_length);
}

Dla takich funkcji jak wyżej jest zmierzony czas 15s podany wcześniej. Stwórzmy teraz wersję wyspecjalizowaną dla liczb rzeczywistych:

public static Complex[] fft(double[] a_ar, int a_length = -1)
{
    if (a_length == -1)
        a_length = a_ar.Length;

    Complex[] result = new Complex[a_length];

    for (int k = 0; k < result.Length; k++)
    {
        Complex step = Complex.FromPolar(-2 * MathExtensions.PI * k / a_ar.Length);
        Complex coeff = Complex.ONE;

        for (int n = 0; n < a_ar.Length; n++)
        {
            result[k] += a_ar[n] * coeff;
            coeff *= step;
        }
    }

    return result;
}

Jest ona identyczna z poprzednią. Czas wykonania: 14.3s.

Wersje 2D możemy zrównoleglić. Najłatwiej podział danych przeprowadzić po rzędach i kolumnach:

public static Complex[,] fft2(double[,] a)
{
    int w = a.GetLength(0);
    int h = a.GetLength(1);

    Complex[,] result = new Complex[w, h];

    Parallel.For(0, h, (y) =>
    {
        double[] row = new double[w];

        for (int x = 0; x < w; x++)
            row[x] = a[x, y];

        Complex[] frow = fft(row);

        for (int x = 0; x < w; x++)
            result[x, y] = frow[x];
    });

    Parallel.For(0, w, (x) =>
    {
        Complex[] col = new Complex[h];

        for (int y = 0; y < h; y++)
            col[y] = result[x, y];

        Complex[] fcol = fft(col);

        for (int y = 0; y < h; y++)
            result[x, y] = fcol[y];
    });

    return result;
}

Czas wykonania: 8s. Czyli skalowanie jest całkiem dobre (mam dwa rdzenie u siebie).

Dla wersji dwuwymierowej spóbujmy wstępnie policzyć nasze współczynniki coeff.

private static Complex[,] PrecalculateCoeffs(int a_in_length, int a_out_length)
{
    Complex[,] result = new Complex[a_out_length, a_in_length];

    for (int k = 0; k < a_out_length; k++)
    {
        Complex step = Complex.FromPolar(-2 * MathExtensions.PI * k / a_in_length);
        Complex coeff = Complex.ONE;

        for (int n = 0; n < a_in_length; n++)
        {
            result[k, n] = coeff;
            coeff *= step;
        }
    }

    return result;
}

Ważne jest znaczenie indeksów w tej tablicy.

Najpierw poprawmy nasze funkcje fft2 w następujący sposób, by wielkość tablicy wyjścia mogła się różnić od wejścia:

public static Complex[,] fft2(double[,] a_in, int a_width = -1, int a_height = -1)
{
    if (a_width == -1)
        a_width = a_in.GetLength(0);
    if (a_height == -1)
        a_height = a_in.GetLength(1);

    Complex[,] temp = new Complex[a_width, a_in.GetLength(1)];

    Parallel.For(0, a_in.GetLength(1), (y) =>
    {
        double[] row = new double[a_width];

        for (int x = 0; x < a_width; x++)
            row[x] = a_in[x, y];

        Complex[] frow = fft(row);

        for (int x = 0; x < a_width; x++)
            temp[x, y] = frow[x];
    });

    Complex[,] result2 = new Complex[a_width, a_height];

    Parallel.For(0, a_width, (x) =>
    {
        Complex[] col = new Complex[a_height];

        for (int y = 0; y < a_height; y++)
            col[y] = temp[x, y];

        Complex[] fcol = fft(col);

        for (int y = 0; y < a_height; y++)
            result2[x, y] = fcol[y];
    });

    return result2;
}

Wersja fft() korzystająca z prekalkulowanej tablicy:

private static Complex[] fft(Complex[] a_ar, Complex[,] a_coeffs)
{
    Complex[] result = new Complex[a_coeffs.GetLength(0)];

    for (int k = 0; k < result.Length; k++)
    {
        for (int n = 0; n < a_ar.Length; n++)
        {
            result[k] += a_ar[n] * a_coeffs[k, n];
        }
    }

    return result;
}

No i wreście nowa wersja fft2():

public static Complex[,] fft2(double[,] a_in, int a_width = -1, int a_height = -1)
{
    if (a_width == -1)
        a_width = a_in.GetLength(0);
    if (a_height == -1)
        a_height = a_in.GetLength(1);

    Complex[,] temp = new Complex[a_width, a_in.GetLength(1)];
    Complex[,] coeffs = PrecalculateCoeffs(a_in.GetLength(0), a_width);

    Parallel.For(0, a_in.GetLength(1), (y) =>
    {
        double[] row = new double[a_width];

        for (int x = 0; x < a_width; x++)
            row[x] = a_in[x, y];

        Complex[] frow = fft(row, coeffs);

        for (int x = 0; x < a_width; x++)
            temp[x, y] = frow[x];
    });

    coeffs = PrecalculateCoeffs(a_in.GetLength(1), a_height);

    Complex[,] result2 = new Complex[a_width, a_height];

    Parallel.For(0, a_width, (x) =>
    {
        Complex[] col = new Complex[a_height];

        for (int y = 0; y < a_height; y++)
            col[y] = temp[x, y];

        Complex[] fcol = fft(col, coeffs);

        for (int y = 0; y < a_height; y++)
            result2[x, y] = fcol[y];
    });

    return result2;
}

Czas wykonania: 6.42.

Po pierwsze z dużo wcześniej przeprowadzonych doświadczeń nie wynika, że tablicy Complex[][] są szybsze od Complex[,]. Po drugie zmiena prekalkulacji współczynników z rekurencyjnej na bezpośrednią nic nie zmienia w dokładności wyników. Przynajmniej w moim wypadku, kiedy FFT jest konwertowana na bitmapę odcieni szarości. Ciągłe kopiowania i tworzenie tablic w fft2() też nie stanowi problemu. W naszym przypadku to około 33ms.

Dalsze przyspieszenie to już implementacja tego co właściwie nazywamy FFT. Samo FFT ma złożonośc obliczeniową $O(n\log n)$, czyli nasz testowy przypadek powinien przyspieszyć do około 0.9s. Przełączając się na alternatywną wersję poskładaną z kodu z biblioteki AlgLib rzeczywiście tak się dzieje. Ich wersja, która stara się zachować złożoność $O(n\log n)$ niezależnie od wielkości danych wejściowych dokładnie tyle się wykonuje.

Istnieje wiele wersji FFT. Głównym problemem jest tutaj wielkość tablicy wejściowej. Samo FFT polega na rekurencyjnym podziale tablicy. Najlepiej jeśli wielkość tej tablicy to wielokrotność potęgi 2 (albo 4, 8, ...). Wtedy podział przebiega gładko, wszystko ładnie mieści się cache. Dla pozostałych przypków wielkości tablicy także istnieją metody gwarantujące zachowanie złożoności $O(n\log n)$, także dla liczb pierwszych. Istnieją także dedykowane algorytmy działające od razy w dwóch albo w większej ilości wymiarów. Ich złożoność obliczeniowa nie ulega zmianie.

My zaimplementujemy tylko algorytm wymagający by wielkość tablicy była rzędu $2^N$. Pozostałe zostawimy. Nie są one dla mnie istotne. Juź i tak implemtuje FFT z czystej ciekawości, gdyż na razie wszystko mi działa z AlgLib.

Algorytm ten nazywamy RADIX-2. Z definicji DFT:

$X_k = \sum_{n=0}^{N-1} x_n \cdot e^{\textstyle -i 2 \pi \frac{k}{N} n}$

$X_k = \sum_{m=0}^{N/2-1} x_{2m} \cdot e^{\textstyle -i \frac{2 \pi}{N} (2m+k)}$

$X_k = \sum_{m=0}^{N/2-1} x_{2m} \cdot e^{\textstyle -i \frac{2 \pi}{N} 2mk} + \sum_{m=0}^{N/2-1} x_{2m+1} \cdot e^{\textstyle -i \frac{2 \pi}{N} (2m+1)k}$

$X_k = \sum_{m=0}^{N/2-1} x_{2m} \cdot e^{\textstyle -i \frac{2 \pi}{N/2} mk} + e^{\textstyle -i \frac{2 \pi}{N} k} \sum_{m=0}^{N/2-1} x_{2m+1} \cdot e^{\textstyle -i \frac{2 \pi}{N/2} mk}$

$X_k = E_k + e^{\textstyle -i \frac{2 \pi}{N} k}O_k$

$E_k$ to DFT parzystych wyrazów ciągu, $O_k$ to DFT nieparzystych wyrazów ciągów.

Rozbiliśmy więc jedną dużą DFT na dwie małe. Cały proces możemy powtarzać, aż długość DFT osiągnie 1. Jednak z punktu prędkości obliczeń może to nie być najlepszy moment na przerwanie rekursji.

Dodatkowo zauważmy, że:

$X_{\textstyle k+\frac{N}{2}} = E_{\textstyle k+\frac{N}{2}} + e^{\textstyle -i \frac{2 \pi}{N} (k+\frac{N}{2})}O_{\textstyle k+\frac{N}{2}}$

$\begin{align*} E_{k+N/2} &= \sum_{m=0}^{N/2-1} x_{2m} \cdot e^{\textstyle -i \frac{2 \pi}{N/2} m(k+\frac{N}{2})} \\ &= \sum_{m=0}^{N/2-1} x_{2m} \cdot e^{\textstyle -i \frac{2 \pi}{N/2} mk} \cdot e^{\textstyle -i \cdot {2 \pi} m} \\ &= \sum_{m=0}^{N/2-1} x_{2m} \cdot e^{\textstyle -i \frac{2 \pi}{N/2} mk} \\ &= E_k \end{align*}$

$\begin{align*} O_{k+N/2} &= e^{\textstyle -i \frac{2 \pi}{N} (k+N/2)} \sum_{m=0}^{N/2-1} x_{2m+1} \cdot e^{\textstyle -i \frac{2 \pi}{N/2} m(k+N/2))} \\ &= e^{\textstyle-i \cdot \pi } \cdot e^{\textstyle -i \frac{2 \pi}{N} k} \sum_{m=0}^{N/2-1} x_{2m+1} \cdot e^{\textstyle -i \cdot 2\pi} \cdot e^{\textstyle -i \frac{2 \pi}{N/2} mk)} \\ &= -e^{\textstyle -i \frac{2 \pi}{N} k} \sum_{m=0}^{N/2-1} x_{2m+1} \cdot e^{\textstyle -i \frac{2 \pi}{N/2} mk)} \\ &= -e^{\textstyle -i \frac{2 \pi}{N} k} \cdot O_k \end{align*}$

$X_{\textstyle k+\frac{N}{2}} = E_k - e^{\textstyle -i \frac{2 \pi}{N} k} \cdot O_k$

Czyli tak naprawdę druga połówka DFT jest pewnym przekształceniem pierwszej i nie trzeba jej liczyć wprost.

Musiałem także zmienić rozmiar obrazka z 600x600 na 512x512. Nowy czas wykonania to 4.8s.

Wywaliłem także kod odpowiedzialny za możliwość tworzenia DFT o różnej długości. DFT jest długości tablicy wejściowej.

Korzystając z tych informacji możemy napisać rekurencyjną nową wersję fft():

public static Complex[] fft_radix2(Complex[] a_ar, int a_length, int a_from, int a_step, Complex[,] a_coeffs = null)
{
    if (a_length == 1)
        return new Complex[] { a_ar[a_from] };

    int len2 = a_length / 2;

    Complex[] odds = fft_radix2(a_ar, len2, a_from, a_step * 2, a_coeffs);
    Complex[] evens = fft_radix2(a_ar, len2, a_from + a_step, a_step * 2, a_coeffs);

    Complex[] result = new Complex[a_length];

    for (int k = 0; k < result.Length / 2; k++)
    {
        Complex even = evens[k];
        Complex odd = odds[k];

        Complex temp = Complex.FromPolar(
            -2 * MathExtensions.PI * k / a_length) * odd;
        result[k] = even + temp;
        result[k + len2] = even - temp;
    }

    return result;
}

Czas wykonania to 750ms. Zwiększamy obraz do 1024x1024. Czas wykonania to 2.8s.

Wprowadźmy rekurencyjne liczenie współczynnika temp:

public static Complex[] fft_radix2(
    Complex[] a_ar, int a_length, int a_from, int a_step, 
    Complex[,] a_coeffs = null)
{
    if (a_length == 1)
        return new Complex[] { a_ar[a_from] };

    int len2 = a_length / 2;

    Complex[] odds = fft_radix2(a_ar, len2, a_from, a_step * 2, a_coeffs);
    Complex[] evens = fft_radix2(a_ar, len2, a_from + a_step, a_step * 2, a_coeffs);

    Complex[] result = new Complex[a_length];

    Complex step = Complex.FromPolar(-2 * MathExtensions.PI / a_length);
    Complex coeff = Complex.ONE;

    for (int k = 0; k < result.Length / 2; k++)
    {
        Complex even = evens[k];
        Complex odd = odds[k];

        Complex temp = coeff * odd;
        result[k] = even + temp;
        result[k + len2] = even - temp;
        coeff *= step;
    }

    return result;
}

Czas wykonania to 2.3s.

Teraz spóbujmy ustalić rozmiar tablicy przy którym powinniśmy przerwać rekurencje. Po testach okazało się, że nic to nie daje. Zarówno jeśli wykorzystamy już istniejącą wersję fft(), lub napiszemy szybką wersję dla tablicy o wielkości 2. Tak więc zostaje jak jest.

W porównaniu do AlgLib moja wersja jest ciut szybsza. Na pewno można ją próbować przyspieszyć. Mnie to co jest zadawala i wystarcza. Nie potrzebuje szybciej ani bardziej uniwersalnie.

Pełny kod:

public static class FFT
{
    #region FFT

    public static Complex[] fft(Complex[] a_ar)
    {
        if (Bits.IsBase2(a_ar.Length))
            return fft_radix2(a_ar, a_ar.Length, 0, 1);
        else
            return fft_slow(a_ar);
    }

    public static Complex[] fft(double[] a_ar)
    {
        if (Bits.IsBase2(a_ar.Length))
            return fft_radix2(a_ar, a_ar.Length, 0, 1);
        else
            return fft_slow(a_ar);
    }

    private static Complex[] fft_radix2(Complex[] a_ar, int a_length, 
        int a_from, int a_step)
    {
        if (a_length == 1)
            return new Complex[] { a_ar[a_from] };

        int len2 = a_length / 2;

        Complex[] odds = fft_radix2(a_ar, len2, a_from, a_step * 2);
        Complex[] evens = fft_radix2(a_ar, len2, a_from + a_step, a_step * 2);

        Complex[] result = new Complex[a_length];

        Complex step = Complex.FromPolar(-2 * MathExtensions.PI / a_length);
        Complex coeff = Complex.ONE;

        for (int k = 0; k < result.Length / 2; k++)
        {
            Complex even = evens[k];

            Complex temp = coeff * odds[k];
            result[k] = even + temp;
            result[k + len2] = even - temp;
            coeff *= step;
        }

        return result;
    }

    private static Complex[] fft_radix2(double[] a_ar, int a_length, 
        int a_from, int a_step)
    {
        if (a_length == 1)
            return new Complex[] { new Complex(a_ar[a_from]) };

        int len2 = a_length / 2;

        Complex[] odds = fft_radix2(a_ar, len2, a_from, a_step * 2);
        Complex[] evens = fft_radix2(a_ar, len2, a_from + a_step, a_step * 2);

        Complex[] result = new Complex[a_length];

        Complex step = Complex.FromPolar(-2 * MathExtensions.PI / a_length);
        Complex coeff = Complex.ONE;

        for (int k = 0; k < result.Length/2; k++)
        {
            Complex even = evens[k];

            Complex temp = coeff * odds[k];
            result[k] = even + temp;
            result[k + len2] = even - temp;
            coeff *= step;
        }

        return result;
    }

    private static Complex[] fft_slow(Complex[] a_ar, Complex[,] a_coeffs = null)
    {
        if (a_coeffs == null)
        {
            Complex[] result = new Complex[a_ar.Length];

            for (int k = 0; k < result.Length; k++)
            {
                Complex step = Complex.FromPolar(
                    -2 * MathExtensions.PI * k / a_ar.Length);
                Complex coeff = Complex.ONE;

                for (int n = 0; n < a_ar.Length; n++)
                {
                    result[k] += a_ar[n] * coeff;
                    coeff *= step;
                }
            }

            return result;
        }
        else
        {
            Complex[] result = new Complex[a_ar.GetLength(0)];

            for (int k = 0; k < result.Length; k++)
            {
                for (int n = 0; n < a_ar.Length; n++)
                {
                    result[k] += a_ar[n] * a_coeffs[k, n];
                }
            }

            return result;
        }
    }

    private static Complex[] fft_slow(double[] a_ar, Complex[,] a_coeffs = null)
    {
        if (a_coeffs == null)
        {
            Complex[] result = new Complex[a_ar.Length];

            for (int k = 0; k < result.Length; k++)
            {
                Complex step = Complex.FromPolar(
                    -2 * MathExtensions.PI * k / a_ar.Length);
                Complex coeff = Complex.ONE;

                for (int n = 0; n < a_ar.Length; n++)
                {
                    result[k] += a_ar[n] * coeff;
                    coeff *= step;
                }
            }

            return result;
        }
        else
        {
            Complex[] result = new Complex[a_ar.Length];

            for (int k = 0; k < result.Length; k++)
            {
                for (int n = 0; n < a_ar.Length; n++)
                {
                    result[k] += a_ar[n] * a_coeffs[k, n];
                }
            }

            return result;
        }
    }

    #endregion

    #region FFT2

    public static Complex[,] fft2(double[,] a_in)
    {
        int w = a_in.GetLength(0);
        int h = a_in.GetLength(1);

        Complex[,] temp = new Complex[w, a_in.GetLength(1)];
        Complex[,] coeffs = null;

        if (Bits.IsBase2(w))
            PrecalculateCoeffs(a_in.GetLength(0), w);

        Parallel.For(0, a_in.GetLength(1), (y) =>
        {
            double[] row = new double[w];

            for (int x = 0; x < w; x++)
                row[x] = a_in[x, y];

            Complex[] frow;
            if (Bits.IsBase2(row.Length))
                frow = fft_radix2(row, row.Length, 0, 1);
            else
                frow = fft_slow(row, coeffs);

            for (int x = 0; x < w; x++)
                temp[x, y] = frow[x];

        });

        if (Bits.IsBase2(h))
            coeffs = PrecalculateCoeffs(a_in.GetLength(1), h);

        Complex[,] result2 = new Complex[w, h];

        Parallel.For(0, w, (x) =>
        {
            Complex[] col = new Complex[h];

            for (int y = 0; y < h; y++)
                col[y] = temp[x, y];

            Complex[] fcol;
            if (Bits.IsBase2(col.Length))
                fcol = fft_radix2(col, col.Length, 0, 1);
            else
                fcol = fft_slow(col, coeffs);

            for (int y = 0; y < h; y++)
                result2[x, y] = fcol[y];
        });

        return result2;
    }

    public static Complex[,] fft2(Complex[,] a_in)
    {
        int w = a_in.GetLength(0);
        int h = a_in.GetLength(1);

        Complex[,] temp = new Complex[w, a_in.GetLength(1)];
        Complex[,] coeffs = null;

        if (Bits.IsBase2(w))
            PrecalculateCoeffs(a_in.GetLength(0), w);

        Parallel.For(0, a_in.GetLength(1), (y) =>
        {
            Complex[] row = new Complex[w];

            for (int x = 0; x < w; x++)
                row[x] = a_in[x, y];

            Complex[] frow;
            if (Bits.IsBase2(row.Length))
                frow = fft_radix2(row, row.Length, 0, 1);
            else
                frow = fft_slow(row, coeffs);

            for (int x = 0; x < w; x++)
                temp[x, y] = frow[x];

        });

        if (Bits.IsBase2(h))
            coeffs = PrecalculateCoeffs(a_in.GetLength(1), h);

        Complex[,] result2 = new Complex[w, h];

        Parallel.For(0, w, (x) =>
        {
            Complex[] col = new Complex[h];

            for (int y = 0; y < h; y++)
                col[y] = temp[x, y];

            Complex[] fcol;
            if (Bits.IsBase2(col.Length))
                fcol = fft_radix2(col, col.Length, 0, 1);
            else
                fcol = fft_slow(col, coeffs);

            for (int y = 0; y < h; y++)
                result2[x, y] = fcol[y];
        });

        return result2;
    }

    private static Complex[,] PrecalculateCoeffs(int a_in_length, int a_out_length)
    {
        Complex[,] result = new Complex[a_out_length, a_in_length];

        for (int k = 0; k < a_out_length; k++)
        {
            Complex step = Complex.FromPolar(
                -2 * MathExtensions.PI * k / a_in_length);
            Complex coeff = Complex.ONE;

            for (int n = 0; n < a_in_length; n++)
            {
                result[k, n] = coeff;
                coeff *= step;
            }
        }

        return result;
    }

    #endregion

    #region FFT SHIFT

    public static T[] fftshift<T>(T[] a_ar)
    {
        int delta = (a_ar.Length / 2.0).Floor();
        return generic_shift(a_ar, delta);
    }

    public static T[] ifftshift<T>(T[] a_ar)
    {
        int delta = (a_ar.Length / 2.0).Ceiling();
        return generic_shift(a_ar, delta);
    }

    private static T[] generic_shift<T>(T[] a_ar, int delta)
    {
        T[] result = new T[a_ar.Length];

        for (int i = 0; i < a_ar.Length; i++)
        {
            int ci = OverlapOverlayCorrector.Correct_XY(i + delta, a_ar.Length);
            result[ci] = a_ar[i];
        }

        return result;
    }

    public static T[,] fftshift<T>(T[,] a_ar)
    {
        int deltax = (a_ar.GetLength(0) / 2.0).Floor();
        int deltay = (a_ar.GetLength(1) / 2.0).Floor();
        return generic_shift(a_ar, deltax, deltay);
    }

    public static T[,] ifftshift<T>(T[,] a_ar)
    {
        int deltax = (a_ar.GetLength(0) / 2.0).Ceiling();
        int deltay = (a_ar.GetLength(1) / 2.0).Ceiling();
        return generic_shift(a_ar, deltax, deltay);
    }

    private static T[,] generic_shift<T>(T[,] a_ar, int deltax, int deltay)
    {
        T[,] result = new T[a_ar.GetLength(0), a_ar.GetLength(1)];

        for (int x = 0; x < a_ar.GetLength(0); x++)
        {
            int cx = OverlapOverlayCorrector.Correct_XY(
                x + deltax, a_ar.GetLength(0));

            for (int y = 0; y < a_ar.GetLength(1); y++)
            {
                int cy = OverlapOverlayCorrector.Correct_XY(
                    y + deltay, a_ar.GetLength(1));
                result[cx, cy] = a_ar[x, y];
            }
        }

        return result;
    }
    #endregion
}

Bisection Root Finder

Jedna z prostszych metod do znajdowania miejsc zerowych funkcji. Miejsce zerowe poszukiwane jest w zadanym obszarze (a,b). Wymagane jest (jak zresztą dla większości metod znajdowania miejsc zerowych) by $sign(f(a)) \ne sign(f(b))$. Wtedy też mamy pewność, że o ile funkcja jest ciągła (różniczkowalna) to istnieje miejsce zerowe (jedno lub więcej). Metoda ta podobnie jak inne znajduje tylko jedno miejsce zerowe. Ponieważ jest to metoda numeryczna miejsce zerowe znajdowane jest tylko z pewnym przybliżeniem - precyzją, która jest parametrem metody. Metoda ta potrafi radzić sobie z lokalnymi ekstremami w zakresie (a,b), czego nie potrafią metody oparte na podążaniu wzdłuż pochodnej.

Metoda ta to nic innego jak zmodyfikowane poszukiwanie binarne. Dzięki temu możemy z góry określić ile maksymalnie potrzebujemy kroków by dojść do zakładanej precyzji. W dowolnym momencie mamy trzy punkty $(y_1,y_2=\frac{y_1+y+2}{2},y_3)$. Znaki $y_1,y_3$ są przeciwne. W zależności od znaku $y_2$ wybieramy taki nowy zakres poszukiwań, by granice poszukiwań miały przeciwne znaki.

Jeśli stworzymy klasę która jako parametr będzie wymagała funktora (wyrażenie lambda) to otrzymamy dość uniwersalne narzędzie. Oprócz znajdowania pierwiastków wielomianów w zakresach wstępnie przybliżonych np. metodą Sturma możemy ją na przykład wykorzystać do znajdowania maksymalnego stopnia kompresji by wielkość pliku była poniżej np. 100KB.

Wszystkie metody mają różne tempo zbieżności tzn. ile iteracji jest potrzebnych by przybliżyć miejsce zerowe z zadaną precyzją. Może to przekładać się na prędkość danej metody ale nie musi (drugi czynnik to czas spędzany w pojedynczym kroku)

Przykład implementacji w C#:

public class BisectionRootFinder : RootFinder
{
    public override double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func)
    {
        Debug.Assert(a_a < a_b);

        double ya = a_func(a_a);
        double yb = a_func(a_b);

        if (Math.Abs(ya) < Constants.ROOT_FINDER_PRECISION)
            return a_a;

        if (Math.Abs(yb) < Constants.ROOT_FINDER_PRECISION)
            return a_b;

        if (Math.Sign(ya) * Math.Sign(yb) > 0.0)
            return Double.PositiveInfinity;

        double x = 0;

        for (int i = 0; i < Constants.ROOT_FINDER_MAXIMUM_ITERATIONS; i++)
        {
            x = (a_a + a_b) / 2;
            double y = a_func(x);

            if (Math.Abs(x) > Constants.ROOT_FINDER_PRECISION)
            {
                if (Math.Abs(y / x) < Constants.ROOT_FINDER_PRECISION)
                    return x;
            }
            else if (Math.Abs(y) < Constants.ROOT_FINDER_PRECISION)
                return x;

            if ((Math.Sign(ya) * Math.Sign(y)) < 0)
            {
                a_b = x;
                yb = y;
            }
            else
            {
                a_a = x;
                ya = y;
            }
        }

        return x;
    }
}

Implementacja klasy wielomian

Oto moja implementacja klasy wielomian, umożliwia znajdowanie miejsc zerowych metodą bezpośrednią i metodą Sturma.

public struct Polynomial
{
    public static readonly Polynomial ZERO = new Polynomial(0.0);
    public static readonly Polynomial ONE = new Polynomial(1.0);

    private int m_order;
    private readonly double[] m_coeffs;

    public Polynomial(Polynomial a_poly)
    {
        m_coeffs = new double[a_poly.m_order + 1];
        m_order = a_poly.m_order;

        Array.Copy(a_poly.m_coeffs, m_coeffs, m_coeffs.Length);
    }

    private Polynomial(int a_order)
    {
        Debug.Assert(a_order >= 0);

        m_order = a_order;
        m_coeffs = new double[m_order + 1];
                
        Check();
    }

    public Polynomial(double a_c0)
    {
        m_coeffs = new double[1];
        m_order = 0;
        m_coeffs[0] = a_c0;

        Check();
    }

    public Polynomial(double a_c0, double a_c1)
    {
        m_coeffs = new double[2];
        m_order = 1;
        m_coeffs[0] = a_c0;
        m_coeffs[1] = a_c1;

        OptimizeOrder();

        Check();
    }

    public Polynomial(double a_c0, double a_c1, double a_c2)
    {
        m_coeffs = new double[3];
        m_order = 2;
        m_coeffs[0] = a_c0;
        m_coeffs[1] = a_c1;
        m_coeffs[2] = a_c2;

        OptimizeOrder();

        Check();
    }

    public Polynomial(double a_c0, double a_c1, double a_c2, double a_c3)
    {
        m_coeffs = new double[4];
        m_order = 3;
        m_coeffs[0] = a_c0;
        m_coeffs[1] = a_c1;
        m_coeffs[2] = a_c2;
        m_coeffs[3] = a_c3;

        OptimizeOrder();

        Check();
    }

    public Polynomial(double a_c0, double a_c1, double a_c2, double a_c3, 
        double a_c4)
    {
        m_coeffs = new double[5];
        m_order = 4;
        m_coeffs[0] = a_c0;
        m_coeffs[1] = a_c1;
        m_coeffs[2] = a_c2;
        m_coeffs[3] = a_c3;
        m_coeffs[4] = a_c4;

        OptimizeOrder();

        Check();
    }

    public Polynomial(double a_c0, double a_c1, double a_c2, double a_c3, 
        double a_c4, double a_c5)
    {
        m_coeffs = new double[6];
        m_order = 5;
        m_coeffs[0] = a_c0;
        m_coeffs[1] = a_c1;
        m_coeffs[2] = a_c2;
        m_coeffs[3] = a_c3;
        m_coeffs[4] = a_c4;
        m_coeffs[5] = a_c5;

        OptimizeOrder();

        Check();
    }

    public Polynomial(double a_c0, double a_c1, double a_c2, double a_c3, 
        double a_c4, double a_c5, double a_c6)
    {
        m_coeffs = new double[7];
        m_order = 6;
        m_coeffs[0] = a_c0;
        m_coeffs[1] = a_c1;
        m_coeffs[2] = a_c2;
        m_coeffs[3] = a_c3;
        m_coeffs[4] = a_c4;
        m_coeffs[5] = a_c5;
        m_coeffs[6] = a_c6;

        OptimizeOrder();

        Check();
    }

    public Polynomial(params double[] a_coeffs)
    {
        Debug.Assert(a_coeffs.Length > 0);

        m_coeffs = new double[a_coeffs.Length];
        m_order = m_coeffs.Length - 1;
        Array.Copy(a_coeffs, m_coeffs, m_coeffs.Length);

        OptimizeOrder();

        Check();
    }

    public override bool Equals(object a_obj)
    {
        if (a_obj == null)
            return false;

        if (!typeof(Polynomial).Equals(a_obj.GetType()))
            return false;
        Polynomial poly = (Polynomial)a_obj;

        if (m_order != poly.m_order)
            return false;

        for (int i = 0; i <= m_order; i++)
        {
            if (m_coeffs[i] != poly.m_coeffs[i])
                return false;
        }

        return true;
    }

    public bool Equals(Polynomial a_poly)
    {
        if (m_order != a_poly.m_order)
            return false;

        for (int i = 0; i <= m_order; i++)
        {
            if (m_coeffs[i] != a_poly.m_coeffs[i])
                return false;
        }

        return true;
    }

    public bool IsAlmostRelativeEquals(Polynomial a_poly, 
        double a_precision = Constants.DOUBLE_PRECISION)
    {
        if (a_poly.m_order > m_order)
        {
            for (int i = m_order + 1; i <= a_poly.m_order; i++)
            {
                if (!a_poly.m_coeffs[i].IsAlmostRelativeEquals(0, a_precision))
                    return false;
            }
        }

        if (m_order > a_poly.m_order)
        {
            for (int i = a_poly.m_order + 1; i <= m_order; i++)
            {
                if (!m_coeffs[i].IsAlmostRelativeEquals(0, a_precision))
                    return false;
            }
        }

        for (int i = 0; i <= Math.Min(m_order, a_poly.m_order); i++)
        {
            if (!m_coeffs[i].IsAlmostRelativeEquals(a_poly.m_coeffs[i],
                a_precision))
            {
                return false;
            }
        }

        return true;
    }

    public bool IsAlmostEquals(Polynomial a_poly, 
        double a_precision = Constants.DOUBLE_PRECISION)
    {
        if (a_poly.m_order > m_order)
        {
            for (int i = m_order + 1; i <= a_poly.m_order; i++)
            {
                if (!a_poly.m_coeffs[i].IsAlmostEquals(0, a_precision))
                    return false;
            }
        }

        if (m_order > a_poly.m_order)
        {
            for (int i = a_poly.m_order + 1; i <= m_order; i++)
            {
                if (!m_coeffs[i].IsAlmostEquals(0, a_precision))
                    return false;
            }
        }

        for (int i = 0; i <= Math.Min(m_order, a_poly.m_order); i++)
        {
            if (!m_coeffs[i].IsAlmostEquals(a_poly.m_coeffs[i], a_precision))
                return false;
        }

        return true;
    }

    //

    public override int GetHashCode()
    {
        return ArrayExtensions.GetHashCode(m_coeffs);
    }

    public static bool operator !=(Polynomial a_a, Polynomial a_b)
    {
        if (a_a.m_order != a_b.m_order)
            return true;

        for (int i = 0; i <= a_a.m_order; i++)
        {
            if (a_a.m_coeffs[i] != a_b.m_coeffs[i])
                return true;
        }

        return false;
    }

    public static bool operator ==(Polynomial a_a, Polynomial a_b)
    {
        if (a_a.m_order != a_b.m_order)
            return false;

        for (int i = 0; i <= a_a.m_order; i++)
        {
            if (a_a.m_coeffs[i] != a_b.m_coeffs[i])
                return false;
        }

        return true;
    }

    [Conditional("DEBUG")]
    private void Check()
    {
        for (int i = 0; i <= m_order; i++)
            Debug.Assert(m_coeffs[i].IsNumber());

        Debug.Assert(m_order >= 0);
        Debug.Assert(m_coeffs.Length >= 1);
    }

    public int Order
    {
        get
        {
            return m_order;
        }
    }

    public double this[int a_index]
    {
        get
        {
            if (a_index > m_order)
                throw new IndexOutOfRangeException();

            return m_coeffs[a_index];
        }
        private set
        {
            if (a_index > m_order)
                throw new IndexOutOfRangeException();

            m_coeffs[a_index] = value;

            OptimizeOrder();

            Check();
        }
    }

    public double LeadCoefficient
    {
        get
        {
            return m_coeffs[m_order];
        }
    }

    public double LastCoefficient
    {
        get
        {
            return m_coeffs[0];
        }
    }

    public Polynomial Normalize()
    {
        if (m_order == 0)
            return Polynomial.ONE;

        double div = m_coeffs[m_order];

        if (div == 1)
            return this;

        Polynomial result = new Polynomial(this);

        for (int i = 0; i < result.m_order; i++)
            result.m_coeffs[i] /= div;

        result.m_coeffs[m_order] = 1;

        return result;  
    }

    public double Evaluate(double a_x)
    {
        double r = m_coeffs[m_order];

        for (int i = m_order - 1; i >= 0; i--)
            r = m_coeffs[i] + r * a_x;

        return r;
    }

    private void OptimizeOrder()
    {
        while (m_order > 0 && (m_coeffs[m_order] == 0.0))
            m_order--;
    }

    public static Polynomial operator * (Polynomial a_left, double a_right)
    {
        Polynomial result = new Polynomial(a_left);

        for (int i = 0; i <= result.m_order; i++)
            result.m_coeffs[i] *= a_right;

        result.OptimizeOrder();

        return result;
    }

    public static Polynomial operator * (double a_left, Polynomial a_right)
    {
        Polynomial result = new Polynomial(a_right);

        for (int i = 0; i <= result.m_order; i++)
            result.m_coeffs[i] *= a_left;

        result.OptimizeOrder();

        return result;
    }

    public static Polynomial operator /(Polynomial a_left, double a_right)
    {
        Polynomial result = new Polynomial(a_left);

        for (int i = 0; i <= result.m_order; i++)
            result.m_coeffs[i] /= a_right;

        result.OptimizeOrder();

        return result;
    }

    public Polynomial Differentiate()
    {
        if (m_order == 0)
            return Polynomial.ZERO;

        Polynomial result = new Polynomial(m_order - 1);

        for (int i = 1; i <= m_order; i++)
            result.m_coeffs[i - 1] += i * m_coeffs[i];

        return result;
    }

    public static Polynomial operator *(Polynomial a_u, Polynomial a_v)
    {
        Polynomial result = new Polynomial(a_u.m_order + a_v.m_order);

        for (int i = 0; i <= a_u.m_order; i++)
            for (int j = 0; j <= a_v.m_order; j++)
                result.m_coeffs[i + j] += a_u.m_coeffs[i] * a_v.m_coeffs[j];

        result.OptimizeOrder();

        return result;
    }

    public static Polynomial operator +(Polynomial a_u, Polynomial a_v)
    {
        if (a_u.m_order == a_v.m_order)
        {
            Polynomial result = new Polynomial(a_u.m_order);

            for (int i = 0; i <= a_u.m_order; i++)
                result.m_coeffs[i] = a_u.m_coeffs[i] + a_v.m_coeffs[i];

            result.OptimizeOrder();

            return result;
        }
        else if (a_u.m_order < a_v.m_order)
        {
            Polynomial result = new Polynomial(a_v.m_order);

            for (int i = a_u.m_order + 1; i <= a_v.m_order; i++)
                result.m_coeffs[i] = a_v.m_coeffs[i];

            for (int i = 0; i <= a_u.m_order; i++)
                result.m_coeffs[i] = a_u.m_coeffs[i] + a_v.m_coeffs[i];

            result.OptimizeOrder();

            return result;
        }
        else
        {
            Polynomial result = new Polynomial(a_u.m_order);

            for (int i = a_v.m_order + 1; i <= a_u.m_order; i++)
                result.m_coeffs[i] = a_u.m_coeffs[i];

            for (int i = 0; i <= a_v.m_order; i++)
                result.m_coeffs[i] = a_u.m_coeffs[i] + a_v.m_coeffs[i];

            result.OptimizeOrder();

            return result;
        }
    }

    public static Polynomial operator -(Polynomial a_poly)
    {
        Polynomial result = new Polynomial(a_poly.m_order);

        for (int i = 0; i <= result.m_order; i++)
            result.m_coeffs[i] = -a_poly.m_coeffs[i];

        return result;
    }

    public static Polynomial operator -(Polynomial a_u, Polynomial a_v)
    {
        if (a_u.m_order == a_v.m_order)
        {
            Polynomial result = new Polynomial(a_u.m_order);

            for (int i = 0; i <= a_u.m_order; i++)
                result.m_coeffs[i] = a_u.m_coeffs[i] - a_v.m_coeffs[i];

            result.OptimizeOrder();

            return result;
        }
        else if (a_u.m_order < a_v.m_order)
        {
            Polynomial result = new Polynomial(a_v.m_order);

            for (int i = a_u.m_order + 1; i <= a_v.m_order; i++)
                result.m_coeffs[i] = -a_v.m_coeffs[i];

            for (int i = 0; i <= a_u.m_order; i++)
                result.m_coeffs[i] = a_u.m_coeffs[i] - a_v.m_coeffs[i];

            result.OptimizeOrder();

            return result;
        }
        else
        {
            Polynomial result = new Polynomial(a_u.m_order);

            for (int i = a_v.m_order + 1; i <= a_u.m_order; i++)
                result.m_coeffs[i] = a_u.m_coeffs[i];

            for (int i = 0; i <= a_v.m_order; i++)
                result.m_coeffs[i] = a_u.m_coeffs[i] - a_v.m_coeffs[i];

            result.OptimizeOrder();

            return result;
        }
    }

    public static Polynomial operator %(Polynomial a_u, Polynomial a_v)
    {
        if (a_u.m_order < a_v.m_order)
            return a_u;

        Polynomial reminder = new Polynomial(a_u);

        while (reminder.m_order >= a_v.m_order)
        {
            double d = reminder.LeadCoefficient / a_v.LeadCoefficient;

            for (int j = 0; j <= a_v.m_order; j++)
            {
                reminder.m_coeffs[reminder.m_order - a_v.m_order + j] -=
                    a_v.m_coeffs[j] * d;
            }

            if (reminder.m_order > 0)
                reminder.m_order--;
            else
                break;
        }

        reminder.OptimizeOrder();
        reminder.Check();

        return reminder;
    }

    public static Polynomial operator /(Polynomial a_u, Polynomial a_v)
    {
        Polynomial reminder = new Polynomial(a_u);

        if (a_u.m_order < a_v.m_order)
            return Polynomial.ZERO;

        Polynomial result = new Polynomial(a_u.m_order - a_v.m_order);

        while (reminder.m_order >= a_v.m_order)
        {
            double d = reminder.LeadCoefficient / a_v.LeadCoefficient;
            result.m_coeffs[reminder.m_order - a_v.m_order] = d;

            for (int j = 0; j <= a_v.m_order; j++)
            {
                reminder.m_coeffs[reminder.m_order - a_v.m_order + j] -=
                    a_v.m_coeffs[j] * d;
            }

            reminder.m_order--;
        }

        result.OptimizeOrder();
        result.Check();

        return result;
    }

    public static Polynomial GCD(Polynomial a_a, Polynomial a_b, 
        double a_precision = Constants.DOUBLE_PRECISION)
    {
        if (a_a.Order < a_b.Order)
        {
            Polynomial t = a_b;
            a_b = a_a;
            a_a = t;
        }

        while (!a_b.IsAlmostEquals(Polynomial.ZERO, a_precision))
        {
            Polynomial t = a_b;

            a_b = a_a % a_b;
            a_a = t;
        }

        if (a_a.LeadCoefficient < 0)
            a_a = -a_a;

        return a_a.Normalize();
    }

    public static Polynomial DivRem(Polynomial a_divider, Polynomial a_divisor, 
        out Polynomial a_reminder)
    {
        a_reminder = new Polynomial(a_divider);

        if (a_divider.m_order < a_divisor.m_order)
            return Polynomial.ZERO;

        Polynomial result = new Polynomial(a_reminder.m_order - a_divisor.m_order);

        while (a_reminder.m_order >= a_divisor.m_order)
        {
            double d = a_reminder.LeadCoefficient / a_divisor.LeadCoefficient;
            result.m_coeffs[a_reminder.m_order - a_divisor.m_order] = d;

            for (int j = 0; j <= a_divisor.m_order; j++)
            {
                a_reminder.m_coeffs[a_reminder.m_order - a_divisor.m_order + j] -=
                    a_divisor.m_coeffs[j] * d;
            }

            if (a_reminder.m_order > 0)
                a_reminder.m_order--;
            else
                break;
        }

        result.OptimizeOrder();
        result.Check();

        a_reminder.OptimizeOrder();
        a_reminder.Check();

        return result;
    }

    public override string ToString()
    {
        string coeffs = "";

        for (int i = 0; i <= m_order; i++)
        {
            coeffs += m_coeffs[i].ToString(CultureInfo.InvariantCulture);

            if (i != m_order)
                coeffs += ", ";
        }

        return String.Format("Order: {0}, Coeffs: {1}", m_order, coeffs);
    }

    public static Polynomial operator ^(Polynomial a_left, int a_right)
    {
        if (a_right < 0)
        {
            Debug.Fail("");
            return Polynomial.ZERO;
        }

        if (a_right == 0)
            return Polynomial.ONE;

        Polynomial result = new Polynomial(a_left);

        while (a_right > 1)
        {
            a_right--;
            result *= a_left;
        }

        return result;
    }

    public Polynomial RemoveZeroes()
    {
        Polynomial result = new Polynomial(this);

        while (result.m_order > 0 && (result.m_coeffs[0].IsAlmostEquals(0)))
        {
            for (int i = 1; i <= result.m_order; i++)
                result[i - 1] = result[i];

            result.m_order--;
        }

        return result;
    }

    public IEnumerable<double> Coeffs
    {
        get
        {
            return m_coeffs;
        }
    }
}

I testy:

[TestMethod]
public void Polynomial_Test()
{
    Polynomial p_1 = new Polynomial(-1.0, 1.0) * 2;
    Polynomial p_2 = new Polynomial(-2.0, 1.0) * -0.3;
    Polynomial p_3 = new Polynomial(-3.0, 1.0) * -3;
    Polynomial p_4 = new Polynomial(-4.0, 1.0) * 0.2;
    Polynomial p_5 = new Polynomial(-5.0, 1.0) * 1.2;
    Polynomial p_6 = new Polynomial(-6.0, 1.0) * -3.2;

    // Constructors.

    {
        Polynomial p = new Polynomial();
        Assert.IsTrue(p.Order == 0);
    }

    {
        Polynomial p = new Polynomial(1.2);
        Assert.IsTrue(p.Order == 0);
        Assert.IsTrue(p[0].IsAlmostEquals(1.2));
        Assert.IsTrue(p.LeadCoefficient.IsAlmostEquals(1.2));
        Assert.IsTrue(p.LastCoefficient.IsAlmostEquals(1.2));
    }

    {
        Polynomial p = new Polynomial(1.3, 1.2);
        Assert.IsTrue(p.Order == 1);
        Assert.IsTrue(p[1].IsAlmostEquals(1.2));
        Assert.IsTrue(p[0].IsAlmostEquals(1.3));
        Assert.IsTrue(p.LeadCoefficient.IsAlmostEquals(1.2));
        Assert.IsTrue(p.LastCoefficient.IsAlmostEquals(1.3));
    }

    {
        Polynomial p = new Polynomial(1.4, 1.3, 1.2);
        Assert.IsTrue(p.Order == 2);
        Assert.IsTrue(p[2].IsAlmostEquals(1.2));
        Assert.IsTrue(p[1].IsAlmostEquals(1.3));
        Assert.IsTrue(p[0].IsAlmostEquals(1.4));
        Assert.IsTrue(p.LeadCoefficient.IsAlmostEquals(1.2));
        Assert.IsTrue(p.LastCoefficient.IsAlmostEquals(1.4));
    }

    {
        Polynomial p = new Polynomial(1.5, 1.4, 1.3, 1.2);
        Assert.IsTrue(p.Order == 3);
        Assert.IsTrue(p[3].IsAlmostEquals(1.2));
        Assert.IsTrue(p[2].IsAlmostEquals(1.3));
        Assert.IsTrue(p[1].IsAlmostEquals(1.4));
        Assert.IsTrue(p[0].IsAlmostEquals(1.5));
        Assert.IsTrue(p.LeadCoefficient.IsAlmostEquals(1.2));
        Assert.IsTrue(p.LastCoefficient.IsAlmostEquals(1.5));
    }

    {
        Polynomial p = new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2);
        Assert.IsTrue(p.Order == 4);
        Assert.IsTrue(p[4].IsAlmostEquals(1.2));
        Assert.IsTrue(p[3].IsAlmostEquals(1.3));
        Assert.IsTrue(p[2].IsAlmostEquals(1.4));
        Assert.IsTrue(p[1].IsAlmostEquals(1.5));
        Assert.IsTrue(p[0].IsAlmostEquals(1.6));
        Assert.IsTrue(p.LeadCoefficient.IsAlmostEquals(1.2));
        Assert.IsTrue(p.LastCoefficient.IsAlmostEquals(1.6));
    }

    {
        Polynomial p = new Polynomial(1.7, 1.6, 1.5, 1.4, 1.3, 1.2, 1.21, 1.22);
        Assert.IsTrue(p.Order == 7);
        Assert.IsTrue(p[7].IsAlmostEquals(1.22));
        Assert.IsTrue(p[6].IsAlmostEquals(1.21));
        Assert.IsTrue(p[5].IsAlmostEquals(1.2));
        Assert.IsTrue(p[4].IsAlmostEquals(1.3));
        Assert.IsTrue(p[3].IsAlmostEquals(1.4));
        Assert.IsTrue(p[2].IsAlmostEquals(1.5));
        Assert.IsTrue(p[1].IsAlmostEquals(1.6));
        Assert.IsTrue(p[0].IsAlmostEquals(1.7));
        Assert.IsTrue(p.LeadCoefficient.IsAlmostEquals(1.22));
        Assert.IsTrue(p.LastCoefficient.IsAlmostEquals(1.7));
    }

    {
        Polynomial p = new Polynomial(new double[] { 1.7, 1.6, 1.5, 1.4, 1.3, 1.2 });
        Assert.IsTrue(p.Order == 5);
        Assert.IsTrue(p[5].IsAlmostEquals(1.2));
        Assert.IsTrue(p[4].IsAlmostEquals(1.3));
        Assert.IsTrue(p[3].IsAlmostEquals(1.4));
        Assert.IsTrue(p[2].IsAlmostEquals(1.5));
        Assert.IsTrue(p[1].IsAlmostEquals(1.6));
        Assert.IsTrue(p[0].IsAlmostEquals(1.7));
        Assert.IsTrue(p.LeadCoefficient.IsAlmostEquals(1.2));
        Assert.IsTrue(p.LastCoefficient.IsAlmostEquals(1.7));
    }

    {
        Polynomial p = new Polynomial(1.6, 1.5, 1.4, 1.3, 0);
        Assert.IsTrue(p.Order == 3);
        Assert.IsTrue(p[3].IsAlmostEquals(1.3));
        Assert.IsTrue(p[2].IsAlmostEquals(1.4));
        Assert.IsTrue(p[1].IsAlmostEquals(1.5));
        Assert.IsTrue(p[0].IsAlmostEquals(1.6));
        Assert.IsTrue(p.LeadCoefficient.IsAlmostEquals(1.3));
        Assert.IsTrue(p.LastCoefficient.IsAlmostEquals(1.6));
    }

    {
        Polynomial p = new Polynomial(0.0);
        Assert.IsTrue(p.Order == 0);
        Assert.IsTrue(p[0].IsAlmostEquals(0.0));
        Assert.IsTrue(p.LeadCoefficient.IsAlmostEquals(0.0));
        Assert.IsTrue(p.LastCoefficient.IsAlmostEquals(0.0));
    }

    // Equals(obj), Equals(Polynomial), ==, !=

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, -1.4, 1.3, 1.2);
        Polynomial p2 = new Polynomial(1.6, 1.5, -1.4, 1.3, 1.2);

        Assert.IsTrue(p1 == p2);
        Assert.IsTrue(p2 == p1);
        Assert.IsFalse(p1 != p2);
        Assert.IsFalse(p2 != p1);
        Assert.AreEqual(p1, p2);
        Assert.AreEqual(p2, p1);
        Assert.AreEqual((object)p1, (object)p2);
        Assert.AreEqual((object)p2, (object)p1);

        Assert.AreEqual(new Polynomial(1.6, 1.5, -1.4, 1.3, 1.2), p1);
        Assert.AreEqual(new Polynomial(1.6, 1.5, -1.4, 1.3, 1.2), p2);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, 1.4, 1.3, -1.21);
        Polynomial p2 = new Polynomial(1.6, 1.5, 1.4, 1.3, -1.2);

        Assert.IsFalse(p1 == p2);
        Assert.IsFalse(p2 == p1);
        Assert.IsTrue(p1 != p2);
        Assert.IsTrue(p2 != p1);
        Assert.AreNotEqual(p1, p2);
        Assert.AreNotEqual(p2, p1);
        Assert.AreNotEqual((object)p1, (object)p2);
        Assert.AreNotEqual((object)p2, (object)p1);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, 1.4, 1.3);
        Polynomial p2 = new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2);

        Assert.IsFalse(p1 == p2);
        Assert.IsFalse(p2 == p1);
        Assert.IsTrue(p1 != p2);
        Assert.IsTrue(p2 != p1);
        Assert.AreNotEqual(p1, p2);
        Assert.AreNotEqual(p2, p1);
        Assert.AreNotEqual((object)p1, (object)p2);
        Assert.AreNotEqual((object)p2, (object)p1);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, -1.4, -1.3);
        Polynomial p2 = new Polynomial(1.6, 1.5, -1.4, -1.3, 0);

        Assert.IsTrue(p1 == p2);
        Assert.IsTrue(p2 == p1);
        Assert.IsFalse(p1 != p2);
        Assert.IsFalse(p2 != p1);
        Assert.AreEqual(p1, p2);
        Assert.AreEqual(p2, p1);
        Assert.AreEqual((object)p1, (object)p2);
        Assert.AreEqual((object)p2, (object)p1);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2);
        Polynomial p2 = new Polynomial(1.6 / 1.2, 1.5 / 1.2, 1.4 / 1.2, 1.3 / 1.2, 1);

        Assert.IsFalse(p1 == p2);
        Assert.IsFalse(p2 == p1);
        Assert.IsTrue(p1 != p2);
        Assert.IsTrue(p2 != p1);
        Assert.AreNotEqual(p1, p2);
        Assert.AreNotEqual(p2, p1);
        Assert.AreNotEqual((object)p1, (object)p2);
        Assert.AreNotEqual((object)p2, (object)p1);
    }

    // Copying constructor.

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, -1.4, 1.3, 1.2);
        Polynomial p2 = new Polynomial(p1);

        Assert.AreEqual(p1, p2);
        Assert.AreEqual(new Polynomial(1.6, 1.5, -1.4, 1.3, 1.2), p1);
    }

    // IsAlmostRelativeEquals

    {
        Polynomial p1 = new Polynomial(0);
        Polynomial p2 = new Polynomial(1);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == false);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(0);
        Polynomial p2 = new Polynomial(0);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(0);
        Polynomial p2 = new Polynomial(Constants.DOUBLE_PRECISION / 2);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(0);
        Polynomial p2 = new Polynomial(Constants.DOUBLE_PRECISION * 2);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == false);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1);
        Polynomial p2 = new Polynomial(Constants.DOUBLE_PRECISION / 2);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == false);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1);
        Polynomial p2 = new Polynomial(Constants.DOUBLE_PRECISION * 2);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == false);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1);
        Polynomial p2 = new Polynomial(1 + Constants.DOUBLE_PRECISION / 2);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1);
        Polynomial p2 = new Polynomial(1 + Constants.DOUBLE_PRECISION * 2);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == false);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1, 0);
        Polynomial p2 = new Polynomial(1);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1, 0, 0);
        Polynomial p2 = new Polynomial(1);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1, 0, 0);
        Polynomial p2 = new Polynomial(1, 0);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1, 0, 0);
        Polynomial p2 = new Polynomial(1, 0, 0);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(2, 3);
        Polynomial p2 = new Polynomial(2, 3);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(2, 3, 0);
        Polynomial p2 = new Polynomial(2, 3);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(2, 3, 0, 0);
        Polynomial p2 = new Polynomial(2, 3);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(2, 3, 0, 0);
        Polynomial p2 = new Polynomial(2, 3);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(2, 3, 0, Constants.DOUBLE_PRECISION / 2);
        Polynomial p2 = new Polynomial(2, 3);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(2, 3, 0, Constants.DOUBLE_PRECISION * 2);
        Polynomial p2 = new Polynomial(2, 3);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == false);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1 / Constants.DOUBLE_PRECISION, 3);
        Polynomial p2 = new Polynomial(1 / Constants.DOUBLE_PRECISION - 0.5, 3);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1 / Constants.DOUBLE_PRECISION, 3);
        Polynomial p2 = new Polynomial(1 / Constants.DOUBLE_PRECISION - 2, 3);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == false);
        Assert.IsTrue(p1.IsAlmostRelativeEquals(p2) == p2.IsAlmostRelativeEquals(p1));
    }

    // IsAlmostEquals.

    {
        Polynomial p1 = new Polynomial(0);
        Polynomial p2 = new Polynomial(1);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == false);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(0);
        Polynomial p2 = new Polynomial(0);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(0);
        Polynomial p2 = new Polynomial(Constants.DOUBLE_PRECISION / 2);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(0);
        Polynomial p2 = new Polynomial(Constants.DOUBLE_PRECISION * 2);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == false);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1);
        Polynomial p2 = new Polynomial(Constants.DOUBLE_PRECISION / 2);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == false);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1);
        Polynomial p2 = new Polynomial(Constants.DOUBLE_PRECISION * 2);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == false);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1);
        Polynomial p2 = new Polynomial(1 + Constants.DOUBLE_PRECISION / 2);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1);
        Polynomial p2 = new Polynomial(1 + Constants.DOUBLE_PRECISION * 2);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == false);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1, 0);
        Polynomial p2 = new Polynomial(1);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1, 0, 0);
        Polynomial p2 = new Polynomial(1);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1, 0, 0);
        Polynomial p2 = new Polynomial(1, 0);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1, 0, 0);
        Polynomial p2 = new Polynomial(1, 0, 0);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(2, 3);
        Polynomial p2 = new Polynomial(2, 3);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(2, 3, 0);
        Polynomial p2 = new Polynomial(2, 3);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(2, 3, 0, 0);
        Polynomial p2 = new Polynomial(2, 3);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(2, 3, 0, 0);
        Polynomial p2 = new Polynomial(2, 3);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(2, 3, 0, Constants.DOUBLE_PRECISION / 2);
        Polynomial p2 = new Polynomial(2, 3);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == true);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(2, 3, 0, Constants.DOUBLE_PRECISION * 2);
        Polynomial p2 = new Polynomial(2, 3);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == false);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1 / Constants.DOUBLE_PRECISION, 3);
        Polynomial p2 = new Polynomial(1 / Constants.DOUBLE_PRECISION - 0.5, 3);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == false);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1 / Constants.DOUBLE_PRECISION, 3);
        Polynomial p2 = new Polynomial(1 / Constants.DOUBLE_PRECISION - 2, 3);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == false);
        Assert.IsTrue(p1.IsAlmostEquals(p2) == p2.IsAlmostEquals(p1));
    }

    {
        Polynomial p1 = new Polynomial(1 / 1e-6, 3);
        Polynomial p2 = new Polynomial(1 / 1e-6 - 0.5 * 1e-6, 3);
        Assert.IsTrue(p1.IsAlmostEquals(p2, 1e-6) == true);
        Assert.IsTrue(p1.IsAlmostEquals(p2, 1e-6) == p2.IsAlmostEquals(p1, 1e-6));
    }

    {
        Polynomial p1 = new Polynomial(1 / 1e-6, 3);
        Polynomial p2 = new Polynomial(1 / 1e-6 - 2 * 1e-6, 3);
        Assert.IsTrue(p1.IsAlmostEquals(p2, 1e-6) == false);
        Assert.IsTrue(p1.IsAlmostEquals(p2, 1e-6) == p2.IsAlmostEquals(p1, 1e-6));
    }

    // Evaluate.

    {
        AreAlmostRelativeEquals(2.0, new Polynomial(2).Evaluate(2.1));
        AreAlmostRelativeEquals(7.2, new Polynomial(3, 2).Evaluate(2.1));
        AreAlmostRelativeEquals(19.12, new Polynomial(4, 3, 2).Evaluate(2.1));
        AreAlmostRelativeEquals(28.352000000000004, new Polynomial(5, -4, 3, 2).Evaluate(2.1));
        AreAlmostRelativeEquals(100.81920000000001, new Polynomial(6, 5, 4, 3, 2).Evaluate(2.1));
        AreAlmostRelativeEquals(-6.0, new Polynomial(-6, 5, 4, -3, 2).Evaluate(0));
        AreAlmostRelativeEquals(24.253200000000003, new Polynomial(6, 5, 4, 3, 2).Evaluate(-2.1));
        AreAlmostRelativeEquals(17.807100000000005, new Polynomial(6 / 2, 5 / 2, 4 / 2, 3 / 2, 1).Evaluate(-2.1));
        AreAlmostRelativeEquals(45.152, new Polynomial(5, 4, 3, 2, 0).Evaluate(2.1));
        AreAlmostRelativeEquals(19.12, new Polynomial(4, 3, 2, 0, 0).Evaluate(2.1));

        Polynomial p = new Polynomial(4, 3, 2, 0, 0);
        double x = p.Evaluate(4.4);
        Assert.AreEqual(new Polynomial(4, 3, 2, 0, 0), p);
    }

    // Diffrentatie.

    {
        AreAlmostRelativeEquals(0.0, new Polynomial().Differentiate().Evaluate(2.1));
        AreAlmostRelativeEquals(0, new Polynomial(2).Differentiate().Evaluate(2.1));
        AreAlmostRelativeEquals(2.0, new Polynomial(-3, 2).Differentiate().Evaluate(2.1));
        AreAlmostRelativeEquals(11.4, new Polynomial(4, 3, 2).Differentiate().Evaluate(2.1));
        AreAlmostRelativeEquals(17.860000000000003, new Polynomial(5, 4, -3, 2).Differentiate().Evaluate(2.1));
        AreAlmostRelativeEquals(135.57800000000003, new Polynomial(6, 5, 4, 3, 2).Differentiate().Evaluate(2.1));
        AreAlmostRelativeEquals(5, new Polynomial(6, 5, 4, -3, 2).Differentiate().Evaluate(0));
        AreAlmostRelativeEquals(-46.198000000000008, new Polynomial(6, 5, 4, 3, 2).Differentiate().Evaluate(-2.1));
        AreAlmostRelativeEquals(-30.214000000000006, new Polynomial(6 / 2, 5 / 2, 4 / 2, 3 / 2, 1).Differentiate().Evaluate(-2.1));
        AreAlmostRelativeEquals(43.06, new Polynomial(5, 4, 3, 2, 0).Differentiate().Evaluate(2.1));
        AreAlmostRelativeEquals(11.4, new Polynomial(4, 3, 2, 0, 0).Differentiate().Evaluate(2.1));

        Polynomial p = new Polynomial(6, 5, 4, -3, 2);
        Polynomial d = p.Differentiate();
        Assert.AreEqual(new Polynomial(6, 5, 4, -3, 2), p);
    }

    // Normalize.

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2).Normalize();
        Polynomial p2 = new Polynomial(1.6 / 1.2, 1.5 / 1.2, 1.4 / 1.2, 1.3 / 1.2, 1);

        Assert.AreEqual(p1, p2);
    }

    {
        Polynomial p = new Polynomial(1.6, 1.5, 1.4, 1.3, -1.2);
        Polynomial p1 = p.Normalize();
        Polynomial p2 = new Polynomial(1.6 / -1.2, 1.5 / -1.2, 1.4 / -1.2, 1.3 / -1.2, 1);

        Assert.AreEqual(p1, p2);
        Assert.AreEqual(p, new Polynomial(1.6, 1.5, 1.4, 1.3, -1.2));
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, 1.4, 1.3, 0).Normalize();
        Polynomial p2 = new Polynomial(1.6, 1.5, 1.4, 1.3, 0);

        Assert.AreNotEqual(p1, p2);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, -1.4, 1.3, 0).Normalize();
        Polynomial p2 = new Polynomial(1.6 / 1.3, 1.5 / 1.3, -1.4 / 1.3, 1);

        Assert.AreEqual(p1, p2);
    }

    // operator / , DivRem(Rem)

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2);
        Polynomial p2 = new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2);

        Polynomial p3 = p1 / p2;
        Polynomial rem;
        Polynomial p4 = Polynomial.DivRem(p1, p2, out rem);

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(new Polynomial(1.0), p3);

        Assert.AreEqual(new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2), p1);
        Assert.AreEqual(new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2), p2);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2);
        Polynomial p2 = new Polynomial(1.6, 1.5, 1.4, 1.3);

        Polynomial p3 = p1 / p2;
        Polynomial rem;
        Polynomial p4 = Polynomial.DivRem(p1, p2, out rem);

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(new Polynomial(0.0059171597633137594, 0.923076923076923), p3);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2);
        Polynomial p2 = new Polynomial(1.6, 1.5, 1.4);

        Polynomial p3 = p1 / p2;
        Polynomial rem;
        Polynomial p4 = Polynomial.DivRem(p1, p2, out rem);

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(new Polynomial(0.00947521865889193, 0.010204081632653, 0.857142857142857), p3);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, 1.4, -1.3, 1.2);
        Polynomial p2 = new Polynomial(1.6, 1.5);

        Polynomial p3 = p1 / p2;
        Polynomial rem;
        Polynomial p4 = Polynomial.DivRem(p1, p2, out rem);

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(new Polynomial(-1.95253333333333, 2.768, -1.72, 0.8), p3);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, -1.4, 1.3, 1.2);
        Polynomial p2 = new Polynomial(1.6);

        Polynomial p3 = p1 / p2;
        Polynomial rem;
        Polynomial p4 = Polynomial.DivRem(p1, p2, out rem);

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(new Polynomial(1, 0.9375, -0.875, 0.8125, 0.75), p3);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, 1.4, 1.3);
        Polynomial p2 = new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2);

        Polynomial p3 = p1 / p2;
        Polynomial rem;
        Polynomial p4 = Polynomial.DivRem(p1, p2, out rem);

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(Polynomial.ZERO, p3);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, 1.4, 1.3, 0);
        Polynomial p2 = new Polynomial(1.6, 1.5, 1.4, 1.3);

        Polynomial p3 = p1 / p2;
        Polynomial rem;
        Polynomial p4 = Polynomial.DivRem(p1, p2, out rem);

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(new Polynomial(1.0), p3);
    }

    // operator %, DivRem(div)

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2);
        Polynomial p2 = new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2);

        Polynomial p3 = p1 % p2;
        Polynomial rem;
        Polynomial.DivRem(p1, p2, out rem);

        Assert.AreEqual(p3, rem);
        AreAlmostRelativeEquals(Polynomial.ZERO, p3);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2);
        Polynomial p2 = new Polynomial(1.6, 1.5, 1.4, 1.3);

        Polynomial p3 = p1 % p2;
        Polynomial rem;
        Polynomial.DivRem(p1, p2, out rem);

        Assert.AreEqual(p3, rem);
        AreAlmostRelativeEquals(new Polynomial(1.5905325443787, 0.0142011834319524, 0.00710059171597607), p3);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, -1.4, 1.3, 1.2);
        Polynomial p2 = new Polynomial(1.6, 1.5, 1.4);

        Polynomial p3 = p1 % p2;
        Polynomial rem;
        Polynomial.DivRem(p1, p2, out rem);

        Assert.AreEqual(p3, rem);
        AreAlmostRelativeEquals(new Polynomial(4.78483965014577, 4.46946064139942), p3);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, -1.4, 1.3, 1.2);
        Polynomial p2 = new Polynomial(1.6, 1.5);

        Polynomial p3 = p1 % p2;
        Polynomial rem;
        Polynomial.DivRem(p1, p2, out rem);

        Assert.AreEqual(p3, rem);
        AreAlmostRelativeEquals(new Polynomial(-1.61716148148148), p3);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2);
        Polynomial p2 = new Polynomial(1.6);

        Polynomial p3 = p1 % p2;
        Polynomial rem;
        Polynomial.DivRem(p1, p2, out rem);

        Assert.AreEqual(p3, rem);
        AreAlmostRelativeEquals(Polynomial.ZERO, p3);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, 1.4, 1.3);
        Polynomial p2 = new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2);

        Polynomial p3 = p1 % p2;
        Polynomial rem;
        Polynomial.DivRem(p1, p2, out rem);

        Assert.AreEqual(p3, rem);
        AreAlmostRelativeEquals(new Polynomial(1.6, 1.5, 1.4, 1.3), p3);

        Assert.AreEqual(new Polynomial(1.6, 1.5, 1.4, 1.3), p1);
        Assert.AreEqual(new Polynomial(1.6, 1.5, 1.4, 1.3, 1.2), p2);
    }

    {
        Polynomial p1 = new Polynomial(1.6, 1.5, 1.4, 1.3, 0);
        Polynomial p2 = new Polynomial(1.6, 1.5, 1.4, 1.3);

        Polynomial p3 = p1 % p2;
        Polynomial rem;
        Polynomial.DivRem(p1, p2, out rem);

        Assert.AreEqual(p3, rem);
        AreAlmostRelativeEquals(Polynomial.ZERO, p3);
    }

    // operator * polynomials
    {
        Polynomial p1 = Polynomial.ZERO;
        Polynomial p2 = Polynomial.ZERO;

        Polynomial p3 = p1 * p2;
        Polynomial p4 = p2 * p1;

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(Polynomial.ZERO, p3);
    }

    {
        Polynomial p1 = Polynomial.ZERO;
        Polynomial p2 = new Polynomial(1.6, 1.5, 1.4, 1.3);

        Polynomial p3 = p1 * p2;
        Polynomial p4 = p2 * p1;

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(Polynomial.ZERO, p3);
    }

    {
        Polynomial p1 = new Polynomial(1.5);
        Polynomial p2 = new Polynomial(1.3);

        Polynomial p3 = p1 * p2;
        Polynomial p4 = p2 * p1;

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(new Polynomial(1.95), p3);
    }

    {
        Polynomial p1 = new Polynomial(1.5);
        Polynomial p2 = new Polynomial(1.4, 1.3);

        Polynomial p3 = p1 * p2;
        Polynomial p4 = p2 * p1;

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(new Polynomial(2.1, 1.95), p3);
    }

    {
        Polynomial p1 = new Polynomial(1.4, 1.5);
        Polynomial p2 = new Polynomial(1.6, 1.3);

        Polynomial p3 = p1 * p2;
        Polynomial p4 = p2 * p1;

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(new Polynomial(2.24, 4.22, 1.95), p3);
    }

    {
        Polynomial p1 = new Polynomial(-0.9, 1.3, 1.4, 1.5);
        Polynomial p2 = new Polynomial(1.6, 1.3);

        Polynomial p3 = p1 * p2;
        Polynomial p4 = p2 * p1;

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(new Polynomial(-1.44, 0.91, 3.93, 4.22, 1.95), p3);

        Assert.AreEqual(p1, new Polynomial(-0.9, 1.3, 1.4, 1.5));
        Assert.AreEqual(p2, new Polynomial(1.6, 1.3));
    }

    {
        Polynomial p1 = new Polynomial(0.9, 1.3, 1.4, 1.5);
        Polynomial p2 = new Polynomial(1.6, 1.3);

        Polynomial p3 = p1 * p2;
        Polynomial p4 = p2 * p1;

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(new Polynomial(1.44, 3.25, 3.93, 4.22, 1.95), p3);
    }

    {
        Polynomial p1 = new Polynomial(-1.3, 1.4, 1.5);
        Polynomial p2 = new Polynomial(1.2, -1.6, 1.3);

        Polynomial p3 = p1 * p2;
        Polynomial p4 = p2 * p1;

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(new Polynomial(-1.56, 3.76, -2.13, -0.580000000000001, 1.95), p3);
    }

    // operator +  
    {
        Polynomial p1 = Polynomial.ZERO;
        Polynomial p2 = Polynomial.ZERO;

        Polynomial p3 = p1 + p2;
        Polynomial p4 = p2 + p1;

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(Polynomial.ZERO, p3);
    }

    {
        Polynomial p1 = Polynomial.ZERO;
        Polynomial p2 = new Polynomial(1.6, 1.5, 1.4, 1.3);

        Polynomial p3 = p1 + p2;
        Polynomial p4 = p2 + p1;

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(new Polynomial(1.6, 1.5, 1.4, 1.3), p3);
    }

    {
        Polynomial p1 = new Polynomial(1.5);
        Polynomial p2 = new Polynomial(-1.3);

        Polynomial p3 = p1 + p2;
        Polynomial p4 = p2 + p1;

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(new Polynomial(0.19999999999999996), p3);
    }

    {
        Polynomial p1 = new Polynomial(1.7, -1.6, 1.5);
        Polynomial p2 = new Polynomial(0.3, -0.2, 0.1, 1.4, 1.3);

        Polynomial p3 = p1 + p2;
        Polynomial p4 = p2 + p1;

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(new Polynomial(2, -1.8, 1.6, 1.4, 1.3), p3);

        Assert.AreEqual(new Polynomial(1.7, -1.6, 1.5), p1);
        Assert.AreEqual(new Polynomial(0.3, -0.2, 0.1, 1.4, 1.3), p2);
    }

    {
        Polynomial p1 = new Polynomial(0.3, 0.2, -0.1, 1.4, 1.3);
        Polynomial p2 = new Polynomial(-0.5, 0.4, 0.3, 2.4, 2.3);

        Polynomial p3 = p1 + p2;
        Polynomial p4 = p2 + p1;

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(new Polynomial(-0.2, 0.6, 0.2, 3.8, 3.6), p3);
    }

    {
        Polynomial p1 = new Polynomial(0.3, 0.2, -0.1, 1.4, 1.3);
        Polynomial p2 = -new Polynomial(0.3, 0.2, -0.1, 1.4, 1.3);

        Polynomial p3 = p1 + p2;
        Polynomial p4 = p2 + p1;

        Assert.AreEqual(p3, p4);
        AreAlmostRelativeEquals(Polynomial.ZERO, p3);
    }

    // operator -
    {
        Polynomial p1 = Polynomial.ZERO;
        Polynomial p2 = Polynomial.ZERO;

        Polynomial p3 = p1 - p2;
        Polynomial p4 = p2 - p1;

        AreAlmostRelativeEquals(Polynomial.ZERO, p3);
        AreAlmostRelativeEquals(Polynomial.ZERO, p4);
    }

    {
        Polynomial p1 = Polynomial.ZERO;
        Polynomial p2 = new Polynomial(1.6, 1.5, 1.4, 1.3);

        Polynomial p3 = p1 - p2;
        Polynomial p4 = p2 - p1;

        AreAlmostRelativeEquals(new Polynomial(-1.6, -1.5, -1.4, -1.3), p3);
        AreAlmostRelativeEquals(-new Polynomial(-1.6, -1.5, -1.4, -1.3), p4);
    }

    {
        Polynomial p1 = new Polynomial(1.5);
        Polynomial p2 = new Polynomial(-1.3);

        Polynomial p3 = p1 - p2;
        Polynomial p4 = p2 - p1;

        AreAlmostRelativeEquals(new Polynomial(2.8), p3);
        AreAlmostRelativeEquals(-new Polynomial(2.8), p4);
    }

    {
        Polynomial p1 = new Polynomial(1.7, -1.6, 1.5);
        Polynomial p2 = new Polynomial(0.3, -0.2, 0.1, 1.4, 1.3);

        Polynomial p3 = p1 - p2;
        Polynomial p4 = p2 - p1;

        AreAlmostRelativeEquals(new Polynomial(1.4, -1.4, 1.4, -1.4, -1.3), p3);
        AreAlmostRelativeEquals(-new Polynomial(1.4, -1.4, 1.4, -1.4, -1.3), p4);

        Assert.AreEqual(new Polynomial(1.7, -1.6, 1.5), p1);
        Assert.AreEqual(new Polynomial(0.3, -0.2, 0.1, 1.4, 1.3), p2);
    }

    {
        Polynomial p1 = new Polynomial(0.3, 0.2, -0.1, 1.4, 1.3);
        Polynomial p2 = new Polynomial(-0.5, 0.4, 0.3, 2.4, 2.3);

        Polynomial p3 = p1 - p2;
        Polynomial p4 = p2 - p1;

        AreAlmostRelativeEquals(new Polynomial(0.8, -0.2, -0.4, -1, -1), p3);
        AreAlmostRelativeEquals(-new Polynomial(0.8, -0.2, -0.4, -1, -1), p4);
    }

    {
        Polynomial p1 = new Polynomial(1.3, 1.4, -0.1, 0.2, 0.3);
        Polynomial p2 = new Polynomial(1.3, 1.4, -0.1, 0.2, 0.3);

        Polynomial p3 = p1 - p2;
        Polynomial p4 = p2 - p1;

        AreAlmostRelativeEquals(Polynomial.ZERO, p3);
        AreAlmostRelativeEquals(Polynomial.ZERO, p4);
    }

    // operator - unary

    {
        Polynomial p = Polynomial.ZERO;
        AreAlmostRelativeEquals(p, - -p);
        AreAlmostRelativeEquals(-p, Polynomial.ZERO - p);
    }

    {
        Polynomial p = new Polynomial(1.0);
        AreAlmostRelativeEquals(p, - -p);
        AreAlmostRelativeEquals(-p, Polynomial.ZERO - p);
    }

    {
        Polynomial p = new Polynomial(6.7, -4.5, 3.4, -2.0);
        AreAlmostRelativeEquals(p, - -p);
        AreAlmostRelativeEquals(-p, Polynomial.ZERO - p);

        Polynomial pp = -p;
        Assert.AreEqual(new Polynomial(6.7, -4.5, 3.4, -2.0), p);
    }

    // operator * const

    {
        Polynomial p1 = Polynomial.ZERO * 6;
        Polynomial p2 = 6 * Polynomial.ZERO;

        AreAlmostRelativeEquals(Polynomial.ZERO, p1);
        AreAlmostRelativeEquals(p1, p2);
    }

    {
        Polynomial p1 = new Polynomial(1.0) * 5;
        Polynomial p2 = 5 * new Polynomial(1.0);

        AreAlmostRelativeEquals(new Polynomial(5.0), p1);
        AreAlmostRelativeEquals(p1, p2);
    }

    {
        Polynomial p = new Polynomial(6.7, -4.5, 3.4, -2.0);
        Polynomial p1 = p * -2;
        Polynomial p2 = -2 * new Polynomial(6.7, -4.5, 3.4, -2.0);

        AreAlmostRelativeEquals(new Polynomial(-13.4, 9, -6.8, 4.0), p1);
        AreAlmostRelativeEquals(p1, p2);
        Assert.AreEqual(new Polynomial(6.7, -4.5, 3.4, -2.0), p);
    }

    {
        Polynomial p = new Polynomial(6.7, -4.5, 3.4, -2.0) * 0;

        AreAlmostRelativeEquals(Polynomial.ZERO, p);
    }

    // operator / const

    {
        Polynomial p = Polynomial.ZERO / 6;

        AreAlmostRelativeEquals(Polynomial.ZERO, p);
    }

    {
        Polynomial p = new Polynomial(5.0) / 5;

        AreAlmostRelativeEquals(new Polynomial(1.0), p);
    }

    {
        Polynomial pp = new Polynomial(-13.4, 9, -6.8, 4.0);
        Polynomial p = pp / -2;

        AreAlmostRelativeEquals(new Polynomial(6.7, -4.5, 3.4, -2.0), p);
        Assert.AreEqual(new Polynomial(-13.4, 9, -6.8, 4.0), pp);
    }

    {
        Polynomial p = new Polynomial(-13.4, 9, -6.8, 4.0) / Double.PositiveInfinity;

        AreAlmostRelativeEquals(Polynomial.ZERO, p);
    }

    // Pow 

    {
        Polynomial p1 = p_1;
        Polynomial p2 = p1 ^ 0;

        AreAlmostRelativeEquals(Polynomial.ONE, p2);
    }

    {
        Polynomial p1 = p_1;
        Polynomial p2 = p1 ^ 1;

        AreAlmostRelativeEquals(p_1, p2);
    }

    {
        Polynomial p1 = p_1 * p_1;
        Polynomial p2 = p1 ^ 1;

        AreAlmostRelativeEquals(p_1 * p_1, p2);
    }

    {
        Polynomial p1 = p_1 * p_1;
        Polynomial p2 = p1 ^ 2;

        AreAlmostRelativeEquals(p_1 * p_1 * p_1 * p_1, p2);
    }

    {
        Polynomial p1 = p_1;
        Polynomial p2 = p1 ^ 4;

        AreAlmostRelativeEquals(p_1 * p_1 * p_1 * p_1, p2);
    }

    // GCD

    {
        Polynomial A = p_1;
        Polynomial B = p_2;
        Polynomial GCD = Polynomial.GCD(A, B);
        Polynomial expected = Polynomial.ONE;

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = p_2;
        Polynomial B = p_1;
        Polynomial GCD = Polynomial.GCD(A, B);
        Polynomial expected = Polynomial.ONE;

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = p_1;
        Polynomial B = p_1;
        Polynomial GCD = Polynomial.GCD(A, B);
        Polynomial expected = p_1.Normalize();

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = p_1;
        Polynomial B = Polynomial.ZERO;
        Polynomial GCD = Polynomial.GCD(A, B);
        Polynomial expected = p_1.Normalize();

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = p_1;
        Polynomial B = Polynomial.ONE;
        Polynomial GCD = Polynomial.GCD(A, B);
        Polynomial expected = Polynomial.ONE;

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = Polynomial.ZERO;
        Polynomial B = p_2;
        Polynomial GCD = Polynomial.GCD(A, B);
        Polynomial expected = p_2.Normalize();

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = Polynomial.ONE;
        Polynomial B = p_2;
        Polynomial GCD = Polynomial.GCD(A, B);
        Polynomial expected = Polynomial.ONE;

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = p_1 * p_1;
        Polynomial B = p_1;
        Polynomial GCD = Polynomial.GCD(A, B);
        Polynomial expected = p_1.Normalize();

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = 2 * p_1 * p_1;
        Polynomial B = p_1;
        Polynomial GCD = Polynomial.GCD(A, B);
        Polynomial expected = p_1.Normalize();

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = 2 * p_1 * p_1;
        Polynomial B = 2 * p_1;
        Polynomial GCD = Polynomial.GCD(A, B);
        Polynomial expected = p_1.Normalize();

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = p_1 ^ 3;
        Polynomial B = p_1 ^ 4;
        Polynomial GCD = Polynomial.GCD(A, B);
        Polynomial expected = (p_1 ^ 3).Normalize();

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = p_1 * p_2 * p_2 * p_2;
        Polynomial B = A.Differentiate();
        Polynomial GCD = Polynomial.GCD(A, B);
        Polynomial expected = (p_2 * p_2).Normalize();

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = p_1 * p_1 * p_2 * p_2;
        Polynomial B = A.Differentiate();
        Polynomial GCD = Polynomial.GCD(A, B);
        Polynomial expected = (p_1 * p_2).Normalize();

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = p_1 * p_2 * (p_4 ^ 3);
        Polynomial B = p_1 * (p_2 ^ 2) * (p_3 ^ 3) * (p_4 ^ 4);
        Polynomial GCD = Polynomial.GCD(A, B, 1e-9);
        Polynomial expected = (p_1 * p_2 * (p_4 ^ 3)).Normalize();

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = p_1 * p_2 * (p_4 ^ 3);
        Polynomial B = A.Differentiate();
        Polynomial GCD = Polynomial.GCD(A, B);
        Polynomial expected = (p_4 ^ 2).Normalize();

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = p_1 * p_2 * (p_4 ^ 2) * (p_5 ^ 3) * (p_6 ^ 2);
        Polynomial B = p_1 * (p_2 ^ 2) * (p_4 ^ 3) * (p_5 ^ 2) * (p_6 ^ 2);
        Polynomial GCD = Polynomial.GCD(A, B, 1e-7);
        Polynomial expected = (p_1 * p_2 * (p_4 ^ 2) * (p_5 ^ 2) * (p_6 ^ 2)).Normalize();

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = p_1 * p_2 * (p_4 ^ 2) * (p_5 ^ 2);
        Polynomial B = A.Differentiate();
        Polynomial GCD = Polynomial.GCD(A, B, 1e-11);
        Polynomial expected = (p_4 * p_5).Normalize();

        AreAlmostRelativeEquals(expected, GCD);
    }

    {
        Polynomial A = p_1 * p_2 * (p_4 ^ 2) * (p_5 ^ 3);
        Polynomial B = A.Differentiate();
        Polynomial GCD = Polynomial.GCD(A, B, 1e-9);
        Polynomial expected = (p_4 * (p_5 ^ 2)).Normalize();

        AreAlmostRelativeEquals(expected, GCD, 1e-11);
    }

    // RemoveZeroes.

    {
        Polynomial p1 = new Polynomial(Constants.DOUBLE_PRECISION / 10);
        Polynomial p2 = p1.RemoveZeroes();

        AreAlmostRelativeEquals(p1, p2);
    }

    {
        Polynomial p1 = new Polynomial(Constants.DOUBLE_PRECISION / 10, 1);
        Polynomial p2 = p1.RemoveZeroes();

        AreAlmostRelativeEquals(new Polynomial(1), p2);
    }

    {
        Polynomial p1 = new Polynomial(Constants.DOUBLE_PRECISION / 10, Constants.DOUBLE_PRECISION / 10, 2, 4);
        Polynomial p2 = p1.RemoveZeroes();

        AreAlmostRelativeEquals(new Polynomial(2, 4), p2);
    }

    {
        Polynomial p1 = new Polynomial(Constants.DOUBLE_PRECISION / 10, Constants.DOUBLE_PRECISION * 10, 2, 4);
        Polynomial p2 = p1.RemoveZeroes();

        AreAlmostRelativeEquals(new Polynomial(Constants.DOUBLE_PRECISION * 10, 2, 4), p2);
    }
}