diff --git a/.idea/runConfigurations/bench_ratf.xml b/.idea/runConfigurations/bench_ratf.xml
new file mode 100644
index 0000000..383f929
--- /dev/null
+++ b/.idea/runConfigurations/bench_ratf.xml
@@ -0,0 +1,7 @@
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/ALFI/ALFI.h b/ALFI/ALFI.h
index 46649e2..9fba159 100644
--- a/ALFI/ALFI.h
+++ b/ALFI/ALFI.h
@@ -5,4 +5,5 @@
#include "ALFI/dist.h"
#include "ALFI/misc.h"
#include "ALFI/poly.h"
+#include "ALFI/ratf.h"
#include "ALFI/spline.h"
\ No newline at end of file
diff --git a/ALFI/ALFI/ratf.h b/ALFI/ALFI/ratf.h
new file mode 100644
index 0000000..c13102c
--- /dev/null
+++ b/ALFI/ALFI/ratf.h
@@ -0,0 +1,225 @@
+#pragma once
+
+#include
+
+#include "config.h"
+#include "util/numeric.h"
+#include "util/poly.h"
+
+/**
+ @namespace alfi::ratf
+ @brief Namespace providing support for rational functions.
+
+ This namespace provides types and functions for representing, computing and evaluating rational functions of the form
+ \f[
+ f(x) = \frac{A(x)}{B(x)}
+ \f]
+ where \f(A(x)\f) and \f(B(x)\f) are polynomials.
+*/
+namespace alfi::ratf {
+ /**
+ @brief Represents a rational function \f(\displaystyle f(x) = \frac{A(x)}{B(x)}\f), where \f(A(x)\f) and \f(B(x)\f) are polynomials.
+
+ A pair (`std::pair`) of polynomials `{.first as numerator, .second as denominator}` stored as containers of coefficients in descending degree order.
+ */
+ template class Container = DefaultContainer>
+ using RationalFunction = std::pair, Container>;
+
+ /**
+ @brief Evaluates the rational function at a scalar point using Horner's method.
+
+ Computes \f(f(x) = \frac{A(x)}{B(x)}\f) by evaluating the values of the numerator \f(A\f) and the denominator \f(B\f)
+ using Horner's method, and then dividing the results.
+
+ @ref val utilizes this function for \f(|x| \leq 1\f).
+ */
+ template class Container = DefaultContainer>
+ std::enable_if_t::value, Number>
+ val_mul(const RationalFunction& rf, Number x) {
+ Number n = 0;
+ for (const auto& c : rf.first) {
+ n = n * x + c;
+ }
+ Number d = 0;
+ for (const auto& c : rf.second) {
+ d = d * x + c;
+ }
+ return n / d;
+ }
+
+ /**
+ @brief Evaluates the rational function at each point in the container using @ref val_mul for scalar values.
+ */
+ template class Container = DefaultContainer>
+ std::enable_if_t>::value, Container>
+ val_mul(const RationalFunction& rf, const Container& xx) {
+ Container result(xx.size());
+#if defined(_OPENMP) && !defined(ALFI_DISABLE_OPENMP)
+#pragma omp parallel for
+#endif
+ for (SizeT i = 0; i < xx.size(); ++i) {
+ result[i] = val_mul(rf, xx[i]);
+ }
+ return result;
+ }
+
+ /**
+ @brief Evaluates the rational function at a scalar point by factoring out powers of x.
+
+ Computes \f(f(x) = \frac{A(x)}{B(x)}\f) by evaluating both numerator and denominator in reverse order,
+ effectively factoring out the dominant power of \f(x\f) to improve numerical stability for large \f(|x|\f).
+
+ @ref val utilizes this function for \f(|x| > 1\f).
+ */
+ template class Container = DefaultContainer>
+ std::enable_if_t::value, Number>
+ val_div(const RationalFunction& rf, Number x) {
+ const auto& numerator = rf.first;
+ const auto& denominator = rf.second;
+ Number n = 0;
+ for (auto i = numerator.rbegin(); i != numerator.rend(); ++i) {
+ n = n / x + *i;
+ }
+ Number d = 0;
+ for (auto i = denominator.rbegin(); i != denominator.rend(); ++i) {
+ d = d / x + *i;
+ }
+ const auto numerator_degree = numerator.empty() ? 0 : numerator.size() - 1;
+ const auto denominator_degree = denominator.empty() ? 0 : denominator.size() - 1;
+ if (numerator_degree >= denominator_degree) {
+ return n / d * util::numeric::binpow(x, numerator_degree - denominator_degree);
+ } else {
+ return n / d / util::numeric::binpow(x, denominator_degree - numerator_degree);
+ }
+ }
+
+ /**
+ @brief Evaluates the rational function at each point in the container using @ref val_div for scalar values.
+ */
+ template class Container = DefaultContainer>
+ std::enable_if_t>::value, Container>
+ val_div(const RationalFunction& rf, const Container& xx) {
+ Container result(xx.size());
+#if defined(_OPENMP) && !defined(ALFI_DISABLE_OPENMP)
+#pragma omp parallel for
+#endif
+ for (SizeT i = 0; i < xx.size(); ++i) {
+ result[i] = val_div(rf, xx[i]);
+ }
+ return result;
+ }
+
+ /**
+ @brief Evaluates the rational function at a scalar point.
+
+ Calls @ref val_mul for \f(|x| \leq 1\f) and @ref val_div otherwise for the sake of numerical stability.
+ */
+ template class Container = DefaultContainer>
+ std::enable_if_t::value, Number>
+ val(const RationalFunction& rf, Number x) {
+ if (std::abs(x) <= 1) {
+ return val_mul(rf, x);
+ } else {
+ return val_div(rf, x);
+ }
+ }
+
+ /**
+ @brief Evaluates the rational function at each point in the container using @ref val for scalar values.
+ */
+ template class Container = DefaultContainer>
+ std::enable_if_t>::value, Container>
+ val(const RationalFunction& rf, const Container& xx) {
+ Container result(xx.size());
+#if defined(_OPENMP) && !defined(ALFI_DISABLE_OPENMP)
+#pragma omp parallel for
+#endif
+ for (SizeT i = 0; i < xx.size(); ++i) {
+ result[i] = val(rf, xx[i]);
+ }
+ return result;
+ }
+
+ /**
+ @brief Computes the [@p n / @p m] Pade approximant of the polynomial @p P about the point \f(x = 0\f).
+
+ The Pade approximant is given by the formula
+ \f[
+ P(x) = \frac{A(x)}{B(x)} = \frac{\sum_{i=0}^{n}{a_ix^i}}{b_0+\sum_{j=0}^{m}{b_jx^j}}
+ \f]
+ where \f(A(x)\f) and \f(B(x)\f) are the numerator and denominator polynomials, respectively; \f(b_0 \neq 0\f).
+
+ This function computes the Pade approximant of a polynomial by applying a modified extended Euclidean algorithm for polynomials.
+
+ The modification consists in that:
+ - the algorithm may terminate early if the numerator's degree already meets the requirements,
+ - or perform an extra iteration involving a division by a zero polynomial in special cases.
+
+ The latter is necessary to avoid false negatives, for example, when computing the `[2/2]` approximant of the function \f(x^5\f).
+
+ Without the additional check, this also may lead to false positives, as in the case of computing the `[2/2]` approximant of \f(x^4\f).@n
+ This is prevented by verifying that the constant term of the denominator is non-zero after the algorithm completes.
+
+ @param P the polynomial to approximate (a container of coefficients in descending degree order)
+ @param n the maximum degree for the numerator
+ @param m the maximum degree for the denominator
+ @param epsilon the tolerance used to determine whether a coefficient is considered zero (default is machine epsilon)
+ @return a pair `{numerator, denominator}` representing the Pade approximant; if an approximant does not exist, an empty pair is returned
+ */
+ template class Container = DefaultContainer>
+ RationalFunction pade(Container P, SizeT n, SizeT m, Number epsilon = std::numeric_limits::epsilon()) {
+ if constexpr (std::is_signed_v) {
+ if (n < 0 || m < 0) {
+ return {{}, {}};
+ }
+ }
+
+ util::poly::normalize(P);
+
+ Container Xnm1(n + m + 2, 0);
+ Xnm1[0] = 1;
+
+ // Modified extended Euclidean algorithm `gcd(a,b)=as+bt` without s variable
+ // a = Xnm1
+ // b = P
+ Container old_r, r, old_t, t;
+ if (Xnm1.size() >= P.size()) {
+ old_r = std::move(Xnm1), r = std::move(P);
+ old_t = {0}, t = {1};
+ } else {
+ old_r = std::move(P), r = std::move(Xnm1);
+ old_t = {1}, t = {0};
+ }
+
+ // `old_r.size()` strictly decreases, except maybe the first iteration
+ // ReSharper disable once CppDFALoopConditionNotUpdated
+ while (old_r.size() > n + 1) {
+ auto [q, new_r] = util::poly::div(old_r, r, epsilon);
+
+ const auto qt = util::poly::mul(q, t);
+ const auto new_t_size = std::max(old_t.size(), qt.size());
+ Container new_t(new_t_size, 0);
+ for (SizeT i = 0, offset = new_t_size - old_t.size(); i < old_t.size(); ++i) {
+ new_t[offset+i] = old_t[i];
+ }
+ for (SizeT i = 0, offset = new_t_size - qt.size(); i < qt.size(); ++i) {
+ new_t[offset+i] -= qt[i];
+ }
+
+ old_r = std::move(r);
+ r = std::move(new_r);
+ old_t = std::move(t);
+ t = std::move(new_t);
+
+ util::poly::normalize(old_r, epsilon);
+ }
+
+ util::poly::normalize(old_t, epsilon);
+
+ if (old_t.size() > m + 1 || std::abs(old_t[old_t.size()-1]) <= epsilon) {
+ return {{}, {}};
+ }
+
+ return std::make_pair(old_r, old_t);
+ }
+}
\ No newline at end of file
diff --git a/ALFI/ALFI/util/numeric.h b/ALFI/ALFI/util/numeric.h
index 24b5965..033bae3 100644
--- a/ALFI/ALFI/util/numeric.h
+++ b/ALFI/ALFI/util/numeric.h
@@ -9,4 +9,34 @@ namespace alfi::util::numeric {
bool are_equal(Number a, Number b, Number epsilon = std::numeric_limits::epsilon()) {
return std::abs(a - b) <= epsilon || std::abs(a - b) <= std::max(std::abs(a), std::abs(b)) * epsilon;
}
+
+ /**
+ @brief Computes the power of a number using binary exponentiation.
+
+ Calculates \f(x^n\f) in \f(O(\log{n})\f) operations using the binary (exponentiation by squaring) method.
+
+ It supports both signed and unsigned exponent types (@p ExponentType).@n
+ If the exponent is negative, the function computes the reciprocal of the positive exponentiation.
+
+ @param x the base
+ @param n the exponent
+ @return \f(x^n\f)
+ */
+ template
+ Number binpow(Number x, ExponentType n) {
+ if constexpr (std::is_signed_v) {
+ if (n < 0) {
+ return 1 / binpow(x, -n);
+ }
+ }
+ Number result = 1;
+ while (n > 0) {
+ if ((n & 1) == 1) {
+ result *= x;
+ }
+ x *= x;
+ n >>= 1;
+ }
+ return result;
+ }
}
\ No newline at end of file
diff --git a/ALFI/ALFI/util/poly.h b/ALFI/ALFI/util/poly.h
new file mode 100644
index 0000000..0aed7bb
--- /dev/null
+++ b/ALFI/ALFI/util/poly.h
@@ -0,0 +1,116 @@
+#pragma once
+
+#include
+#include
+
+#include "../config.h"
+
+/**
+ @namespace alfi::util::poly
+ @brief Namespace with utility functions for polynomial operations.
+
+ This namespace provides a set of functions for manipulating polynomials, such as normalization, multiplication, and division.
+ All polynomials are represented as containers of coefficients in descending degree order.
+*/
+namespace alfi::util::poly {
+ /**
+ @brief Normalizes a polynomial by removing leading zero coefficients.
+
+ Removes leading coefficients that are considered zero within a given tolerance @p epsilon.
+ If all coefficients are close to zero, the last one is preserved.
+ If the container is empty, a single zero coefficient is inserted.
+
+ @param p the polynomial to normalize (a container of coefficients in descending degree order)
+ @param epsilon the tolerance used to determine whether a coefficient is considered zero (default is machine epsilon)
+ */
+ template class Container = DefaultContainer>
+ void normalize(Container& p, Number epsilon = std::numeric_limits::epsilon()) {
+ if (p.empty()) {
+ return p.push_back(0);
+ }
+ auto p_start = std::find_if(p.begin(), p.end(), [&epsilon](Number v) { return std::abs(v) > epsilon; });
+ if (p_start == p.end()) {
+ --p_start;
+ }
+ if (p_start > p.begin()) {
+ p.erase(p.begin(), p_start);
+ }
+ }
+
+ /**
+ @brief Multiplies two polynomials.
+
+ Given two polynomials \f(p_1\f) and \f(p_2\f) represented as containers of coefficients,
+ this function computes their product using the convolution formula:
+ \f[
+ (p_1 \cdot p_2)[k] = \sum_{i+j=k}{p_1[i] \cdot p_2[j]}
+ \f]
+
+ If either polynomial is empty, the function returns an empty container.
+
+ @param p1 the first polynomial
+ @param p2 the second polynomial
+ @return the product polynomial (either empty or of size `p1.size() + p2.size() - 1`)
+ */
+ template class Container = DefaultContainer>
+ Container mul(const Container& p1, const Container& p2) {
+ if (p1.empty() || p2.empty()) {
+ return {};
+ }
+ Container result(p1.size() + p2.size() - 1);
+ for (SizeT i = 0; i < p1.size(); ++i) {
+ for (SizeT j = 0; j < p2.size(); ++j) {
+ result[i+j] += p1[i] * p2[j];
+ }
+ }
+ return result;
+ }
+
+ /**
+ @brief Divides one polynomial by another.
+
+ This function divides the @p dividend polynomial \f(A\f) by the @p divisor polynomial \f(B\f),
+ and returns the @p quotient \f(Q\f) and the @p remainder \f(R\f) such that:
+ \f[
+ A(x) = B(X) \cdot Q(x) + R(x)
+ \f]
+ with either \f(R\f) being effectively zero or the degree of \f(R\f) being lower than the degree of \f(B\f).
+
+ The division is performed using a tolerance @p epsilon to determine when a coefficient is considered zero.
+
+ If the divisor is effectively zero or if the dividend has a lower degree than the divisor,
+ the function returns an empty quotient and the dividend as the remainder.
+
+ @param dividend the polynomial to be divided
+ @param divisor the polynomial to divide by
+ @param epsilon the tolerance used to determine whether a coefficient is considered zero (default is machine epsilon)
+ @return a pair `{quotient, remainder}`
+ */
+ template class Container = DefaultContainer>
+ std::pair,Container> div(const Container& dividend, const Container& divisor, Number epsilon = std::numeric_limits::epsilon()) {
+ const auto divisor_start = std::find_if(divisor.begin(), divisor.end(), [&epsilon](Number v) { return std::abs(v) > epsilon; });
+
+ if (divisor_start == divisor.end() || dividend.size() < divisor.size()) {
+ return {{}, dividend};
+ }
+
+ const auto divisor_start_idx = divisor_start - divisor.begin();
+
+ const auto n = dividend.size();
+ const auto m = divisor.size() - divisor_start_idx;
+
+ Container quotient(n - m + 1, 0);
+ Container remainder = dividend;
+
+ for (SizeT i = 0; i <= n - m; ++i) {
+ const Number factor = remainder[i] / divisor[divisor_start_idx];
+ quotient[i] = factor;
+ for (SizeT j = 0; j < m; ++j) {
+ remainder[i+j] -= factor * divisor[divisor_start_idx+j];
+ }
+ }
+
+ remainder.erase(remainder.begin(), remainder.end() - m + 1);
+ return {quotient, remainder};
+ }
+}
\ No newline at end of file
diff --git a/benches/CMakeLists.txt b/benches/CMakeLists.txt
index 749b26b..032d214 100644
--- a/benches/CMakeLists.txt
+++ b/benches/CMakeLists.txt
@@ -14,6 +14,7 @@ macro(add_bench_executable bench_name)
endmacro()
add_bench_executable(bench_poly poly/bench_poly.cpp)
+add_bench_executable(bench_ratf ratf/bench_ratf.cpp)
add_bench_executable(bench_barycentric misc/bench_barycentric.cpp)
add_bench_executable(bench_polyeqv spline/bench_polyeqv.cpp)
add_bench_executable(bench_step spline/bench_step.cpp)
diff --git a/benches/ratf/bench_ratf.cpp b/benches/ratf/bench_ratf.cpp
new file mode 100644
index 0000000..7da0a40
--- /dev/null
+++ b/benches/ratf/bench_ratf.cpp
@@ -0,0 +1,94 @@
+#include
+
+#include
+
+#include "../bench_utils.h"
+
+static const std::vector> intervals = {
+ {-1, 1},
+ {5, 10},
+ {0, 2},
+};
+
+static void BM_val_scalar(benchmark::State& state) {
+ const auto& n = state.range(0);
+ const auto& interval = intervals[state.range(1)];
+ alfi::ratf::RationalFunction<> rf;
+ rf.first.resize(n);
+ rf.second.resize(n);
+ for (auto& c : rf.first) {
+ c = random_value(-1, 1);
+ }
+ for (auto& c : rf.second) {
+ c = random_value(-1, 1);
+ }
+ const auto x = random_value(interval.first, interval.second);
+ double result;
+ for (auto _ : state) {
+ result = alfi::ratf::val(rf, x);
+ benchmark::DoNotOptimize(result);
+ benchmark::ClobberMemory();
+ }
+}
+BENCHMARK(BM_val_scalar)
+ ->ArgsProduct({
+ benchmark::CreateRange(1 << 3, 1 << 12, 1 << 3),
+ benchmark::CreateDenseRange(0, static_cast(intervals.size()) - 1, 1),
+ })
+ ->ArgNames({"n", "interval"});
+
+static void BM_val_vector(benchmark::State& state) {
+ const auto& n = state.range(0);
+ const auto& nn = state.range(1);
+ const auto& interval = intervals[state.range(2)];
+ alfi::ratf::RationalFunction<> rf;
+ rf.first.resize(n);
+ rf.second.resize(n);
+ for (auto& c : rf.first) {
+ c = random_value(-1, 1);
+ }
+ for (auto& c : rf.second) {
+ c = random_value(-1, 1);
+ }
+ std::vector xx(nn);
+ for (auto& x : xx) {
+ x = random_value(interval.first, interval.second);
+ }
+ std::vector result;
+ for (auto _ : state) {
+ result = alfi::ratf::val(rf, xx);
+ benchmark::DoNotOptimize(result);
+ benchmark::ClobberMemory();
+ }
+}
+BENCHMARK(BM_val_vector)
+ ->ArgsProduct({
+ benchmark::CreateRange(1 << 3, 1 << 12, 1 << 3),
+ benchmark::CreateRange(1 << 3, 1 << 12, 1 << 3),
+ benchmark::CreateDenseRange(0, static_cast(intervals.size()) - 1, 1),
+ })
+ ->ArgNames({"n", "nn", "interval"});
+
+static void BM_pade(benchmark::State& state) {
+ const auto& size = state.range(0);
+ const auto degree = size - 1;
+ const auto n = degree/2, m = n;
+ std::vector p(size);
+ for (size_t i = 0; i < p.size(); ++i) {
+ p[i] = random_value(-2, 2) * static_cast(i + 1);
+ }
+ alfi::ratf::RationalFunction<> result;
+ for (auto _ : state) {
+ result = alfi::ratf::pade(p, n, m);
+ benchmark::DoNotOptimize(result);
+ benchmark::ClobberMemory();
+ }
+ state.SetBytesProcessed(state.iterations() * size * static_cast(sizeof(double)));
+ state.SetComplexityN(size);
+}
+BENCHMARK(BM_pade)
+ ->ArgsProduct({benchmark::CreateRange(1, 1 << 16, 2)})
+ ->ArgName("n")
+ ->Complexity();
+
+BENCHMARK_MAIN();
\ No newline at end of file
diff --git a/examples/interpolation/interpolation.cpp b/examples/interpolation/interpolation.cpp
index 1a1cd7e..ef6d406 100644
--- a/examples/interpolation/interpolation.cpp
+++ b/examples/interpolation/interpolation.cpp
@@ -116,6 +116,7 @@ class PlotWindow final : public QWidget {
_function_checkbox = new QCheckBox("Function");
_points_checkbox = new QCheckBox("Points");
_poly_checkbox = new QCheckBox("Polynomial");
+ _pade_checkbox = new QCheckBox("Pade Approximant");
_barycentric_checkbox = new QCheckBox("Barycentric Formula");
_poly_eqv_spline_checkbox = new QCheckBox("Poly. Eqv. Spline");
_step_spline_checkbox = new QCheckBox("Step Spline");
@@ -163,7 +164,7 @@ class PlotWindow final : public QWidget {
connect(_b_spin_box, QOverload::of(&QDoubleSpinBox::valueChanged), this, &PlotWindow::update_plot);
for (const auto checkbox : {
- _function_checkbox, _points_checkbox, _poly_checkbox, _barycentric_checkbox, _poly_eqv_spline_checkbox,
+ _function_checkbox, _points_checkbox, _poly_checkbox, _pade_checkbox, _barycentric_checkbox, _poly_eqv_spline_checkbox,
_step_spline_checkbox, _linear_spline_checkbox, _quadratic_spline_checkbox, _cubic_spline_checkbox}) {
connect(checkbox, &QCheckBox::toggled, this, &PlotWindow::update_plot);
}
@@ -211,6 +212,7 @@ class PlotWindow final : public QWidget {
interpolation_layout->addWidget(_points_checkbox);
interpolation_layout->addWidget(_poly_checkbox);
interpolation_layout->addWidget(_poly_combo);
+ interpolation_layout->addWidget(_pade_checkbox);
interpolation_layout->addWidget(_barycentric_checkbox);
interpolation_layout->addWidget(_barycentric_combo);
interpolation_layout->addWidget(_poly_eqv_spline_checkbox);
@@ -340,6 +342,9 @@ class PlotWindow final : public QWidget {
if (_poly_checkbox->isChecked()) {
add_graph("Polynomial", xx, poly_types[_poly_combo->currentIndex()].second(X, Y, xx));
}
+ if (_pade_checkbox->isChecked()) {
+ add_graph("Pade Approximation", xx, alfi::ratf::val(alfi::ratf::pade(alfi::poly::newton(X, Y), N / 2, N / 2), xx));
+ }
if (_barycentric_checkbox->isChecked()) {
int barycentric_dist_type = _barycentric_combo->currentIndex();
if (barycentric_dist_type == 0) {
@@ -387,6 +392,7 @@ class PlotWindow final : public QWidget {
QCheckBox* _function_checkbox;
QCheckBox* _points_checkbox;
QCheckBox* _poly_checkbox;
+ QCheckBox* _pade_checkbox;
QCheckBox* _barycentric_checkbox;
QCheckBox* _poly_eqv_spline_checkbox;
QCheckBox* _step_spline_checkbox;
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index b74a0c0..9f1ae98 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -19,6 +19,8 @@ add_test_executable(test_dist dist/test_dist.cpp)
add_test_executable(test_poly poly/test_poly.cpp)
+add_test_executable(test_ratf ratf/test_ratf.cpp)
+
add_test_executable(test_barycentric misc/test_barycentric.cpp)
add_test_executable(test_polyeqv spline/test_polyeqv.cpp)
diff --git a/tests/ratf/test_ratf.cpp b/tests/ratf/test_ratf.cpp
new file mode 100644
index 0000000..2164e94
--- /dev/null
+++ b/tests/ratf/test_ratf.cpp
@@ -0,0 +1,60 @@
+#include
+#include
+
+#include "../test_utils.h"
+
+TEST(RationalFunctionsTest, Evaluation) {
+ alfi::ratf::RationalFunction<> rf({}, {});
+ for (const auto x : alfi::dist::uniform(11, -2.0, 2.0)) {
+ EXPECT_TRUE(std::isnan(alfi::ratf::val(rf, x)));
+ }
+ rf = {{1, 0}, {1}};
+ for (const auto x : alfi::dist::uniform(11, -2.0, 2.0)) {
+ EXPECT_EQ(alfi::ratf::val(rf, x), x);
+ }
+ rf = {{2, 3, 4}, {1, 2, 2}};
+ EXPECT_EQ(alfi::ratf::val(rf, -2.0), 3.0);
+ EXPECT_EQ(alfi::ratf::val(rf, -1.0), 3.0);
+ EXPECT_EQ(alfi::ratf::val(rf, 0.0), 2.0);
+ EXPECT_EQ(alfi::ratf::val(rf, 1.0), 9.0/5.0);
+ EXPECT_EQ(alfi::ratf::val(rf, 2.0), 9.0/5.0);
+}
+
+TEST(RationalFunctionsTest, Pade) {
+ // [0/0] empty array (effectively 0)
+ auto rf = alfi::ratf::pade({}, 0, 0);
+ expect_eq(rf.first, {0});
+ expect_eq(rf.second, {1});
+ // [1/1] empty array (effectively 0)
+ rf = alfi::ratf::pade({}, 1, 1);
+ expect_eq(rf.first, {0});
+ expect_eq(rf.second, {1});
+ // [1/1] 10
+ rf = alfi::ratf::pade({10}, 1, 1);
+ expect_eq(rf.first, {10});
+ expect_eq(rf.second, {1});
+ // [2/2] x^4 + 2x^3 + 3x^2 + 4x + 5
+ rf = alfi::ratf::pade({1, 2, 3, 4, 5}, 2, 2);
+ expect_eq(rf.first, {-6, 5});
+ expect_eq(rf.second, {1, -2, 1});
+ // [2/2] x^4 - 2x^3 + 3x^2 - 4x + 5
+ rf = alfi::ratf::pade({1, -2, 3, -4, 5}, 2, 2);
+ expect_eq(rf.first, {6, 5});
+ expect_eq(rf.second, {1, 2, 1});
+ // [2/2] exp(x)
+ rf = alfi::ratf::pade({1.0/2/3/4, 1.0/2/3, 1.0/2, 1, 1}, 2, 2);
+ expect_eq(rf.first, {0.25, 1.5, 3});
+ expect_eq(rf.second, {0.25, -1.5, 3});
+ // [2/2] sin(x)
+ rf = alfi::ratf::pade({1.0/2/3/4, 0, -1.0/2, 0, 1}, 2, 2);
+ expect_eq(rf.first, {-5.0/12, 0, 1});
+ expect_eq(rf.second, {1.0/12, 0, 1});
+ // [2/2] x^5
+ rf = alfi::ratf::pade({1, 0, 0, 0, 0, 0}, 2, 2);
+ expect_eq(rf.first, {0});
+ expect_eq(rf.second, {-1});
+ // [2/2] x^4
+ rf = alfi::ratf::pade({1, 0, 0, 0, 0}, 2, 2);
+ expect_eq(rf.first, {});
+ expect_eq(rf.second, {});
+}
\ No newline at end of file