diff --git a/src/math/map.hpp b/src/math/map.hpp index b040519..94a9c1d 100644 --- a/src/math/map.hpp +++ b/src/math/map.hpp @@ -25,7 +25,7 @@ namespace math { /** * Remaps a number from one range to another. * - * @param value Value to remap. + * @param x Value to remap. * @param from_min Start of the first range. * @param from_max End of the first range. * @param to_min Start of the second range. @@ -33,9 +33,9 @@ namespace math { * @return Unclamped remapped value. */ template -T map(T value, T from_min, T from_max, T to_min, T to_max) +T map(T x, T from_min, T from_max, T to_min, T to_max) { - return to_min + (to_max - to_min) * ((value - from_min) / (from_max - from_min)); + return to_min + (x - from_min) * (to_max - to_min) / (from_max - from_min); } } // namespace math diff --git a/src/math/polynomial.hpp b/src/math/polynomial.hpp index 6ce660e..50a3753 100644 --- a/src/math/polynomial.hpp +++ b/src/math/polynomial.hpp @@ -20,9 +20,12 @@ #ifndef ANTKEEPER_MATH_POLYNOMIAL_HPP #define ANTKEEPER_MATH_POLYNOMIAL_HPP +#include "math/constants.hpp" +#include "math/map.hpp" + namespace math { -/// Polynomial evaluation functions. +/// Polynomial functions. namespace polynomial { /** @@ -37,12 +40,91 @@ namespace polynomial { template T horner(InputIt first, InputIt last, T x) { - T sum = *first; + T y = *first; for (++first; first != last; ++first) - sum = sum * x + *first; - return sum; + y = y * x + *first; + return y; } +/** Chebychev polynomials. + * + * @see https://en.wikipedia.org/wiki/Chebyshev_polynomials + */ +namespace chebyshev { + + /** + * Generates a Chebyshev approximation of a function. + * + * @param[out] first,last Range of Chebyshev polynomial coefficients. + * @param[in] f Unary function to approximate. + * @param[in] min,max Domain of @p f. + */ + template + void approximate(OutputIt first, OutputIt last, UnaryOp f, T min, T max) + { + std::size_t n = last - first; + const T two_over_n = T(2) / static_cast(n); + const T pi_over_n = math::pi / static_cast(n); + + last = first; + for (std::size_t i = 0; i < n; ++i) + *(last++) = T(0); + + for (std::size_t i = 0; i < n; ++i) + { + const T y = pi_over_n * (static_cast(i) + T(0.5)); + + T x = f(math::map(std::cos(y), T(-1), T(1), min, max)) * two_over_n; + + *first += x; + last = first; + for (std::size_t j = 1; j < n; ++j) + { + *(++last) += x * std::cos(y * static_cast(j)); + } + } + } + + /** + * Evaluates a Chebyshev polynomial. + * + * @param[in] first,last Range of Chebychev polynomial coefficients. + * @param[in] x Value on the interval `[-1, 1]`. + * + * @return Evaluated value. + */ + template + T evaluate(InputIt first, InputIt last, T x) + { + T y = *(first++) * T(0.5) + *(first++) * x; + + const T x2 = x * T(2); + for (T n2 = T(1), n1 = x, n0; first != last; n2 = n1, n1 = n0) + { + n0 = x2 * n1 - n2; + y += *(first++) * n0; + } + + return y; + } + + /** + * Evaluates a Chebyshev polynomial. + * + * @param first,last Range of Chebychev polynomial coefficients. + * @param min,max Domain of the approximated function. + * @param x Value on the interval `[min, max]`. + * + * @return Evaluated value. + */ + template + T evaluate(InputIt first, InputIt last, T min, T max, T x) + { + return evaluate(first, last, math::map(x, min, max, T(-1), T(1))); + } + +} // namespace chebyshev + } // namespace polynomial } // namespace math