Shinya Kitaoka 810553
#pragma once
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
#ifndef TCG_POLY_OPS_H
Toshihiro Shimizu 890ddd
#define TCG_POLY_OPS_H
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
// tcg includes
Toshihiro Shimizu 890ddd
#include "macros.h"
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
/*!
Toshihiro Shimizu 890ddd
  \file     tcg_poly_ops.h
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
  \brief    This file contains useful functions for operating with polynomials.
Toshihiro Shimizu 890ddd
*/
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
//*******************************************************************************
Toshihiro Shimizu 890ddd
//    Polynomial Operations
Toshihiro Shimizu 890ddd
//*******************************************************************************
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
namespace tcg {
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
//! Contains useful functions for operating with polynomials.
Shinya Kitaoka 120a6e
namespace poly_ops {
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
/*!
Toshihiro Shimizu 890ddd
  \brief    Evaluates a polynomial using Horner's algorithm
Toshihiro Shimizu 890ddd
  \return   The value of the input polynomial at the specified parameter
Toshihiro Shimizu 890ddd
*/
Toshihiro Shimizu 890ddd
template <typename scalar=""></typename>
Shinya Kitaoka 120a6e
Scalar evaluate(const Scalar poly[],  //!< Coefficients of the input polynomial,
Shinya Kitaoka 38fd86
                                      //! indexed by degree
Shinya Kitaoka 120a6e
                int deg,              //!< Degree of the polynomial function
Shinya Kitaoka 120a6e
                Scalar x)  //!< Parameter the polynomial will be evaluated on
Toshihiro Shimizu 890ddd
{
Shinya Kitaoka 120a6e
  // ((poly[deg] * x) + poly[deg-1]) * x + poly[deg - 2] + ...
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  Scalar value = poly[deg];
Shinya Kitaoka 120a6e
  for (int d = deg - 1; d >= 0; --d) value = value * x + poly[d];
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  return value;
Toshihiro Shimizu 890ddd
}
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
//-------------------------------------------------------------------------------------------
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
/*!
Shinya Kitaoka 120a6e
  \brief  Reduces the degree of the input polynomial, discarding all leading
Shinya Kitaoka 120a6e
  coefficients
Toshihiro Shimizu 890ddd
          whose absolute value is below the specified tolerance threshold.
Toshihiro Shimizu 890ddd
*/
Toshihiro Shimizu 890ddd
template <typename scalar=""></typename>
Toshihiro Shimizu 890ddd
void reduceDegree(
Shinya Kitaoka 120a6e
    const Scalar poly[],  //!< Input polynomial to be reduced.
Shinya Kitaoka 120a6e
    int °,             //!< Input/output polynomial degree.
Shinya Kitaoka 120a6e
    Scalar tolerance  //!< Coefficients threshold to reduce the polynomial with.
Shinya Kitaoka 120a6e
    ) {
Shinya Kitaoka 120a6e
  while ((deg > 0) && (std::abs(poly[deg]) < tolerance)) --deg;
Toshihiro Shimizu 890ddd
}
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
//-------------------------------------------------------------------------------------------
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
/*!
Toshihiro Shimizu 890ddd
  \brief    Adds two polynomials and returns the sum.
Toshihiro Shimizu 890ddd
  \remark   The supplied polynomials can actually be the same.
Toshihiro Shimizu 890ddd
*/
Toshihiro Shimizu 890ddd
template <typename a,="" b,="" c,="" deg="" int="" typename=""></typename>
Shinya Kitaoka 120a6e
void add(const A (&poly1)[deg],  //!< First polynomial addendum.
Shinya Kitaoka 120a6e
         const B (&poly2)[deg],  //!< Second polynomial addendum.
Shinya Kitaoka 120a6e
         C (&result)[deg])       //!< Resulting sum.
Toshihiro Shimizu 890ddd
{
Shinya Kitaoka 120a6e
  for (int d = 0; d != deg; ++d) result[d] = poly1[d] + poly2[d];
Toshihiro Shimizu 890ddd
}
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
//-------------------------------------------------------------------------------------------
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
/*!
Toshihiro Shimizu 890ddd
  \brief  Subtracts two polynomials /p poly1 and \p poly2 and returns the
Toshihiro Shimizu 890ddd
          difference <tt>poly1 - poly2</tt>.
Toshihiro Shimizu 890ddd
*/
Toshihiro Shimizu 890ddd
template <typename a,="" b,="" c,="" deg="" int="" typename=""></typename>
Shinya Kitaoka 120a6e
void sub(const A (&poly1)[deg],  //!< First polynomial addendum.
Shinya Kitaoka 120a6e
         const B (&poly2)[deg],  //!< Second polynomial addendum.
Shinya Kitaoka 120a6e
         C (&result)[deg])       //!< Resulting difference.
Toshihiro Shimizu 890ddd
{
Shinya Kitaoka 120a6e
  for (int d = 0; d != deg; ++d) result[d] = poly1[d] - poly2[d];
Toshihiro Shimizu 890ddd
}
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
//-------------------------------------------------------------------------------------------
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
/*!
Toshihiro Shimizu 890ddd
  \brief    Multiplies two polynomials into a polynomial whose degree is the
Toshihiro Shimizu 890ddd
            \a sum of the multiplicands' degrees.
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
  \warning  The resulting polynomial is currently not allowed to be one
Toshihiro Shimizu 890ddd
            of the multiplicands.
Toshihiro Shimizu 890ddd
*/
Toshihiro Shimizu 890ddd
template <typename a,="" b,="" c,="" deg1,="" deg2,="" degr="" int="" typename=""></typename>
Shinya Kitaoka 120a6e
void mul(const A (&poly1)[deg1],  //!< First multiplicand.
Shinya Kitaoka 120a6e
         const B (&poly2)[deg2],  //!< Second multiplicand.
Shinya Kitaoka 120a6e
         C (&result)[degR])       //!< Resulting polynomial.
Toshihiro Shimizu 890ddd
{
Shinya Kitaoka 120a6e
  TCG_STATIC_ASSERT(degR == deg1 + deg2 - 1);
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  for (int a = 0; a != deg1; ++a) {
Shinya Kitaoka 120a6e
    for (int b = 0; b != deg2; ++b) result[a + b] += poly1[a] * poly2[b];
Shinya Kitaoka 120a6e
  }
Toshihiro Shimizu 890ddd
}
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
//-------------------------------------------------------------------------------------------
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
/*!
Shinya Kitaoka 120a6e
  \brief    Calculates the chaining <tt>poly1 o poly2</tt> of two given
Shinya Kitaoka 120a6e
  polynomial
Toshihiro Shimizu 890ddd
            \p poly1 and \p poly2.
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
  \warning  The resulting polynomial is currently not allowed to be one
Toshihiro Shimizu 890ddd
            of the multiplicands.
Toshihiro Shimizu 890ddd
  \warning  This function is still \b untested.
Toshihiro Shimizu 890ddd
*/
Toshihiro Shimizu 890ddd
template <typename a,="" b,="" c,="" deg1,="" deg2,="" degr="" int="" typename=""></typename>
Shinya Kitaoka 120a6e
void chain(const A (&poly1)[deg1],  //!< First polynomial.
Shinya Kitaoka 120a6e
           const B (&poly2)[deg2],  //!< Second polynomial.
Shinya Kitaoka 120a6e
           C (&result)[degR])       //!< Resulting polynomial.
Toshihiro Shimizu 890ddd
{
Shinya Kitaoka 120a6e
  for (int a = 0; a != deg1; ++a) {
Shinya Kitaoka 120a6e
    B pow[degR][2] = {{}};
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
    // Build poly2 powers
Shinya Kitaoka 120a6e
    {
Shinya Kitaoka 120a6e
      std::copy(poly2, poly2 + deg2, pow[1]);
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
      for (int p = 1; p < a; ++p) poly_mul(pow[p % 2], poly2, pow[(p + 1) % 2]);
Shinya Kitaoka 120a6e
    }
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
    B(&pow_add)[degR] = pow[a % 2];
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
    // multiply by poly1[a]
Shinya Kitaoka 120a6e
    C addendum[degR];
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
    for (int c = 0; c != degR; ++c) addendum[c] = poly1[c] * pow_add[c];
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
    poly_add(addendum, result);
Shinya Kitaoka 120a6e
  }
Toshihiro Shimizu 890ddd
}
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
//-------------------------------------------------------------------------------------------
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
template <typename a,="" b,="" deg,="" degr="" int="" typename=""></typename>
Shinya Kitaoka 120a6e
void derivative(const A (&poly)[deg],  //!< Polynomial to be derived.
Shinya Kitaoka 120a6e
                B (&result)[degR])     //!< Resulting derivative polynomial.
Toshihiro Shimizu 890ddd
{
Shinya Kitaoka 120a6e
  TCG_STATIC_ASSERT(degR == deg - 1);
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  for (int c = 1; c != deg; ++c) result[c - 1] = c * poly[c];
Toshihiro Shimizu 890ddd
}
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
//-------------------------------------------------------------------------------------------
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
/*!
Toshihiro Shimizu 890ddd
  \brief    Solves the 1st degree equation: $c[1] t + c[0] = 0$
Shinya Kitaoka 120a6e
  \return   The number of solutions found under the specified divide-by
Shinya Kitaoka 120a6e
  tolerance
Toshihiro Shimizu 890ddd
*/
Toshihiro Shimizu 890ddd
template <typename scalar=""></typename>
Toshihiro Shimizu 890ddd
inline unsigned int solve_1(
Shinya Kitaoka 120a6e
    Scalar c[2],  //!< Polynomial coefficients array
Shinya Kitaoka 120a6e
    Scalar s[1],  //!< Solutions array
Shinya Kitaoka 120a6e
    Scalar tol =
Shinya Kitaoka 120a6e
        0)  //!< Leading coefficient tolerance, the equation has no solution
Shinya Kitaoka 120a6e
            //!  if the leading coefficient is below this threshold
Toshihiro Shimizu 890ddd
{
Shinya Kitaoka 120a6e
  if (std::abs(c[1]) <= tol) return 0;
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  s[0] = -c[0] / c[1];
Shinya Kitaoka 120a6e
  return 1;
Toshihiro Shimizu 890ddd
}
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
//-------------------------------------------------------------------------------------------
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
/*!
Toshihiro Shimizu 890ddd
  \brief    Solves the 2nd degree equation: $c[2] t^2 + c[1] t + c[0] = 0$
Shinya Kitaoka 120a6e
  \return   The number of #real# solutions found under the specified divide-by
Shinya Kitaoka 120a6e
  tolerance
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
  \remark   The returned solutions are sorted, with $s[0] <= s[1]$
Toshihiro Shimizu 890ddd
*/
Toshihiro Shimizu 890ddd
template <typename scalar=""></typename>
Toshihiro Shimizu 890ddd
unsigned int solve_2(
Shinya Kitaoka 120a6e
    Scalar c[3],  //!< Polynomial coefficients array
Shinya Kitaoka 120a6e
    Scalar s[2],  //!< Solutions array
Shinya Kitaoka 120a6e
    Scalar tol =
Shinya Kitaoka 120a6e
        0)  //!< Leading coefficient tolerance, the equation has no solution
Shinya Kitaoka 120a6e
            //!  if the leading coefficient is below this threshold
Toshihiro Shimizu 890ddd
{
Shinya Kitaoka 120a6e
  if (std::abs(c[2]) <= tol)
Shinya Kitaoka 120a6e
    return solve_1(c, s, tol);  // Reduces to first degree
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  Scalar nc[2] = {
Shinya Kitaoka 120a6e
      c[0] / c[2],
Shinya Kitaoka 120a6e
      c[1] / (c[2] + c[2])};  // NOTE: nc[1] gets further divided by 2
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  Scalar delta = nc[1] * nc[1] - nc[0];
Shinya Kitaoka 120a6e
  if (delta < 0) return 0;
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  delta = sqrt(delta);
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  s[0] = -delta - nc[1];
Shinya Kitaoka 120a6e
  s[1] = delta - nc[1];
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  return 2;
Toshihiro Shimizu 890ddd
}
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
//-------------------------------------------------------------------------------------------
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
/*!
Shinya Kitaoka 120a6e
  \brief    Solves the 3rd degree equation: $c[3]t^3 + c[2] t^2 + c[1] t + c[0]
Shinya Kitaoka 120a6e
  = 0$
Shinya Kitaoka 120a6e
  \return   The number of #real# solutions found under the specified divide-by
Shinya Kitaoka 120a6e
  tolerance
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
  \remark   The returned solutions are sorted, with $s[0] <= s[1] <= s[2]$
Toshihiro Shimizu 890ddd
*/
Toshihiro Shimizu 890ddd
template <typename scalar=""></typename>
Shinya Kitaoka 120a6e
unsigned int solve_3(Scalar c[4],     //!< Polynomial coefficients array
Shinya Kitaoka 120a6e
                     Scalar s[3],     //!< Solutions array
Shinya Kitaoka 120a6e
                     Scalar tol = 0)  //!< Leading coefficient tolerance, the
Shinya Kitaoka 38fd86
                                      //! equation is reduced to 2nd degree
Shinya Kitaoka 120a6e
//!  if the leading coefficient is below this threshold
Toshihiro Shimizu 890ddd
{
Shinya Kitaoka 120a6e
  if (std::abs(c[3]) <= tol)
Shinya Kitaoka 120a6e
    return solve_2(c, s, tol);  // Reduces to second degree
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  Scalar nc[3] = {c[0] / c[3], c[1] / c[3], c[2] / c[3]};
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  Scalar b2     = nc[2] * nc[2];
Shinya Kitaoka 120a6e
  Scalar p      = nc[1] - b2 / 3;
Shinya Kitaoka 120a6e
  Scalar q      = nc[2] * (b2 + b2 - 9 * nc[1]) / 27 + nc[0];
Shinya Kitaoka 120a6e
  Scalar p3     = p * p * p;
Shinya Kitaoka 120a6e
  Scalar d      = q * q + 4 * p3 / 27;
Shinya Kitaoka 120a6e
  Scalar offset = -nc[2] / 3;
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  if (d >= 0) {
Shinya Kitaoka 120a6e
    // Single solution
Shinya Kitaoka 120a6e
    Scalar z = sqrt(d);
Shinya Kitaoka 120a6e
    Scalar u = (-q + z) / 2;
Shinya Kitaoka 120a6e
    Scalar v = (-q - z) / 2;
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
    u = (u < 0) ? -pow(-u, 1 / Scalar(3)) : pow(u, 1 / Scalar(3));
Shinya Kitaoka 120a6e
    v = (v < 0) ? -pow(-v, 1 / Scalar(3)) : pow(v, 1 / Scalar(3));
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
    s[0] = offset + u + v;
Shinya Kitaoka 120a6e
    return 1;
Shinya Kitaoka 120a6e
  }
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  assert(p3 < 0);
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  Scalar u = sqrt(-p / Scalar(3));
Shinya Kitaoka 120a6e
  Scalar v = acos(-sqrt(-27 / p3) * q / Scalar(2)) / Scalar(3);
Shinya Kitaoka 120a6e
  Scalar m = cos(v), n = sin(v) * 1.7320508075688772935274463415059;  // sqrt(3)
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  s[0] = offset - u * (n + m);
Shinya Kitaoka 120a6e
  s[1] = offset + u * (n - m);
Shinya Kitaoka 120a6e
  s[2] = offset + u * (m + m);
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  assert(s[0] <= s[1] && s[1] <= s[2]);
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
  return 3;
Toshihiro Shimizu 890ddd
}
Toshihiro Shimizu 890ddd
}
Shinya Kitaoka 120a6e
}  // namespace tcg::poly_ops
Toshihiro Shimizu 890ddd
Shinya Kitaoka 120a6e
#endif  // TCG_POLY_OPS_H