前言
在数值分析中,拉格朗日插值法是一种多项式插值方法。
拉格朗日插值法
下文,本蒟蒻可能有说的不妥甚至不对的地方, 欢迎大神来打脸。
定义多项式函数
f(x)
f
(
x
)
,则点值运算就是给定
x
x
,求
插值是点值的逆运算,就是给定
n
n
次多项式的个点值表达
(xi,yi)
(
x
i
,
y
i
)
,要求出
f(x)
f
(
x
)
各次项系数
插值法有许多种,当然有名的有拉格朗日插值法与牛顿插值法。
有一个结论,n次的函数可以由n+1个点确定。例如两个点确定一次函数;三个点确定二次函数。因为待定系数包括常数项有n+1项,所以要有n+1个方程来解。
拉格朗日插值基函数
对于每个
(xk,yk)
(
x
k
,
y
k
)
,定义
这个函数有一个好处,就是 ∀i≤n ∀ i ≤ n ,当且仅当 i=k i = k 时, lk(xi)=1 l k ( x i ) = 1 ,否则 lk(xi)=0 l k ( x i ) = 0
于是把每个点值处的基函数求和
就能符合每个点值了, L(x) L ( x ) 即为所求多项式函数
应用
自然数幂和
引理: fk(n)= f k ( n ) = ∑ni=1ik ∑ i = 1 n i k 必定能用 k+1 k + 1 次多项式函数表示
我们只要取
k+2
k
+
2
个点值表达,即可插值出这个多项式函数
为了简便直接取
(0,fk(0)),(1,fk(1))...(k+1,fk(k+1))
(
0
,
f
k
(
0
)
)
,
(
1
,
f
k
(
1
)
)
.
.
.
(
k
+
1
,
f
k
(
k
+
1
)
)
代入
L(x)
L
(
x
)
,得
但是它的局限性在于,需要预处理出阶乘的逆元,因此对模数有要求
时间复杂度可以做到 O(k) O ( k )
例题 51nod1258 Code
#include<cstdio>
#include<cstring>
#include<algorithm>
#define fo(i,a,b) for(ll i=a;i<=b;i++)
#define fd(i,b,a) for(ll i=b;i>=a;i--)
#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)<(y)?(x):(y))
#define mset(a,x) memset(a,x,sizeof(a))
using namespace std;
typedef long long ll;
const ll mo=1e9+7,K=5e4+10;
ll k,fk[K],fac[K],inv[K],pre[K],suf[K],pr[K];
bool bz[K];
ll qmi(ll x,ll n)
{
ll t=1;
for(x%=mo;n;n>>=1,x=x*x%mo)
if(n&1) t=t*x%mo;
return t;
}
void sieve(ll n)
{
mset(fk,0);mset(inv,0);mset(fac,0);
fk[0]=0,fk[1]=fac[0]=inv[0]=fac[1]=inv[1]=1;
pr[0]=0;
fo(i,2,n)
{
fac[i]=fac[i-1]*i%mo,fk[i]=(fk[i]+fk[i-1])%mo;
if(!bz[i]) pr[++pr[0]]=i,fk[i]=(fk[i-1]+qmi(i,k))%mo,inv[i]=qmi(i,mo-2);
inv[i]=inv[i]*inv[i-1]%mo;
fo(j,1,pr[0])
{
ll x=i*pr[j];
if(x>n) break;
bz[x]=1;
inv[x]=fac[i-1]*inv[i]%mo*fac[pr[j]-1]%mo*inv[pr[j]]%mo;
fk[x]=(fk[x]+(fk[i]-fk[i-1]+mo)*(fk[pr[j]]-fk[pr[j]-1]+mo)%mo)%mo;
if(i%pr[j]==0) break;
}
}
}
ll calc(ll n)
{
if(!n) return 0;
if(!k) return n;
if(n<=k+1) return fk[n];
ll ans=0;
int t=0;
pre[0]=1;
fo(i,n-k-1,n) ++t,pre[t]=pre[t-1]*(i%mo)%mo;
t=0;
suf[0]=1;
fd(i,n,n-k-1) ++t,suf[t]=suf[t-1]*(i%mo)%mo;
fo(j,0,k+1)
{
ll t=fk[j]*pre[k-j+1]%mo*suf[j]%mo*inv[j]%mo*inv[k+1-j]%mo;
if((k+1-j)&1) ans=(ans-t+mo)%mo;
else ans=(ans+t)%mo;
}
return ans;
}
int main()
{
int T;ll n;
scanf("%d",&T);
while(T--)
{
scanf("%lld %lld",&n,&k);
sieve(k+1);
printf("%lld\n",calc(n));
}
return 0;
}