lagrange_interpolation.hpp¶
Evaluate the degree < values.size() polynomial known at consecutive points 0,1,... in O(n) over a field.
\[
\displaystyle f(x)=\sum_{i=0}^{n} y_i\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}
\]
Implementation¶
#ifndef NOYA_LAGRANGE_INTERPOLATION_HPP
#define NOYA_LAGRANGE_INTERPOLATION_HPP 1
/// @complexity Time: O(n) per evaluation at consecutive points.
/// Space: O(n).
#include <cassert>
#include <cstdint>
#include <type_traits>
#include <vector>
namespace noya {
/// @brief Evaluate the degree < values.size() polynomial known at consecutive
/// points 0,1,... in O(n) over a field.
template <class T, class Integer>
T lagrange_consecutive(const std::vector<T> &values, Integer x) {
static_assert(std::is_integral_v<Integer>);
assert(!values.empty());
if constexpr (std::is_signed_v<Integer>) {
if (x >= 0 && std::uint64_t(x) < values.size()) {
return values[std::size_t(x)];
}
} else if (x < values.size()) {
return values[std::size_t(x)];
}
int n = int(values.size());
T point = T(x);
std::vector<T> prefix(n + 1, T(1));
std::vector<T> suffix(n + 1, T(1));
for (int index = 0; index < n; index++) {
prefix[index + 1] = prefix[index] * (point - T(index));
}
for (int index = n - 1; index >= 0; index--) {
suffix[index] = suffix[index + 1] * (point - T(index));
}
std::vector<T> inverse_factorial(n, T(1));
T factorial = T(1);
for (int value = 1; value < n; value++) {
factorial *= T(value);
}
inverse_factorial[n - 1] = T(1) / factorial;
for (int value = n - 1; value >= 1; value--) {
inverse_factorial[value - 1] = inverse_factorial[value] * T(value);
}
T result{};
for (int index = 0; index < n; index++) {
T coefficient = prefix[index] * suffix[index + 1] *
inverse_factorial[index] * inverse_factorial[n - 1 - index];
if ((n - 1 - index) & 1) {
coefficient = -coefficient;
}
result += values[index] * coefficient;
}
return result;
}
/// @brief Return sum_{i=1}^n i^exponent over a field in O(exponent).
template <class T> T power_sum(std::uint64_t n, int exponent) {
assert(exponent >= 0);
auto power = [&](T value, int degree) {
T result = T(1);
while (degree > 0) {
if (degree & 1) {
result *= value;
}
value *= value;
degree >>= 1;
}
return result;
};
std::vector<T> values(exponent + 2);
for (int point = 1; point < int(values.size()); point++) {
values[point] = values[point - 1] + power(T(point), exponent);
}
return lagrange_consecutive(values, n);
}
} // namespace noya
#endif // NOYA_LAGRANGE_INTERPOLATION_HPP
#include <cassert>
#include <cstdint>
#include <type_traits>
#include <vector>
/// @complexity Time: O(n) per evaluation at consecutive points.
/// Space: O(n).
namespace noya {
/// @brief Evaluate the degree < values.size() polynomial known at consecutive
/// points 0,1,... in O(n) over a field.
template <class T, class Integer>
T lagrange_consecutive(const std::vector<T> &values, Integer x) {
static_assert(std::is_integral_v<Integer>);
assert(!values.empty());
if constexpr (std::is_signed_v<Integer>) {
if (x >= 0 && std::uint64_t(x) < values.size()) {
return values[std::size_t(x)];
}
} else if (x < values.size()) {
return values[std::size_t(x)];
}
int n = int(values.size());
T point = T(x);
std::vector<T> prefix(n + 1, T(1));
std::vector<T> suffix(n + 1, T(1));
for (int index = 0; index < n; index++) {
prefix[index + 1] = prefix[index] * (point - T(index));
}
for (int index = n - 1; index >= 0; index--) {
suffix[index] = suffix[index + 1] * (point - T(index));
}
std::vector<T> inverse_factorial(n, T(1));
T factorial = T(1);
for (int value = 1; value < n; value++) {
factorial *= T(value);
}
inverse_factorial[n - 1] = T(1) / factorial;
for (int value = n - 1; value >= 1; value--) {
inverse_factorial[value - 1] = inverse_factorial[value] * T(value);
}
T result{};
for (int index = 0; index < n; index++) {
T coefficient = prefix[index] * suffix[index + 1] *
inverse_factorial[index] * inverse_factorial[n - 1 - index];
if ((n - 1 - index) & 1) {
coefficient = -coefficient;
}
result += values[index] * coefficient;
}
return result;
}
/// @brief Return sum_{i=1}^n i^exponent over a field in O(exponent).
template <class T> T power_sum(std::uint64_t n, int exponent) {
assert(exponent >= 0);
auto power = [&](T value, int degree) {
T result = T(1);
while (degree > 0) {
if (degree & 1) {
result *= value;
}
value *= value;
degree >>= 1;
}
return result;
};
std::vector<T> values(exponent + 2);
for (int point = 1; point < int(values.size()); point++) {
values[point] = values[point - 1] + power(T(point), exponent);
}
return lagrange_consecutive(values, n);
}
} // namespace noya