[GDKOI2016]小学生数学题

题面

分析

设$$F(n,k)equivsum_{i=1}^{nfrac{1}{i}pmod{p}k}$$

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

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

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

则有：

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

题目已保证除(p)必定整除。但为了除到(p^k)范围内，将其模数从(p^k)改为(p^{k+1})：$$F(n,k)equivfrac{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})：

[egin{align*} frac{1}{a+bp}&equivfrac{a^{-1}}{1+a^{-1}bp}\ &equiv a^{-1}sum_{i=0}^infty left(-frac{b p}{a} ight)^i\ &equiv frac{1}{a}sum_{i=0}^infty (-1)^ifrac{b^i}{a^i}p^ipmod{p^k}\ &equiv frac{1}{a}sum_{i=0}^{k-1} (-1)^ifrac{b^i}{a^i}p^i\ end{align*}]

故

[egin{align*} G(n,k)&equivsum_{i=1}^{p-1}sum_{j=0}^{frac{n}{p}-1}frac{1}{i}sum_{w=0}^{k-1} (-1)^wfrac{j^w}{i^w}p^w\ &equivsum_{i=1}^{p-1}sum_{w=0}^{k-1}(-1)^w p^wfrac{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 plog 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
",f(n,k,m));
}