[CodeForces622F]The Sum of the k-th Powers

题意:这是一个模板题 ni=1im ∑ i = 1 n i m

可以知道这个解显然是一个 d=m+1 d = m + 1 次多项式

考虑拉格朗日插值求这个多项式的解

假设有 d+1 d + 1 组点值 (x0,y0),(x1,y1),,(xd,yd) ( x 0 , y 0 ) , ( x 1 , y 1 ) , … , ( x d , y d )

那么这个多项式可以表示成 f(n)=di=0li(n)yi f ( n ) = ∑ i = 0 d l i ( n ) y i

li(n)=dj=0,jinxjxixj l i ( n ) = ∏ j = 0 , j ≠ i d n − x j x i − x j

一般的我们取 xi=i,yi=f(i) x i = i , y i = f ( i )

li=dj=0,jinjij ⇒ l i = ∏ j = 0 , j ≠ i d n − j i − j

=n(n1)(nd)(ni)(i1j=0ij)(dj=i+1ji)(1)di = n ( n − 1 ) … ( n − d ) ( n − i ) ( ∏ j = 0 i − 1 i − j ) ( ∏ j = i + 1 d j − i ) ( − 1 ) d − i

=(1)din(n1)(nd)(ni)i!(di)! = ( − 1 ) d − i n ( n − 1 ) … ( n − d ) ( n − i ) i ! ( d − i ) !

=(1)din(n1)(ni+1)(ni1)(nd)i!(di)! = ( − 1 ) d − i n ( n − 1 ) … ( n − i + 1 ) ( n − i − 1 ) … ( n − d ) i ! ( d − i ) !

f(n)=di=0(1)dif(i)n(n1)(ni+1)(ni1)(nd)i!(di)! ⇒ f ( n ) = ∑ i = 0 d ( − 1 ) d − i f ( i ) n ( n − 1 ) … ( n − i + 1 ) ( n − i − 1 ) … ( n − d ) i ! ( d − i ) !

考虑到 g(x)=xm g ( x ) = x m 是积性函数,可以在 O(mlnmlog(P2)) O ( m ln ⁡ m log ⁡ ( P − 2 ) ) 的时间预处理

然后预处理阶乘逆元,前缀积和后缀积就可以在 O(m) O ( m ) 的时间内完成插值,或者你嫌麻烦也可以 O(mlogm) O ( m log ⁡ m )

#include<bits/stdc++.h>
#define fp(i,a,b) for(int i=a,I=b;i<=I;++i)
#define fd(i,a,b) for(int i=a,I=b;i>=I;--i)
#define file(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
using namespace std;
const int N=1e6+5,P=1e9+7;
typedef int arr[N];
typedef long long ll;
int n,m,d;arr is,pr,g,f,pre,suf,fac,ifac;ll ans;
inline int fpm(int a,int b){int x=1;for(;b;b>>=1,a=1ll*a*a%P)if(b&1)x=1ll*a*x%P;return x;}
int main(){
    #ifndef ONLINE_JUDGE
        file("s");
    #endif
    scanf("%d%d",&n,&m);d=m+1;
    g[1]=fac[0]=1;
    fp(i,2,d){
        if(!is[i])pr[++pr[0]]=i,g[i]=fpm(i,m);
        for(int j=1,x;j<=pr[0]&&(x=i*pr[j])<=d;++j){
            g[x]=1ll*g[i]*g[pr[j]]%P;is[x]=1;
            if(i%pr[j]==0)break;
        }
    }
    fp(i,1,d)f[i]=(f[i-1]+g[i])%P;
    if(n<=d)return printf("%d",f[n]),0;
    fp(i,1,d)fac[i]=1ll*fac[i-1]*i%P;ifac[d]=fpm(fac[d],P-2);
    fd(i,d,1)ifac[i-1]=1ll*ifac[i]*i%P;
    pre[0]=1;fp(i,1,d)pre[i]=1ll*pre[i-1]*(n-i+1)%P;
    suf[d]=1;fd(i,d,1)suf[i-1]=1ll*suf[i]*(n-i)%P;
    fp(i,0,d)
        ans+=1ll*((d-i)&1?-1:1)*f[i]*pre[i]%P*suf[i]%P*ifac[d-i]%P*ifac[i]%P;
    printf("%lld",(ans%P+P)%P);
return 0;
}

其实我们还可以 O(mlogm) O ( m log ⁡ m ) 多项式求逆求出伯努利级数,然后 O(m) O ( m ) 扫一遍

不过伯努利级数貌似比较局限(其实就是蒟蒻码力有限)

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值