/*
  *  MoMEMta: a modular implementation of the Matrix Element Method
  *  Copyright (C) 2016  Universite catholique de Louvain (UCL), Belgium
  *
  *  This program 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.
  *
  *  This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
 #pragma once
 
 #include <cmath>
 #include <vector>
 
 #define SQ(x) ((x) * (x))
 #define CB(x) ((x) * (x) * (x))
 #define QU(x) ((x) * (x) * (x) * (x))
 
 template <typename T> T sign(const T& x) {
     if (x > 0)
         return 1;
     else if (!x)
         return 0;
     else
         return -1;
 }
 
 template <typename T> inline double dP_over_dE(const T& v) {
     const double rad = SQ(v.E()) - SQ(v.M());
     if (rad <= 0)
         return 0.;
     else
         return v.E() / std::sqrt(rad);
 }
 
 // If the NWA is used for a particle, a multiplication factor has to be introduced
 // because of the integrated-out delta function
 inline double jacobianNWA(const double mass, const double width) {
     return (M_PI / 2. + std::atan(mass / width)) * mass * width;
 }
 
 // Compute cos(x +- 2*pi/3) in a more "analytical" way (pm = +- 1)
 // Useful for solveCubic
 inline double cosXpm2PI3(const double x, const double pm) {
     return -0.5 * (std::cos(x) + pm * std::sin(x) * std::sqrt(3.));
 }
 
 bool solveQuadratic(const double a, const double b, const double c, std::vector<double>& roots,
                     bool verbose = false);
 
 bool solveCubic(const double a, const double b, const double c, const double d,
                 std::vector<double>& roots, bool verbose = false);
 
 // Finds the real solutions to a*x^4 + b*x^3 + c*x^2 + d*x + e = 0
 // Handles special case a=0.
 // Appends the solutions to the std::vector roots, making no attempt to check whether the vector is
 // empty.
 // Multiple roots appear multiple times.
 // If verbose is true (default is false), the solutions are printed,
 // as well as the polynomial evaluated on these solutions
 //
 // See https://en.wikipedia.org/wiki/Quartic_function#Solving_by_factoring_into_quadratics
 //    https://fr.wikipedia.org/wiki/Méthode_de_Descartes
 //
 // .
 
 bool solveQuartic(const double a, const double b, const double c, const double d, const double e,
                   std::vector<double>& roots, bool verbose = false);
 
 bool solve2Quads(const double a20, const double a02, const double a11, const double a10,
                  const double a01, const double a00, const double b20, const double b02,
                  const double b11, const double b10, const double b01, const double b00,
                  std::vector<double>& E1, std::vector<double>& E2, bool verbose = false);
 
 bool solve2QuadsDeg(const double a11, const double a10, const double a01, const double a00,
                     const double b11, const double b10, const double b01, const double b00,
                     std::vector<double>& E1, std::vector<double>& E2, bool verbose = false);
 
 bool solve2Linear(const double a10, const double a01, const double a00, const double b10,
                   const double b01, const double b00, std::vector<double>& E1,
                   std::vector<double>& E2, bool verbose = false);
 
 double BreitWigner(const double s, const double m, const double g);
 
 bool ApproxComparison(double value, double expected);