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
Toshihiro Shimizu 890ddd
namespace tcg
Toshihiro Shimizu 890ddd
{
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
//! Contains useful functions for operating with polynomials.
Toshihiro Shimizu 890ddd
namespace poly_ops
Toshihiro Shimizu 890ddd
{
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>
Toshihiro Shimizu 890ddd
Scalar evaluate(
Toshihiro Shimizu 890ddd
	const Scalar poly[], //!< Coefficients of the input polynomial, indexed by degree
Toshihiro Shimizu 890ddd
	int deg,			 //!< Degree of the polynomial function
Toshihiro Shimizu 890ddd
	Scalar x)			 //!< Parameter the polynomial will be evaluated on
Toshihiro Shimizu 890ddd
{
Toshihiro Shimizu 890ddd
	// ((poly[deg] * x) + poly[deg-1]) * x + poly[deg - 2] + ...
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	Scalar value = poly[deg];
Toshihiro Shimizu 890ddd
	for (int d = deg - 1; d >= 0; --d)
Toshihiro Shimizu 890ddd
		value = value * x + poly[d];
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	return value;
Toshihiro Shimizu 890ddd
}
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
//-------------------------------------------------------------------------------------------
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
/*!
Toshihiro Shimizu 890ddd
  \brief  Reduces the degree of the input polynomial, discarding all leading 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(
Toshihiro Shimizu 890ddd
	const Scalar poly[], //!< Input polynomial to be reduced.
Toshihiro Shimizu 890ddd
	int °,			 //!< Input/output polynomial degree.
Toshihiro Shimizu 890ddd
	Scalar tolerance	 //!< Coefficients threshold to reduce the polynomial with.
Toshihiro Shimizu 890ddd
	)
Toshihiro Shimizu 890ddd
{
Toshihiro Shimizu 890ddd
	while ((deg > 0) && (std::abs(poly[deg]) < tolerance))
Toshihiro Shimizu 890ddd
		--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>
Toshihiro Shimizu 890ddd
void add(
Toshihiro Shimizu 890ddd
	const A(&poly1)[deg], //!< First polynomial addendum.
Toshihiro Shimizu 890ddd
	const B(&poly2)[deg], //!< Second polynomial addendum.
Toshihiro Shimizu 890ddd
	C(&result)[deg])	  //!< Resulting sum.
Toshihiro Shimizu 890ddd
{
Toshihiro Shimizu 890ddd
	for (int d = 0; d != deg; ++d)
Toshihiro Shimizu 890ddd
		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>
Toshihiro Shimizu 890ddd
void sub(
Toshihiro Shimizu 890ddd
	const A(&poly1)[deg], //!< First polynomial addendum.
Toshihiro Shimizu 890ddd
	const B(&poly2)[deg], //!< Second polynomial addendum.
Toshihiro Shimizu 890ddd
	C(&result)[deg])	  //!< Resulting difference.
Toshihiro Shimizu 890ddd
{
Toshihiro Shimizu 890ddd
	for (int d = 0; d != deg; ++d)
Toshihiro Shimizu 890ddd
		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>
Toshihiro Shimizu 890ddd
void mul(
Toshihiro Shimizu 890ddd
	const A(&poly1)[deg1], //!< First multiplicand.
Toshihiro Shimizu 890ddd
	const B(&poly2)[deg2], //!< Second multiplicand.
Toshihiro Shimizu 890ddd
	C(&result)[degR])	  //!< Resulting polynomial.
Toshihiro Shimizu 890ddd
{
Toshihiro Shimizu 890ddd
	TCG_STATIC_ASSERT(degR == deg1 + deg2 - 1);
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	for (int a = 0; a != deg1; ++a) {
Toshihiro Shimizu 890ddd
		for (int b = 0; b != deg2; ++b)
Toshihiro Shimizu 890ddd
			result[a + b] += poly1[a] * poly2[b];
Toshihiro Shimizu 890ddd
	}
Toshihiro Shimizu 890ddd
}
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
//-------------------------------------------------------------------------------------------
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
/*!
Toshihiro Shimizu 890ddd
  \brief    Calculates the chaining <tt>poly1 o poly2</tt> of two given 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>
Toshihiro Shimizu 890ddd
void chain(
Toshihiro Shimizu 890ddd
	const A(&poly1)[deg1], //!< First polynomial.
Toshihiro Shimizu 890ddd
	const B(&poly2)[deg2], //!< Second polynomial.
Toshihiro Shimizu 890ddd
	C(&result)[degR])	  //!< Resulting polynomial.
Toshihiro Shimizu 890ddd
{
Toshihiro Shimizu 890ddd
	for (int a = 0; a != deg1; ++a) {
Toshihiro Shimizu 890ddd
		B pow[degR][2] = {{}};
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
		// Build poly2 powers
Toshihiro Shimizu 890ddd
		{
Toshihiro Shimizu 890ddd
			std::copy(poly2, poly2 + deg2, pow[1]);
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
			for (int p = 1; p < a; ++p)
Toshihiro Shimizu 890ddd
				poly_mul(pow[p % 2], poly2, pow[(p + 1) % 2]);
Toshihiro Shimizu 890ddd
		}
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
		B(&pow_add)
Toshihiro Shimizu 890ddd
		[degR] = pow[a % 2];
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
		// multiply by poly1[a]
Toshihiro Shimizu 890ddd
		C addendum[degR];
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
		for (int c = 0; c != degR; ++c)
Toshihiro Shimizu 890ddd
			addendum[c] = poly1[c] * pow_add[c];
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
		poly_add(addendum, result);
Toshihiro Shimizu 890ddd
	}
Toshihiro Shimizu 890ddd
}
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
//-------------------------------------------------------------------------------------------
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
template <typename a,="" b,="" deg,="" degr="" int="" typename=""></typename>
Toshihiro Shimizu 890ddd
void derivative(
Toshihiro Shimizu 890ddd
	const A(&poly)[deg], //!< Polynomial to be derived.
Toshihiro Shimizu 890ddd
	B(&result)[degR])	//!< Resulting derivative polynomial.
Toshihiro Shimizu 890ddd
{
Toshihiro Shimizu 890ddd
	TCG_STATIC_ASSERT(degR == deg - 1);
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	for (int c = 1; c != deg; ++c)
Toshihiro Shimizu 890ddd
		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$
Toshihiro Shimizu 890ddd
  \return   The number of solutions found under the specified divide-by tolerance
Toshihiro Shimizu 890ddd
*/
Toshihiro Shimizu 890ddd
template <typename scalar=""></typename>
Toshihiro Shimizu 890ddd
inline unsigned int solve_1(
Toshihiro Shimizu 890ddd
	Scalar c[2],	//!< Polynomial coefficients array
Toshihiro Shimizu 890ddd
	Scalar s[1],	//!< Solutions array
Toshihiro Shimizu 890ddd
	Scalar tol = 0) //!< Leading coefficient tolerance, the equation has no solution
Toshihiro Shimizu 890ddd
					//!  if the leading coefficient is below this threshold
Toshihiro Shimizu 890ddd
{
Toshihiro Shimizu 890ddd
	if (std::abs(c[1]) <= tol)
Toshihiro Shimizu 890ddd
		return 0;
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	s[0] = -c[0] / c[1];
Toshihiro Shimizu 890ddd
	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$
Toshihiro Shimizu 890ddd
  \return   The number of #real# solutions found under the specified divide-by 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(
Toshihiro Shimizu 890ddd
	Scalar c[3],	//!< Polynomial coefficients array
Toshihiro Shimizu 890ddd
	Scalar s[2],	//!< Solutions array
Toshihiro Shimizu 890ddd
	Scalar tol = 0) //!< Leading coefficient tolerance, the equation has no solution
Toshihiro Shimizu 890ddd
					//!  if the leading coefficient is below this threshold
Toshihiro Shimizu 890ddd
{
Toshihiro Shimizu 890ddd
	if (std::abs(c[2]) <= tol)
Toshihiro Shimizu 890ddd
		return solve_1(c, s, tol); // Reduces to first degree
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	Scalar nc[2] = {c[0] / c[2], c[1] / (c[2] + c[2])}; // NOTE: nc[1] gets further divided by 2
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	Scalar delta = nc[1] * nc[1] - nc[0];
Toshihiro Shimizu 890ddd
	if (delta < 0)
Toshihiro Shimizu 890ddd
		return 0;
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	delta = sqrt(delta);
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	s[0] = -delta - nc[1];
Toshihiro Shimizu 890ddd
	s[1] = delta - nc[1];
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	return 2;
Toshihiro Shimizu 890ddd
}
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
//-------------------------------------------------------------------------------------------
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
/*!
Toshihiro Shimizu 890ddd
  \brief    Solves the 3rd degree equation: $c[3]t^3 + c[2] t^2 + c[1] t + c[0] = 0$
Toshihiro Shimizu 890ddd
  \return   The number of #real# solutions found under the specified divide-by 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>
Toshihiro Shimizu 890ddd
unsigned int solve_3(
Toshihiro Shimizu 890ddd
	Scalar c[4],	//!< Polynomial coefficients array
Toshihiro Shimizu 890ddd
	Scalar s[3],	//!< Solutions array
Toshihiro Shimizu 890ddd
	Scalar tol = 0) //!< Leading coefficient tolerance, the equation is reduced to 2nd degree
Toshihiro Shimizu 890ddd
					//!  if the leading coefficient is below this threshold
Toshihiro Shimizu 890ddd
{
Toshihiro Shimizu 890ddd
	if (std::abs(c[3]) <= tol)
Toshihiro Shimizu 890ddd
		return solve_2(c, s, tol); // Reduces to second degree
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	Scalar nc[3] = {c[0] / c[3], c[1] / c[3], c[2] / c[3]};
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	Scalar b2 = nc[2] * nc[2];
Toshihiro Shimizu 890ddd
	Scalar p = nc[1] - b2 / 3;
Toshihiro Shimizu 890ddd
	Scalar q = nc[2] * (b2 + b2 - 9 * nc[1]) / 27 + nc[0];
Toshihiro Shimizu 890ddd
	Scalar p3 = p * p * p;
Toshihiro Shimizu 890ddd
	Scalar d = q * q + 4 * p3 / 27;
Toshihiro Shimizu 890ddd
	Scalar offset = -nc[2] / 3;
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	if (d >= 0) {
Toshihiro Shimizu 890ddd
		// Single solution
Toshihiro Shimizu 890ddd
		Scalar z = sqrt(d);
Toshihiro Shimizu 890ddd
		Scalar u = (-q + z) / 2;
Toshihiro Shimizu 890ddd
		Scalar v = (-q - z) / 2;
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
		u = (u < 0) ? -pow(-u, 1 / Scalar(3)) : pow(u, 1 / Scalar(3));
Toshihiro Shimizu 890ddd
		v = (v < 0) ? -pow(-v, 1 / Scalar(3)) : pow(v, 1 / Scalar(3));
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
		s[0] = offset + u + v;
Toshihiro Shimizu 890ddd
		return 1;
Toshihiro Shimizu 890ddd
	}
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	assert(p3 < 0);
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	Scalar u = sqrt(-p / Scalar(3));
Toshihiro Shimizu 890ddd
	Scalar v = acos(-sqrt(-27 / p3) * q / Scalar(2)) / Scalar(3);
Toshihiro Shimizu 890ddd
	Scalar m = cos(v), n = sin(v) * 1.7320508075688772935274463415059; // sqrt(3)
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	s[0] = offset - u * (n + m);
Toshihiro Shimizu 890ddd
	s[1] = offset + u * (n - m);
Toshihiro Shimizu 890ddd
	s[2] = offset + u * (m + m);
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	assert(s[0] <= s[1] && s[1] <= s[2]);
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
	return 3;
Toshihiro Shimizu 890ddd
}
Toshihiro Shimizu 890ddd
}
Toshihiro Shimizu 890ddd
} // namespace tcg::poly_ops
Toshihiro Shimizu 890ddd
Toshihiro Shimizu 890ddd
#endif // TCG_POLY_OPS_H