【GDKOI2016】小学生数学题 【自然数幂求和】【斯特林数】

题意：求(sum_{i=1}^ni^{-1} mod p^k)的值。保证(p)为质数，(n imes p^k≤10^{18})
题解：神题！细节巨多！
设(f(n,k)=sum_{i=1}^{n}i^{-1} mod p^k)，(g(n,k)=sum_{i=1,i eq jp}^{n}i^{-1} mod p^k)
则(f(n,k)=g(n,k)+frac{f(lfloorfrac{n}{p} floor,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}^{lfloorfrac{n-a}{p} floor}frac{1}{a+bp})
(=sum_{a=1}^{p-1}sum_{b=0}^{lfloorfrac{n-a}{p} floor}(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)^isum_{a=1}^{p-1}frac{1}{a^{i+1}}sum_{b=0}^{lfloorfrac{n-a}{p} floor}b^i)
我们令(S_k(n))表示(sum_{i=0}^{n}i^k)
则(g(n,k)=sum_{i=0}^{k-1}(-p)^isum_{a=1}^{p-1}frac{1}{a^{i+1}}S_i(lfloorfrac{n-a}{p} floor))
如何求(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 imes 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=n!C_{x}^{n}-sum_{j=1}^{n-1}s(n,j)x^j ]

其实挺显然的，上面那个证明的式子(x^n)的项系数为1，所以一减就可以得到这个式子。
我们还要证一个东西：

[sum_{i=k}^nC_{i}^{k}=C_{n+1}^{k+1} ]

怎么证？

[sum_{i=k}^nC_{i}^{k}=sum_{i=k}^{n-1}C_{i}^{k}+C_{n}^{k}\=C_{n}^{k}+C_{n}^{k+1}\=C_{n+1}^{k+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)^isum_{a=1}^{p-1}frac{1}{a^{i+1}}S_i(lfloorfrac{n-a}{p} floor))
对于一些零碎的部分，我是直接暴力计算的。
这一题数域可能很大，我写了(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
",f(n,k));
	return 0;
}