/* Try to solve a polynomial, using home-grown techniques. * By: Dirk Mittler * Date: January 29, 2018 */ #include #include #include #include #include #define MIN_A0 0.01 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; 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, \ 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; float 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; } // Do the actual work. bool search_root(Polynomial *a, CC *root, double epsilon = 1.0E-12) { int N = a->deg; if (N < 1) return false; CC xa = -1.0; CC xb = cexp(M_PI * 1i / (N + 1)); CC ya = apply_poly(a, xa); CC yb = apply_poly(a, xb); double spread = abs(yb - ya); 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; double epsilond = epsilon; CC t = 0; CC xt = 0; // Really bad form below! while (true) { if (N > 1) { epsilond = std::max(N * (pow(abs(xa), N - 1)), 1.0 / N) * epsilon; } if (abs(ya) <= epsilond || abs(yb - ya) <= epsilond) { *root = xa; return true; } t = (-ya) / (yb - ya); xt = (t * xb) + ((1.0 - t) * xa); xb = xt + (CC(0.0, 0.5) * (xa - xt)); xa = xt; ya = apply_poly(a, xa); yb = apply_poly(a, xb); if (abs(ya) >= spread) { misses++; } if (misses > N + 2) { return false; } spread = abs(ya); } } // Return number of derivatives... int create_deriv_list(Polynomial a, Polynomial **l) { CC coeff = 0.0; *l = new Polynomial[std::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 std::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, \ 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); } // Finally, the driver: bool poly_solve(Polynomial *a, Polynomial **PL, bool polish) { create_deriv_list(*a, PL); int s = 0; double epsilon = 1.0E-12; Polynomial *L = *PL; CC current_root = 0.0; CC current_conjugate = 0.0; 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 (imag(current_root) >= epsilon || \ 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 (imag(current_root) < epsilon && \ imag(current_root) > -epsilon) { current_root = CC(real(current_root), 0.0); } if (real(current_root) < epsilon && \ 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; cout << setprecision(14); A.offset = 0; A.deg = 3; A.coeffs = new CC[4]; A.coeffs[0] = CC(+1.0, 0.0); A.coeffs[1] = CC(-3.0, 0.0); A.coeffs[2] = CC(-2.0, 0.0); A.coeffs[3] = CC(+6.0, 0.0); cout << "x^3 -3x^2 -2x + 6 = 0 :" << endl; poly_solve(&A, &PL, true); destroy_deriv_list(PL); delete[] A.coeffs; A.offset = 0; A.deg = 3; A.coeffs = new CC[4]; A.coeffs[0] = CC(+1.0, 0.0); A.coeffs[1] = CC(-3.0, 0.0); A.coeffs[2] = CC(-2.0, 0.0); A.coeffs[3] = CC(+15.0, 0.0); cout << "x^3 -3x^2 -2x + 15 = 0 :" << endl; poly_solve(&A, &PL, true); destroy_deriv_list(PL); delete[] A.coeffs; A.offset = 0; A.deg = 4; A.coeffs = new CC[5]; A.coeffs[0] = CC(+1.0, 0.0); A.coeffs[1] = CC(0.0, 0.0); A.coeffs[2] = CC(+2.0, 0.0); A.coeffs[3] = CC(0.0, 0.0); A.coeffs[4] = CC(+1.0, 0.0); cout << "x^4 +2x^2 + 1 = 0 :" << endl; poly_solve(&A, &PL, true); destroy_deriv_list(PL); delete[] A.coeffs; A.offset = 0; A.deg = 5; A.coeffs = new CC[6]; A.coeffs[0] = CC(+1.0, 0.0); A.coeffs[1] = CC(0.0, 0.0); A.coeffs[2] = CC(0.0, 0.0); A.coeffs[3] = CC(0.0, 0.0); A.coeffs[4] = CC(0.0, 0.0); A.coeffs[5] = CC(0.0, 0.0); cout << "x^5 = 0 :" << endl; poly_solve(&A, &PL, true); destroy_deriv_list(PL); delete[] A.coeffs; A.offset = 0; A.deg = 8; A.coeffs = new CC[9]; A.coeffs[0] = CC(+1.0, 0.0); A.coeffs[1] = CC(0.0, 0.0); A.coeffs[2] = CC(+3.0, 0.0); A.coeffs[3] = CC(-3.0, 0.0); A.coeffs[4] = CC(+1.0, 0.0); A.coeffs[5] = CC(+1.0, 0.0); A.coeffs[6] = CC(+2.0, 0.0); A.coeffs[7] = CC(0.0, 0.0); A.coeffs[8] = CC(1.0, 0.0); cout << "x^8 + 3x^6 -3x^5 +x^4" << endl; cout << " +x^3 +2x^2 + 1 = 0 :" << endl; poly_solve(&A, &PL, true); destroy_deriv_list(PL); delete[] A.coeffs; return 0; }