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.
Pokazywanie postów oznaczonych etykietą Math. Pokaż wszystkie posty
Pokazywanie postów oznaczonych etykietą Math. Pokaż wszystkie posty
2012-01-27
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:
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.
$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:
$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:
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;
}
2012-01-22
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#:
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.
I testy:
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);
}
}
2012-01-20
Klasa liczby zespolonej
Moja implementacja klasy liczby zespolonej, nie jest może wszystko-mająca, ale do FFT wystarczy. Do tego testy. Publikuje to dlatego, że przechodzę na implementację klasy Complex z System.Numerics. Jedyną różnicę jaką zauważyłem to zachowanie się w przypadku liczenia przesunięcia fazowego z (0,0). Zgodnie z definicją z wikipedii i nazwijmy to zdrowym rozsądkiem wartość ta powinna być nieokreślona, tymczasem w System.Numerics jest to zero.
Testy:
public struct Complex
{
public static readonly Complex ZERO = new Complex(0, 0);
public static readonly Complex ONE = new Complex(1, 0);
public static readonly Complex i = new Complex(0, 1);
public static readonly Complex UNDEFINED = new Complex(Double.NaN, Double.NaN);
public readonly double Real;
public readonly double Image;
public Complex(double a_real, double a_image = 0)
{
Real = a_real;
Image = a_image;
}
public static bool operator ==(Complex a_a, Complex a_b)
{
return (a_a.Real == a_b.Real) & (a_a.Image == a_b.Image);
}
public static bool operator !=(Complex a_a, Complex a_b)
{
return (a_a.Real != a_b.Real) | (a_a.Image != a_b.Image);
}
public static Complex operator +(Complex a_a)
{
return a_a;
}
public static Complex FromPolar(double a_magnitude, double a_phase)
{
return new Complex(a_magnitude * Math.Cos(a_phase),
a_magnitude * Math.Sin(a_phase));
}
public static Complex FromPolar(double a_phase)
{
return new Complex(Math.Cos(a_phase), Math.Sin(a_phase));
}
public static Complex operator -(Complex a_a)
{
return new Complex(-a_a.Real, -a_a.Image);
}
public static Complex operator +(Complex a_a, Complex a_b)
{
return new Complex(
a_a.Real + a_b.Real,
a_a.Image + a_b.Image);
}
public static Complex operator -(Complex a_a, Complex a_b)
{
return new Complex(
a_a.Real - a_b.Real,
a_a.Image - a_b.Image);
}
public static Complex operator *(Complex a_a, Complex a_b)
{
return new Complex(
a_a.Real * a_b.Real - a_a.Image * a_b.Image,
a_a.Real * a_b.Image + a_a.Image * a_b.Real);
}
public static Complex operator *(Complex a_a, double a_b)
{
return new Complex(
a_a.Real * a_b,
a_a.Image * a_b);
}
public static Complex operator *(double a_a, Complex a_b)
{
return new Complex(
a_b.Real * a_a,
a_b.Image * a_a);
}
public static Complex operator /(Complex a_a, Complex a_b)
{
double den = a_b.Real * a_b.Real + a_b.Image * a_b.Image;
return new Complex(
(a_a.Real * a_b.Real + a_a.Image * a_b.Image) / den,
(a_a.Image * a_b.Real - a_a.Real * a_b.Image) / den);
}
public Complex Scale(double a_scale)
{
return new Complex(
Real * a_scale,
Image * a_scale);
}
public bool IsAlmostEquals(Complex a_a,
double a_precision = Constants.DOUBLE_PRECISION)
{
return Real.IsAlmostEquals(a_a.Real, a_precision) &&
Image.IsAlmostEquals(a_a.Image, a_precision);
}
public bool IsAlmostRelativeEquals(Complex a_a,
double a_precision = Constants.DOUBLE_PRECISION)
{
return Real.IsAlmostRelativeEquals(a_a.Real, a_precision) &&
Image.IsAlmostRelativeEquals(a_a.Image, a_precision);
}
public override int GetHashCode()
{
return Real.GetHashCode() ^ Image.GetHashCode();
}
public override bool Equals(object a_obj)
{
if (a_obj == null)
return false;
if (!typeof(Complex).Equals(a_obj.GetType()))
return false;
Complex complex = (Complex)a_obj;
return (complex.Real == Real) && (complex.Image == Image);
}
public bool Equals(Complex a_complex)
{
return (a_complex.Real == Real) && (a_complex.Image == Image);
}
[DebuggerHidden]
public Complex Conjugate
{
get
{
return new Complex(Real, -Image);
}
}
public double Abs
{
get
{
return MathExtensions.Hypot(Real, Image);
}
}
public double Phase
{
get
{
if (IsAlmostRelativeEquals(ZERO))
return Double.NaN;
return Math.Atan2(Image, Real);
}
}
public override string ToString()
{
if ((Real != 0) && (Image != 0))
return String.Format("({0} + {1}i)", Real, Image);
else if (Real != 0)
return Real.ToString();
else if (Image != 0)
return String.Format("{1}i", Image);
else
return "0";
}
}Testy:
[TestMethod]
public void Complex_Test()
{
// Constructor
{
Assert.AreEqual(4, new Complex(4).Real);
Assert.AreEqual(0, new Complex(4).Image);
Assert.AreEqual(5, new Complex(5, 6).Real);
Assert.AreEqual(6, new Complex(5, 6).Image);
}
// ==, !=
{
Assert.IsTrue(new Complex(5, 4) == new Complex(5, 4));
Assert.IsTrue(new Complex(5.3, 4.7) == new Complex(5.3, 4.7));
Assert.IsFalse(new Complex(5, 4) == new Complex(5, 5));
Assert.IsFalse(new Complex(5, 4) == new Complex(4, 4));
Assert.IsFalse(new Complex(5, 5) == new Complex(5, 4));
Assert.IsFalse(new Complex(4, 4) == new Complex(5, 4));
Assert.IsFalse(new Complex(5, 4) != new Complex(5, 4));
Assert.IsTrue(new Complex(5, 4) != new Complex(5, 5));
Assert.IsTrue(new Complex(5, 4) != new Complex(4, 4));
Assert.IsTrue(new Complex(5, 5) != new Complex(5, 4));
Assert.IsTrue(new Complex(4, 4) != new Complex(5, 4));
}
// unary -
{
Assert.IsTrue(new Complex(-1, 1) == -new Complex(1, -1));
Assert.IsTrue(new Complex(-1.6, 1.5) == -new Complex(1.6, -1.5));
Assert.IsTrue(new Complex(-1, -1) == -new Complex(1, 1));
Assert.IsTrue(new Complex(1, -1) == -new Complex(-1, 1));
Assert.IsTrue(new Complex(1, 0) == -new Complex(-1, 0));
Assert.IsTrue(new Complex(0, 0) == -new Complex(0, 0));
Assert.IsTrue(new Complex(0, 1) == -new Complex(0, -1));
Assert.IsTrue(-new Complex(-1, 1) == - -new Complex(1, -1));
Assert.IsTrue(-new Complex(-1, -1) == - -new Complex(1, 1));
Assert.IsTrue(-new Complex(1, -1) == - -new Complex(-1, 1));
Assert.IsTrue(-new Complex(1, 0) == - -new Complex(-1, 0));
Assert.IsTrue(-new Complex(0, 0) == - -new Complex(0, 0));
Assert.IsTrue(-new Complex(0, 1) == - -new Complex(0, -1));
}
// unary +
{
Assert.IsTrue(-new Complex(-1, 1) == +new Complex(1, -1));
Assert.IsTrue(-new Complex(-1, -1) == +new Complex(1, 1));
Assert.IsTrue(-new Complex(-1.2, -1.5) == +new Complex(1.2, 1.5));
Assert.IsTrue(-new Complex(1, -1) == +new Complex(-1, 1));
Assert.IsTrue(-new Complex(1, 0) == +new Complex(-1, 0));
Assert.IsTrue(-new Complex(0, 0) == +new Complex(0, 0));
Assert.IsTrue(-new Complex(0, 1) == +new Complex(0, -1));
}
// +
{
Assert.IsTrue(new Complex(0, 0) == new Complex(0, 0) + new Complex(0, 0));
Assert.IsTrue(new Complex(0, 1) == new Complex(0, 1) + new Complex(0, 0));
Assert.IsTrue(new Complex(0, 1) == new Complex(0, 0) + new Complex(0, 1));
Assert.IsTrue(new Complex(1, 0) == new Complex(1, 0) + new Complex(0, 0));
Assert.IsTrue(new Complex(1, 0) == new Complex(0, 0) + new Complex(1, 0));
Assert.IsTrue(new Complex(1, 1) == new Complex(1, 0) + new Complex(0, 1));
Assert.IsTrue(new Complex(1, 1) == new Complex(0, 1) + new Complex(1, 0));
Assert.IsTrue(new Complex(3, 3) == new Complex(1, 2) + new Complex(2, 1));
Assert.IsTrue(new Complex(4.1, 3.7) == new Complex(1.4, 2.6) +
new Complex(2.7, 1.1));
}
// -
{
Assert.IsTrue(new Complex(0, 0) == new Complex(0, 0) - -new Complex(0, 0));
Assert.IsTrue(new Complex(0, 1) == new Complex(0, 1) - -new Complex(0, 0));
Assert.IsTrue(new Complex(0, 1) == new Complex(0, 0) - -new Complex(0, 1));
Assert.IsTrue(new Complex(1, 0) == new Complex(1, 0) - -new Complex(0, 0));
Assert.IsTrue(new Complex(1, 0) == new Complex(0, 0) - -new Complex(1, 0));
Assert.IsTrue(new Complex(1, 1) == new Complex(1, 0) - -new Complex(0, 1));
Assert.IsTrue(new Complex(1, 1) == new Complex(0, 1) - -new Complex(1, 0));
Assert.IsTrue(new Complex(3, 3) == new Complex(1, 2) - -new Complex(2, 1));
Assert.IsTrue(new Complex(3.8, 3.8) == new Complex(1.7, 2.6) -
-new Complex(2.1, 1.2));
}
// * Complex
{
Assert.IsTrue(new Complex(0, 0) == new Complex(0, 0) * new Complex(0, 0));
Assert.IsTrue(new Complex(0, 0) == new Complex(0, 1) * new Complex(0, 0));
Assert.IsTrue(new Complex(0, 0) == new Complex(0, 0) * new Complex(0, 1));
Assert.IsTrue(new Complex(0, 0) == new Complex(1, 0) * new Complex(0, 0));
Assert.IsTrue(new Complex(0, 0) == new Complex(0, 0) * new Complex(1, 0));
Assert.IsTrue(new Complex(0, 1) == new Complex(1, 0) * new Complex(0, 1));
Assert.IsTrue(new Complex(0, 1) == new Complex(0, 1) * new Complex(1, 0));
Assert.IsTrue(new Complex(0, 5) == new Complex(1, 2) * new Complex(2, 1));
Assert.IsTrue(new Complex(-14, 31) == new Complex(2, 3) *
new Complex(5, 8));
AreAlmostRelativeEquals(new Complex(-14.32, 37.47),
new Complex(2.3, 3.2) * new Complex(5.6, 8.5));
}
// Complex * double, double * Complex, Scale
{
Action<Complex, Complex, double> scale_test = (c1, c2, s) =>
{
Assert.IsTrue(c1 == c2 * s);
Assert.IsTrue(c1 == s * c2);
Assert.IsTrue(c1 == c2.Scale(s));
};
scale_test(new Complex(0, 0), new Complex(0, 0), 0);
scale_test(new Complex(0, 0), new Complex(4, 4), 0);
scale_test(new Complex(0, 0), new Complex(4, 0), 0);
scale_test(new Complex(0, 0), new Complex(0, 4), 0);
scale_test(new Complex(0, 0), new Complex(0, 0), 4);
scale_test(new Complex(16, 16), new Complex(4, 4), 4);
scale_test(new Complex(16, 0), new Complex(4, 0), 4);
scale_test(new Complex(0, 16), new Complex(0, 4), 4);
scale_test(new Complex(0, 0), new Complex(0, 0), -4);
scale_test(new Complex(-16, -16), new Complex(4, 4), -4);
scale_test(new Complex(-16, 0), new Complex(4, 0), -4);
scale_test(new Complex(0, -16), new Complex(0, 4), -4);
scale_test(new Complex(-16, -12), new Complex(4, 3), -4);
scale_test(new Complex(-16, 12), new Complex(4, -3), -4);
scale_test(new Complex(16, -12), new Complex(-4, 3), -4);
scale_test(new Complex(17.63, -13.12), new Complex(-4.3, 3.2), -4.1);
}
// /
{
AreAlmostRelativeEquals(new Complex(23.0 / 41.0, 2.0 / 41.0),
new Complex(2, 3) / new Complex(4, 5));
AreAlmostRelativeEquals(new Complex(23.0 / 13.0, -2.0 / 13.0),
new Complex(4, 5) / new Complex(2, 3));
AreAlmostRelativeEquals(
new Complex(1.4896226415094338, -0.13867924528301892),
new Complex(4.4, 5.3) / new Complex(2.6, 3.8));
}
// IsAlmostEquals
{
Complex c = new Complex(0, 0);
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(2, 2) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(2, -2) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-2, 2) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-2, -2) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(2, 0.5) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-2, 0.5) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(2, -0.5) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-2, -0.5) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(0.5, 2) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(0.5, -2) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-0.5, 2) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-0.5, -2) *
Constants.DOUBLE_PRECISION));
Assert.IsTrue(c.IsAlmostEquals(c + new Complex(0.5, 0.5) *
Constants.DOUBLE_PRECISION));
Assert.IsTrue(c.IsAlmostEquals(c + new Complex(0.5, -0.5) *
Constants.DOUBLE_PRECISION));
Assert.IsTrue(c.IsAlmostEquals(c + new Complex(-0.5, 0.5) *
Constants.DOUBLE_PRECISION));
Assert.IsTrue(c.IsAlmostEquals(c + new Complex(-0.5, -0.5) *
Constants.DOUBLE_PRECISION));
c = new Complex(40, -40);
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(2, 2) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(2, -2) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-2, 2) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-2, -2) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(2, 0.5) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-2, 0.5) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(2, -0.5) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-2, -0.5) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(0.5, 2) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(0.5, -2) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-0.5, 2) *
Constants.DOUBLE_PRECISION));
Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-0.5, -2) *
Constants.DOUBLE_PRECISION));
Assert.IsTrue(c.IsAlmostEquals(c + new Complex(0.5, 0.5) *
Constants.DOUBLE_PRECISION));
Assert.IsTrue(c.IsAlmostEquals(c + new Complex(0.5, -0.5) *
Constants.DOUBLE_PRECISION));
Assert.IsTrue(c.IsAlmostEquals(c + new Complex(-0.5, 0.5) *
Constants.DOUBLE_PRECISION));
Assert.IsTrue(c.IsAlmostEquals(c + new Complex(-0.5, -0.5) *
Constants.DOUBLE_PRECISION));
}
// IsAlmostRelativeEquals
{
Complex c = new Complex(1 / Constants.DOUBLE_PRECISION,
1 / Constants.DOUBLE_PRECISION);
Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(2, 2)));
Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(2, -2)));
Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(-2, 2)));
Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(-2, -2)));
Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(2, 0.5)));
Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(-2, 0.5)));
Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(2, -0.5)));
Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(-2, -0.5)));
Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(0.5, 2)));
Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(0.5, -2)));
Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(-0.5, 2)));
Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(-0.5, -2)));
Assert.IsTrue(c.IsAlmostRelativeEquals(c + new Complex(0.5, 0.5)));
Assert.IsTrue(c.IsAlmostRelativeEquals(c + new Complex(0.5, -0.5)));
Assert.IsTrue(c.IsAlmostRelativeEquals(c + new Complex(-0.5, 0.5)));
Assert.IsTrue(c.IsAlmostRelativeEquals(c + new Complex(-0.5, -0.5)));
}
// Equals(obj), Equals(complex)
{
Assert.AreEqual(new Complex(5, 4), new Complex(5, 4));
Assert.AreEqual(new Complex(5.4, 4.4), new Complex(5.4, 4.4));
Assert.AreNotEqual(new Complex(5, 4), new Complex(5, 5));
Assert.AreNotEqual(new Complex(5, 4), new Complex(4, 4));
Assert.AreNotEqual(new Complex(5, 5), new Complex(5, 4));
Assert.AreNotEqual(new Complex(4, 4), new Complex(5, 4));
Assert.IsTrue(new Complex(5, 4).Equals(new Complex(5, 4)));
Assert.IsTrue(new Complex(5.3, 4.3).Equals(new Complex(5.3, 4.3)));
Assert.IsFalse(new Complex(5, 4).Equals(new Complex(5, 5)));
Assert.IsFalse(new Complex(5, 4).Equals(new Complex(4, 4)));
Assert.IsFalse(new Complex(5, 5).Equals(new Complex(5, 4)));
Assert.IsFalse(new Complex(4, 4).Equals(new Complex(5, 4)));
}
// Conjugate
{
Assert.IsTrue(new Complex(0, 0) == new Complex(0, 0).Conjugate);
Assert.IsTrue(new Complex(1, 0) == new Complex(1, 0).Conjugate);
Assert.IsTrue(new Complex(0, -1) == new Complex(0, 1).Conjugate);
Assert.IsTrue(new Complex(0, 1) == new Complex(0, -1).Conjugate);
Assert.IsTrue(new Complex(2, -3) == new Complex(2, 3).Conjugate);
Assert.IsTrue(new Complex(2, 3) == new Complex(2, -3).Conjugate);
Assert.IsTrue(new Complex(-2, -3) == new Complex(-2, 3).Conjugate);
Assert.IsTrue(new Complex(-2, 3) == new Complex(-2, -3).Conjugate);
Assert.IsTrue(new Complex(-2.2, 3.2) == new Complex(-2.2, -3.2).Conjugate);
}
// Abs
{
Assert.IsTrue(0 == new Complex(0, 0).Abs);
Assert.IsTrue(1 == new Complex(0, 1).Abs);
Assert.IsTrue(1 == new Complex(0, -1).Abs);
Assert.IsTrue(1 == new Complex(1, 0).Abs);
Assert.IsTrue(1 == new Complex(-1, 0).Abs);
Assert.IsTrue(5 == new Complex(3, 4).Abs);
Assert.IsTrue(5 == new Complex(-3, 4).Abs);
Assert.IsTrue(5 == new Complex(3, -4).Abs);
Assert.IsTrue(5 == new Complex(-3, -4).Abs);
Assert.IsTrue(2*MathExtensions.SQRT2 == new Complex(2, 2).Abs);
}
// Phase
{
Assert.AreEqual(Double.NaN,
MathExtensions.ToDeg(new Complex(0, 0).Phase));
Assert.AreEqual(Double.NaN,
MathExtensions.ToDeg(new Complex(Constants.DOUBLE_PRECISION / 2,
Constants.DOUBLE_PRECISION / 2).Phase));
Assert.AreEqual(0,
MathExtensions.ToDeg(new Complex(1, 0).Phase));
Assert.AreEqual(45,
MathExtensions.ToDeg(new Complex(1, 1).Phase));
Assert.AreEqual(90,
MathExtensions.ToDeg(new Complex(0, 1).Phase));
Assert.AreEqual(135,
MathExtensions.ToDeg(new Complex(-1, 1).Phase));
Assert.AreEqual(180,
MathExtensions.ToDeg(new Complex(-1, 0).Phase));
Assert.AreEqual(225 - 360,
MathExtensions.ToDeg(new Complex(-1, -1).Phase));
Assert.AreEqual(270 - 360,
MathExtensions.ToDeg(new Complex(0, -1).Phase));
Assert.AreEqual(315 - 360,
MathExtensions.ToDeg(new Complex(1, -1).Phase));
}
}
Subskrybuj:
Posty (Atom)