lagrange_interpolation.hpp¶
已知次数受限多项式在 0,1,… 的值，以线性时间求任意一点的值。
\[
\displaystyle f(x)=\sum_{i=0}^{n} y_i\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}
\]
Complexity: Time: O(n) per evaluation at consecutive points. Space: O(n).
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @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 < vs.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> &vs, Integer x) {
static_assert(std::is_integral_v<Integer>);
assert(!vs.empty());
if constexpr (std::is_signed_v<Integer>) {
if (x >= 0 && std::uint64_t(x) < vs.size()) {
return vs[std::size_t(x)];
}
} else if (x < vs.size()) {
return vs[std::size_t(x)];
}
int n = int(vs.size());
T poi = T(x);
std::vector<T> pre(n + 1, T(1));
std::vector<T> suf(n + 1, T(1));
for (int idx = 0; idx < n; idx++) {
pre[idx + 1] = pre[idx] * (poi - T(idx));
}
for (int idx = n - 1; idx >= 0; idx--) {
suf[idx] = suf[idx + 1] * (poi - T(idx));
}
std::vector<T> ifc(n, T(1));
T fac = T(1);
for (int val = 1; val < n; val++) {
fac *= T(val);
}
ifc[n - 1] = T(1) / fac;
for (int val = n - 1; val >= 1; val--) {
ifc[val - 1] = ifc[val] * T(val);
}
T res{};
for (int idx = 0; idx < n; idx++) {
T cf = pre[idx] * suf[idx + 1] * ifc[idx] * ifc[n - 1 - idx];
if ((n - 1 - idx) & 1) {
cf = -cf;
}
res += vs[idx] * cf;
}
return res;
}
/// @brief Return sum_{i=1}^n i^exp over a field in O(exp).
template <class T> T power_sum(std::uint64_t n, int exp) {
assert(exp >= 0);
auto pw = [&](T val, int deg) {
T res = T(1);
while (deg > 0) {
if (deg & 1) {
res *= val;
}
val *= val;
deg >>= 1;
}
return res;
};
std::vector<T> vs(exp + 2);
for (int poi = 1; poi < int(vs.size()); poi++) {
vs[poi] = vs[poi - 1] + pw(T(poi), exp);
}
return lagrange_consecutive(vs, n);
}
} // namespace noya
#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 < vs.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> &vs, Integer x) {
static_assert(std::is_integral_v<Integer>);
assert(!vs.empty());
if constexpr (std::is_signed_v<Integer>) {
if (x >= 0 && std::uint64_t(x) < vs.size()) {
return vs[std::size_t(x)];
}
} else if (x < vs.size()) {
return vs[std::size_t(x)];
}
int n = int(vs.size());
T poi = T(x);
std::vector<T> pre(n + 1, T(1));
std::vector<T> suf(n + 1, T(1));
for (int idx = 0; idx < n; idx++) {
pre[idx + 1] = pre[idx] * (poi - T(idx));
}
for (int idx = n - 1; idx >= 0; idx--) {
suf[idx] = suf[idx + 1] * (poi - T(idx));
}
std::vector<T> ifc(n, T(1));
T fac = T(1);
for (int val = 1; val < n; val++) {
fac *= T(val);
}
ifc[n - 1] = T(1) / fac;
for (int val = n - 1; val >= 1; val--) {
ifc[val - 1] = ifc[val] * T(val);
}
T res{};
for (int idx = 0; idx < n; idx++) {
T cf = pre[idx] * suf[idx + 1] * ifc[idx] * ifc[n - 1 - idx];
if ((n - 1 - idx) & 1) {
cf = -cf;
}
res += vs[idx] * cf;
}
return res;
}
/// @brief Return sum_{i=1}^n i^exp over a field in O(exp).
template <class T> T power_sum(std::uint64_t n, int exp) {
assert(exp >= 0);
auto pw = [&](T val, int deg) {
T res = T(1);
while (deg > 0) {
if (deg & 1) {
res *= val;
}
val *= val;
deg >>= 1;
}
return res;
};
std::vector<T> vs(exp + 2);
for (int poi = 1; poi < int(vs.size()); poi++) {
vs[poi] = vs[poi - 1] + pw(T(poi), exp);
}
return lagrange_consecutive(vs, 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 < vs.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> &vs, Integer x) {
static_assert(std::is_integral_v<Integer>);
assert(!vs.empty());
if constexpr (std::is_signed_v<Integer>) {
if (x >= 0 && std::uint64_t(x) < vs.size()) {
return vs[std::size_t(x)];
}
} else if (x < vs.size()) {
return vs[std::size_t(x)];
}
int n = int(vs.size());
T poi = T(x);
std::vector<T> pre(n + 1, T(1));
std::vector<T> suf(n + 1, T(1));
for (int idx = 0; idx < n; idx++) {
pre[idx + 1] = pre[idx] * (poi - T(idx));
}
for (int idx = n - 1; idx >= 0; idx--) {
suf[idx] = suf[idx + 1] * (poi - T(idx));
}
std::vector<T> ifc(n, T(1));
T fac = T(1);
for (int val = 1; val < n; val++) {
fac *= T(val);
}
ifc[n - 1] = T(1) / fac;
for (int val = n - 1; val >= 1; val--) {
ifc[val - 1] = ifc[val] * T(val);
}
T res{};
for (int idx = 0; idx < n; idx++) {
T cf = pre[idx] * suf[idx + 1] * ifc[idx] * ifc[n - 1 - idx];
if ((n - 1 - idx) & 1) {
cf = -cf;
}
res += vs[idx] * cf;
}
return res;
}
/// @brief Return sum_{i=1}^n i^exp over a field in O(exp).
template <class T> T power_sum(std::uint64_t n, int exp) {
assert(exp >= 0);
auto pw = [&](T val, int deg) {
T res = T(1);
while (deg > 0) {
if (deg & 1) {
res *= val;
}
val *= val;
deg >>= 1;
}
return res;
};
std::vector<T> vs(exp + 2);
for (int poi = 1; poi < int(vs.size()); poi++) {
vs[poi] = vs[poi - 1] + pw(T(poi), exp);
}
return lagrange_consecutive(vs, n);
}
} // namespace noya