JZOJ 5746 和

Description

给定$m$和$k$，共$T$次询问，每次询问一个$n$，输出$\sum_{i=1}^ni^k\ mod\ m$
数据保证$m$的最大值质因子不超过$3*10^5$。

Data Constraints

$2\leq n,m,k\leq 10^{18}\ \\ 1\leq T\leq 3*10^3$

Solution

首先这题用到了中国剩余定理，将$m$分解质因数，分解成$m=\Pi_{i=1}^v p_i^{a_i}$的形式。
那对于某一询问，我们只需求出对于每个$i\in [1,v]$，求出$\sum_{j=1}^nj^k\ mod\ p_i^{a_i}$，最后用中国剩余定理组合起来即可。

现在考虑怎么求$\sum_{j=1}^nj^k\ mod\ p_i^{a_i}$，可以观察到${a_i}$不会超过$59$我们可以将每一个数拆成$xp_i+y$的形式，于是我们有

$$\sum_{j=1}^nj^k\ mod\ p_i^{a_i}=\sum_x\sum_{0\leq y<p_i}[x*p_i+y\leq n](xp_i+y)^k\ mod\ p_i^{a_i}$$
为了方便，我们将答案写成
$$\sum_{j=1}^n(x_jp_i+y_j)^k\ mod\ p_i^{a_i}$$
我们用二次项定理将 $k$次方展开，得
$$\sum_{j=1}^n\sum_{v=0}^k C_{k}^v(x_jp_i)^vy_j^{k-v}\ mod\ p_i^{a_i}$$
可以看到，当 $v\geq a_i$时，值在模意义下就变成 $0$了，因此我们又能把式子写成
$$\sum_{j=1}^n\sum_{v=0}^{min(a_i,k)} C_{k}^v(x_jp_i)^vy_j^{k-v}\ mod\ p_i^{a_i}$$
对于每个可能的 $v$的取值分开求解，同时预处理出前缀和数组$H$,
$$H_{i,u,v}=\sum_{j=1}^u j^{k-v}\ mod\ p_i^{a_i}(u\leq p_i)$$
这样我们就能继续化简式子，得
$$\sum_{v=0}^{min(a_i,k)}C_{k}^v[({\lfloor\frac{n}{p_i}\rfloor}p_i)^vH_{i,(n-1)\ mod\ p_i+1,v}+H_{i,\ p_i,v}p_i^v\sum_{j=0}^{\lfloor\frac{n}{p_i}\rfloor-1}j^v]$$

上式中$C_{k}^v \ mod \ p_i^{a_i}$可以通过预处理出$F_v=C_{k}^v \ mod\ m$，再用$F_v \ mod \ p_i^{a_i}$得到。

至于$\sum_{j=0}^{\lfloor\frac{n}{p_i}\rfloor-1}j^v$，考虑$v$并不会超过$59$，所以可以用第二类斯特林数在$O(v^2)$的时间复杂度内求解自然数幂和。

乘法要用快速乘。

时间复杂度我没去算（本蒟蒻太懒了），反正是$O($能过 ) ) <script type="math/tex" id="MathJax-Element-3993">)</script>。

Code

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#pragma GCC optimize(3)

#define fo(i,j,l) for(int i=j;i<=l;++i)
#define fd(i,j,l) for(int i=j;i>=l;--i)

using namespace std;
typedef long long ll;
const ll N=12e6,K=1e3,sj=1e7,NN=1000,P=12e2,W=3e5;

ll x[N],oo,ss[NN],gg[NN];
int t,gs[NN];
ll n,m,k;
ll kt[120][305000];
int be[K],en[K];

ll S[P][P],c[P];
ll jj[P];

inline ll mod(ll x,ll y,ll m)
{
    ll tmp=(ll)((long double)x*y/m+1e-8)*m;
    return (x*y-tmp+m)%m;
}

inline int min(int a,ll b)
{return a<b?a:b;}

inline ll mmo(ll p,ll mo)
{return (p>=mo)?p-mo:p;}

inline ll mmo_m(ll p)
{return (p>=0)?p:(p+m);}

