[GDKOI2016]小学生数学题

题面

分析

设$$F(n,k)\equiv\sum_{i=1}^{n\frac{1}{i}\pmod{p}k}$$

\(p\)的倍数都没有逆元，因此一定是把\(p\)的倍数的倒数里的\(p\)提出后剩余部分为\(p\)的倍数。因为题目保证有解，因此我们按这个思路分开求就好了。

不妨先设\(n\)为\(p\)的倍数，剩余暴力即可。

再设出不为\(p\)的倍数的部分$$G(n,k)\equiv\sum_{i=1}^{{p-1}\sum_{j=0}}-1}\frac{1}{i+jp}\pmod{p^k}$$

则有：

\[\begin{align*} F(n,k)&\equiv\sum_{i=1}^n\frac{1}{i}\\ &\equiv\sum_{i=1}^\frac{n}{p}\frac{1}{i p}+G(n,k)\\ &\equiv\frac{1}{p}\sum_{i=1}^\frac{n}{p}\frac{1}{i}+G(n,k)\\ &\equiv\frac{F(\frac{n}{p},k)}{p}+G(n,k)\\ \end{align*}\]

题目已保证除\(p\)必定整除。但为了除到\(p^k\)范围内，将其模数从\(p^k\)改为\(p^{k+1}\)：$$F(n,k)\equiv\frac{F(\frac{n}{p},k+1)}{p}+G(n,k)$$

第一部分递归即可。则问题转换为如何快速求\(G(n,k)\)。

根据$$\frac{1}{1-x}=\sum_{k=0}^\infty x^k$$

考虑展开\(\frac{1}{a+bp}\)：

\[\begin{align*} \frac{1}{a+bp}&\equiv\frac{a^{-1}}{1+a^{-1}bp}\\ &\equiv a^{-1}\sum_{i=0}^\infty \left(-\frac{b p}{a}\right)^i\\ &\equiv \frac{1}{a}\sum_{i=0}^\infty (-1)^i\frac{b^i}{a^i}p^i\pmod{p^k}\\ &\equiv \frac{1}{a}\sum_{i=0}^{k-1} (-1)^i\frac{b^i}{a^i}p^i\\ \end{align*}\]

故

\[\begin{align*} G(n,k)&\equiv\sum_{i=1}^{p-1}\sum_{j=0}^{\frac{n}{p}-1}\frac{1}{i}\sum_{w=0}^{k-1} (-1)^w\frac{j^w}{i^w}p^w\\ &\equiv\sum_{i=1}^{p-1}\sum_{w=0}^{k-1}(-1)^w p^w\frac{1}{i^{w+1}}\sum_{j=0}^{\frac{n}{p}-1}j^w \end{align*}\]

实验后我们发现应该约定\(0^0=1\)。

那么我们预处理逆元和\(p^w\)和自然数幂求和，就可以\(O(k p)\)计算了。

自然数幂求和可以\(O(k^2)\)求出，具体可以看我的博客自然数幂求和(注意边界)

时间复杂度\(O(k p\log n)\)

代码

#include<iostream>
#include<cstdio>
using namespace std;
typedef long long ll;
ll n,p,k,m=1,S2[101][101],inv[100001],Snm[101];
inline ll add(ll a,ll b,ll m){return a+b>=m?a+b-m:a+b;}
inline ll mul(ll a,ll b,ll m){return ((a*b-(ll)((double)a/m*b+0.5)*m)%m+m)%m;}
inline ll fp(ll a,ll m){return a&1?m-1:1;}
void init(ll k,ll m){
	S2[0][0]=1;
	for(int i=1;i<=k;i++)for(int j=1;j<=k;j++)S2[i][j]=add(S2[i-1][j-1],mul(S2[i-1][j],j,m),m);
}
ll S(ll w,ll n,ll m){
	if(!w)return n+1;
	ll ans=0,facpw=n;
	for(int i=1;i<=w;i++){
		ans=add(ans,mul(S2[w][i],mul(facpw,inv[i+1],m),m),m);
		facpw=mul(facpw,n+m-i,m);
	}
	return mul(ans,n+1,m);
}
ll g(ll n,ll k,ll m){
	ll r=n%p,c=n/p,ans=0;
	for(int i=2;i<p;i++)inv[i]=mul(m-m/i,inv[m%i],m);
	for(int a=1;a<=r;a++){
		ll cnt=0,tpow=1,cc=mul(mul(inv[a],c,m),p,m);
		for(int i=0;i<k;i++){
			cnt=add(cnt,mul(fp(i,m),tpow,m),m);
			tpow=mul(tpow,cc,m);
		}
		ans=add(ans,mul(cnt,inv[a],m),m);
	}
	if(c){
		init(k-1,m);
		for(int w=0;w<k;w++)Snm[w]=S(w,c-1,m);
		for(int i=1;i<p;i++){
			ll tpow=inv[i],cc=mul(inv[i],p,m);
			for(int w=0;w<k;w++){
				ans=add(ans,mul(fp(w,m),mul(tpow,Snm[w],m),m),m);
				tpow=mul(tpow,cc,m);
			}
		}
	}
	return ans;
}
ll f(ll n,ll k,ll m){
	if(!n)return 0;
	else return add(g(n,k,m),f(n/p,k+1,m*p)/p,m);
}
int main(){
	scanf("%lld%lld%lld",&p,&k,&n);
	for(int i=1;i<=k;i++)m*=p;
	inv[1]=1;
	printf("%lld\n",f(n,k,m));
}