[JZOJ5746]和 自然数幂和+中国剩余定理

首先考虑计算模质数下的自然数幂和，通过stirling数转化成下降幂，

$$\sum_{i=0}^n i^k=\sum_{i=0}^n\sum_{j=0}^k \begin{Bmatrix} k\\j \end{Bmatrix} i^{\underline j} =\sum_{j=0}^k\begin{Bmatrix} k\\j \end{Bmatrix}\sum_{i=0}^n i^{\underline j}$$
然后考虑对后面那个下降幂的和进行积分，我们注意到 $(x+1)^{\underline k}-x^{\underline k}=kx^{\underline {k-1}}$，变一下就得到 $x^{\underline {k}}=\frac{(x+1)^{\underline {k+1}}-x^{\underline {k+1}}}{k+1}$，设 $g(x)=\frac{x^{\underline {k+1}}}{k+1}$，那么 $\Delta g(x)=g(x+1)-g(x)=x^{\underline {k}}$，那么和式就是值就是 $g(n+1)$了（和普通幂函数的求导好像形式一样），
$$=\sum_{j=0}^k \begin{Bmatrix} k\\j \end{Bmatrix} \frac{(n+1)^{\underline {j+1}}}{j+1}$$
于是可以 $O(k^2)$解决问题。自然数幂和还有更优的做法，但这题有这个就可以解决。

将模数$m$分解质因数$=p_1^{\alpha _1}p_2^{\alpha _2}... (p_i\le 300000)$，假设所有$\alpha\le k$，因为$\alpha$最多$59$，当$k\le 100$可以直接用上述算法解决，考虑如何计算模$p^\alpha$的值。
设$i=ap+b$，那么有

