Loading...
Searching...
No Matches
cardano.cpp
Go to the documentation of this file.
/**
* \file pappsomspp/vendors/tims/mzcalibration/cardano.cpp
* \date 17/12/2022
* \brief cubic solver adapted from
* https://www.codeproject.com/articles/798474/to-solve-a-cubic-equation
thanks
* to "Sergey Bochkanov" <sergey.bochkanov@alglib.net> for his advise
*/
/*******************************************************************************
* Copyright (c) 2022 Olivier Langella <Olivier.Langella@u-psud.fr>.
*
* This file is part of the PAPPSOms++ library.
*
* PAPPSOms++ is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* PAPPSOms++ is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with PAPPSOms++. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#include "cardano.h"
#include <cmath>
#include <QDebug>
const double BUFFER_SQRT3{std::sqrt(3.0)};
const double BUFFER_inv27{1.0 / 27.0};
const double BUFFER_pow11{std::pow(10.0, -11.0)};
void
cubic_solver(
InHousePolynomialSolverResult &res, double a1, double b, double c, double d)
{
/**
* @todo Cardaono cubic solver
*
* adapted in c++ from
* https://www.codeproject.com/articles/798474/to-solve-a-cubic-equation
thanks
* to "Sergey Bochkanov" <sergey.bochkanov@alglib.net> for his advise
*
* Cubic Equation
https://github.com/harveytriana/CubicEquation
Quartic Equation
https://github.com/harveytriana/QuarticEcuation
*/
double a, p, q, u, v;
double r, alpha;
a = b / a1;
b = c / a1;
c = d / a1;
double aover3(a / 3.0);
p = -(a * aover3) + b;
q = (2.0 * BUFFER_inv27 * a * a * a) - (b * aover3) + c;
d = q * q / 4.0 + p * p * p * BUFFER_inv27;
if(std::abs(d) < BUFFER_pow11)
d = 0;
// 3 cases D > 0, D == 0 and D < 0
if(d > 1e-20)
{
double dsqrt = std::sqrt(d);
u = std::cbrt(-q / 2.0 + dsqrt);
v = std::cbrt(-q / 2.0 - dsqrt);
res.x1 = u + v - aover3;
/*
res.x2.real(-(u + v) / 2.0 - aover3);
res.x2.imag(BUFFER_SQRT3 / 2.0 * (u - v));
res.x3.real(res.x2.real());
res.x3.imag(-res.x2.imag());
*/
res.type = CardanoResultCase::positived;
}
if(std::abs(d) <= 1e-20)
{
u = std::cbrt(-q / 2.0);
v = u;
res.x1 = u + v - aover3;
// res.x2.real(-(u + v) / 2.0 - aover3);
res.type = CardanoResultCase::zerod;
}
if(d < -1e-20)
{
r = std::cbrt(std::sqrt(-p * p * p * BUFFER_inv27));
alpha = std::atan(std::sqrt(-d) / -q * 2.0);
if(q > 0) // if q > 0 the angle becomes 2 * PI - alpha
alpha = 2.0 * M_PI - alpha;
res.x1 =
r * (std::cos((6.0 * M_PI - alpha) / 3.0) + std::cos(alpha / 3.0)) -
aover3;
/*
res.x2.real(r * (std::cos((2.0 * M_PI + alpha) / 3.0) +
std::cos((4.0 * M_PI - alpha) / 3.0)) -
aover3);
res.x3.real(r * (std::cos((4.0 * M_PI + alpha) / 3.0) +
std::cos((2.0 * M_PI - alpha) / 3.0)) -
aover3);
*/
res.type = CardanoResultCase::negatived;
}
}
InHousePolynomialSolverResult
inHousePolynomialSolve(const std::vector<double> &polynome)
{
InHousePolynomialSolverResult res;
res.type = CardanoResultCase::notvalid;
std::size_t polynome_size = polynome.size();
double a, b, c, delta;
switch(polynome_size)
{
case 0:
break;
case 1:
break;
case 2: // linear// equation ax + b = 0
a = polynome[1];
b = polynome[0];
if(a != 0)
{
res.x1 = -b / a;
res.type = CardanoResultCase::line;
}
break;
case 3:
// quadratic equation ax**2 + bx + c = 0
a = polynome[2];
if(a == 0)
return res;
b = polynome[1];
c = polynome[0];
if(a == 0)
return res;
// calculate the discriminant
delta = (b * b) - (a * c * 4);
// qDebug() << delta;
if(delta < 0)
{
}
else if(std::abs(delta) <= 1e-20)
{
res.x1 = -b / (a * 2);
// qDebug() << x1.real() << " " << x2.real();
res.type = CardanoResultCase::quadratic;
}
else
{
// find two solutions
delta = std::sqrt(delta);
res.x1 = (-b + delta) / (a * 2);
res.x2 = (-b - delta) / (a * 2);
// qDebug() << x1.real() << " " << x2.real();
res.type = CardanoResultCase::quadratic;
}
break;
case 4:
cubic_solver(res, polynome[3], polynome[2], polynome[1], polynome[0]);
}
return res;
}
BUFFER_SQRT3
const double BUFFER_SQRT3
Definition cardano.cpp:37
cubic_solver
void cubic_solver(InHousePolynomialSolverResult &res, double a1, double b, double c, double d)
Definition cardano.cpp:42
BUFFER_inv27
const double BUFFER_inv27
Definition cardano.cpp:38
BUFFER_pow11
const double BUFFER_pow11
Definition cardano.cpp:39
inHousePolynomialSolve
InHousePolynomialSolverResult inHousePolynomialSolve(const std::vector< double > &polynome)
Definition cardano.cpp:118
cardano.h
cubic solver adapted from https://www.codeproject.com/articles/798474/to-solve-a-cubic-equation thank...
CardanoResultCase::positived
@ positived
CardanoResultCase::line
@ line
CardanoResultCase::notvalid
@ notvalid
CardanoResultCase::zerod
@ zerod
CardanoResultCase::quadratic
@ quadratic
CardanoResultCase::negatived
@ negatived
InHousePolynomialSolverResult
Definition cardano.h:47