题目描述
给你(n,k,p)（(p)为质数），求
[sum_{i=1}^nfrac{1}{i}mod p^k
]
保证有解。
(pleq {10}^5,np^kleq {10}^{18})
题解
为什么会有解？
可能会发生这样的情况：
[egin{align}
&sum_{i=1}^6frac{1}{i}mod 3\
=&frac{1}{1}+frac{1}{2}+frac{1}{4}+frac{1}{5}+frac{1}{3}+frac{1}{6}mod 3\
=&frac{1}{1}+frac{1}{2}+frac{1}{4}+frac{1}{5}+frac{1}{2}mod 3\
=&cdots
end{align}
]
因为几个分母是(p)的倍数的数的逆元加起来后分子可能会变成(p)的倍数然后就约掉了。
令
[egin{align}
f(n)&=sum_{i=1}^nfrac{1}{i}\
g(n)&=sum_{i=1,i
eq jp}^nfrac{1}{i}
end{align}
]
那么
[f(n)=g(n)+frac{f(lfloorfrac{n}{p}
floor)}{p}
]
因为我们的答案要对(p^k)取模，所以(f(lfloorfrac{n}{p} floor))在计算的时候要对(p^{k+1})取模，才能得到正确答案。
[egin{align}
g(n)&=sum_{i=a+bp}^{n}frac{1}{i}\
&=sum_{a=1}^{p-1}sum_{b=0}^{lfloorfrac{n-a}{p}
floor}frac{1}{a+bp}\
&=sum_{a=1}^{p-1}sum_{b=0}^{lfloorfrac{n-a}{p}
floor}a^{-1}sum_{i=0}^{k-1}{(-1)}^ifrac{b^i}{a^i}p^i\
&=sum_{i=0}^{k-1}{(-1)}^ip^isum_{a=1}^{p-1}frac{1}{a^{i+1}}sum_{b=0}^{lfloorfrac{n-a}{p}
floor}b^i
end{align}
]
这时候要算自然数幂和，但我们不能线性插值。
那就大力推一波公式吧。
(s_s(n,m)) 为带符号第一类斯特林数。
[egin{align}
S_k(n)&=sum_{i=1}^ni^k\
x^underline{n}&=sum_{k=0}^ns_s(n,k)x^k\
k!inom{n}{k}&=n^underline{k}\
i^k&=s_s(k,k)i^k\
i^k&=k!inom{i}{k}-i^underline{k}+s_s(k,k)i^k\
i^k&=k!inom{i}{k}-(i^underline{k}-s_s(k,k)i^k)\
i^underline{k}-s_s(k,k)i^k&=sum_{j=0}^ks_s(k,j)i^j-s_s(k,k)i^k\
&=sum_{j=0}^{k-1}s_s(k,j)i^j\
i^k&=k!inom{i}{k}-sum_{j=0}^{k-1}s_s(k,j)i^j\
S_k(n)&=sum_{i=0}^{n}(k!inom{i}{k}-sum_{j=0}^{k-1}s_s(k,j)i^j)\
&=sum_{i=0}^nk!inom{i}{k}-sum_{j=0}^{k-1}s_s(k,j)sum_{i=0}^ni^j\
&=k!inom{n+1}{k+1}-sum_{i=0}^{k-1}s_s(k,i)S_j(n)\
&=frac{{(n+1)}^underline{k+1}}{k+1}-sum_{i=0}^{k-1}s_s(k,i)S_j(n)
end{align}
]
还有一种方法
[egin{align}
S_m(n)&=sum_{i=1}^ni^m\
&=sum_{i=1}^nsum_{j=1}^megin{Bmatrix}m\jend{Bmatrix}i^underline{j}\
&=sum_{j=1}^megin{Bmatrix}m\jend{Bmatrix}sum_{i=1}^ni^underline{j}\
&=sum_{j=1}^megin{Bmatrix}m\jend{Bmatrix}frac{{(n+1)}^underline{j+1}}{j+1}
end{align}
]
预处理出斯特林数就可以快速算(S_k(n))了。
你可能会注意到，(frac{{(n+1)}^underline{k+1}}{k+1})有除法。
但是分子是(k+1)个连续的自然数幂的乘积，一定有一个是(k+1)的倍数，这样就可以消除分母的影响。
算(g)的话，每层是(O(kp))的，一共有(O(log_pn))层，所以总时间复杂度是(O(kplog_pn))。
乘法取模可以用快速乘（多一个(log)），也可以用黑科技，我偷懒用了int128。
代码
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
typedef __int128 lll;
ll p;
ll mul(ll a,ll b,const ll c)
{
return (lll)a*b%c;
}
ll fp(ll a,ll b){ll s=1;for(;b;b>>=1,a*=a)if(b&1)s*=a;return s;}
ll fp(ll a,ll b,const ll c){ll s=1;for(;b;b>>=1,a=mul(a,a,c))if(b&1)s=mul(s,a,c);return s;}
int pri[1000010];
int b[1000010];
ll pw[1000010];
int cnt;
ll ss[100][100];
void getstirling(ll k,ll md)
{
ss[0][0]=1;
for(int i=1;i<=k;i++)
for(int j=1;j<=i;j++)
ss[i][j]=(ss[i-1][j-1]-mul(i-1,ss[i-1][j],md))%md;
}
ll gao(ll n,ll k,ll md)
{
ll s=(n+1)/(k+1);
ll v=s*(k+1);
for(ll i=n+1;i>=n-k+1;i--)
s=mul(s,(i==v?1:i),md);
return s;
}
void gets(ll *s,ll n,ll k,ll md)
{
s[0]=n;
for(int i=1;i<=k;i++)
{
s[i]=gao(n,i,md);
for(int j=0;j<i;j++)
s[i]=(s[i]-mul(ss[i][j],s[j],md))%md;
}
s[0]++;
}
ll s1[100];
ll s2[100];
ll g(ll n,ll md,ll k)
{
ll r=n%p;
getstirling(k,md);
gets(s1,n/p-1,k,md);
gets(s2,n/p,k,md);
ll res=0,s;
ll t=fp(p,k)-fp(p,k-1)-1;
for(int i=0;i<k;i++)
{
s=0;
pw[1]=1;
cnt=0;
for(int j=1;j<p;j++)
b[j]=0;
for(int j=2;j<p;j++)
{
if(!b[j])
{
pri[++cnt]=j;
pw[j]=fp(fp(j,i+1),t,md);
}
for(int k=1;k<=cnt&&j*pri[k]<p;k++)
{
b[j*pri[k]]=1;
pw[j*pri[k]]=mul(pw[j],pw[pri[k]],md);
if(j%pri[k]==0)
break;
}
}
if(n%p==0)
for(int j=1;j<p;j++)
s=(s+mul(s1[i],pw[j],md))%md;
else
for(int j=1;j<p;j++)
if(j<=r)
s=(s+mul(s2[i],pw[j],md))%md;
else
s=(s+mul(s1[i],pw[j],md))%md;
s=mul(s*(i&1?-1:1),fp(p,i),md);
res=(res+s)%md;
}
// fprintf(stderr,"G(%lld, %lld) = %lld
",n,md,(res+md)%md);
return res;
}
ll f(ll n,ll md,ll k)
{
if(!n)
return 0;
return (g(n,md,k)+f(n/p,md*p,k+1)/p)%md;
}
int main()
{
#ifndef ONLINE_JUDGE
freopen("c.in","r",stdin);
freopen("c.out","w",stdout);
#endif
ll n,k;
scanf("%lld%lld%lld",&p,&k,&n);
ll md=fp(p,k);
ll ans=f(n,md,k);
ans=(ans+md)%md;
printf("%lld
",ans);
return 0;
}