PolyBase.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include <algorithm>
#include <cmath>
#include <complex>
#include <functional>
#include <gsl/gsl_poly.h>
#include <iterator>
#include <vector>
 
#include "MantidGeometry/Math/PolyBase.h"
#include "MantidKernel/Exception.h"
 
namespace Mantid::mathLevel {
 
std::ostream &operator<<(std::ostream &OX, const PolyBase &A)
{
  A.write(OX);
  return OX;
}
 
PolyBase::PolyBase(const int iD)
    : iDegree((iD > 0) ? iD : 0), afCoeff(iDegree + 1), Eaccuracy(1e-6)
{}
 
PolyBase::PolyBase(const int iD, const double E)
    : iDegree((iD > 0) ? iD : 0), afCoeff(iDegree + 1), Eaccuracy(fabs(E))
{}
 
void PolyBase::setDegree(const int iD)
{
  iDegree = (iD > 0) ? iD : 0;
  afCoeff.resize(iDegree + 1);
}
 
int PolyBase::getDegree() const
{
  return iDegree;
}
 
PolyBase::operator const std::vector<double> &() const
{
  return afCoeff;
}
 
PolyBase::operator std::vector<double> &()
{
  return afCoeff;
}
 
double PolyBase::operator[](const int i) const
{
  if (i > iDegree || i < 0)
    throw Kernel::Exception::IndexError(i, iDegree + 1, "PolyBase::operator[] const");
  return afCoeff[i];
}
 
double &PolyBase::operator[](const int i)
{
  if (i > iDegree || i < 0)
    throw Kernel::Exception::IndexError(i, iDegree + 1, "PolyBase::operator[]");
  return afCoeff[i];
}
 
double PolyBase::operator()(const double X) const
{
  double Result = afCoeff[iDegree];
  for (int i = iDegree - 1; i >= 0; i--) {
    Result *= X;
    Result += afCoeff[i];
  }
  return Result;
}
 
PolyBase &PolyBase::operator+=(const PolyBase &A)
{
  iDegree = (iDegree > A.iDegree) ? iDegree : A.iDegree;
  afCoeff.resize(iDegree + 1);
  for (int i = 0; i <= iDegree && i <= A.iDegree; i++)
    afCoeff[i] += A.afCoeff[i];
  return *this;
}
 
PolyBase &PolyBase::operator-=(const PolyBase &A)
{
  iDegree = (iDegree > A.iDegree) ? iDegree : A.iDegree;
  afCoeff.resize(iDegree + 1);
  for (int i = 0; i <= iDegree && i <= A.iDegree; i++)
    afCoeff[i] -= A.afCoeff[i];
  return *this;
}
 
PolyBase &PolyBase::operator*=(const PolyBase &A)
{
  const int iD = iDegree + A.iDegree;
  std::vector<double> CX(iD + 1, 0.0); // all set to zero
  for (int i = 0; i <= iDegree; i++)
    for (int j = 0; j <= A.iDegree; j++) {
      const int cIndex = i + j;
      CX[cIndex] += afCoeff[i] * A.afCoeff[j];
    }
  afCoeff = CX;
  return *this;
}
 
PolyBase PolyBase::operator+(const PolyBase &A) const
{
  PolyBase kSum(*this);
  return kSum += A;
}
 
PolyBase PolyBase::operator-(const PolyBase &A) const
{
  PolyBase kSum(*this);
  return kSum -= A;
}
 
PolyBase PolyBase::operator*(const PolyBase &A) const
{
  PolyBase kSum(*this);
  return kSum *= A;
}
 
PolyBase PolyBase::operator+(const double V) const
{
  PolyBase kSum(*this);
  return kSum += V;
}
 
PolyBase PolyBase::operator-(const double V) const
{
  PolyBase kSum(*this);
  return kSum -= V;
}
 
PolyBase PolyBase::operator*(const double V) const
{
  PolyBase kSum(*this);
  return kSum *= V;
}
 
PolyBase PolyBase::operator/(const double V) const
{
  PolyBase kSum(*this);
  return kSum /= V;
}
 
// ---------------- SCALARS --------------------------
 
PolyBase &PolyBase::operator+=(const double V)
{
  afCoeff[0] += V; // There is always zero component
  return *this;
}
 
PolyBase &PolyBase::operator-=(const double V)
{
  afCoeff[0] -= V; // There is always zero component
  return *this;
}
 
PolyBase &PolyBase::operator*=(const double V)
{
  using std::placeholders::_1;
  transform(afCoeff.begin(), afCoeff.end(), afCoeff.begin(), std::bind(std::multiplies<double>(), _1, V));
  return *this;
}
 
PolyBase &PolyBase::operator/=(const double V)
{
  using std::placeholders::_1;
  transform(afCoeff.begin(), afCoeff.end(), afCoeff.begin(), std::bind(std::divides<double>(), _1, V));
  return *this;
}
 
PolyBase PolyBase::operator-() const
{
  PolyBase KOut(*this);
  KOut *= -1.0;
  return KOut;
}
 
PolyBase PolyBase::getDerivative() const
{
  PolyBase KOut(*this);
  return KOut.derivative();
}
 
PolyBase &PolyBase::derivative()
{
  if (iDegree < 1) {
    afCoeff[0] = 0.0;
    return *this;
  }
 
  for (int i = 0; i < iDegree; i++)
    afCoeff[i] = afCoeff[i + 1] * (i + 1);
  iDegree--;
 
  return *this;
}
 
PolyBase PolyBase::GetInversion() const
{
  PolyBase InvPoly(iDegree);
  for (int i = 0; i <= iDegree; i++)
    InvPoly.afCoeff[i] = afCoeff[iDegree - i];
  return InvPoly;
}
 
void PolyBase::compress(const double epsilon)
{
  const double eps((epsilon > 0.0) ? epsilon : Eaccuracy);
  for (; iDegree >= 0 && fabs(afCoeff[iDegree]) <= eps; iDegree--)
    ;
 
  if (iDegree >= 0) {
    const double Leading = afCoeff[iDegree];
    afCoeff[iDegree] = 1.0; // avoid rounding issues
    for (int i = 0; i < iDegree; i++)
      afCoeff[i] /= Leading;
  } else
    iDegree = 0;
 
  afCoeff.resize(iDegree + 1);
}
 
void PolyBase::divide(const PolyBase &pD, PolyBase &pQ, PolyBase &pR, const double epsilon) const
{
  const int iQuotDegree = iDegree - pD.iDegree;
  if (iQuotDegree >= 0) {
    pQ.setDegree(iQuotDegree);
    // temporary storage for the remainder
    pR = *this;
 
    // do the division
    double fInv = 1.0 / pD[pD.iDegree];
    for (int iQ = iQuotDegree; iQ >= 0; iQ--) {
      int iR = pD.iDegree + iQ;
      pQ[iQ] = fInv * pR[iR];
      for (iR--; iR >= iQ; iR--)
        pR.afCoeff[iR] -= pQ[iQ] * pD[iR - iQ];
    }
 
    // calculate the correct degree for the remainder
    pR.compress(epsilon);
    return;
  }
 
  pQ.setDegree(0);
  pQ[0] = 0.0;
  pR = *this;
}
 
std::vector<double> PolyBase::realRoots(const double epsilon)
{
  const double eps((epsilon > 0.0) ? epsilon : Eaccuracy);
  std::vector<std::complex<double>> Croots = calcRoots(epsilon);
  std::vector<double> Out;
  std::vector<std::complex<double>>::const_iterator vc;
  for (vc = Croots.begin(); vc != Croots.end(); ++vc) {
    if (fabs(vc->imag()) < eps)
      Out.emplace_back(vc->real());
  }
  return Out;
}
 
std::vector<std::complex<double>> PolyBase::calcRoots(const double epsilon)
{
  compress(epsilon);
  std::vector<std::complex<double>> Out(iDegree);
  // Zero State:
  if (iDegree == 0)
    return Out;
 
  // x+a_0 =0
  if (iDegree == 1) {
    Out[0] = std::complex<double>(-afCoeff[0]);
    return Out;
  }
  // x^2+a_1 x+c = 0
  if (iDegree == 2) {
    solveQuadratic(Out[0], Out[1]);
    return Out;
  }
 
  // x^3+a_2 x^2+ a_1 x+c=0
  if (iDegree == 3) {
    solveCubic(Out[0], Out[1], Out[2]);
    return Out;
  }
  // THERE IS A QUARTIC / QUINTIC Solution availiable but...
  // WS contains the the hessian matrix if required (eigenvalues/vectors)
  //
  gsl_poly_complex_workspace *WS(gsl_poly_complex_workspace_alloc(iDegree + 1));
  auto RZ = new double[2 * (iDegree + 1)];
  gsl_poly_complex_solve(&afCoeff.front(), iDegree + 1, WS, RZ);
  for (int i = 0; i < iDegree; i++)
    Out[i] = std::complex<double>(RZ[2 * i], RZ[2 * i + 1]);
  gsl_poly_complex_workspace_free(WS);
  delete[] RZ;
  return Out;
}
 
int PolyBase::solveQuadratic(std::complex<double> &AnsA, std::complex<double> &AnsB) const
{
  const double b = afCoeff[1];
  const double c = afCoeff[0];
 
  const double cf = b * b - 4.0 * c;
  if (cf >= 0) /* Real Roots */
  {
    const double q = (b >= 0) ? -0.5 * (b + sqrt(cf)) : -0.5 * (b - sqrt(cf));
    AnsA = std::complex<double>(q, 0.0);
    AnsB = std::complex<double>(c / q, 0.0);
    return (cf == 0) ? 1 : 2;
  }
 
  std::complex<double> CQ(-0.5 * b, (b >= 0 ? -0.5 * sqrt(-cf) : 0.5 * sqrt(-cf)));
  AnsA = CQ;
  AnsB = c / CQ;
  return 2;
}
 
int PolyBase::solveCubic(std::complex<double> &AnsA, std::complex<double> &AnsB, std::complex<double> &AnsC) const
{
  double q, r; /* solution parameters */
  double termR, discrim;
  double r13;
  // std::pair<std::complex<double>,std::complex<double> > SQ;
 
  const double b = afCoeff[2];
  const double c = afCoeff[1];
  const double d = afCoeff[0];
 
  q = (3.0 * c - (b * b)) / 9.0;               // -q
  r = -27.0 * d + b * (9.0 * c - 2.0 * b * b); // -r
  r /= 54.0;
 
  discrim = q * q * q + r * r; // r^2-qq^3
  /* The first root is always real. */
  termR = (b / 3.0);
 
  if (discrim > 1e-13) /* one root real, two are complex */
  {
    double s, t, termI;
    s = r + sqrt(discrim);
    s = ((s < 0) ? -pow(-s, (1.0 / 3.0)) : pow(s, (1.0 / 3.0)));
    t = r - sqrt(discrim);
    t = ((t < 0) ? -pow(-t, (1.0 / 3.0)) : pow(t, (1.0 / 3.0)));
    AnsA = std::complex<double>(-termR + s + t, 0.0);
    // second real point.
    termR += (s + t) / 2.0;
    termI = sqrt(3.0) * (-t + s) / 2;
    AnsB = std::complex<double>(-termR, termI);
    AnsC = std::complex<double>(-termR, -termI);
    return 3;
  }
 
  /* The remaining options are all real */
 
  if (discrim < 1e-13) // All roots real
  {
    double q3;
    q = -q;
    q3 = q * q * q;
    q3 = acos(-r / sqrt(q3));
    r13 = -2.0 * sqrt(q);
    AnsA = std::complex<double>(-termR + r13 * cos(q3 / 3.0), 0.0);
    AnsB = std::complex<double>(-termR + r13 * cos((q3 + 2.0 * M_PI) / 3.0), 0.0);
    AnsC = std::complex<double>(-termR + r13 * cos((q3 - 2.0 * M_PI) / 3.0), 0.0);
    return 3;
  }
 
  // Only option left is that all roots are real and unequal
  // (to get here, q*q*q=r*r)
 
  r13 = ((r < 0) ? -pow(-r, (1.0 / 3.0)) : pow(r, (1.0 / 3.0)));
  AnsA = std::complex<double>(-termR + 2.0 * r13, 0.0);
  AnsB = std::complex<double>(-(r13 + termR), 0.0);
  AnsC = std::complex<double>(-(r13 + termR), 0.0);
  return 2;
}
 
void PolyBase::write(std::ostream &OX) const
{
  copy(afCoeff.begin(), afCoeff.end(), std::ostream_iterator<double>(OX, " "));
}
 
} // namespace Mantid::mathLevel