/* polynom.hpp POLYNOM 0.11, simple C++ API for complex polynomial factoring and expansion. Latest version available at: https://ecomaan.nl/cpp/polynom Copyright (c) 2006, 2020 - Pieter Suurmond Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. Any person wishing to distribute modifications to the Software is requested to send the modifications to the original developer so that they can be incorporated into the canonical version. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. Before you include this headerfile, say "using namespace std;" and include and headers. */ int factor (const vector >& coeffs, vector >& roots); int factor_d(const vector >& coeffs, vector >& roots); /* Functions factor() and factor_d() both convert complex polynomials in the expanded form to the factored form. They are identical, except that factor_d() expects coefficients in descending order of degree, whereas factor() expects its' input in ascending order, like this: 2 3 a0 + a1 x + a2 x + a3 x + ... = (x-r0) (x-r1) (x-r2) (... Thus a vector of coefficients (a0, a1, a2, a3, ..) is transformed into a vector of roots (r0, r1, r2, ..). The size of the output vector is usually one less than the input vector's size. However, when one sets coefficient(s) with the highest degree(s) to zero (last coefficient(s) for factor(), or first one(s) for factor_d()), it may happen that we receive less roots than expected. The following kind of equations are solved analytically: a0 + a1 x = 0 (linear equations.) 2 a0 + a1 x + a2 x = 0 (quadratic equations.) 99 a0 + a99 x = 0 (when a1 through a98 are 0, 99 roots are cal- culated easily. See 'de Moivre' and 'roots of unity'.) All other polynomials are first reduced to the second degree using Laguerre's numerical method, synthetically dividing the polynomial after each root found. Input coefficients as well as outputted roots can both be complex. The functions return zero in case of success, and a non-zero error number in case of failure: 0 = Ok. 1 = Less than 2 coefficients supplied (only the first coeff my be zero). 2 = Laguerre's algorithm failed to converge. 3 = Synthetic division left a remainder (>=kEPSILON): Our Laguerre code failed to converge or we reached accuracy-limits. Root-finding example: // Known 4th degree polynomial: 2 4 // -6x - 7x + x = x(x+1)(x+2)(x-3) vector > coeffs(5); // Input real, imag. coeffs[0] = complex( 0.0, 0.0); coeffs[1] = complex(-6.0, 0.0); coeffs[2] = complex(-7.0, 0.0); coeffs[3] = complex( 0.0, 0.0); coeffs[4] = complex( 1.0, 0.0); vector > roots; // Output pre-alloc with 'roots(4)' // is slightly faster. int e = factor(coeffs, // Input in ascending degree. roots); // Output (unordered). // May give: // roots[0] = ( 0.0000000000000, 0.0000000000000) // roots[1] = ( 3.0000000000000, 0.0000000000000) // roots[2] = (-1.0000000000000, 0.0000000000000) // roots[3] = (-2.0000000000000, 0.0000000000000) Because factor() and factor_d() use temporary internal memory, they may fail, throwing a bad_alloc exception. (Perhaps catch this as error #4?) */ int expand (const vector >& roots, vector >& coeffs); int expand_d(const vector >& roots, vector >& coeffs); /* Functions expand() and expand_d() do the opposite of factor() and factor_d() respectively: they convert polynomials in factored form to the expanded form. Function expand() delivers coefficients in ascending order of degree, whereas expand_d() yields coefficients in descending degree. The functions return 1 when not enough roots were supplied; zero on success. Coefficient-calculation example (with the same polynomial): vector > roots(4); // Input. roots[0] = complex( 0.0, 0.0); roots[1] = complex( 3.0, 0.0); roots[2] = complex(-1.0, 0.0); roots[3] = complex(-2.0, 0.0); vector > coeffs; // Output pre-alloc with 'coeffs(5)' // is slightly faster. int e = expand(roots, // Input (unordered). coeffs); // Output in ascending degree. // Gives: // coeffs[0] = ( 0.0000000000000, 0.0000000000000) // coeffs[1] = (-6.0000000000000, 0.0000000000000) // coeffs[2] = (-7.0000000000000, 0.0000000000000) // coeffs[3] = ( 0.0000000000000, 0.0000000000000) // coeffs[4] = ( 1.0000000000000, 0.0000000000000) Sourcefile 'polynom_test.cpp' may serve as another example: it thoroughly tests this 'polynom library' with another million of complex polynomials. Up to the 40th degree or so, factor() and expand() appear to work nicely together (as opposites). */