拉格朗日插值法 在ACM竞赛中的一些应用
什么是拉格朗日插值法
本文转载自 https://riteme.github.io/blog/2017-3-18/lagrange-interpolation.html
大佬写的太好了
当然求自然幂数和有很多种方法 Acdreamer
模板
namespace polysum {
#define rep(i,a,n) for (int i=a;i<n;i++)
#define per(i,a,n) for (int i=n-1;i>=a;i--)
const int D=1e6+10;
ll a[D],f[D],g[D],p[D],p1[D],p2[D],b[D],h[D][2],C[D];
ll powmod(ll a,ll b){ll res=1;a%=mod;assert(b>=0);for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;}
//..........................
// 已知f(0),f(1)..f(d) 求 f(n)
ll calcn(int d,ll *a,ll n) { // a[0].. a[d] a[n]
if (n<=d) return a[n];
p1[0]=p2[0]=1;
rep(i,0,d+1) {
ll t=(n-i+mod)%mod;
p1[i+1]=p1[i]*t%mod;
}
rep(i,0,d+1) {
ll t=(n-d+i+mod)%mod;
p2[i+1]=p2[i]*t%mod;
}
ll ans=0;
rep(i,0,d+1) {
ll t=g[i]*g[d-i]%mod*p1[i]%mod*p2[d-i]%mod*a[i]%mod;
if ((d-i)&1) ans=(ans-t+mod)%mod;
else ans=(ans+t)%mod;
}
return ans;
}
// 初始化，初始化的时候记得将D的值
void init(int M) {
f[0]=f[1]=g[0]=g[1]=1;
rep(i,2,M+5) f[i]=f[i-1]*i%mod;
g[M+4]=powmod(f[M+4],mod-2);
per(i,1,M+4) g[i]=g[i+1]*(i+1)%mod;
}
// 已知 f(0),f(1)...f(m),求 \sum_{i=0 }^{n} f[i]
ll polysum(ll m,ll *a,ll n) { // a[0].. a[m]
ll b[D];
ll b[D];
for(int i=0;i<=m;i++) b[i]=a[i];
b[m+1]=calcn(m,b,m+1);
rep(i,1,m+2) b[i]=(b[i-1]+b[i])%mod;
return calcn(m+1,b,n);// m次多项式的和是m+1 次多项式
}
ll qpolysum(ll R,ll n,ll *a,ll m) {
// a[0].. a[m] \sum_{i=0}^{n-1} a[i]*R^i
if (R==1) return polysum(n,a,m);
a[m+1]=calcn(m,a,m+1);
ll r=powmod(R,mod-2),p3=0,p4=0,c,ans;
h[0][0]=0;h[0][1]=1;
rep(i,1,m+2) {
h[i][0]=(h[i-1][0]+a[i-1])*r%mod;
h[i][1]=h[i-1][1]*r%mod;
}
rep(i,0,m+2) {
ll t=g[i]*g[m+1-i]%mod;
if (i&1) p3=((p3-h[i][0]*t)%mod+mod)%mod,p4=((p4-h[i][1]*t)%mod+mod)%mod;
else p3=(p3+h[i][0]*t)%mod,p4=(p4+h[i][1]*t)%mod;
}
c=powmod(p4,mod-2)*(mod-p3)%mod;
rep(i,0,m+2) h[i][0]=(h[i][0]+h[i][1]*c)%mod;
rep(i,0,m+2) C[i]=h[i][0];
ans=(calcn(m,C,n)*powmod(R,n)-c)%mod;
if (ans<0) ans+=mod;
return ans;
}
} // polysum::init();
例题
1 CF622F
1
LL qpow(LL a,LL b){
LL ans = 1;
while( b > 0){
if(b & 1) ans = ans*a%mod;
a = a*a%mod;
b >>= 1;
}
return ans;
}
namespace polysum {
#define rep(i,a,n) for (int i=a;i<n;i++)
#define per(i,a,n) for (int i=n-1;i>=a;i--)
const int D=1e6+10;
ll a[D],f[D],g[D],p[D],p1[D],p2[D],b[D],h[D][2],C[D];
ll powmod(ll a,ll b){ll res=1;a%=mod;assert(b>=0);for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;}
ll calcn(int d,ll *a,ll n) { // a[0].. a[d] a[n]
if (n<=d) return a[n];
p1[0]=p2[0]=1;
rep(i,0,d+1) {
ll t=(n-i+mod)%mod;
p1[i+1]=p1[i]*t%mod;
}
rep(i,0,d+1) {
ll t=(n-d+i+mod)%mod;
p2[i+1]=p2[i]*t%mod;
}
ll ans=0;
rep(i,0,d+1) {
ll t=g[i]*g[d-i]%mod*p1[i]%mod*p2[d-i]%mod*a[i]%mod;
if ((d-i)&1) ans=(ans-t+mod)%mod;
else ans=(ans+t)%mod;
}
return ans;
}
void init(int M) {
f[0]=f[1]=g[0]=g[1]=1;
rep(i,2,M+5) f[i]=f[i-1]*i%mod;
g[M+4]=powmod(f[M+4],mod-2);
per(i,1,M+4) g[i]=g[i+1]*(i+1)%mod;
}
// a[0].. a[m] \sum_{i=0}^{n-1} a[i]
ll b[D];
for(int i=0;i<=m;i++) b[i]=a[i];
b[m+1]=calcn(m,b,m+1);
rep(i,1,m+2) b[i]=(b[i-1]+b[i])%mod;
return calcn(m+1,b,n-1);
}
ll qpolysum(ll R,ll n,ll *a,ll m) {
// a[0].. a[m] \sum_{i=0}^{n-1} a[i]*R^i
if (R==1) return polysum(n,a,m);
a[m+1]=calcn(m,a,m+1);
ll r=powmod(R,mod-2),p3=0,p4=0,c,ans;
h[0][0]=0;h[0][1]=1;
rep(i,1,m+2) {
h[i][0]=(h[i-1][0]+a[i-1])*r%mod;
h[i][1]=h[i-1][1]*r%mod;
}
rep(i,0,m+2) {
ll t=g[i]*g[m+1-i]%mod;
if (i&1) p3=((p3-h[i][0]*t)%mod+mod)%mod,p4=((p4-h[i][1]*t)%mod+mod)%mod;
else p3=(p3+h[i][0]*t)%mod,p4=(p4+h[i][1]*t)%mod;
}
c=powmod(p4,mod-2)*(mod-p3)%mod;
rep(i,0,m+2) h[i][0]=(h[i][0]+h[i][1]*c)%mod;
rep(i,0,m+2) C[i]=h[i][0];
ans=(calcn(m,C,n)*powmod(R,n)-c)%mod;
if (ans<0) ans+=mod;
return ans;
}
} // polysum::init();
const int maxn = 1e6+100;
LL a[maxn],b[maxn];
int main(void)
{
LL n,k;
cin>>n>>k;
if(k == 0)
return 0*printf("%I64d",n);
polysum::init(k+100);
b[0] = 0;
for(int i = 1; i <= k+1; ++i)
b[i] = qpow(i,k);
LL ans = polysum::polysum(k,b,n+1);
cout<<ans<<endl;
return 0;
}
BZOJ2665
const int maxn = 500+10;
//LL qpow(LL a,LL b){LL s=1;while(b>0){if(b&1)s=s*a%mod;a=a*a%mod;b>>=1;}return s;}
LL gcd(LL a,LL b) {return b?gcd(b,a%b):a;}
LL f[maxn*2][maxn*2];
LL qpow(LL a,LL b,LL m){
LL ans = 1;
while(b > 0){
if(b&1) ans = ans*a%m;
a = a*a%m;
b >>= 1;
}
return ans;
}
int main(void)
{
LL A,n,d;
cin>>A>>n>>d;
me(f);
f[0][0] = 1;
for(int i = 1;i <= 2*n; ++i){
f[i][0] = 1;
for(int j = 1;j <=2 * n; ++j)
f[i][j] = (f[i-1][j-1]*i%d*j+f[i-1][j])%d;
}
LL ans = 0;
if(A <= 2*n)
ans = f[A][n];
else{
for(LL i = 0;i <= 2*n; ++i){
LL tmp1 = 1,tmp2 = 1;
for(LL j = 0;j <= 2*n; ++j){
if(j != i)
tmp1 = tmp1*(A-j)%d,tmp2 = tmp2*(i-j)%d;
}
ans = ans+f[i][n]*(tmp1*qpow(tmp2,d-2,d)%d)%d;
ans = ans%d;
}
}
cout<<(ans%d+d)%d<<endl;
return 0;
}
BZOJ4559
参考博客
https://blog.csdn.net/cdsszjj/article/details/78778488
int R[maxn];
LL U[maxn];
LL polysum[maxn];
LL dp[maxn][maxn],C[maxn][maxn];
LL N,M,K;
LL qpow(LL a,LL b){
LL ans = 1;
while(b > 0){
if(b & 1) ans =ans*a%mod;
a = a*a%mod;
b >>= 1;
}
return ans;
}
LL CC(LL x,LL y){
if( x < 0||y < 0||x < y)
return 0;
return C[x][y];
}
LL poly(LL u, LL r){
LL ans = 0;
int F[maxn];
for(int i = 1;i <= N + 1; ++i){
F[i] = 0;
for(int j = 1;j <= i; ++j){
F[i] = (F[i] + qpow(i-j,r-1)*qpow(j,N-r)%mod)%mod;
}
if(i == u)
return F[i];
}
for(int i = 1;i <= N+1; ++i){
LL tmp1 = 1,tmp2 = 1;
for(int j = 1;j <= N+1; ++j){
if(i != j) tmp1 = tmp1*(u-j)%mod,tmp2 = tmp2*(i-j)%mod;
}
assert(tmp1 != 0&&tmp2 != 0);
ans = (ans + F[i]*tmp1%mod*qpow(tmp2,mod-2)%mod)%mod;
}
return ans;
}
int main(void)
{
C[0][0] = 1;
for(int i = 1;i < maxn; ++i){
C[i][0] = 1;
for(int j = 1;j < maxn; ++j)
C[i][j] = (C[i-1][j-1]+C[i-1][j])%mod;
}
cin>>N>>M>>K;
// cout<<N<<M<<K<<endl;
for(int i = 1;i <= M; ++i)
scanf("%lld",&U[i]);
for(int i = 1;i <= M; ++i)
scanf("%d",&R[i]);
for(int i = 1;i <= M; ++i){
// 拉格朗日算法求解 \sum_{x = 1}{x = U[i]}(U[i]-x)^(R-1)*x^(N-R) = poly[i]
polysum[i] = poly(U[i],R[i]);
}
// for(int i = 1;i <= M; ++i)
// cout<<polysum[i]<<endl;
dp[0][0] = 1;
for(int i = 1;i <= M; ++i){
for(int j = 0;j <= N; ++j){
dp[i][j] = 0;
for(int w = 0;w <= j; ++w){
dp[i][j] = (dp[i][j] + polysum[i]*dp[i-1][w]%mod*CC(w,R[i]-j+w-1)%mod*CC(N-w-1,j-w)%mod)%mod;
}
}
}
cout<<(dp[M][N-K-1]%mod+mod)%mod<<endl;
return 0;
}
腾讯公司福利 465人发布