|
9 | 9 |
|
10 | 10 | using namespace nbl;
|
11 | 11 |
|
12 |
| -#ifdef DEPRECATED_ROOT_FINDER |
13 |
| -#include <boost/math/tools/polynomial.hpp> |
14 |
| -using std::complex; |
15 |
| -using boost::math::tools::polynomial; |
16 |
| - |
17 |
| -template<typename T, size_t Order> |
18 |
| -polynomial<complex<T>> polynomialDerivative(polynomial<complex<T>> values) |
19 |
| -{ |
20 |
| - std::array<complex<T>, Order> output; |
21 |
| - for (uint32_t i = 1; i < Order + 1; i++) |
22 |
| - { |
23 |
| - complex<T> originalValue = values[i]; |
24 |
| - uint32_t originalPosition = i; |
25 |
| - output[i - 1] = T(originalPosition) * originalValue; |
26 |
| - } |
27 |
| - return polynomial<complex<T>>(output.data(), Order - 1); |
28 |
| -} |
29 |
| - |
30 |
| -template<typename T, size_t Order> |
31 |
| -complex<T> evaluatePolynomial(polynomial<complex<T>> values, complex<T> x) |
32 |
| -{ |
33 |
| - if constexpr (Order == -1) |
34 |
| - return complex<T>(0.0, 0.0); |
35 |
| - complex<T> acc = values[Order]; |
36 |
| - for (int i = Order - 1; i >= 0; i--) |
37 |
| - acc = acc * x + values[i]; |
38 |
| - |
39 |
| - return acc; |
40 |
| -} |
41 |
| - |
42 |
| -constexpr double LaguerrePolynomialThreshhold = 2e-300; |
43 |
| -constexpr double LaguerreAThreshhold = 1e-40; |
44 |
| -constexpr uint32_t MaxIterations = 256; |
45 |
| -constexpr double LaguerreInitialGuess = 0.0; |
46 |
| -constexpr double EPS = 2e-6; |
47 |
| - |
48 |
| -// https://en.m.wikipedia.org/wiki/Laguerre's_method |
49 |
| -template<typename float_t, size_t Order> |
50 |
| -complex<float_t> solveLaguerreRoot(polynomial<complex<float_t>> p, complex<float_t> initialGuess) |
51 |
| -{ |
52 |
| - using complex_t = complex<float_t>; |
53 |
| - using polynomial_t = polynomial<complex_t>; |
54 |
| - |
55 |
| - polynomial_t pd = polynomialDerivative<float_t, Order>(p); |
56 |
| - polynomial_t pdd = polynomialDerivative <float_t, Order - 1>(pd); |
57 |
| - |
58 |
| - assert(p.size() == Order + 1); |
59 |
| - assert(pd.size() == Order); |
60 |
| - assert(pdd.size() == Order - 1); |
61 |
| - |
62 |
| - const uint32_t n = Order; |
63 |
| - uint32_t k = 0; |
64 |
| - complex_t x = initialGuess; |
65 |
| - complex_t a; |
66 |
| - |
67 |
| - do { |
68 |
| - complex_t f = evaluatePolynomial<float_t, Order>(p, x); // p(x) |
69 |
| - complex_t fd = evaluatePolynomial<float_t, Order - 1>(pd, x); // p'(x) |
70 |
| - complex_t fdd = evaluatePolynomial<float_t, Order - 2>(pdd, x); // p''(x) |
71 |
| - // If p(x) is very small, exit the loop |
72 |
| - if (abs(f) < LaguerrePolynomialThreshhold) break; |
73 |
| - |
74 |
| - // TODO: Using polynomial math (division, multiplication & subtraction) we can turn G and H into |
75 |
| - // polynomials and evaluate them here |
76 |
| - complex_t G = fd / f; // p'(x) / p(x) |
77 |
| - complex_t H = G * G - fdd / f; // G² - p''(x) / p(x) |
78 |
| - |
79 |
| - // this value fed to sqrt could end up negative, meaning we can get |
80 |
| - // an imaginary value with complex sqrt |
81 |
| - complex_t det = std::sqrt(complex_t(n - 1, 0.0) * (complex_t(n, 0.0) * H - G * G)); |
82 |
| - |
83 |
| - // n / (G +- sqrt((n - 1) * (nH - G²))) |
84 |
| - complex_t aSign0 = G + det; |
85 |
| - complex_t aSign1 = G - det; |
86 |
| - a = complex_t(n, 0.0) / (abs(aSign0) > abs(aSign1) ? aSign0 : aSign1); |
87 |
| - //if (core::isnan(a) || isinf(a)) return nbl::core::nan<double>(); |
88 |
| - x -= a; |
89 |
| - k++; |
90 |
| - // Repeat until a is small enough or if the maximum number of iterations has been reached. |
91 |
| - } while (abs(a) >= LaguerreAThreshhold && k < MaxIterations); |
92 |
| - |
93 |
| - // Collapse complex into real variable |
94 |
| - if (abs(x.imag()) < 2.0 * EPS * abs(x.real())) |
95 |
| - return complex_t(x.real(), 0); |
96 |
| - return x; |
97 |
| -} |
98 |
| - |
99 |
| -template<typename T> |
100 |
| -polynomial<complex<T>> deflatePolynomial(polynomial<complex<T>> poly, complex<T> root) |
101 |
| -{ |
102 |
| - // divide by x - root |
103 |
| - // (-root, 1) |
104 |
| - std::array<complex<T>, 2> divisionPoly = { -root, complex<T>(1.0, 0.0) }; |
105 |
| - return poly / polynomial<complex<T>>(divisionPoly.data(), 1); |
106 |
| -} |
107 |
| - |
108 |
| -// https://stackoverflow.com/a/12683219 |
109 |
| -template<size_t Order, typename T> |
110 |
| -class SolvePolynomialRoots |
111 |
| -{ |
112 |
| -public: |
113 |
| - static std::array<complex<T>, Order> solveRoots(polynomial<complex<T>> p) |
114 |
| - { |
115 |
| - std::array<complex<T>, Order> roots = {}; |
116 |
| - polynomial <complex<T>> poly = p; |
117 |
| - |
118 |
| - // Values of 0 may reduce the order of the polynomial, so we always check before performing |
119 |
| - // solveLaguerreRoot for that order |
120 |
| - if (p.size() == Order + 1) |
121 |
| - { |
122 |
| - roots[0] = solveLaguerreRoot<T, Order>(p, LaguerreInitialGuess); |
123 |
| - // If a root has been found, the corresponding linear factor can be removed from p. |
124 |
| - // This deflation step reduces the degree of the polynomial by one, so that eventually, |
125 |
| - // approximations for all roots of p can be found.Note however that deflation can lead |
126 |
| - // to approximate factors that differ significantly from the corresponding exact factors. |
127 |
| - // This error is least if the roots are found in the order of increasing magnitude. |
128 |
| - poly = deflatePolynomial<T>(poly, roots[0]); |
129 |
| - } |
130 |
| - |
131 |
| - auto otherRoots = SolvePolynomialRoots<Order - 1, T>::solveRoots(poly); |
132 |
| - std::copy(otherRoots.begin(), otherRoots.end(), roots.begin() + 1); |
133 |
| - return roots; |
134 |
| - } |
135 |
| -}; |
136 |
| - |
137 |
| -template<typename T> |
138 |
| -class SolvePolynomialRoots<1, T> |
139 |
| -{ |
140 |
| -public: |
141 |
| - static std::array<complex<T>, 1> solveRoots(polynomial<complex<T>> p) |
142 |
| - { |
143 |
| - if (p.size() == 2) |
144 |
| - return { solveLaguerreRoot<T, 1>(p, LaguerreInitialGuess) }; |
145 |
| - return { complex<T>(core::nan<double>()) }; |
146 |
| - } |
147 |
| -}; |
148 |
| - |
149 |
| -template<typename T> |
150 |
| -std::array<complex<T>, 4> solveQuarticRootsLaguerre(std::array<T, 5> p) |
151 |
| -{ |
152 |
| - std::array<complex<T>, 5> polynomialValues = { |
153 |
| - complex<T>(p[0]),complex<T>(p[1]),complex<T>(p[2]),complex<T>(p[3]),complex<T>(p[4]) |
154 |
| - }; |
155 |
| - std::array<complex<T>, 4> roots = SolvePolynomialRoots<4, T>::solveRoots(polynomial <complex<T>>(polynomialValues.data(), 4)); |
156 |
| - |
157 |
| - // Improve initial guess (polish numerical error) |
158 |
| - // TODO: use Halley for this? |
159 |
| - // TODO: find better way of checking for the size |
160 |
| - polynomial<complex<T>> poly(polynomialValues.data(), 4); |
161 |
| - if (poly.size() == 5) |
162 |
| - { |
163 |
| - for (uint32_t i = 0; i < roots.size(); i++) |
164 |
| - { |
165 |
| - roots[i] = solveLaguerreRoot<T, 4>(poly, roots[i]); |
166 |
| - } |
167 |
| - } |
168 |
| - |
169 |
| - return roots; |
170 |
| -} |
171 |
| - |
172 |
| -std::array<double, 4> Hatch::solveQuarticRootsLaguerre(double a, double b, double c, double d, double e, double t_start, double t_end) |
173 |
| -{ |
174 |
| - auto laguerreRoots = solveQuarticRootsLaguerre<double>({ e, d, c, b, a }); |
175 |
| - std::array<double, 4> roots = { -1.0, -1.0, -1.0, -1.0 }; // only two candidates in range, ever |
176 |
| - uint32_t realRootCount = 0; |
177 |
| - |
178 |
| - for (uint32_t i = 0; i < laguerreRoots.size(); i++) |
179 |
| - { |
180 |
| - if (laguerreRoots[i].imag() == 0.0) |
181 |
| - { |
182 |
| - bool duplicatedRoot = false; |
183 |
| - |
184 |
| - for (uint32_t realRoot = 0; realRoot < realRootCount; realRoot++) |
185 |
| - if (roots[realRoot] == laguerreRoots[i].real()) |
186 |
| - duplicatedRoot = true; |
187 |
| - |
188 |
| - if (duplicatedRoot) |
189 |
| - continue; |
190 |
| - |
191 |
| - roots[realRootCount++] = laguerreRoots[i].real(); |
192 |
| - } |
193 |
| - } |
194 |
| - |
195 |
| - return roots; |
196 |
| -} |
197 |
| -#endif |
198 |
| - |
199 | 12 | bool Hatch::Segment::isStraightLineConstantMajor() const
|
200 | 13 | {
|
201 | 14 | auto major = (uint32_t)SelectedMajorAxis;
|
|
0 commit comments