首先考虑计算模质数下的自然数幂和,通过stirling数转化成下降幂,
∑i=0nik=∑i=0n∑j=0k{kj}ij–=∑j=0k{kj}∑i=0nij–
∑
i
=
0
n
i
k
=
∑
i
=
0
n
∑
j
=
0
k
{
k
j
}
i
j
_
=
∑
j
=
0
k
{
k
j
}
∑
i
=
0
n
i
j
_
然后考虑对后面那个下降幂的和进行积分,我们注意到 (x+1)k––−xk––=kxk−1––––– ( x + 1 ) k _ − x k _ = k x k − 1 _ ,变一下就得到 xk––=(x+1)k+1–––––−xk+1–––––k+1 x k _ = ( x + 1 ) k + 1 _ − x k + 1 _ k + 1 ,设 g(x)=xk+1–––––k+1 g ( x ) = x k + 1 _ k + 1 ,那么 Δg(x)=g(x+1)−g(x)=xk–– Δ g ( x ) = g ( x + 1 ) − g ( x ) = x k _ ,那么和式就是值就是 g(n+1) g ( n + 1 ) 了(和普通幂函数的求导好像形式一样),
=∑j=0k{kj}(n+1)j+1––––j+1
=
∑
j
=
0
k
{
k
j
}
(
n
+
1
)
j
+
1
_
j
+
1
于是可以 O(k2) O ( k 2 ) 解决问题。自然数幂和还有更优的做法,但这题有这个就可以解决。
将模数
m
m
分解质因数,假设所有
α≤k
α
≤
k
,因为
α
α
最多
59
59
,当
k≤100
k
≤
100
可以直接用上述算法解决,考虑如何计算模
pα
p
α
的值。
设
i=ap+b
i
=
a
p
+
b
,那么有
∑i=0nik=∑i=0n(ap+b)k=∑i=0n∑j=0k(kj)ajpjbk−j
∑
i
=
0
n
i
k
=
∑
i
=
0
n
(
a
p
+
b
)
k
=
∑
i
=
0
n
∑
j
=
0
k
(
k
j
)
a
j
p
j
b
k
−
j
我们注意到当 j≥α j ≥ α 的时候值就为 0 0 了,上界可以直接改成,
∑i=0n∑j=0α−1(kj)ajpjbk−j=∑j=0α−1∑a=0⌊np⌋aj∑b=0nmodpbk−j
∑
i
=
0
n
∑
j
=
0
α
−
1
(
k
j
)
a
j
p
j
b
k
−
j
=
∑
j
=
0
α
−
1
∑
a
=
0
⌊
n
p
⌋
a
j
∑
b
=
0
n
mod
p
b
k
−
j
这个式子不是严格的, a a 可能只能到,然后剩下不完整的一块 b b ,的幂和直接使用之前的 O(k2) O ( k 2 ) 做法, b b 的可以预处理对于每一个的前缀和(因为还要解决不完整的一块),那就可以在 O(∑α−1j=0j2) O ( ∑ j = 0 α − 1 j 2 ) 的复杂度搞定。
最后CRT合并就行了。注意因为有可能没有逆元,stirling数做法里的那个除 j+1 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;
}