inline ll ksm(ll o,ll t,ll m)
{
    ll y=1; o%=m;
    for(;t;t>>=1,o=mod(o,o,m))
    if(t&1)y=mod(y,o,m);
    return y;
}

inline ll ksm_mod(ll o,ll t,ll mo)
{
    ll y=1;
    for(;t;t>>=1,o=o*o%mo)
    if(t&1)y=y*o%mo;
    return y;
}

inline ll iabs(ll o)
{return (o>0)?o:(-o);}

inline ll gcd(ll a,ll b)
{return (b==0)?a:gcd(b,a%b);}

inline void prepare()
{
    S[1][1]=1; S[1][0]=0;
    fo(i,2,K){
        S[i][1]=1; S[i][0]=0;
        fo(l,2,i)S[i][l]=mmo(S[i-1][l-1]+mod(S[i-1][l],l,m),m);
    }
}

inline ll js(ll p,int k,ll m)
{
    if(k==0)return (p+1)%m;
    ll ans=0;
    fo(i,1,k){
        ll lj=1;
        fo(l,-1,i-1)if((p-l)%(i+1)==0)lj=mod(lj,(p-l)/(i+1),m);
        else lj=mod(lj,(p-l),m);
        ans=mmo(ans+mod(S[k][i],lj,m),m);
    }
    return ans;
}

inline ll C(ll a,ll b,ll mo)
{
    fo(i,1,b)jj[i]=a-i+1;
    fo(i,2,b){
        ll dq=i;
        fo(i,1,b)if(gcd(jj[i],dq)!=1){
            ll up=gcd(jj[i],dq);
            dq/=up;
            jj[i]/=up;
            if(dq==1)break;
        }
    }
    ll op=1;
    fo(i,1,b)op=mod(op,jj[i],mo);
    return op;
}

inline ll countt(ll x,ll len,int gg,int bh)
{
    ll num=(x/len)%len;
    ll mo=1,op=0;
    fo(i,1,gg)mo=mo*len;
    ll kk=(x/len); ll w=1;
    fo(i,0,min(gg,k))
    {
        ll lj=mod(kt[be[bh]+i][x%len],w,mo);
        lj=mod(lj,c[i]%mo,mo);
        op=mmo(op+lj,mo);
        w=mod(w,len,mo);
        w=mod(w,kk,mo);
    }
    w=1;
    if(kk)
    fo(i,0,min(gg,k)){
        ll lj=js(kk-1,i,mo);
        lj=mod(lj,w,mo);
        lj=mod(lj,kt[be[bh]+i][len],mo);
        lj=mod(lj,c[i]%mo,mo);
        op=mmo(op+lj,mo);
        w=mod(w,len,mo);
    }
    return op;
}

int main()
{
    cin>>m>>k>>t;
    ll op=m; prepare();
    fo(i,0,60)c[i]=C(k,i,m);
    fo(i,2,W)if(op%i==0){
        ss[++oo]=i;
        while(op%i==0)++gs[oo],op/=i;
    }
    sort(ss+1,ss+oo+1); int ok=1;
    fo(i,1,oo){
        ll mo=1;
        be[i]=en[i-1]+1; en[i]=be[i]+gs[i];
        fo(l,1,gs[i])mo=mo*ss[i];
        fo(j,0,min(gs[i],k))
        fo(l,1,ss[i])
        kt[j+be[i]][l]=mmo(kt[j+be[i]][l-1]+ksm(l,k-j,mo),mo);
    }
    fo(i,1,t){
        scanf("%lld",&n);
        ll ans=0;
        fo(i,1,oo){
            ll p=countt(n,ss[i],gs[i],i);
            ll p2=m,mo=1;
            fo(j,1,gs[i])mo=mo*ss[i],p2=p2/ss[i];
            ll p1=ksm(p2%mo,(mo/ss[i])*(ss[i]-1)-1,mo);
            p=mod(p,p1,m); p=mod(p,p2,m);
            ans=mmo(ans+p,m);
        }
        printf("%lld\n",ans);
    }
}