/* Try to solve a polynomial, using home-grown techniques. * By: Dirk Mittler * Date: September 22, 2019 * * First Windows-patch added: May 16, 2020 * Second Windows-patch added: May 18, 2020 * */ #include #include #include #include #include #include #define MIN_A0 0.01 #ifdef _WIN32 #define M_PI 3.14159265358979323846264338327950288 #endif #ifndef max #define max(a, b) ((a > b) ? a : b) #endif using std::complex; using namespace std::complex_literals; using std::endl; using std::cout; using std::setprecision; typedef complex CC; struct Polynomial { CC* coeffs; int deg; int offset; CC* roots; int solved; }; void normalize(Polynomial* a) { int array_end = a->offset + a->deg; CC first_coeff = a->coeffs[a->offset]; a->coeffs[a->offset] = 1.0; for (int i = a->offset + 1; i <= array_end; i++) { a->coeffs[i] /= first_coeff; } return; } // Apply polynomial to parameter. CC apply_poly(Polynomial* a, CC param) { int array_end = a->offset + a->deg; CC power = 1.0; CC result = (a->coeffs)[array_end]; for (int i = 1; i <= a->deg; i++) { power = 1.0; for (int j = 1; j <= i; j++) { power *= param; } result += power * (a->coeffs)[array_end - i]; } return result; } // In-place division of polynomial. // Returns degree of remaining plynomial in same array. int divide(Polynomial* a, CC root, CC* remaind) { CC multiple = 0; int array_end = a->offset + a->deg; for (int i = a->offset; i < array_end; i++) { multiple = a->coeffs[i]; a->coeffs[i + 1] += multiple * root; } *remaind = a->coeffs[array_end]; a->deg--; return a->deg; } // Augmented division of polynomial, // may or may not polish. bool try_root(Polynomial* a, CC* root, bool polish, \ long double epsilon = 1.0E-12) { CC adjustment = 0.0; CC remaind = 0.0; CC deriv = 0.0; if (a->deg < 1) return false; int N = a->deg; double deriv2_a0 = 0.0; int passes = 1; if (N == 1) polish = false; if (polish) passes = 3; int array_end = a->offset + N; Polynomial quot; quot.coeffs = new CC[array_end + 1]; while (passes-- > 0) { *root += adjustment; quot.offset = a->offset; quot.deg = N; for (int i = 0; i <= array_end; i++) { quot.coeffs[i] = a->coeffs[i]; } // (1.0)^0 is possible here... deriv2_a0 = N * N * pow(norm(*root) + 1.0, N - 1); if (norm(apply_poly(", *root)) > \ epsilon * epsilon * deriv2_a0) { delete[] quot.coeffs; return false; } divide(", *root, &remaind); adjustment = CC(0.0, 0.0); deriv = -apply_poly(", *root); if (polish && norm(deriv) > abs(remaind)) { adjustment = remaind / deriv; } } for (int i = 0; i < array_end; i++) { a->coeffs[i] = quot.coeffs[i]; } a->deg = quot.deg; a->offset = quot.offset; delete[] quot.coeffs; return true; } // Solver for linear case. // Assumes that the equation has been normalized. inline bool solve_lin(Polynomial* a, CC* root) { if (a->deg != 1) return false; *root = -a->coeffs[a->offset + 1]; return true; } // Do the actual work. bool search_root(Polynomial* a, CC* root, long double epsilon = 1.0E-12) { int N = a->deg; if (N < 1) return false; // !! if (N == 1) return solve_lin(a, root); CC xa = -1.0; // CC xb = exp(CC(0.0, M_PI / (N + 1))); CC xb = exp((M_PI / (N + 1)) * 1i); CC ya = apply_poly(a, xa); CC yb = apply_poly(a, xb); double spread = abs(yb - ya); double spread2 = 0.0; if (N <= 4 && pow(spread, N) < epsilon) { return false; } else if (N > 4 && pow(spread, N - 2) < N * N * epsilon) { return false; } int misses = 0; long double epsilond = epsilon; CC t = 0.0; CC xt = 0.0; // Really bad form: while (true) while (true) { epsilond = N * (pow(abs(xa), N - 1)) * epsilon; if (epsilond < epsilon / N) epsilond = epsilon / N; if (norm(ya) <= epsilond * epsilond) { *root = xa; return true; } if (norm(ya - yb) <= epsilon * epsilon * epsilon) { return false; } t = (-ya) / (yb - ya); xt = (t * xb) + ((CC(1.0, 0.0) - t) * xa); xb = xt + (CC(0.0, 0.5) * (xa - xt)); xa = xt; ya = apply_poly(a, xa); yb = apply_poly(a, xb); spread2 = abs(ya); if (spread2 >= spread) { misses++; } if (misses > N + 2) { return false; } spread = spread2; } } // Return number of derivatives... int create_deriv_list(Polynomial a, Polynomial** l) { CC coeff = 0.0; *l = new Polynomial[(max(a.deg, 1))]; Polynomial* L = *l; int g_offset = a.offset; int array_end = a.deg + g_offset; L[0].deg = a.deg; L[0].offset = g_offset; L[0].coeffs = new CC[array_end + 1]; L[0].solved = 0; L[0].roots = new CC[a.deg]; for (int i = 0; i <= array_end; i++) { L[0].coeffs[i] = a.coeffs[i]; } if (abs(L[0].coeffs[g_offset]) < MIN_A0) { L[0].coeffs[g_offset] = CC(1.0, 0.0); } normalize(L); for (int i = 1; i < a.deg; i++) { L[i].deg = L[i - 1].deg - 1; L[i].offset = g_offset; L[i].coeffs = new CC[array_end + 1]; L[i].solved = 0; L[i].roots = new CC[L[i].deg]; for (int j = 0; j < L[i - 1].deg; j++) { coeff = L[i - 1].coeffs[L[i].deg + g_offset - j]; L[i].coeffs[L[i].deg + g_offset - j] = coeff * CC(1.0 + j, 0.0); } normalize(L + i); } return max(a.deg - 1, 0); } void destroy_deriv_list(Polynomial* l) { int array_num = l[0].deg; for (int i = 0; i < array_num; i++) { delete[] l[i].coeffs; delete[] l[i].roots; } delete[] l; return; } bool climb_root(Polynomial* L, CC* root, int i, int mult, \ long double epsilon = 1.0E-12) { int s = 0; int mult_temp = 0; bool last_attempt = false; if (i < 0) return true; while ((last_attempt = try_root(L + i, root, false, epsilon))) { s = L[i].solved; L[i].roots[s] = *root; L[i].solved++; if (++mult_temp >= mult) break; } if (!last_attempt) return false; return climb_root(L, root, i - 1, mult + 1, epsilon); } // Define absolute of real number. // Should not be necessary. template T abs_d(T x) { if (x < 0) return -x; return x; } // Finally, the driver: bool poly_solve(Polynomial* a, Polynomial** PL, bool polish) { create_deriv_list(*a, PL); int s = 0; long double epsilon = 1.0E-12; Polynomial* L = *PL; CC current_root = 0.0; CC current_conjugate = 0.0; if (sizeof(CC) == 32) { epsilon = 1.0L; int limits_10 = std::numeric_limits::digits10; for (int i = 3; i < limits_10; i++) epsilon /= 10.0L; } for (int i = a->deg - 1; i >= 0; i--) { while (search_root(L + i, ¤t_root, epsilon) && \ try_root(L + i, ¤t_root, polish, epsilon)) { s = L[i].solved; L[i].roots[s] = current_root; L[i].solved++; climb_root(L, ¤t_root, i - 1, 2, epsilon); if (abs_d(imag(current_root)) > epsilon) { current_conjugate = conj(current_root); if (try_root(L + i, ¤t_conjugate, false, epsilon)) { s = L[i].solved; L[i].roots[s] = current_conjugate; L[i].solved++; climb_root(L, ¤t_conjugate, i - 1, 2, epsilon); } } } } // Improvised feature: // Print solutions. // This wouldn't be done in actual software. // I.e., Edit out. cout << endl; for (int i = 0; i < L[0].solved; i++) { current_root = L[0].roots[i]; // Secondary Polishing... if (abs_d(imag(current_root)) <= epsilon) { current_root = CC(real(current_root), 0.0); } if (abs_d(real(current_root)) <= epsilon) { current_root = CC(0.0, imag(current_root)); } cout << real(current_root) << " + "; cout << imag(current_root) << "I"; cout << endl; } cout << endl; cout << endl; return true; } int main(int argc, char* argv[]) { Polynomial* PL = NULL; Polynomial A; CC temp = 0.0; char* comma = NULL; int arg_offset = 0; bool polish = true; if (argc > 1 && *argv[1] == '!') { arg_offset = 1; polish = false; } cout << setprecision(14); if (sizeof(CC) == 32) { cout << setprecision(std::numeric_limits::digits10 - 2); } if (argc - arg_offset < 3) { cout << "Usage:" << endl; cout << "./roots_2 [!] coeff1 coeff2 coeff3 coeff4 ..." << endl; cout << "In the above case, prints the solutions to:" << endl; cout << "coeff1*x^3 + coeff2*x^2 + coeff3*x + coeff4 = 0.0" << endl; cout << "All command-line arguments must be numbers," << endl; cout << "separated by spaces. The only way to omit a term," << endl; cout << "is to put a zero." << endl; cout << "Minimum degree: 1 (a linear equation)." << endl; cout << endl; cout << "Because all the polynomials need to be normalized," << endl; cout << "i.e., the coefficient of the highest exponent" << endl; cout << "made equal to 1.0, values deemed too low as the first" << endl; cout << "argument, will arbitrarily be set to +1.0 ." << endl; cout << endl; cout << "Precede coefficients-list with an exclamation-mark," << endl; cout << "i.e., with a ! , to turn polishing off." << endl; cout << "Sometimes, polishing can break finding roots." << endl; cout << endl; cout << "Imaginary components can be added to coefficients, by typing" << endl; cout << "a comma after the real component, followed without space" << endl; cout << "by the imaginary component(s)." << endl; cout << endl; return 1; } A.offset = 0; A.deg = argc - arg_offset - 2; A.coeffs = new CC[A.deg + 1]; for (int c = 0; c <= A.deg; c++) { temp = CC(strtod(argv[c + arg_offset + 1], &comma), 0.0); if (*comma == ',') { temp = CC(real(temp), strtod(comma + 1, NULL)); } A.coeffs[c] = temp; } poly_solve(&A, &PL, polish); destroy_deriv_list(PL); delete[] A.coeffs; return 0; }