#pragma once

#ifndef TMATH_UTIL_H

#define TMATH_UTIL_H

//#include "tcommon.h"

#include "texception.h"

#include <numeric></numeric>

#include <complex></complex>

#ifndef __sgi

#include <cmath></cmath>

#else

#include <math.h></math.h>

#endif

#undef DVAPI

#undef DVVAR

#ifdef TNZCORE_EXPORTS

#define DVAPI DV_EXPORT_API

#define DVVAR DV_EXPORT_VAR

#else

#define DVAPI DV_IMPORT_API

#define DVVAR DV_IMPORT_VAR

#endif

enum TMathError { INFINITE_SOLUTIONS = -1 };

//-----------------------------------------------------------------------------

class DVAPI TMathException final : public TException {

  TString m_msg;

public:

  TMathException(std::string msg);

  virtual ~TMathException() {}

  TString getMessage() const override { return m_msg; }

};

//-----------------------------------------------------------------------------

/*!

  Decompose a square matrix [A] dim(A)=[n x n]  using the LU method.

  A is a matrix stored for row.

  \par A array of coeff

  \par n size of matrix

  \par indx array with permutation from pivoting

  \par d +1/-1 and depend whether the number of row interchanges (even or odd)

  \note If error a TException is thrown.

  \include mathutils_ex2.cpp

 */

DVAPI void tLUDecomposition(double *A, int n, int *indx, double &d);

//-----------------------------------------------------------------------------

/*!

  Use the back-substitution method to solve the linear

  system A*x=b.

  \par A matrix exited from tLUDecomposition

  \par n size of matrix

  \par indx array exited from tLUDecomposition

  \par b in input is array of element

  \ret b in output contains result

  \note A is the square matrix [n x n] returned from tLUDecomposition routine.

  \sa tLUDecomposition

 */

DVAPI void tbackSubstitution(double *A, int n, int *indx, double *b);

//-----------------------------------------------------------------------------

/*!

  Solve the system A*x=b.

  It`s needed to insert b in res.

  A and res are modifyed in the routine.

  \par A matrix

  \par n size of matrix

  \note A is square [n x n]

  \include mathutils_ex1.cpp

 */

DVAPI void tsolveSistem(double *A, int n, double *res);

//-----------------------------------------------------------------------------

/*!

  Solve the system A*x=b.

  It`s needed to insert b in res.

  A and res are modifyed in the routine.

  \note A is square [n x n]

  \include mathutils_ex1.cpp

 */

inline void tsolveSistem(std::vector<double> &A, std::vector<double> &res) {</double></double>

  assert(res.size() * res.size() == A.size());

  tsolveSistem(&A[0], res.size(), &res[0]);

}

//-----------------------------------------------------------------------------

/*!

  Find determinant of square matrix A

  \par A is matrix to test

  \par n size of matrix

 */

DVAPI double tdet(double *A, int n);

//-----------------------------------------------------------------------------

/*!

  Find determinant of square matrix A [n x n], using information

  from tLUDecomposition decomposition.

  \par LUa output matrix from tLUDecomposition

  \par n size of matrix

  \par d parameter from tLUDecomposition

  \sa tLUDecomposition

 */

DVAPI double tdet(double *LUa, int n, double d);

//-----------------------------------------------------------------------------

/*!

  Find real root of a polynomious

  \par poly is a vector with coeff of polynomious in crescent order.

  \ret sol is the vector of solution and contain real found solution.

 */

DVAPI int rootFinding(const std::vector<double> &poly,</double>

                      std::vector<double> &sol);</double>

//-----------------------------------------------------------------------------

/*!

  Find complex roots of a quadratic equation

    a*x^2 + b*x + c = 0

  \ret count of roots found:

    0 if a and b are almost zero

    1 if a is almost zero

    2 if a is not almost zero

    roots may have duplicates, count depends only from type of equation

 */

DVAPI int solveEquation2(std::complex<double> *roots, double a, double b, double c);</double>

//-----------------------------------------------------------------------------

/*!

  Find complex roots of a cubic equation:

    a*x^3 + b*x^2 + c*x + d = 0

  \ret count of roots found:

    0 if a, b and c are almost zero

    1 if a and b are almost zero

    2 if a is almost zero

    3 if a is not almost zero

    roots may have duplicates, count depends only from type of equation

 */

DVAPI int solveEquation3(std::complex<double> *roots, double a, double b, double c, double d);</double>

//-----------------------------------------------------------------------------

/*!

  Find complex roots of a power of four equation:

    a*x^4 + b*x^3 + c*x^2 + d*x + e = 0

  \ret count of roots found:

    0 if a, b, c and d are almost zero

    1 if a, b and c are almost zero

    2 if a and b are almost zero

    3 if a is almost zero

    4 if a is not almost zero

    roots may have duplicates, count depends only from type of equation

 */

DVAPI int solveEquation4(std::complex<double> *roots, double a, double b, double c, double d, double e);</double>

//-----------------------------------------------------------------------------

/*!

  Check if val is nearest to epsilon.

  \par val value to test

  \par eps tolerance required

  \ret true if value is nearest to zero

 */

inline bool isAlmostZero(double val, double eps = TConsts::epsilon) {

  return -eps < val && val < eps;

}

//-----------------------------------------------------------------------------

/*!

  Check if complex val is nearest to epsilon.

  \par val value to test

  \par eps tolerance required

  \ret true if value is nearest to zero

 */

inline bool isAlmostZero(std::complex<double> val, double eps = TConsts::epsilon) {</double>

  return val.real()*val.real() + val.imag()*val.imag() < eps*eps;

}

//-----------------------------------------------------------------------------

/*!

  Find number of roots of a poly in interval min max.

  \par order degree of polynomious

  \par coeff of poly in crescent order

  \par lesser extreme

  \par greater extreme

  \ret number of root

 */

DVAPI int numberOfRootsInInterval(int order, const double *polyH, double min,

                                  double max);

//-----------------------------------------------------------------------------

/*!

  Find the real root of the poly a*x^3+b*x^2+c*x+d in [0,1]

  \par a coeff of x^3

  \par b coeff of x^2

  \par c coeff of x^1

  \par d coeff of x^0

  \ret the real root of the poly, or 1.0 if the root is not in [0,1]

*/

DVAPI double cubicRoot(double a, double b, double c, double d);

//-----------------------------------------------------------------------------

/*!

  Find the root of the poly a*x^2+b*x+c

  \par a coeff of x^2

  \par b coeff of x^1

  \par c coeff of x^0

  \ret the root of the poly, or 1.0 if it does not exist

*/

DVAPI double quadraticRoot(double a, double b, double c);

//-----------------------------------------------------------------------------

#endif  // #ifndef TMATH_UTIL_H

//-----------------------------------------------------------------------------

// End Of File

//-----------------------------------------------------------------------------