73 lines · plain
1/*2 * Copyright Thomas Dybdahl Ahle, Nick Thompson, Matt Borland, John Maddock, 20233 * Use, modification and distribution are subject to the4 * Boost Software License, Version 1.0. (See accompanying file5 * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 */7#ifndef BOOST_MATH_TOOLS_ESTRIN_HPP8#define BOOST_MATH_TOOLS_ESTRIN_HPP9 10#include <array>11#include <vector>12#include <type_traits>13#include <boost/math/tools/assert.hpp>14 15namespace boost {16namespace math {17namespace tools {18 19template <typename RandomAccessContainer1, typename RandomAccessContainer2, typename RealOrComplex>20inline RealOrComplex evaluate_polynomial_estrin(RandomAccessContainer1 const &coeffs, RandomAccessContainer2 &scratch, RealOrComplex z) {21 // Does anyone care about the complex coefficients, real argument case?22 // I've never seen it used, and this static assert makes the error messages much better:23 static_assert(std::is_same<typename RandomAccessContainer2::value_type, RealOrComplex>::value,24 "The value type of the scratch space must be the same as the abscissa.");25 auto n = coeffs.size();26 BOOST_MATH_ASSERT_MSG(scratch.size() >= (n + 1) / 2, "The scratch space must be at least N+1/2");27 28 if (n == 0) {29 return static_cast<RealOrComplex>(0);30 }31 for (decltype(n) i = 0; i < n / 2; i++) {32 scratch[i] = coeffs[2 * i] + coeffs[2 * i + 1] * z;33 }34 if (n & 1) {35 scratch[n / 2] = coeffs[n - 1];36 }37 auto m = (n + 1) / 2;38 39 while (m != 1) {40 z = z * z;41 for (decltype(n) i = 0; i < m / 2; i++) {42 scratch[i] = scratch[2 * i] + scratch[2 * i + 1] * z;43 }44 if (m & 1) {45 scratch[m / 2] = scratch[m - 1];46 }47 m = (m + 1) / 2;48 }49 return scratch[0];50}51 52// The std::array template specialization doesn't need to allocate:53template <typename RealOrComplex1, size_t n, typename RealOrComplex2>54inline RealOrComplex2 evaluate_polynomial_estrin(const std::array<RealOrComplex1, n> &coeffs, RealOrComplex2 z) {55 std::array<RealOrComplex2, (n + 1) / 2> ds;56 return evaluate_polynomial_estrin(coeffs, ds, z);57}58 59template <typename RandomAccessContainer, typename RealOrComplex>60inline RealOrComplex evaluate_polynomial_estrin(const RandomAccessContainer &coeffs, RealOrComplex z) {61 auto n = coeffs.size();62 // Normally, I'd make `ds` a RandomAccessContainer, but its value type needs to be RealOrComplex,63 // and the value_type of the passed RandomAccessContainer can just be Real.64 // Allocation of the std::vector is not ideal, but I have no other ideas at the moment:65 std::vector<RealOrComplex> ds((n + 1) / 2);66 return evaluate_polynomial_estrin(coeffs, ds, z);67}68 69} // namespace tools70} // namespace math71} // namespace boost72#endif73