libpappsomspp
Library for mass spectrometry
Loading...
Searching...
No Matches
cardano.cpp
Go to the documentation of this file.
1/**
2 * \file pappsomspp/vendors/tims/mzcalibration/cardano.cpp
3 * \date 17/12/2022
4 * \brief cubic solver adapted from
5 * https://www.codeproject.com/articles/798474/to-solve-a-cubic-equation
6 thanks
7 * to "Sergey Bochkanov" <sergey.bochkanov@alglib.net> for his advise
8 */
9
10/*******************************************************************************
11 * Copyright (c) 2022 Olivier Langella <Olivier.Langella@u-psud.fr>.
12 *
13 * This file is part of the PAPPSOms++ library.
14 *
15 * PAPPSOms++ is free software: you can redistribute it and/or modify
16 * it under the terms of the GNU General Public License as published by
17 * the Free Software Foundation, either version 3 of the License, or
18 * (at your option) any later version.
19 *
20 * PAPPSOms++ is distributed in the hope that it will be useful,
21 * but WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 * GNU General Public License for more details.
24 *
25 * You should have received a copy of the GNU General Public License
26 * along with PAPPSOms++. If not, see <http://www.gnu.org/licenses/>.
27 *
28 ******************************************************************************/
29
30#include "cardano.h"
31
32
33#include <cmath>
34#include <QDebug>
35
36
37const double BUFFER_SQRT3{std::sqrt(3.0)};
38const double BUFFER_inv27{1.0 / 27.0};
39const double BUFFER_pow11{std::pow(10.0, -11.0)};
40
41void
43 InHousePolynomialSolverResult &res, double a1, double b, double c, double d)
44{
45
46 /**
47 * @todo Cardaono cubic solver
48 *
49 * adapted in c++ from
50 * https://www.codeproject.com/articles/798474/to-solve-a-cubic-equation
51 thanks
52 * to "Sergey Bochkanov" <sergey.bochkanov@alglib.net> for his advise
53 *
54 * Cubic Equation
55 https://github.com/harveytriana/CubicEquation
56 Quartic Equation
57 https://github.com/harveytriana/QuarticEcuation
58 */
59 double a, p, q, u, v;
60 double r, alpha;
61 a = b / a1;
62 b = c / a1;
63 c = d / a1;
64
65 double aover3(a / 3.0);
66 p = -(a * aover3) + b;
67 q = (2.0 * BUFFER_inv27 * a * a * a) - (b * aover3) + c;
68 d = q * q / 4.0 + p * p * p * BUFFER_inv27;
69 if(std::abs(d) < BUFFER_pow11)
70 d = 0;
71 // 3 cases D > 0, D == 0 and D < 0
72 if(d > 1e-20)
73 {
74 double dsqrt = std::sqrt(d);
75 u = std::cbrt(-q / 2.0 + dsqrt);
76 v = std::cbrt(-q / 2.0 - dsqrt);
77 res.x1 = u + v - aover3;
78 /*
79 res.x2.real(-(u + v) / 2.0 - aover3);
80 res.x2.imag(BUFFER_SQRT3 / 2.0 * (u - v));
81 res.x3.real(res.x2.real());
82 res.x3.imag(-res.x2.imag());
83 */
85 }
86 if(std::abs(d) <= 1e-20)
87 {
88 u = std::cbrt(-q / 2.0);
89 v = u;
90 res.x1 = u + v - aover3;
91 // res.x2.real(-(u + v) / 2.0 - aover3);
92 res.type = CardanoResultCase::zerod;
93 }
94 if(d < -1e-20)
95 {
96 r = std::cbrt(std::sqrt(-p * p * p * BUFFER_inv27));
97 alpha = std::atan(std::sqrt(-d) / -q * 2.0);
98 if(q > 0) // if q > 0 the angle becomes 2 * PI - alpha
99 alpha = 2.0 * M_PI - alpha;
100
101 res.x1 =
102 r * (std::cos((6.0 * M_PI - alpha) / 3.0) + std::cos(alpha / 3.0)) -
103 aover3;
104 /*
105 res.x2.real(r * (std::cos((2.0 * M_PI + alpha) / 3.0) +
106 std::cos((4.0 * M_PI - alpha) / 3.0)) -
107 aover3);
108 res.x3.real(r * (std::cos((4.0 * M_PI + alpha) / 3.0) +
109 std::cos((2.0 * M_PI - alpha) / 3.0)) -
110 aover3);
111 */
113 }
114}
115
116
118inHousePolynomialSolve(const std::vector<double> &polynome)
119{
122 std::size_t polynome_size = polynome.size();
123
124 double a, b, c, delta;
125
126 switch(polynome_size)
127 {
128 case 0:
129 break;
130 case 1:
131 break;
132 case 2: // linear// equation ax + b = 0
133 a = polynome[1];
134 b = polynome[0];
135 if(a != 0)
136 {
137 res.x1 = -b / a;
138 res.type = CardanoResultCase::line;
139 }
140 break;
141 case 3:
142 // quadratic equation ax**2 + bx + c = 0
143 a = polynome[2];
144 if(a == 0)
145 return res;
146 b = polynome[1];
147 c = polynome[0];
148
149 if(a == 0)
150 return res;
151 // calculate the discriminant
152 delta = (b * b) - (a * c * 4);
153 // qDebug() << delta;
154 if(delta < 0)
155 {
156 }
157 else if(std::abs(delta) <= 1e-20)
158 {
159 res.x1 = -b / (a * 2);
160 // qDebug() << x1.real() << " " << x2.real();
162 }
163 else
164 {
165 // find two solutions
166 delta = std::sqrt(delta);
167 res.x1 = (-b + delta) / (a * 2);
168 res.x2 = (-b - delta) / (a * 2);
169
170 // qDebug() << x1.real() << " " << x2.real();
172 }
173 break;
174 case 4:
175 cubic_solver(res, polynome[3], polynome[2], polynome[1], polynome[0]);
176 }
177 return res;
178}
const double BUFFER_SQRT3
Definition cardano.cpp:37
void cubic_solver(InHousePolynomialSolverResult &res, double a1, double b, double c, double d)
Definition cardano.cpp:42
const double BUFFER_inv27
Definition cardano.cpp:38
const double BUFFER_pow11
Definition cardano.cpp:39
InHousePolynomialSolverResult inHousePolynomialSolve(const std::vector< double > &polynome)
Definition cardano.cpp:118
cubic solver adapted from https://www.codeproject.com/articles/798474/to-solve-a-cubic-equation thank...