文章目录
STEP1 问题引入
现在给定在平面直角坐标系中的 k + 1 k+1 k+1 个点,记 P i ( x i , y i ) P_i(x_i,y_i) Pi(xi,yi)。现在要求一个 k k k 次多项式函数 y = ∑ i = 0 k a i x i y=\sum_{i=0}^ka_ix^i y=∑i=0kaixi,使得它经过这 k + 1 k+1 k+1 个点,即 ∀ i ∣ 0 ≤ i ≤ k , y ∣ x = x i = y i \forall_{i|0\le i\le k},y|_{x=x_i}=y_i ∀i∣0≤i≤k,y∣x=xi=yi。输出 x = n x=n x=n 时 y y y 的取值。
我们首先想到的是把 k + 1 k+1 k+1 个点的坐标分别代入 y = ∑ i = 0 k a i x i y=\sum_{i=0}^ka_ix^i y=∑i=0kaixi,则可以得到 k + 1 k+1 k+1 个方程。于是便可以用高斯消元法求出所有的 a i a_i ai,最后把 n n n 代进去算就好了时间复杂度为 O ( k 3 ) O(k^3) O(k3)。这种方法效率非常低,需要用更高效的方法代替。
STEP2 拉格朗日插值公式
存在性
先考察多项式
l
i
=
∏
j
≠
i
x
−
x
j
x
i
−
x
j
l_i=\prod_{j\not=i}\frac{x-x_j}{x_i-x_j}
li=∏j=ixi−xjx−xj。不难发现,如果
x
=
x
i
x=x_i
x=xi,那么每个分数分子与分母相等,即每个分数都为
1
1
1,则此时
l
i
=
1
l_i=1
li=1;如果
x
≠
x
i
x\not=x_i
x=xi且
x
x
x 为某个点
P
j
P_j
Pj 的横坐标,那么分子有一项为
x
−
x
j
=
0
x-x_j=0
x−xj=0,那么此时
l
i
=
0
l_i=0
li=0。像这样的多项式,我们称之为拉格朗日基本多项式。
那么如果将
l
i
l_i
li 乘上
y
i
y_i
yi ,就可以得到在
x
i
x_i
xi 上取值为
y
i
y_i
yi,在
x
j
x_j
xj (
j
≠
i
j\not=i
j=i)上取值为
0
0
0 的一系列多项式。我们把所有这种多项式累加起来,就可以得到要求的多项式
L
(
x
)
=
∑
i
=
0
k
y
i
l
i
=
∑
i
=
0
k
y
i
∑
0
≤
j
≤
k
,
j
≠
i
x
−
x
j
x
i
−
x
j
L(x)=\sum_{i=0}^ky_il_i=\sum_{i=0}^ky_i\sum_{0\le j\le k,j\not=i}\frac{x-x_j}{x_i-x_j}
L(x)=i=0∑kyili=i=0∑kyi0≤j≤k,j=i∑xi−xjx−xj
时间复杂度是
O
(
k
2
)
O(k^2)
O(k2)。
由于每个点的横坐标不相等,
i
≠
j
i\not=j
i=j,故
x
i
≠
x
j
x_i\not=x_j
xi=xj,即
x
i
−
x
j
≠
0
x_i-x_j\not=0
xi−xj=0。也就是说,上式中的所有分母不等于零,故这样的多项式是存在的。
唯一性
假设有两个不同的 k k k 次多项式 L 1 ≠ L 2 L_1\not=L_2 L1=L2,记 L = L 1 − L 2 L=L_1-L_2 L=L1−L2,则 L L L 的次数不大于 k k k。又 L 1 L_1 L1 和 L 2 L_2 L2 在 x 0 , x 1 , … , x k x_0,x_1,\dots,x_k x0,x1,…,xk 上的取值相同,则 L L L 在这些点上取值都为零,也就是说 L = a ∏ i = 0 k ( x − x i ) L=a\prod_{i=0}^k(x-x_i) L=a∏i=0k(x−xi)。若 a ≠ 0 a\not=0 a=0,则这是一个 k + 1 k+1 k+1次多项式,与它次数不大于 k k k 相矛盾,故 a = 0 a=0 a=0,那么 L = 0 L=0 L=0 即 L 1 = L 2 L_1=L_2 L1=L2。这又与 L 1 ≠ L 2 L_1\not=L_2 L1=L2 相矛盾。故有且只有一个 k k k 次多项式经过平面上的 k + 1 k+1 k+1 个点。
代码模板
#include<iostream>
#include<cstdio>
#define ll long long
using namespace std;
int n;
ll x[2001],y[2001],mod=1e9+7;
ll k,ans,w,l;
inline ll pw(ll a,int b){ //快速幂
if(b==0) return 1ll;
if(b==1) return a;
ll res=pw(a,b>>1);
(res*=res)%=mod;
if(b&1) (res*=a)%=mod;
return res;
}
inline ll sub(ll a,ll b){ //模下的减法
return (a-b>0)?(a-b):(a-b+mod);
}
inline ll inv(ll a){ //乘法逆元
return pw(a,998244351);
}
inline void lagrange(){
for(register int i=0;i<=k;i+=1){
l=1ll; w=1ll;
for(register int j=0;j<=k;j+=1){
if(i==j) continue;
(l*=sub(n,x[j]))%=mod; //计算分子
(w*=sub(x[i],x[j]))%=mod; //计算分母
}
(l*=inv(w))%=mod; //拉格朗日基本多项式的值
(ans+=l*y[i])%=mod; //累加答案
}
return;
}
int main(){
scanf("%d%lld",&n,&k);
for(register int i=0;i<=k;i+=1){ //下标从0开始
scanf("%lld%lld",&x[i],&y[i]);
}
lagrange();
printf("%lld\n",ans);
return 0;
}
STEP3 算法优化
重心拉格朗日插值
定义
l
(
x
)
=
∏
i
=
0
k
(
x
−
x
i
)
l(x)=\prod_{i=0}^k(x-x_i)
l(x)=∏i=0k(x−xi),那么有
l
i
(
x
)
=
l
x
−
x
i
l_i(x)=\frac{l}{x-x_i}
li(x)=x−xil。
原式变成
L
(
x
)
=
∑
i
=
0
k
y
i
l
(
x
)
(
x
−
x
i
)
∏
0
≤
j
≤
k
,
j
≠
i
(
x
i
−
x
j
)
L(x)=\sum_{i=0}^k\frac{y_il(x)}{(x-x_i)\prod_{0\le j\le k,j\not=i}(x_i-x_j)}
L(x)=∑i=0k(x−xi)∏0≤j≤k,j=i(xi−xj)yil(x)。定义重心权
w
i
=
1
∏
0
≤
j
≤
k
,
j
≠
i
(
x
i
−
x
j
)
w_i=\frac{1}{\prod_{0\le j\le k,j\not=i}{(x_i-x_j)}}
wi=∏0≤j≤k,j=i(xi−xj)1,则
L
(
x
)
=
l
(
x
)
∑
i
=
0
k
y
i
w
i
x
−
x
i
L(x)=l(x)\sum_{i=0}^k\frac{y_iw_i}{x-x_i}
L(x)=l(x)i=0∑kx−xiyiwi
这样
l
i
(
x
)
l_i(x)
li(x) 就能以
O
(
1
)
O(1)
O(1) 的时间计算出来,总共耗时
O
(
n
)
O(n)
O(n)。
x x x连续取值的情况
当所有
x
x
x 的取值为连续整数时,基本多项式的分母被
x
i
x_i
xi 分成两段。因此我们与处理好阶乘的乘法逆元,就可以以
O
(
n
)
O(n)
O(n) 的复杂度计算分母了。
上述两种方法结合起来,就能把算法的时间复杂度从
O
(
n
2
)
O(n^2)
O(n2) 变为
O
(
n
)
O(n)
O(n)。
STEP4 经典例题
The Sum of the k-th Powers(CF622F)
There are well-known formulas: ∑ i = 1 n i = 1 + 2 + ⋯ + n = n ( n + 1 ) 2 \sum_{i=1}^ni=1+2+\dots+n=\frac{n(n+1)}{2} ∑i=1ni=1+2+⋯+n=2n(n+1), ∑ i = 1 n i 2 = 1 2 + 2 2 + ⋯ + n 2 = n ( 2 n + 1 ) ( n + 1 ) 6 \sum_{i=1}^ni^2=1^2+2^2+\dots+n^2=\frac{n(2n+1)(n+1)}{6} ∑i=1ni2=12+22+⋯+n2=6n(2n+1)(n+1) , ∑ i = 1 n i 3 = 1 3 + 2 3 + … n 3 = ( n ( n + 1 ) 2 ) 2 \sum_{i=1}^ni^3=1^3+2^3+\dots_n^3=(\frac{n(n+1)}{2})^2 ∑i=1ni3=13+23+…n3=(2n(n+1))2 . Also mathematicians found similar formulas for higher degrees.
Find the value of the sum ∑ i = 1 n i k = 1 k + 2 k + ⋯ + n k \sum_{i=1}^ni^k=1^k+2^k+\dots+n^k ∑i=1nik=1k+2k+⋯+nk modulo 1 0 9 + 7 10^9 + 7 109 + 7 (so you should find the remainder after dividing the answer by the value 1 0 9 + 7 10^9 + 7 109 + 7).
题目大意: 求 ∑ i = 1 n i k ( m o d 1 0 9 + 7 ) \sum_{i=1}^ni^k(mod 10^9+7) ∑i=1nik(mod109+7),其中 1 ≤ n ≤ 1 0 9 1\le n\le 10^9 1≤n≤109, 0 ≤ k ≤ 1 0 6 0\le k\le 10^6 0≤k≤106。
记 F ( x ) = ∑ i = 1 x i k F(x)=\sum_{i=1}^xi^k F(x)=∑i=1xik 跟据题目的提示,这是一个 k + 1 k+1 k+1 次多项式函数。因此只要计算 F ( 1 ) , F ( 2 ) , … , F ( k + 2 ) F(1),F(2),\dots,F(k+2) F(1),F(2),…,F(k+2),再用拉格朗日插值法就能求得 F ( n ) F(n) F(n)。注意到选取的 x x x 为连续整数,可以用到上面提到的两个优化,总时间复杂度为 O ( k ) O(k) O(k)。
#include<iostream>
#include<cstdio>
#define ll long long
using namespace std;
int n,k;
ll y[1000005],jc[1000005],l=1ll,w,ans,mod=1e9+7;
ll qp(ll a,int b){
if(b==0) return 1ll;
if(b==1) return a;
ll res=qp(a,b>>1);
(res*=res)%=mod;
if(b&1) (res*=a)%=mod;
return res;
}
ll sub(ll a,ll b){
return (a>b)?(a-b):(a-b+mod);
}
ll inv(ll a){
return qp(a,mod-2);
}
void lagrange(){
for(int i=1;i<=k+2;i+=1){
(w*=sub(n,i))%=mod;
}
jc[0]=1ll;
for(int i=1;i<=k+2;i+=1){
jc[i]=(jc[i-1]*i*1ll)%mod; //阶乘预处理
}
for(int i=1;i<=k+2;i+=1){ //拉格朗日插值
l=(jc[i-1]*jc[k+2-i])%mod;
if((k+2-i)&1) l=sub(0,w);
l=inv(l);
(l*=inv(n-i))%=mod;
(l*=w)%=mod;
(l*=y[i])%=mod;
(ans+=l)%=mod;
}
return;
}
int main(){
scanf("%d%d",&n,&k);
for(int i=1;i<=k+2;i+=1){
y[i]=qp(i*1ll,k);
(y[i]+=y[i-1])%=mod; //计算y
}
if(n<=k+2){
ans=y[n]; //这种情况直接可以获得答案
}
else lagrange();
printf("%lld\n",ans);
return 0;
}
谢谢观看,记得点赞