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;
}