适用范围:给出n个点$x_i,y_i)$,过这n个点能够确定一个最高n-1次的多项式$fx)$,求$fk)$
做法:如图所示,我们将每一个点$x_i,y_i)$在x轴上的投影$x_i,0)$记为$H_i$。对每一个i,我们选择一个点集$lbrace P_ibrace cup lbrace H_j vert 1 le ile n, j
eq ibrace$,求过这n个点的最高n-1次的多项式$g_ix)$。这样我们就得到了n个$g_ix)$,它们都在各自对应的$x_i$处的值为$y_i$,并且在其它$x_ji
eq j)$处值为0。
很容易就能够构造出$g_ix)$的表达式:
$$g_ix)=y_i*prod_{j
eq i}frac{x-x_j}{x_i-x_j}$$
显然最后有:
$$fx)=sum_{i=1}^{n}g_ix)=sum _{i=1}^{n}y_i*prod_{j
eq i}frac{x-x_j}{x_i-x_j}$$
由于只用求$fk)$的值,代入得$fk)=sum _{i=1}^{n}y_i*prod_{j
eq i}frac{k-x_j}{x_i-x_j}$
例题:
#include <iostream> #include <algorithm> #include <cstring> #include <cstdio> using namespace std; typedef long long ll; const ll mod = 998244353; const int N = 2010; int n; ll k, x[N], y[N]; ll powerll a, ll n) { ll res = 1; while n) { if n & 1) res = res * a % mod; a = a * a % mod; n >>= 1; } return res; } int main) { scanf"%d%lld", &n, &k); for int i = 1; i <= n; i++) scanf"%lld%lld", &x[i], &y[i]); ll res = 0; for int i = 1; i <= n; i++) { ll ta = y[i] % mod, tb = 1; for int h = 1; h <= n; h++) { if h == i) continue; ta = ta * k - x[h]) % mod + mod) % mod) % mod; tb = tb * x[i] - x[h]) % mod + mod) % mod) % mod; } res = res + ta * powertb, mod - 2) % mod) % mod; } printf"%lld ", res); return 0; }
[模板]拉格朗日插值
Codeforces – 622F The Sum of the k-th Powers
题意:求$sum_{i=1}^{n}i^k,1leq nleq 10^9,0leq k leq 10^6$
思路:首先有一个结论,$sum_{i=1}^{n}i^k$为k+1阶多项式,所以我们只需要暴力算出$fn)=sum_{i=1}^{n}i^k$的前k+2项,然后用拉格朗日插值法求第n项即可
下面简单证明一下$sum_{i=1}^{n}i^k$为k+1阶多项式
对于一个数列${a_n}$来说,把数列${a_n}$的元素两两做差得到数列${b_n}$,我们称数列${b_n}$为数列${a_n}$的一阶阶差数列,如果再将数列${b_n}$的元素两两做差得到数列${c_n}$,我们称数列${c_n}$为数列${a_n}$的一阶阶差数列,依次类推定义p阶阶差数列
定理:数列${a_n}$是一个p阶等差数列的充要条件是数列的通项$a_n$为n的一个p次多项式
证明:设$fx)=sum_{i=0}^{n}u_i x^i$,令${b_n}$为${a_n}$的一阶阶差数列
那么$Delta fx)=fx+1)-fx)=sum_{i=0}^{n}u_i x+1)^i-sum_{i=0}^{n}u_i x^i$
我们只考虑$x^n$有$Delta fx)=u_nx+1)^n-u_n x^n$
将$u_nx+1)^n$二项式展开,仅考虑$x^n$有$Delta fx)=u_n x^n-u_n x^n=0$
所以每次差分后多项式的最高次数会减1,p次差分后,变为常数项,此时多项式的最高次数为0,定理成立
显然,如果一个多项式k次差分之后变为0,那么这个多项式的最高次数为k-1
我们将$sum_{i=1}^{n}i^k$差分后得$1^k,2^k,3^kdots n^k$
显然$1^k,2^k,3^kdots n^k$为k阶多项式,所以$sum_{i=1}^{n}i^k$为k+1阶多项式
#include <iostream> #include <algorithm> #include <cstring> #include <cstdio> using namespace std; typedef long long ll; const int N = 2000010; const ll mod = 1000000007; int n, k; ll a[N], f[N], nf[N], now; ll powerll a, ll n) { ll res = 1; while n) { if n & 1) res = res * a % mod; a = a * a % mod; n >>= 1; } return res; } ll solve) { now = f[0] = nf[0] = 1; for int i = 1; i <= k + 2; i++) { now = now * n - i) % mod; f[i] = f[i - 1] * i % mod; nf[i] = -nf[i - 1] * i % mod; } ll res = 0; for int i = 1; i <= k + 2; i++) { ll t = powerf[i - 1] * nf[k + 2 - i] % mod, mod - 2); res = res + a[i] * now % mod * powern - i, mod - 2) % mod * t % mod) % mod; res = res + mod) % mod; } return res; } int main) { // freopen"in.txt", "r", stdin); // freopen"out.txt", "w", stdout); scanf"%d%d", &n, &k); for int i = 1; i <= mink + 2, n); i++) a[i] = a[i - 1] + poweri, k)) % mod; if n <= k + 2) printf"%lld ", a[n]); else printf"%lld ", solve)); return 0; }
The Sum of the k-th Powers
题意:求$sum_{i=0}^m fn-a_i)-sum_{i=0}^{m-1} sum_{j=i+1}^{m}a_j-a_i)^m+1$,其中$fn)=sum_{i=1}^{n} i^{m+1}$
思路:后半部分$sum_{i=0}^{m-1} sum_{j=i+1}^{m}a_j-a_i)^m+1$可以直接暴力求解
$fn)=sum_{i=1}^{n} i^{m+1}$为m+2阶多项式,考虑到m不是很大,所以我们可以插入m+3个值后,暴力对每一个$a_i$求一次fn-a_i),最后相加即可
#include <iostream> #include <algorithm> #include <cstring> #include <cstdio> using namespace std; typedef long long ll; const int N = 60; const ll mod = 1000000007; int T; ll n, m, now; ll a[N], c[N], f[N], nf[N]; ll powerll a, ll n) { ll res = 1; while n) { if n & 1) res = res * a % mod; a = a * a % mod; n >>= 1; } return res; } ll solvell n) { now = f[0] = nf[0] = 1; for int i = 1; i <= m + 3; i++) { now = now * n - i) % mod; f[i] = f[i - 1] * i % mod; nf[i] = -nf[i - 1] * i % mod; } ll res = 0; for int i = 1; i <= m + 3; i++) { ll t = powerf[i - 1] * nf[m + 3 - i] % mod, mod - 2); res = res + c[i] * now % mod * powern - i, mod - 2) % mod * t % mod) % mod; res = res + mod) % mod; } return res; } ll calcll x) { if x <= m + 3) return c[x]; return solvex); } int main) { // freopen"in.txt", "r", stdin); // freopen"out.txt", "w", stdout); scanf"%d", &T); while T--) { scanf"%lld%lld", &n, &m); for int i = 1; i <= m; i++) scanf"%lld", &a[i]); sorta + 1, a + m + 1); for int i = 1; i <= minm + 3, n); i++) c[i] = c[i - 1] + poweri, m + 1)) % mod; ll res = 0; for int i = 0; i <= m; i++) res = res + calcn - a[i])) % mod; for int i = 0; i <= m - 1; i++) { for int h = i + 1; h <= m; h++) { ll t = powera[h] - a[i], m + 1); res = res - t) % mod + mod) % mod; } } printf"%lld ", res); } return 0; }
[TJOI2018]教科书般的亵渎
题意:求$sum_{i=0}^n sum_{j=1}^{a+i*d} sum_{x=1}^{j}x^k$
思路:$fn)=sum_{x=1}^{n}x^k$,为k+1阶多项式
$gn)=sum_{i=1}^{n} sum_{j=1}^{i} j^k$,gn)差分后的结果为fn),所以gn)为k+2阶多项式
$resn)=sum_{i=0}^{n}ga+i*d)$,resn)进行k+5次差分后为0,所以resn)为k+4阶多项式
然后插值求解即可
#include <iostream> #include <algorithm> #include <cstring> #include <cstdio> using namespace std; typedef long long ll; const int N = 200; const ll mod = 1234567891; int T, k; ll a, n, d, now, fc[N], nfc[N]; ll f[N], g[N]; ll powerll a, ll n) { ll res = 1; while n) { if n & 1) res = res * a % mod; a = a * a % mod; n >>= 1; } return res; } ll solvell n, int m, ll *c) { now = fc[0] = nfc[0] = 1; for int i = 1; i <= m; i++) { now = now * n - i) % mod; fc[i] = fc[i - 1] * i % mod; nfc[i] = -nfc[i - 1] * i % mod; } ll res = 0; for int i = 1; i <= m; i++) { ll t = powerfc[i - 1] * nfc[m - i] % mod, mod - 2); res = res + c[i] * now % mod * powern - i, mod - 2) % mod * t % mod) % mod; res = res + mod) % mod; } return res; } ll calcll n, int m, ll *c) { if n <= m) return c[n]; return solven, m, c); } int main) { // freopen"in.txt", "r", stdin); // freopen"out.txt", "w", stdout); scanf"%d", &T); while T--) { scanf"%d%lld%lld%lld", &k, &a, &n, &d); for int i = 0; i <= k + 3; i++) f[i] = f[i - 1] + poweri, k)) % mod; for int i = 1; i <= k + 3; i++) f[i] = f[i] + f[i - 1]) % mod; for int i = 0; i <= k + 5; i++) g[i] = calca + i * d) % mod, k + 3, f); for int i = 1; i <= k + 5; i++) g[i] = g[i] + g[i - 1]) % mod; printf"%lld ", calcn, k + 5, g)); } return 0; }
tyvj 1858 XLkxc