$$\sum_{i=0}^n i^k=\sum_{i=0}^n(ap+b)^k=\sum_{i=0}^n\sum_{j=0}^k{k \choose j}a^jp^jb^{k-j}$$
我们注意到当 $j\ge \alpha$的时候值就为 $0$了，上界可以直接改成$\alpha-1$，
$$\sum_{i=0}^n\sum_{j=0}^{\alpha -1}{k \choose j}a^jp^jb^{k-j}=\sum_{j=0}^{\alpha -1}\sum_{a=0}^{\lfloor \frac{n}{p}\rfloor }a^j\sum_{b=0}^{n \bmod p}b^{k-j}$$
这个式子不是严格的， $a$可能只能到$\lfloor \frac{n}{p}\rfloor-1$，然后剩下不完整的一块 $b$，$a$的幂和直接使用之前的 $O(k^2)$做法， $b$的可以预处理对于每一个$j$的前缀和（因为还要解决不完整的一块），那就可以在 $O(\sum_{j=0}^{\alpha -1}j^2)$的复杂度搞定。
最后CRT合并就行了。注意因为有可能没有逆元，stirling数做法里的那个除 $j+1$和后来的组合数要特殊算。
还有，常数真难卡。
贴一个没卡常但比较正常的代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<map>
#define ll long long
#define fs first
#define sc second
using namespace std;
ll m,K,pri[300010],num,S[110][110],f[62][300010];
int ca,nd,top,cnt[18],vis[62];
bool flag[300010];
pair<int,ll> d[18];
ll mul(ll a,ll b,ll mod)
{
    a%=mod;b%=mod;
    ll tmp=(a*b-(ll)((long double)a/mod*b+1e-8)*mod)%mod;
    return (tmp<0?tmp+mod:tmp);
}
ll ksm(ll a,ll b,ll mod)
{
    ll r=1;
    for(;b;b>>=1)
    {
        if(b&1) r=mul(r,a,mod);
        a=mul(a,a,mod);
    }
    return r;
}
ll gcd(ll x,ll y)
{
    if(!y) return x;
    return gcd(y,x%y);
}
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 x,ll mod)
{
    ll xx,yy;
    exgcd(x,mod,xx,yy);
    return (xx%mod+mod)%mod;
}
void getpri(int n)
{
    memset(flag,1,sizeof(flag));
    flag[1]=0;
    for(int i=2;i<=n;i++)
    {
        if(flag[i]) pri[++num]=i;
        for(int j=1;j<=num&&i*pri[j]<=n;j++)
        {
            flag[i*pri[j]]=0;
            if(i%pri[j]==0) break;
        }
    }
}
ll C(ll a,ll b,ll mod)
{
    int mxp=0;
    for(ll i=2;i<=b;i++)
    {
        ll t=i;int j=1;
        for(;j<=num&&t>1;j++)while(t%pri[j]==0) t/=pri[j],vis[j]++;
        mxp=max(mxp,j); 
    }
    ll re=1;
    for(ll i=a;i>a-b;i--)
    {
        ll t=i;
        for(int j=1;j<=mxp;j++)
            while(vis[j]&&t%pri[j]==0) t/=pri[j],vis[j]--;
        re=mul(re,t,mod);
    }
    return re;
}
ll CRT(ll ans[])
{
    ll re=0;
    for(int t=1;t<=nd;t++)
    {
        ll Mi=m/d[t].sc;
        re=(re+mul(mul(ans[t],Mi,m),inv(Mi,d[t].sc),m))%m;
    }
    return re;
}
ll solve_K_100(ll n,ll K)
{
    ll ans[18];
    memset(ans,0,sizeof(ans));
    for(int t=1;t<=nd;t++)
    {
        ll mod=d[t].sc;
        for(ll j=0;j<=min(n,K);j++)
        {
            ll mi=1,fm=j+1;
            for(ll i=0;i<=j;i++)
            {
                ll g=gcd(fm,n+1-i);fm/=g;
                mi=mul(mi,(n+1-i)/g,mod);
            }
            ans[t]=(ans[t]+mul(S[K][j],mi,mod))%mod;
        }
    }   
    return CRT(ans);        
}
ll solve(ll n)
{
    ll ans[18];
    memset(ans,0,sizeof(ans));
    for(int t=1;t<=nd;t++)
    {
        ll q=n/d[t].fs,r=n%d[t].fs,mod=d[t].sc;
        for(int j=0;j<cnt[t];j++)
        {
            ll cal=mul(solve_K_100(q-1,j),f[j][d[t].fs-1],mod);
            cal=(cal+mul(ksm(q,j,mod),f[j][r],mod))%mod;
            cal=mul(ksm(d[t].fs,j,mod),cal,mod);
            cal=mul(cal,C(K,j,mod),mod);
            ans[t]=(ans[t]+cal)%mod;

        }
    }
    return CRT(ans);
}
int main()
{
    scanf("%lld%lld%d",&m,&K,&ca);
    getpri(300000);
    bool fff=1;
    for(int i=1;i<=num;i++)
    {
        ll t=m;
        if(t%pri[i]==0)
        {
            nd++;d[nd].fs=pri[i];d[nd].sc=1;
            while(t%pri[i]==0) t/=pri[i],d[nd].sc=d[nd].sc*pri[i],cnt[nd]++;
            if(cnt[nd]>1) fff=0;
        }
    }   

    S[0][0]=1;
    for(int i=1;i<=100;i++)
        for(int j=1;j<=i;j++)
            S[i][j]=(mul(S[i-1][j],j,m)+S[i-1][j-1])%m;

    if(K>100)
    {
        for(int i=1;i<=300000;i++)
            f[60][i]=ksm(i,K-60,m);
        for(int j=59;j>=0;j--)
            for(int i=1;i<=300000;i++)
                f[j][i]=mul(f[j+1][i],i,m);
        for(int j=0;j<=60;j++)
            for(int i=1;i<=300000;i++)
                f[j][i]=(f[j][i-1]+f[j][i])%m;  
    }   
    while(ca--)
    {
        ll n;
        scanf("%lld",&n);
        if(K<=100) printf("%lld\n",solve_K_100(n,K));
        else printf("%lld\n",solve(n));
    }
    return 0;
}