Description
给定
m,k
m
,
k
,共
T
T
次询问,每次输入一个,求
∑ni=1ik
∑
i
=
1
n
i
k
在
modm
mod
m
意义下的值
m
m
的最大质因子
2≤n,m,k≤1018,1≤T≤3∗103
2
≤
n
,
m
,
k
≤
10
18
,
1
≤
T
≤
3
∗
10
3
Analysis
首先分解质因数
m=∏peii
m
=
∏
p
i
e
i
,考虑对每个不同质因子分别求出
modpe
mod
p
e
意义下的值再用中国剩余定理求答案
将每一个
ik
i
k
拆成
(ip+j)k,i≤⌊n/p⌋,j<p
(
i
p
+
j
)
k
,
i
≤
⌊
n
/
p
⌋
,
j
<
p
,用二项式定理展开以后发现某一项里面如果有
px,x≥e
p
x
,
x
≥
e
那么显然这一项值为
0
0
,通过这个化简式子,综合运用斯特林数求自然数幂和等数学知识,再做些适当(变态)的预处理,就能以的复杂度解决此题了
一个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;
}