lagrange_interpolation.hpp¶
已知次数受限多项式在 0,1,… 的值，以线性时间求任意一点的值。
\[
\displaystyle f(x)=\sum_{i=0}^{n-1} y_i\prod_{\substack{0\le j<n\\j\ne i}}\frac{x-j}{i-j}
\]
Complexity: Time: \(O(n)\) per evaluation at consecutive points. power_sum with exponent \(d\) takes \(O(d + 1)\) field operations. Space: \(O(n)\) for interpolation and \(O(d + 1)\) for power_sum.
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(n) per evaluation at consecutive points.
/// power_sum with exponent d takes O(d+1) field operations.
/// Space: O(n) for interpolation and O(d+1) for power_sum.
#include <cassert>
#include <cstdint>
#include <type_traits>
#include <vector>
namespace noya {
namespace lagrange_internal {
// Sieve all powers 0^d,...,(n-1)^d; 0^0 is taken as 1.
template <class T> std::vector<T> monomial_values(int d, int n) {
std::vector<T> vs;
vs.reserve(n + 1);
if (d == 0) {
vs.assign(n, T(1));
return vs;
}
vs.resize(n);
if (n > 1) vs[1] = T(1);
std::vector<int> lp(n), ps;
for (int x = 2; x < n; ++x) {
if (lp[x] == 0) {
lp[x] = x;
ps.push_back(x);
T a = T(x), res = T(1);
for (int k = d; k > 0; k >>= 1) {
if (k & 1) res *= a;
a *= a;
}
vs[x] = res;
}
for (int p : ps) {
if (p > lp[x] || std::int64_t(p) * x >= n) break;
lp[p * x] = p;
// Each composite is reached once through its smallest prime factor.
vs[p * x] = vs[p] * vs[x];
}
}
return vs;
}
} // namespace lagrange_internal
/// @brief Evaluate the degree < vs.size() polynomial known at consecutive
/// points 0,1,... in O(n) over a field. The sample vector must be nonempty;
/// its sample positions must be distinct in the 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));
for (int idx = 0; idx < n; idx++) {
pre[idx + 1] = pre[idx] * (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{}, suf = T(1);
for (int idx = n - 1; idx >= 0; idx--) {
// suf is the product of (poi-j) over j > idx.
T cf = pre[idx] * suf * ifc[idx] * ifc[n - 1 - idx];
if ((n - 1 - idx) & 1) {
cf = -cf;
}
res += vs[idx] * cf;
suf *= poi - T(idx);
}
return res;
}
/// @brief Return sum_{i=1}^n i^exp over a field in O(exp+1)
/// field operations. Integers 1 through exp+1 must be nonzero in the field.
template <class T> T power_sum(std::uint64_t n, int exp) {
assert(exp >= 0);
auto vs = lagrange_internal::monomial_values<T>(exp, exp + 2);
vs[0] = T(0);
for (int i = 1; i < int(vs.size()); ++i) vs[i] += vs[i - 1];
return lagrange_consecutive(vs, n);
}
} // namespace noya
// NOYA_LAGRANGE_INTERPOLATION_HPP
#ifndef NOYA_LAGRANGE_INTERPOLATION_HPP
#define NOYA_LAGRANGE_INTERPOLATION_HPP 1
/// @complexity Time: O(n) per evaluation at consecutive points.
/// power_sum with exponent d takes O(d+1) field operations.
/// Space: O(n) for interpolation and O(d+1) for power_sum.
#include <cassert>
#include <cstdint>
#include <type_traits>
#include <vector>
namespace noya {
namespace lagrange_internal {
// Sieve all powers 0^d,...,(n-1)^d; 0^0 is taken as 1.
template <class T> std::vector<T> monomial_values(int d, int n) {
std::vector<T> vs;
vs.reserve(n + 1);
if (d == 0) {
vs.assign(n, T(1));
return vs;
}
vs.resize(n);
if (n > 1) vs[1] = T(1);
std::vector<int> lp(n), ps;
for (int x = 2; x < n; ++x) {
if (lp[x] == 0) {
lp[x] = x;
ps.push_back(x);
T a = T(x), res = T(1);
for (int k = d; k > 0; k >>= 1) {
if (k & 1) res *= a;
a *= a;
}
vs[x] = res;
}
for (int p : ps) {
if (p > lp[x] || std::int64_t(p) * x >= n) break;
lp[p * x] = p;
// Each composite is reached once through its smallest prime factor.
vs[p * x] = vs[p] * vs[x];
}
}
return vs;
}
} // namespace lagrange_internal
/// @brief Evaluate the degree < vs.size() polynomial known at consecutive
/// points 0,1,... in O(n) over a field. The sample vector must be nonempty;
/// its sample positions must be distinct in the 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));
for (int idx = 0; idx < n; idx++) {
pre[idx + 1] = pre[idx] * (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{}, suf = T(1);
for (int idx = n - 1; idx >= 0; idx--) {
// suf is the product of (poi-j) over j > idx.
T cf = pre[idx] * suf * ifc[idx] * ifc[n - 1 - idx];
if ((n - 1 - idx) & 1) {
cf = -cf;
}
res += vs[idx] * cf;
suf *= poi - T(idx);
}
return res;
}
/// @brief Return sum_{i=1}^n i^exp over a field in O(exp+1)
/// field operations. Integers 1 through exp+1 must be nonzero in the field.
template <class T> T power_sum(std::uint64_t n, int exp) {
assert(exp >= 0);
auto vs = lagrange_internal::monomial_values<T>(exp, exp + 2);
vs[0] = T(0);
for (int i = 1; i < int(vs.size()); ++i) vs[i] += vs[i - 1];
return lagrange_consecutive(vs, n);
}
} // namespace noya
#endif // NOYA_LAGRANGE_INTERPOLATION_HPP
/// @complexity Time: O(n) per evaluation at consecutive points.
/// power_sum with exponent d takes O(d+1) field operations.
/// Space: O(n) for interpolation and O(d+1) for power_sum.
#include <cassert>
#include <cstdint>
#include <type_traits>
#include <vector>
namespace noya {
namespace lagrange_internal {
// Sieve all powers 0^d,...,(n-1)^d; 0^0 is taken as 1.
template <class T> std::vector<T> monomial_values(int d, int n) {
std::vector<T> vs;
vs.reserve(n + 1);
if (d == 0) {
vs.assign(n, T(1));
return vs;
}
vs.resize(n);
if (n > 1) vs[1] = T(1);
std::vector<int> lp(n), ps;
for (int x = 2; x < n; ++x) {
if (lp[x] == 0) {
lp[x] = x;
ps.push_back(x);
T a = T(x), res = T(1);
for (int k = d; k > 0; k >>= 1) {
if (k & 1) res *= a;
a *= a;
}
vs[x] = res;
}
for (int p : ps) {
if (p > lp[x] || std::int64_t(p) * x >= n) break;
lp[p * x] = p;
// Each composite is reached once through its smallest prime factor.
vs[p * x] = vs[p] * vs[x];
}
}
return vs;
}
} // namespace lagrange_internal
/// @brief Evaluate the degree < vs.size() polynomial known at consecutive
/// points 0,1,... in O(n) over a field. The sample vector must be nonempty;
/// its sample positions must be distinct in the 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));
for (int idx = 0; idx < n; idx++) {
pre[idx + 1] = pre[idx] * (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{}, suf = T(1);
for (int idx = n - 1; idx >= 0; idx--) {
// suf is the product of (poi-j) over j > idx.
T cf = pre[idx] * suf * ifc[idx] * ifc[n - 1 - idx];
if ((n - 1 - idx) & 1) {
cf = -cf;
}
res += vs[idx] * cf;
suf *= poi - T(idx);
}
return res;
}
/// @brief Return sum_{i=1}^n i^exp over a field in O(exp+1)
/// field operations. Integers 1 through exp+1 must be nonzero in the field.
template <class T> T power_sum(std::uint64_t n, int exp) {
assert(exp >= 0);
auto vs = lagrange_internal::monomial_values<T>(exp, exp + 2);
vs[0] = T(0);
for (int i = 1; i < int(vs.size()); ++i) vs[i] += vs[i - 1];
return lagrange_consecutive(vs, n);
}
} // namespace noya
// NOYA_LAGRANGE_INTERPOLATION_HPP