/*****************************************************************************/ /* */ /* Copyright (C) 1999-2023 M. Held */ /* */ /* This code is not in the public domain. All rights reserved! Please make */ /* sure to read the full copyright statement contained in "README.txt" or in */ /* the "main" file of this code, such as "main.cc". */ /* */ /*****************************************************************************/ #ifndef VRONI_ROOTS_H #define VRONI_ROOTS_H int Roots3(double lead, double a, double b, double c, double roots[]); int Roots4(double lead, double a, double b, double c, double d, double roots[]); #define ROOTS_ZERO 1.0e-50 #define ROOTS_ZERO2 1.0e-100 //#define ROOTS_ZERO 1.0e-14 //#define ROOTS_ZERO2 1.0e-28 #define DISCR_ZERO 1.0e-8 #define ROOT_SMALL 0.125 /* ROOT_SMALL = 1/8 >> ZERO */ #define ROOT_INV_SMALL 8.0 /* ROOT_INV_SMALL = 1/ROOT_SMALL */ /* */ /* This macro solves the following second-degree polynomial equation: */ /* */ /* a * x^2 + b * x + c = 0. */ /* */ /* The roots are stored in roots[0] and roots[1]. Note that only real */ /* roots are sought. The number of real roots found is stored in num_roots. */ /* */ #define Roots2abc(a, b, c, roots, num_roots) \ { \ while (eq(a, ROOT_SMALL) && eq(b, ROOT_SMALL) && eq(c, ROOT_SMALL) && ((a != 0.0) || (b != 0.0) || (c != 0.0))) { \ a *= 2.0; \ b *= 2.0; \ c *= 2.0; \ } \ while ((Abs(a) > ROOT_INV_SMALL) && (Abs(b) > ROOT_INV_SMALL) && (Abs(c) > ROOT_INV_SMALL)) { \ a /= 2.0; \ b /= 2.0; \ c /= 2.0; \ } \ if (eq(a, ROOTS_ZERO)) { \ if (eq(b, ROOTS_ZERO)) { \ if (eq(c, ROOTS_ZERO)) { \ num_roots = -1; \ } \ else { \ num_roots = 0; \ } \ } \ else { \ roots[0] = - c / b; \ num_roots = 1; \ } \ } \ else { \ basic_h_local_delta = b * b - 4 * a * c; \ if (basic_h_local_delta > 0.0) { \ if (b > 0) { \ basic_h_local = - 0.5 * (b + sqrt(basic_h_local_delta)); \ } \ else { \ basic_h_local = - 0.5 * (b - sqrt(basic_h_local_delta)); \ } \ if (eq(basic_h_local, ROOTS_ZERO)) { \ roots[0] = basic_h_local / a; \ num_roots = 1; \ } \ else { \ roots[0] = basic_h_local / a; \ roots[1] = c / basic_h_local; \ num_roots = 2; \ } \ } \ else if (eq(basic_h_local_delta, DISCR_ZERO)) { \ roots[0] = - b / (2.0 * a); \ num_roots = 1; \ } \ else { \ num_roots = 0; \ } \ } \ } /* */ /* This macro solves the following second-degree polynomial equation: */ /* */ /* x^2 + p * x + q = 0. */ /* */ /* The roots are stored in roots[0] and roots[1]. Note that only real */ /* roots are sought. The number of real roots found is stored in num_roots. */ /* */ #define Roots2pq(p, q, roots, num_roots) \ { \ basic_h_local_delta = p * p - 4 * q; \ if (basic_h_local_delta > 0.0) { \ if (p > 0) { \ basic_h_local = - 0.5 * (p + sqrt(basic_h_local_delta)); \ } \ else { \ basic_h_local = - 0.5 * (p - sqrt(basic_h_local_delta)); \ } \ if (eq(basic_h_local, ROOTS_ZERO)) { \ roots[0] = basic_h_local; \ num_roots = 1; \ } \ else { \ roots[0] = basic_h_local; \ roots[1] = q / basic_h_local; \ num_roots = 2; \ } \ } \ else if (eq(basic_h_local_delta, ROOTS_ZERO2)) { \ roots[0] = - 0.5 * p; \ num_roots = 1; \ } \ else { \ num_roots = 0; \ } \ } /* */ /* This macro solves the following 2x2 linear system: */ /* */ /* A[0][0] * x + A[0][1] * y = B[0] */ /* A[1][1] * x + A[1][1] * y = B[1] */ /* */ /* If a unique solution exists, then exists := 1, and the solution is stored */ /* in xy[2]. If the solution is not unique, then exists := 2, and a solution */ /* is stored in xy[2]. Otherwise, exists := 0. */ /* */ /* i, j, I, J are dummy integers needed within the macro. */ /* */ #define LinearEqnSolver_2x2(A, B, xy, exists, i, j, I, J) \ { \ /* */ \ /* find a column with a non-zero element */ \ /* */ \ exists = 0; \ if (!eq((A)[0][0], ROOTS_ZERO) || !eq((A)[1][0], ROOTS_ZERO)) { \ I = 0; \ J = 1; \ } \ else if (!eq((A)[0][1], ROOTS_ZERO) || !eq((A)[1][1], ROOTS_ZERO)) { \ I = 1; \ J = 0; \ } \ else { \ if (eq((B)[0], ROOTS_ZERO) && eq((B)[1], ROOTS_ZERO)) { \ (xy)[0] = (xy)[1] = 0.0; \ (exists) = 2; \ } \ I = J = 0; \ } \ \ /* */ \ /* determine i s.t. Abs(A[i][I]) is maximum. */ \ /* */ \ if ((I > 0) || (J > 0)) { \ if (Abs((A)[0][I]) > Abs((A)[1][I])) { \ i = 0; \ j = 1; \ } \ else { \ i = 1; \ j = 0; \ } \ \ basic_h_local_quot = (A)[j][I] / (A)[i][I]; \ basic_h_local_delta = (A)[j][J] - basic_h_local_quot * (A)[i][J]; \ if (!eq(basic_h_local_delta, ROOTS_ZERO)) { \ (xy)[J] = ((B)[j] - basic_h_local_quot * (B)[i]) / basic_h_local_delta; \ (xy)[I] = ((B)[i] - (xy)[J] * (A)[i][J]) / (A)[i][I]; \ (exists) = 1; \ } \ else { \ basic_h_local_delta = (B)[j] - basic_h_local_quot * (B)[i]; \ if (eq(basic_h_local_delta, ROOTS_ZERO)) { \ (xy)[J] = 0.0; \ (xy)[I] = (B)[i] / (A)[i][I]; \ (exists) = 2; \ } \ else { \ (exists) = 0; \ } \ } \ } \ } /* */ /* This macro solves the following 2x2 linear system: */ /* */ /* A[0][0] * x + A[0][1] * y = B[0] */ /* A[1][1] * x + A[1][1] * y = B[1] */ /* */ /* If a unique solution exists, then exists := 1, and the solution is stored */ /* in xy[2]. If the solution is not unique, then exists := 2, and a solution */ /* is stored in xy[2]. Otherwise, exists := 0. */ /* */ /* i, j, I, J are dummy integers needed within the macro. */ /* */ #define LinearEqnSolver_2x2_Zero(A, B, xy, exists, i, j, I, J) \ { \ /* */ \ /* find a column with a non-zero element */ \ /* */ \ exists = 0; \ if (!eq((A)[0][0], ROOTS_ZERO) || !eq((A)[1][0], ROOTS_ZERO)) { \ I = 0; \ J = 1; \ } \ else if (!eq((A)[0][1], ROOTS_ZERO) || !eq((A)[1][1], ROOTS_ZERO)) { \ I = 1; \ J = 0; \ } \ else { \ if (eq((B)[0], ROOTS_ZERO) && eq((B)[1], ROOTS_ZERO)) { \ (xy)[0] = (xy)[1] = 0.0; \ (exists) = 2; \ } \ I = J = 0; \ } \ \ /* */ \ /* determine i s.t. Abs(A[i][I]) is maximum. */ \ /* */ \ if ((I > 0) || (J > 0)) { \ if (Abs((A)[0][I]) > Abs((A)[1][I])) { \ i = 0; \ j = 1; \ } \ else { \ i = 1; \ j = 0; \ } \ \ basic_h_local_quot = (A)[j][I] / (A)[i][I]; \ basic_h_local_delta = (A)[j][J] - basic_h_local_quot * (A)[i][J]; \ if (!eq(basic_h_local_delta, ROOTS_ZERO)) { \ (xy)[J] = ((B)[j] - basic_h_local_quot * (B)[i]) / basic_h_local_delta; \ (xy)[I] = ((B)[i] - (xy)[J] * (A)[i][J]) / (A)[i][I]; \ (exists) = 1; \ } \ else { \ basic_h_local_delta = (B)[j] - basic_h_local_quot * (B)[i]; \ if (eq(basic_h_local_delta, ROOTS_ZERO)) { \ (xy)[J] = 0.0; \ (xy)[I] = (B)[i] / (A)[i][I]; \ (exists) = 2; \ } \ else { \ (exists) = 0; \ } \ } \ } \ } #endif