【JZOJ 5746】一道比较强的 自然数幂和 板题

Description

给定$m,k$，共$T$次询问，每次输入一个$n$，求$\sum_{i=1}^ni^k$在 $\bmod m$意义下的值
$m$的最大质因子$\leq 3*10^5$
$2\leq n,m,k\leq 10^{18},1\leq T\leq 3*10^3$

Analysis

首先分解质因数$m=\prod p_i^{e_i}$，考虑对每个不同质因子分别求出$\bmod p^e$意义下的值再用中国剩余定理求答案
将每一个$i^k$拆成$(ip+j)^k,i\leq \lfloor n/p\rfloor,j<p$，用二项式定理展开以后发现某一项里面如果有$p^x,x\geq e$那么显然这一项值为$0$，通过这个化简式子，综合运用斯特林数求自然数幂和等数学知识，再做些适当（变态）的预处理，就能以$O(Te^2+pe)$的复杂度解决此题了
一个Trick：涉及求组合数时可能出现没有逆元的情况，所以需要用质因数分解来求

Codes

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
#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;
ll gcd(ll a,ll b){return a%b==0?b:gcd(b,a%b);}
ll n,m,mo,k;
ll qmi(ll x,ll n)
{
    ll t=1;
    for(x%=mo;n;n>>=1,x=(__int128)x*x%mo)
        if(n&1) t=(__int128)t*x%mo;
    return t;
}

const int N=3e5+5,E=66;
int now;
ll mxp,mxe,num,a[20],b[E];
vector<ll> Ans[20][N];
ll S[E][E];

void exgcd(ll a,ll b,ll &x,ll &y)
{
    if(!b) {x=1,y=0;return;}
    ll xx,yy;
    exgcd(b,a%b,xx,yy);
    x=yy,y=xx-a/b*yy;
}
ll inv(ll a)
{
    ll x,y;
    exgcd(a,mo,x,y);
    x=(x%mo+mo)%mo;
    return x;
}

ll fz[E],pro[E][E];
ll calc(ll p,ll e)
{
    e=min(e,k);
    S[0][0]=1;
    fo(i,1,e)
        fo(j,1,i) S[i][j]=(S[i-1][j-1]+(__int128)S[i-1][j]*j)%mo;

    ll lim=n/p,r=n%p,ans=0;
    for(ll l=0,_c=1,limp=1;l<=e;++l)
    {
        if(l)
        {
            fz[l]=k-l+1;
            ll _l=l;
            fo(i,1,l) if(_l>1)
            {
                ll d=gcd(fz[i],_l);
                fz[i]/=d,_l/=d;
            }
            _c=1;
            fo(i,1,l) _c=(__int128)_c*fz[i]%mo;
            limp=(__int128)limp*lim%mo*p%mo;
        }
        ll t=(__int128)_c*limp%mo*Ans[now][r][l]%mo;
        ans=(ans+t)%mo;
    }
    if(r==0) ans=qmi((__int128)lim*p,k);

    --lim;
    fo(i,0,e+2)
        fo(j,0,e+2) pro[i][j]=1;
    fo(i,lim-e+1,lim+1)//lim-e+1->1
        for(ll j=i;j<=lim+1;++j)
            pro[i-lim+e][j-lim+e]=(__int128)pro[i-lim+e][j-1-lim+e]*j%mo;
    for(ll l=0,_c=1;l<=e;++l)
    {
        if(l)
        {
            fz[l]=k-l+1;
            ll _l=l;
            fo(i,1,l) if(_l>1)
            {
                ll d=gcd(fz[i],_l);
                fz[i]/=d,_l/=d;
            }
            _c=1;
            fo(i,1,l) _c=(__int128)_c*fz[i]%mo;
        }
        ll x=0;
        fo(j,0,l)
        {
            ll t=S[l][j];
            ll pos=(lim+1)/(j+1)*(j+1);
            t=(__int128)t*pro[(lim-j+1)-lim+e][(pos-1)-lim+e]%mo*pro[(pos+1)-lim+e][(lim+1)-lim+e]%mo;
            t=(__int128)t*(pos/(j+1))%mo;

            x=(x+t)%mo;
        }
        x=(__int128)x*qmi(p,l)%mo*_c%mo*Ans[now][p-1][l]%mo;
        ans=(ans+x)%mo;
    }
    return ans;
}

int main()
{
    freopen("sum.in","r",stdin);
    freopen("sum.out","w",stdout);
    int T;
    scanf("%lld %lld %d",&m,&k,&T);ll _m=m;
    fo(i,2,3e5)
        if(_m%i==0)
            for(a[++num]=i;_m%i==0;_m/=i) ++b[num];//num<=14
    fo(i,1,num) mxp=max(mxp,a[i]),mxe=max(mxe,b[i]);

    fo(id,1,num)
    {
        mo=1;
        fo(j,1,b[id]) mo*=a[id];
        ll p=a[id],e=min(k,b[id]);

        fd(j,e,1) Ans[id][0].push_back(0);
        Ans[id][0].push_back(1);
        fo(i,1,p)
        {
            fo(j,0,e) Ans[id][i].push_back(0);
            for(ll j=e,t=qmi(i,k-e);j>=0;--j,t=(__int128)t*i%mo)
                Ans[id][i][j]=(Ans[id][i-1][j]+t)%mo;
        }
    }

    for(;T;T--)
    {
        scanf("%lld",&n);
        ll ans=0;
        fo(i,1,num)
        {
            now=i;
            mo=1;
            fo(j,1,b[i]) mo*=a[i];
            ll ti=calc(a[i],b[i]);
            ll Mi=m/mo;
            ll pi=inv(Mi);
            ll s=(__int128)ti*Mi%m*pi%m;
            ans=(ans+s)%m;
        }
        printf("%lld\n",ans);
    }
    return 0;
}