Eclipse SUMO - Simulation of Urban MObility
Loading...
Searching...
No Matches
PolySolver.cpp
Go to the documentation of this file.
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/****************************************************************************/
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// ===========================================================================
33std::tuple<int, double, double>
34PolySolver::quadraticSolve(double a, double b, double c) {
35 // ax^2 + bx + c = 0
36 // return only real roots
37 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 if (b == 0) {
40 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 return std::make_tuple(1, -c / b, NAN);
47 }
48 }
49
50 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 return std::make_tuple(2, 0, -b / a);
53 }
54
55 double disc, x1_real, x2_real ;
56 disc = b * b - 4 * a * c;
57
58 if (disc > 0) {
59 // two different real roots
60 x1_real = (-b + sqrt(disc)) / (2 * a);
61 x2_real = (-b - sqrt(disc)) / (2 * a);
62 return std::make_tuple(2, x1_real, x2_real);
63 } else if (disc == 0) {
64 // a real root with multiplicity 2
65 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
73std::tuple<int, double, double, double>
74PolySolver::cubicSolve(double a, double b, double c, double d) {
75 // ax^3 + bx^2 + cx + d = 0
76 // return only real roots
77 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 std::tie(numX, x1, x2) = quadraticSolve(b, c, d);
82 return std::make_tuple(numX, x1, x2, NAN);
83 }
84 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 std::tie(numX, x1, x2) = quadraticSolve(a, b, c);
89 return std::make_tuple(numX + 1, 0, x1, x2);
90 }
91
92 b /= a;
93 c /= a;
94 d /= a;
95
96 double disc, q, r, dum1, s, t, term1, r13;
97 q = (3.0 * c - (b * b)) / 9.0;
98 r = -(27.0 * d) + b * (9.0 * c - 2.0 * (b * b));
99 r /= 54.0;
100 disc = q * q * q + r * r;
101 term1 = (b / 3.0);
102
103 double x1_real, x2_real, x3_real;
104 if (disc > 0) { // One root real, two are complex
105 s = r + sqrt(disc);
106 s = s < 0 ? -cbrt(-s) : cbrt(s);
107 t = r - sqrt(disc);
108 t = t < 0 ? -cbrt(-t) : cbrt(t);
109 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 else if (disc == 0) { // All roots real, at least two are equal.
116 r13 = r < 0 ? -cbrt(-r) : cbrt(r);
117 x1_real = -term1 + 2.0 * r13;
118 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 q = -q;
124 dum1 = q * q * q;
125 dum1 = acos(r / sqrt(dum1));
126 r13 = 2.0 * sqrt(q);
127 x1_real = -term1 + r13 * cos(dum1 / 3.0);
128 x2_real = -term1 + r13 * cos((dum1 + 2.0 * M_PI) / 3.0);
129 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/****************************************************************************/
static std::tuple< int, double, double, double > cubicSolve(double a, double b, double c, double d)
Solver of cubic equation ax^3 + bx^2 + cx + d = 0.
static std::tuple< int, double, double > quadraticSolve(double a, double b, double c)
Solver of quadratic equation ax^2 + bx + c = 0.
#define M_PI
Definition odrSpiral.cpp:45