【GDKOI2016】小学生数学题 【自然数幂求和】【斯特林数】
题意：求\(\sum_{i=1}^ni^{-1}\ mod\ p^k\)的值。保证\(p\)为质数，\(n\times p^k≤10^{18}\)
题解：神题！细节巨多！
设\(f(n,k)=\sum_{i=1}^{n}i^{-1}\ mod\ p^k\)，\(g(n,k)=\sum_{i=1,i\neq jp}^{n}i^{-1}\ mod\ p^k\)
则\(f(n,k)=g(n,k)+\frac{f(\lfloor\frac{n}{p}\rfloor,k+1)}{p}\)
为什么后面那里是k+1呢？
因为若\(a≡b\ mod\ c\)，则\(ak≡bk\ mod\ ck\)
我们把那一部分先整体乘\(p\)，最后再除回来就好。
这样我们就可以递归计算了。
如何求g(n,k)？
\(g(n,k)\)
\(=\sum_{i=a+bp<=n}\frac{1}{i}\)
\(=\sum_{a=1}^{p-1}\sum_{b=0}^{\lfloor\frac{n-a}{p}\rfloor}\frac{1}{a+bp}\)
\(=\sum_{a=1}^{p-1}\sum_{b=0}^{\lfloor\frac{n-a}{p}\rfloor}(a+bp)^{-1}\)
这个东西用广义二项式定理展开一下
\((a+bp)^{-1}=\sum_{i=0}^{\infty}C_{-1}^{i}a^{-1-i}b^ip^i=\sum_{i=0}^{\infty}(-1)^ia^{-1-i}b^ip^i\)
又因为我们模了一个\(p^k\)
所以
\((a+bp)^{-1}=\sum_{i=0}^{k-1}(-1)^ia^{-1-i}b^ip^i\)
所以
\(g(n,k)\)
\(=\sum_{i=0}^{k-1}(-p)^i\sum_{a=1}^{p-1}\frac{1}{a^{i+1}}\sum_{b=0}^{\lfloor\frac{n-a}{p}\rfloor}b^i\)
我们令\(S_k(n)\)表示\(\sum_{i=0}^{n}i^k\)
则\(g(n,k)=\sum_{i=0}^{k-1}(-p)^i\sum_{a=1}^{p-1}\frac{1}{a^{i+1}}S_i(\lfloor\frac{n-a}{p}\rfloor)\)
如何求\(S_k(n)?\)
我们先证一个东西：
\(x^{\overline n}=\Pi_{i=x}^{x+n-1}i=\sum_{k}s(n,k)x^k\)，s是第一类斯特林数。
可以转化为证\([x^p]x^{\overline n}=[x^p]\sum_{k}s(n,k)x^{k}=s(n,p)\)
显然在\(n=1\)的情况下是成立的。
\([x^k]x^{\overline n}\)
\(=[x^k](x+n-1)x^{\overline{n-1}}\)
\(=[x^k](x+n-1)x^{\overline{n-1}}\)
\(=[x^k]x\times x^{\overline{n-1}}+[x^k](n-1)x^{\overline{n-1}}\)
\(=[x^{k-1}]x^{\overline{n-1}}+[x^k](n-1)x^{\overline{n-1}}\)
\(=s(n-1,k-1)+(n-1)s(n-1,k)\)
\(=s(n,k)\)
证明完毕。
我们再证一个东西：
其实挺显然的，上面那个证明的式子\(x^n\)的项系数为1，所以一减就可以得到这个式子。
我们还要证一个东西：
怎么证？
证明完毕。
终于要继续自然数幂求和啦！有了之前推的一堆式子，我们就可以化简式子了。
\(S_{k}(n)\)
\(=\sum_{i=0}^{n}i^k\)
\(=\sum_{i=0}^{n}k!C_{i}^{k}-\sum_{j=1}^{k-1}s(k,j)i^j\)
\(=k!(\sum_{i=0}^{n}C_{i}^{k})-\sum_{i=0}^{n}\sum_{j=1}^{k-1}s(k,j)i^j\)
\(=k!C_{n+1}^{k+1}\sum_{j=1}^{k-1}s(k,j)\sum_{i=0}^{n}i^j\)
对于固定的\(n\)，这个可以\(O(k^2)\)预处理出来。然后代会这个式子里计算。
\(g(n,k)=\sum_{i=0}^{k-1}(-p)^i\sum_{a=1}^{p-1}\frac{1}{a^{i+1}}S_i(\lfloor\frac{n-a}{p}\rfloor)\)
对于一些零碎的部分，我是直接暴力计算的。
这一题数域可能很大，我写了\(O(1)\)快速乘法。
这题好难写啊！QAQ
代码
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define int long long
using namespace std;
int p,k,n,x,y,r[100005],s[105][105];
int mul(int a,int b,int mod){
int c=a*(double)b/mod+0.5;
int res=a*b-c*mod;
if(res<0){
res+=mod;
}
return res;
}
int fastpow(int a,int x,int mod){
a%=mod;
int res=1;
while(x){
if(x&1){
res=mul(res,a,mod);
}
x>>=1;
a=mul(a,a,mod);
}
return res;
}
int exgcd(int a,int b,int &x,int &y){
if(!b){
x=1;
y=0;
return a;
}
int d=exgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}
int getinv(int a,int b){
exgcd(a,b,x,y);
return (x%b+b)%b;
}
void init(int n,int k){
const int mod=pow(p,k);
r[0]=n+1;
s[0][0]=1;
for(int i=1;i<=k;i++){
for(int j=1;j<=k;j++){
s[i][j]=(s[i-1][j-1]+mul(s[i-1][j],i-1,mod))%mod;
}
}
for(int i=1;i<=k;i++){
for(int j=1;j<=i;j++){
if((i+j)&1){
s[i][j]=mod-s[i][j];
}
}
}
for(int i=1;i<=k;i++){
r[i]=1;
bool flag=false;
for(int j=n+1;j>=n-i+1;j--){
if(!flag&&j%(i+1)==0){
flag=true;
r[i]=mul(r[i],j/(i+1),mod);
}else{
r[i]=mul(r[i],j,mod);
}
}
if(!flag){
r[i]=mul(r[i],getinv(i+1,mod),mod);
}
for(int j=1;j<i;j++){
r[i]=(r[i]-mul(s[i][j],r[j],mod)+mod)%mod;
}
}
return;
}
int g(int n,int k){
const int mod=pow(p,k);
int rem=n%p,res=0;
for(int i=n-rem+1;i<=n;i++){
res=(res+getinv(i,mod))%mod;
}
n-=rem;
if(!n){
return res;
}
init((n-1)/p,k);
for(int i=0,base=1;i<k;i++,base=(base*(-p)%mod+mod)%mod){
for(int j=1;j<p;j++){
res=(res+mul(mul(base,getinv(fastpow(j,i+1,mod),mod),mod),r[i],mod))%mod;
}
}
return res;
}
int f(int n,int k){
if(!n){
return 0;
}
return (g(n,k)+f(n/p,k+1)/p)%(int)(pow(p,k));
}
signed main(){
scanf("%lld%lld%lld",&p,&k,&n);
printf("%lld\n",f(n,k));
return 0;
}
