Skip to content

Commit

Permalink
Re-working API for polynomial operations (#527)
Browse files Browse the repository at this point in the history 
* Re-working API for polynomial operations

* Apply clang-format

* Fix benchmarks

---------

Co-authored-by: github-actions[bot] <github-actions[bot]@users.noreply.github.com>
  • Loading branch information
@jatinchowdhury18 @github-actions
jatinchowdhury18 and github-actions[bot] authored Apr 22, 2024
1 parent cf53b11 commit ba945fa
Show file tree
Hide file tree
Showing 5 changed files with 275 additions and 127 deletions.
36 changes: 18 additions & 18 deletions bench/PolynomialBench.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,30 +11,30 @@ auto outVecFloat = std::vector<float> (N, 0.0f);

using namespace chowdsp::Polynomials;

#define POLY_BENCH_7(mode) \
static void bench_Poly7_##mode (benchmark::State& state) \
{ \
for (auto _ : state) \
{ \
for (size_t n = 0; n < N; ++n) \
outVecFloat[n] = mode<7, float> ({ 7.0f, 6.0f, 5.0f, 4.0f, 3.0f, 2.0f, 1.0f, 0.0 }, inputVecFloat[n]); \
} \
} \
#define POLY_BENCH_7(mode) \
static void bench_Poly7_##mode (benchmark::State& state) \
{ \
for (auto _ : state) \
{ \
for (size_t n = 0; n < N; ++n) \
outVecFloat[n] = mode<7, float> ({ { 7.0f, 6.0f, 5.0f, 4.0f, 3.0f, 2.0f, 1.0f, 0.0f }, poly_order_descending {} }, inputVecFloat[n]); \
} \
} \
BENCHMARK (bench_Poly7_##mode)->MinTime (1);

POLY_BENCH_7 (naive)
POLY_BENCH_7 (estrin)
POLY_BENCH_7 (horner)

#define POLY_BENCH_10(mode) \
static void bench_Poly10_##mode (benchmark::State& state) \
{ \
for (auto _ : state) \
{ \
for (size_t n = 0; n < N; ++n) \
outVecFloat[n] = mode<10, float> ({ 10.0f, 9.0f, 8.0f, 7.0f, 6.0f, 5.0f, 4.0f, 3.0f, 2.0f, 1.0f, 0.0 }, inputVecFloat[n]); \
} \
} \
#define POLY_BENCH_10(mode) \
static void bench_Poly10_##mode (benchmark::State& state) \
{ \
for (auto _ : state) \
{ \
for (size_t n = 0; n < N; ++n) \
outVecFloat[n] = mode<10, float> ({ { 10.0f, 9.0f, 8.0f, 7.0f, 6.0f, 5.0f, 4.0f, 3.0f, 2.0f, 1.0f, 0.0 }, poly_order_descending {} }, inputVecFloat[n]); \
} \
} \
BENCHMARK (bench_Poly10_##mode)->MinTime (1);

POLY_BENCH_10 (naive)
Expand Down
8 changes: 4 additions & 4 deletions modules/dsp/chowdsp_math/Math/chowdsp_LogApprox.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,14 +31,14 @@ namespace LogApprox
constexpr auto beta = (T) -1.098865286222744;
constexpr auto gamma = (T) 3.148297929334117;
constexpr auto zeta = (T) -2.213475204444817;
return chowdsp::Polynomials::estrin<3> ({ alpha, beta, gamma, zeta }, x);
return Polynomials::estrin<3> (Polynomial<T, 3> { { alpha, beta, gamma, zeta } }, x);
}
else if constexpr (order == 2)
{
constexpr auto beta = (T) -0.33985;
constexpr auto gamma = (T) 2.01955;
constexpr auto zeta = (T) -1.6797;
return chowdsp::Polynomials::estrin<2> ({ beta, gamma, zeta }, x);
return Polynomials::estrin<2> (Polynomial<T, 2> { { beta, gamma, zeta } }, x);
}
else if constexpr (order == 1)
{
Expand All @@ -59,14 +59,14 @@ namespace LogApprox
constexpr auto beta = (T) -1.098865286222744;
constexpr auto gamma = (T) 3.148297929334117;
constexpr auto zeta = (T) -2.213475204444817;
return chowdsp::Polynomials::estrin<3> ({ alpha, beta, gamma, zeta }, x);
return Polynomials::estrin<3> (Polynomial<T, 3> { { alpha, beta, gamma, zeta } }, x);
}
else if constexpr (order == 2)
{
constexpr auto beta = (T) -0.33985;
constexpr auto gamma = (T) 2.01955;
constexpr auto zeta = (T) -1.6797;
return chowdsp::Polynomials::estrin<2> ({ beta, gamma, zeta }, x);
return Polynomials::estrin<2> (Polynomial<T, 2> { { beta, gamma, zeta } }, x);
}
else if constexpr (order == 1)
{
Expand Down
210 changes: 176 additions & 34 deletions modules/dsp/chowdsp_math/Math/chowdsp_Polynomials.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,91 +3,233 @@
JUCE_BEGIN_IGNORE_WARNINGS_GCC_LIKE ("-Wsign-conversion",
"-Waggressive-loop-optimizations")

namespace chowdsp
{
/**
* Type tag to denote a polynomial stored with coefficients in descending order,
* i.e. { a_n, a_n-1, ..., a_1, a_0 }.
*/
struct poly_order_descending
{
};

/**
* Type tag to denote a polynomial stored with coefficients in ascending order,
* i.e. { a_0, a_1, ..., a_n-1, a_n }.
*/
struct poly_order_ascending
{
};

/** An array wrapper for polynomial coefficients, stored in either ascending or descending order. */
template <typename T, size_t poly_order, typename storage_mode = poly_order_descending>
struct Polynomial
{
static constexpr auto num_coeffs = poly_order + 1;
std::array<T, num_coeffs> coeffs {};

Polynomial() = default;

constexpr Polynomial (const std::array<T, num_coeffs>& init_coeffs, [[maybe_unused]] storage_mode tag = storage_mode {}) // NOLINT(*-explicit-constructor)
{
set_coefficients (init_coeffs);
}

template <typename new_storage_mode>
constexpr Polynomial<T, poly_order, new_storage_mode> convert() const
{
if constexpr (std::is_same_v<storage_mode, new_storage_mode>)
return *this;

Polynomial<T, poly_order, new_storage_mode> new_poly {};
std::reverse_copy (coeffs.begin(), coeffs.end(), new_poly.coeffs.begin());
return new_poly;
}

constexpr void set_coefficients (const std::array<T, num_coeffs>& new_coeffs)
{
coeffs = new_coeffs;
}
};
} // namespace chowdsp

/** Useful methods for working with and evaluating polynomials */
namespace chowdsp::Polynomials
{
using chowdsp::poly_order_ascending;
using chowdsp::poly_order_descending;
using chowdsp::Polynomial;

/** Useful template type representing the product of two other types. */
template <typename T, typename X>
using P = decltype (T {} * X {});

/**
* Evaluates a polynomial of a given order, using a naive method.
* Coefficients should be given in the form { a_n, a_n-1, ..., a_1, a_0 }
* Coefficients should be given in the descending form.
*/
template <int ORDER, typename T, typename X>
constexpr P<T, X> naive (const T (&coeffs)[ORDER + 1], X x)
template <int poly_order, typename T, typename X>
constexpr P<T, X> naive (const Polynomial<T, poly_order, poly_order_descending>& poly, X x)
{
CHOWDSP_USING_XSIMD_STD (pow);

P<T, X> sum = coeffs[ORDER];
for (int n = 0; n < ORDER; ++n)
sum += coeffs[n] * pow (x, X (ORDER - n));
P<T, X> sum = poly.coeffs[poly_order];
for (int n = 0; n < poly_order; ++n)
sum += poly.coeffs[n] * pow (x, X (poly_order - n));

return sum;
}

/**
* Evaluates a polynomial of a given order, using a naive method.
* Coefficients should be given in the ascending form.
*/
template <int poly_order, typename T, typename X>
constexpr P<T, X> naive (const Polynomial<T, poly_order, poly_order_ascending>& poly, X x)
{
CHOWDSP_USING_XSIMD_STD (pow);

P<T, X> sum = poly.coeffs[0];
for (int n = 1; n <= poly_order; ++n)
sum += poly.coeffs[n] * pow (x, static_cast<X> (n));

return sum;
}

/**
* Evaluates a polynomial of a given order, using Horner's method.
* Coefficients should be given in the form { a_n, a_n-1, ..., a_1, a_0 }
* Coefficients should be given in descending form.
* https://en.wikipedia.org/wiki/Horner%27s_method
*/
template <int ORDER, typename T, typename X>
constexpr P<T, X> horner (const T (&coeffs)[ORDER + 1], X x)
template <int poly_order, typename T, typename X>
constexpr P<T, X> horner (const Polynomial<T, poly_order, poly_order_descending>& poly, X x)
{
P<T, X> b = coeffs[0];
for (int n = 1; n <= ORDER; ++n)
b = b * x + coeffs[n];
P<T, X> b = poly.coeffs[0];
for (int n = 1; n <= poly_order; ++n)
b = b * x + poly.coeffs[n];

return b;
}

/**
* Evaluates a polynomial of a given order, using Horner's method.
* Coefficients should be given in ascending form.
* https://en.wikipedia.org/wiki/Horner%27s_method
*/
template <int poly_order, typename T, typename X>
constexpr P<T, X> horner (const Polynomial<T, poly_order, poly_order_ascending>& poly, X x)
{
P<T, X> b = poly.coeffs.back();
for (int n = poly_order; n > 0; --n)
b = b * x + poly.coeffs[n - 1];

return b;
}

/**
* Evaluates a polynomial of a given order, using Estrin's scheme.
* Coefficients should be given in descending form.
* https://en.wikipedia.org/wiki/Estrin%27s_scheme
*/
template <int poly_order, typename T, typename X>
constexpr P<T, X> estrin (const Polynomial<T, poly_order, poly_order_descending>& poly, X x)
{
if constexpr (poly_order <= 1) // base case
{
return poly.coeffs[1] + poly.coeffs[0] * x;
}
else
{
Polynomial<P<T, X>, poly_order / 2, poly_order_descending> temp {};
for (int n = poly_order; n > 0; n -= 2)
temp.coeffs[n / 2] = poly.coeffs[n] + poly.coeffs[n - 1] * x;

if constexpr (poly_order % 2 == 0) // even order polynomial
temp.coeffs[0] = poly.coeffs[0];

return estrin<poly_order / 2> (temp, x * x); // recurse!
}
}

/**
* Evaluates a polynomial of a given order, using Estrin's scheme.
* Coefficients should be given in the form { a_n, a_n-1, ..., a_1, a_0 }
* Coefficients should be given in descending form.
* https://en.wikipedia.org/wiki/Estrin%27s_scheme
*/
template <int ORDER, typename T, typename X>
constexpr P<T, X> estrin (const T (&coeffs)[ORDER + 1], X x)
template <int poly_order, typename T, typename X>
constexpr P<T, X> estrin (const Polynomial<T, poly_order, poly_order_ascending>& poly, X x)
{
if constexpr (ORDER <= 1) // base case
if constexpr (poly_order <= 1) // base case
{
return coeffs[1] + coeffs[0] * x;
return poly.coeffs[poly_order - 1] + poly.coeffs[poly_order] * x;
}
else
{
P<T, X> temp[ORDER / 2 + 1];
for (int n = ORDER; n > 0; n -= 2)
temp[n / 2] = coeffs[n] + coeffs[n - 1] * x;
Polynomial<P<T, X>, poly_order / 2, poly_order_ascending> temp {};
for (int n = 1; n <= poly_order; n += 2)
temp.coeffs[n / 2] = poly.coeffs[n - 1] + poly.coeffs[n] * x;

if constexpr (ORDER % 2 == 0) // even order polynomial
temp[0] = coeffs[0];
if constexpr (poly_order % 2 == 0) // even order polynomial
temp.coeffs[poly_order / 2] = poly.coeffs[poly_order];

return estrin<ORDER / 2> (temp, x * x); // recurse!
return estrin<poly_order / 2> (temp, x * x); // recurse!
}
}

/**
* Computes the coefficients of the antiderivative of a polynomial.
* Coefficients should be given in the form { a_n, a_n-1, ..., a_1, a_0 }
* Coefficients should be given in descending form.
*/
template <int ORDER, typename T>
constexpr void antiderivative (const T (&coeffs)[ORDER + 1], T (&ad_coeffs)[ORDER + 2], T C = (T) 0)
template <int poly_order, typename T>
constexpr Polynomial<T, poly_order + 1> antiderivative (const Polynomial<T, poly_order>& poly, T C = {})
{
for (int n = 0; n <= ORDER; ++n)
ad_coeffs[n] = coeffs[n] / T (ORDER + 1 - n);
Polynomial<T, poly_order + 1> ad_poly {};
for (int n = 0; n <= poly_order; ++n)
ad_poly.coeffs[n] = poly.coeffs[n] / static_cast<T> (poly_order + 1 - n);
ad_poly.coeffs[poly_order + 1] = C;

ad_coeffs[ORDER + 1] = C;
return ad_poly;
}

/**
* Computes the coefficients of the antiderivative of a polynomial.
* Coefficients should be given in ascending form.
*/
template <int poly_order, typename T>
constexpr Polynomial<T, poly_order + 1, poly_order_ascending> antiderivative (const Polynomial<T, poly_order, poly_order_ascending>& poly, T C = {})
{
Polynomial<T, poly_order + 1, poly_order_ascending> ad_poly {};
ad_poly.coeffs[0] = C;
for (int n = 1; n <= poly_order + 1; ++n)
ad_poly.coeffs[n] = poly.coeffs[n - 1] / static_cast<T> (n);

return ad_poly;
}

/**
* Computes the coefficients of the derivative of a polynomial.
* Coefficients should be given in descending form.
*/
template <int poly_order, typename T>
constexpr Polynomial<T, poly_order - 1> derivative (const Polynomial<T, poly_order>& poly)
{
Polynomial<T, poly_order - 1> d_poly {};
for (int n = 0; n < poly_order; ++n)
d_poly.coeffs[n] = poly.coeffs[n] * static_cast<T> (poly_order - n);
return d_poly;
}

/**
* Computes the coefficients of the derivative of a polynomial.
* Coefficients should be given in the form { a_n, a_n-1, ..., a_1, a_0 }
* Coefficients should be given in ascending form.
*/
template <int ORDER, typename T>
constexpr void derivative (const T (&coeffs)[ORDER + 1], T (&d_coeffs)[ORDER])
template <int poly_order, typename T>
constexpr Polynomial<T, poly_order - 1, poly_order_ascending> derivative (const Polynomial<T, poly_order, poly_order_ascending>& poly)
{
for (int n = 0; n < ORDER; ++n)
d_coeffs[n] = coeffs[n] * T (ORDER - n);
Polynomial<T, poly_order - 1, poly_order_ascending> d_poly {};
for (int n = 0; n < poly_order; ++n)
d_poly.coeffs[n] = poly.coeffs[n + 1] * static_cast<T> (n + 1);
return d_poly;
}
} // namespace chowdsp::Polynomials

Expand Down
8 changes: 4 additions & 4 deletions modules/dsp/chowdsp_math/Math/chowdsp_PowApprox.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,14 +31,14 @@ namespace PowApprox
constexpr auto beta = (T) 0.2274112777602189;
constexpr auto gamma = (T) 0.6931471805599453;
constexpr auto zeta = (T) 1.0;
return chowdsp::Polynomials::estrin<3> ({ alpha, beta, gamma, zeta }, x);
return Polynomials::estrin<3> (Polynomial<T, 3> { { alpha, beta, gamma, zeta } }, x);
}
else if constexpr (order == 2)
{
constexpr auto beta = (T) 0.343146;
constexpr auto gamma = (T) 0.656854;
constexpr auto zeta = (T) 1.0;
return chowdsp::Polynomials::estrin<2> ({ beta, gamma, zeta }, x);
return Polynomials::estrin<2> (Polynomial<T, 2> { { beta, gamma, zeta } }, x);
}
else if constexpr (order == 1)
{
Expand All @@ -59,14 +59,14 @@ namespace PowApprox
constexpr auto beta = (T) 0.2274112777602189;
constexpr auto gamma = (T) 0.6931471805599453;
constexpr auto zeta = (T) 1.0;
return chowdsp::Polynomials::estrin<3> ({ alpha, beta, gamma, zeta }, x);
return Polynomials::estrin<3> (Polynomial<T, 3> { { alpha, beta, gamma, zeta } }, x);
}
else if constexpr (order == 2)
{
constexpr auto beta = (T) 0.343146;
constexpr auto gamma = (T) 0.656854;
constexpr auto zeta = (T) 1.0;
return chowdsp::Polynomials::estrin<2> ({ beta, gamma, zeta }, x);
return Polynomials::estrin<2> (Polynomial<T, 2> { { beta, gamma, zeta } }, x);
}
else if constexpr (order == 1)
{
Expand Down
Loading

0 comments on commit ba945fa

Please sign in to comment.