Line data Source code
1 : /****************************************************************************/
2 : // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
3 : // Copyright (C) 2001-2024 German Aerospace Center (DLR) and others.
4 : // This program and the accompanying materials are made available under the
5 : // terms of the Eclipse Public License 2.0 which is available at
6 : // https://www.eclipse.org/legal/epl-2.0/
7 : // This Source Code may also be made available under the following Secondary
8 : // Licenses when the conditions for such availability set forth in the Eclipse
9 : // Public License 2.0 are satisfied: GNU General Public License, version 2
10 : // or later which is available at
11 : // https://www.gnu.org/licenses/old-licenses/gpl-2.0-standalone.html
12 : // SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
13 : /****************************************************************************/
14 : /// @file PolySolver.cpp
15 : /// @author Jakub Sevcik (RICE)
16 : /// @author Jan Prikryl (RICE)
17 : /// @date 2019-12-06
18 : ///
19 : //
20 : /****************************************************************************/
21 : #include <config.h>
22 :
23 : #include <math.h>
24 : #include <cmath>
25 : #include <limits>
26 : #include <utils/geom/GeomHelper.h> // defines M_PI
27 : #include "PolySolver.h"
28 :
29 :
30 : // ===========================================================================
31 : // member method definitions
32 : // ===========================================================================
33 : std::tuple<int, double, double>
34 0 : PolySolver::quadraticSolve(double a, double b, double c) {
35 : // ax^2 + bx + c = 0
36 : // return only real roots
37 0 : if (a == 0) {
38 : //WRITE_WARNING(TL("The coefficient of the quadrat of x is 0. Please use the utility for a FIRST degree linear. No further action taken."));
39 0 : if (b == 0) {
40 0 : if (c == 0) {
41 : return std::make_tuple(2, std::numeric_limits<double>::infinity(), -std::numeric_limits<double>::infinity());
42 : } else {
43 : return std::make_tuple(0, NAN, NAN);
44 : }
45 : } else {
46 0 : return std::make_tuple(1, -c / b, NAN);
47 : }
48 : }
49 :
50 0 : if (c == 0) {
51 : //WRITE_WARNING(TL("One root is 0. Now divide through by x and use the utility for a FIRST degree linear to solve the resulting equation for the other one root. No further action taken."));
52 0 : return std::make_tuple(2, 0, -b / a);
53 : }
54 :
55 : double disc, x1_real, x2_real ;
56 0 : disc = b * b - 4 * a * c;
57 :
58 0 : if (disc > 0) {
59 : // two different real roots
60 0 : x1_real = (-b + sqrt(disc)) / (2 * a);
61 0 : x2_real = (-b - sqrt(disc)) / (2 * a);
62 : return std::make_tuple(2, x1_real, x2_real);
63 0 : } else if (disc == 0) {
64 : // a real root with multiplicity 2
65 0 : x1_real = (-b + sqrt(disc)) / (2 * a);
66 : return std::make_tuple(1, x1_real, NAN);
67 : } else {
68 : // all roots are complex
69 : return std::make_tuple(0, NAN, NAN);
70 : }
71 : }
72 :
73 : std::tuple<int, double, double, double>
74 71 : PolySolver::cubicSolve(double a, double b, double c, double d) {
75 : // ax^3 + bx^2 + cx + d = 0
76 : // return only real roots
77 71 : if (a == 0) {
78 : //WRITE_WARNING(TL("The coefficient of the cube of x is 0. Please use the utility for a SECOND degree quadratic. No further action taken."));
79 : int numX;
80 : double x1, x2;
81 0 : std::tie(numX, x1, x2) = quadraticSolve(b, c, d);
82 : return std::make_tuple(numX, x1, x2, NAN);
83 : }
84 71 : if (d == 0) {
85 : //WRITE_WARNING(TL("One root is 0. Now divide through by x and use the utility for a SECOND degree quadratic to solve the resulting equation for the other two roots. No further action taken."));
86 : int numX;
87 : double x1, x2;
88 0 : std::tie(numX, x1, x2) = quadraticSolve(a, b, c);
89 0 : return std::make_tuple(numX + 1, 0, x1, x2);
90 : }
91 :
92 71 : b /= a;
93 71 : c /= a;
94 71 : d /= a;
95 :
96 : double disc, q, r, dum1, s, t, term1, r13;
97 71 : q = (3.0 * c - (b * b)) / 9.0;
98 71 : r = -(27.0 * d) + b * (9.0 * c - 2.0 * (b * b));
99 71 : r /= 54.0;
100 71 : disc = q * q * q + r * r;
101 71 : term1 = (b / 3.0);
102 :
103 : double x1_real, x2_real, x3_real;
104 71 : if (disc > 0) { // One root real, two are complex
105 0 : s = r + sqrt(disc);
106 0 : s = s < 0 ? -cbrt(-s) : cbrt(s);
107 0 : t = r - sqrt(disc);
108 0 : t = t < 0 ? -cbrt(-t) : cbrt(t);
109 0 : x1_real = -term1 + s + t;
110 : term1 += (s + t) / 2.0;
111 : x3_real = x2_real = -term1;
112 : return std::make_tuple(1, x1_real, NAN, NAN);
113 : }
114 : // The remaining options are all real
115 71 : else if (disc == 0) { // All roots real, at least two are equal.
116 0 : r13 = r < 0 ? -cbrt(-r) : cbrt(r);
117 0 : x1_real = -term1 + 2.0 * r13;
118 0 : x3_real = x2_real = -(r13 + term1);
119 : return std::make_tuple(2, x1_real, x2_real, NAN);
120 : }
121 : // Only option left is that all roots are real and unequal (to get here, q < 0)
122 : else {
123 71 : q = -q;
124 : dum1 = q * q * q;
125 71 : dum1 = acos(r / sqrt(dum1));
126 71 : r13 = 2.0 * sqrt(q);
127 71 : x1_real = -term1 + r13 * cos(dum1 / 3.0);
128 71 : x2_real = -term1 + r13 * cos((dum1 + 2.0 * M_PI) / 3.0);
129 71 : x3_real = -term1 + r13 * cos((dum1 + 4.0 * M_PI) / 3.0);
130 : return std::make_tuple(3, x1_real, x2_real, x3_real);
131 : }
132 : }
133 :
134 :
135 : /****************************************************************************/
|
