JZOJ 5746 和

Description

給定m$m$

和k$k$，共T$T$次詢問，每次詢問一個n$n$，輸出∑ni=1ik mod m$\sum _{i=1}^n{i}^{k}modm$
資料保證m$m$的最大值質因子不超過3∗105$3\ast {10}^5$。

Data Constraints

2≤n,m,k≤1018 1≤T≤3∗103$$\begin{gathered} 2\le n,m,k\le {10}^{18} \\ 1\le T\le 3\ast {10}^3 \end{gathered}$$

Solution

首先這題用到了中國剩餘定理，將m$m$

分解質因數，分解成m=Πvi=1paii$m={\mathrm{\Pi }}_{i=1}^v{p}_i^{a_i}$的形式。
那對於某一詢問，我們只需求出對於每個i∈[1,v]$i\in [1,v]$，求出∑nj=1jk mod paii$\sum _{j=1}^n{j}^{k}mod{p}_i^{a_i}$，最後用中國剩餘定理組合起來即可。

現在考慮怎麼求∑nj=1jk mod paii$\sum _{j=1}^n{j}^{k}mod{p}_i^{a_i}$

，可以觀察到ai$a_i$不會超過59$59$我們可以將每一個數拆成xpi+y$x{p}_i+y$的形式，於是我們有

為了方便，我們將答案寫成

我們用二次項定理將k$k$次方展開，得

可以看到，當v≥ai$v\ge a_i$時，值在模意義下就變成0$0$了，因此我們又能把式子寫成

對於每個可能的v$v$的取值分開求解，同時預處理出字首和陣列H$H$,

這樣我們就能繼續化簡式子，得

上式中Cvk mod paii$C_k^{v}mod{p}_i^{a_i}$

可以通過預處理出Fv=Cvk mod m$F_v=C_k^{v}modm$，再用Fv mod paii$F_{v}mod{p}_i^{a_i}$得到。

至於∑⌊npi⌋−1j=0jv$\sum _{j=0}^{\lfloor \frac{n}{p_i}\rfloor -1}{j}^v$

，考慮v$v$並不會超過59$59$，所以可以用第二類斯特林數在O(v2)$O(v^2)$的時間複雜度內求解自然數冪和。

乘法要用快速乘。

時間複雜度我沒去算（本蒟蒻太懶了），反正是O($O($

能過)$)$。

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);
    }
}