[GDKOI2016]小学生数学题
记
$F(n)=\sum\limits_{i=1}^{n}i^{-1}$
$G(n)=\sum\limits_{i=1,i\neq jp}^{n}i^{-1}$
我们要算$F(n)\%p^k$
那么
$F(n)\%p^k=\frac{F( \left \lfloor \frac{n}{p} \right \rfloor )}{p}\%p^k+G(n)\%p^k$
我们知道$\frac{F( \left \lfloor \frac{n}{p} \right \rfloor )}{p}\%p^k=\frac{F( \left \lfloor \frac{n}{p} \right \rfloor )\%p^{k+1}}{p}$，其中$F( \left \lfloor \frac{n}{p} \right \rfloor )\%p^{k+1}$可以可以递归算，所以我们重点要考虑的是$G(n)\%p^k$
不妨设$p|n$，那么
$G(n)\%p^k=\sum\limits_{a=1}^{p-1}\sum\limits_{b=0}^{\left \lfloor \frac{n}{p} \right \rfloor -1}(a+bp)^{-1}\%p^k$
用广义二项式定理展开：
$(a+bp)^{-1}=\sum\limits_{i=0}^{oo}C_{-1}^{i}a^{-1-i}b^ip^i=\sum\limits_{i=0}^{oo}(-1)^ia^{-1-i}b^ip^i$
又因为是在模$p^k$意义下的，所以
$(a+bp)^{-1}\%p^k=\sum\limits_{i=0}^{k-1}(-1)^ia^{-1-i}b^ip^i\%p^k$
所以
$G(n)\%p^k=\sum\limits_{a=1}^{p-1}\sum\limits_{b=0}^{\left \lfloor \frac{n}{p} \right \rfloor -1}\sum\limits_{i=0}^{k-1}(-1)^ia^{-1-i}b^ip^i\%p^k$
$=\sum\limits_{i=0}^{k-1}(-1)^ip^i\sum\limits_{a=1}^{p-1}a^{-1-i}\sum\limits_{b=0}^{\left \lfloor \frac{n}{p} \right \rfloor -1}b^i\%p^k$
需要注意的一点是，用二项式定理的时候，规定$0^0=1$
我们枚举$i$和$a$，要$O(kp)$的时间复杂度；剩下的$\sum\limits_{b=0}^{\left \lfloor \frac{n}{p} \right \rfloor -1}b^i\%p^k$就是自然数幂和，记$S_{i}(n)=\sum\limits_{b=0}^{n}b^i\%p^k$，可以用$O(k^2)$的时间内预处理出所有的$S_{i}(\left \lfloor \frac{n}{p} \right \rfloor-1)\%p^k(0\leq i\leq k-1)$
如果$n$不是$p$的倍数，剩下的零头乱搞就行了
所以每递归一次的时间复杂度是$O(kp+k^2)$
所以总的时间复杂度是$O(log_{p}n(kp+k^2))$
#include<cstdio> #include<cstdlib> #include<iostream> #include<fstream> #include<algorithm> #include<cstring> #include<string> #include<cmath> #include<queue> #include<stack> #include<map> #include<utility> #include<set> #include<bitset> #include<vector> #include<functional> #include<deque> #include<cctype> #include<climits> #include<complex> #include<cassert> //#include<bits/stdc++.h>适用于CF,UOJ，但不适用于poj using namespace std; typedef long long LL; typedef double DB; typedef pair<int,int> PII; typedef pair<DB,DB> PDD; typedef complex<DB> CP; typedef vector<int> VI; #define mmst(a,v) memset(a,v,sizeof(a)) #define mmcy(a,b) memcpy(a,b,sizeof(a)) #define fill(a,l,r,v) fill(a+l,a+r+1,v) #define re(i,a,b) for(i=(a);i<=(b);i++) #define red(i,a,b) for(i=(a);i>=(b);i--) #define fi first #define se second #define mp(a,b) make_pair(a,b) #define pb(a) push_back(a) #define SF scanf #define PF printf #define two(k) (1<<(k)) #define SZ(x) (int(x.size())) #define all(x) (x).begin(),(x).end() #define ire(i,v,x) for(i=0,v=i<SZ(x)?x[i]:0;i<SZ(x);v=x[++i]) template<class T>inline T sqr(T x){return x*x;} template<class T>inline void upmin(T &t,T tmp){if(t>tmp)t=tmp;} template<class T>inline void upmax(T &t,T tmp){if(t<tmp)t=tmp;} inline int sgn(DB x){if(abs(x)<1e-9)return 0;return(x>0)?1:-1;} const DB Pi=acos(-1.0); int gint() { int res=0;bool neg=0;char z; for(z=getchar();z!=EOF && z!='-' && !isdigit(z);z=getchar()); if(z==EOF)return 0; if(z=='-'){neg=1;z=getchar();} for(;z!=EOF && isdigit(z);res=res*10+z-'0',z=getchar()); return (neg)?-res:res; } LL gll() { LL res=0;bool neg=0;char z; for(z=getchar();z!=EOF && z!='-' && !isdigit(z);z=getchar()); if(z==EOF)return 0; if(z=='-'){neg=1;z=getchar();} for(;z!=EOF && isdigit(z);res=res*10+z-'0',z=getchar()); return (neg)?-res:res; } const LL maxp=100100; const LL maxk=1010; LL multiply(LL a,LL b,LL MOD) { LL t=(LL)((DB)a*b/MOD); return a*b-t*MOD; } LL su[maxk][maxk],sb[maxk]; void calc(LL n,LL d,LL m) { LL i,j; su[0][0]=1; re(i,1,d) { su[i][0]=0; re(j,1,d-1)su[i][j]=(multiply(i-1,su[i-1][j],m)+su[i-1][j-1])%m; su[i][i]=1; } sb[0]=n%m; re(i,1,d) { sb[i]=1; red(j,n+1,n+1-i)if(j%(i+1)==0)sb[i]=multiply(sb[i],j/(i+1),m); else sb[i]=multiply(sb[i],j,m); LL l=(i%2==0)?1:-1; re(j,0,i-1){(sb[i]-=l*multiply(su[i][j],sb[j],m))%=m;l=-l;} } } LL reva[maxp],powreva[maxp]; LL F(LL n,LL p,LL k,LL m) { LL i,j,res=0; if(n<p) { reva[1]=1;re(j,2,n)reva[j]=multiply(reva[m%j],(m-m/j),m); re(j,1,n)(res+=reva[j])%=m; return res; } calc(n/p-1,k-1,m); (sb[0]+=1)%=m; reva[1]=1;re(j,2,p-1)reva[j]=multiply(reva[m%j],(m-m/j),m); re(j,1,p-1)powreva[j]=1; LL sp=1,sa; re(i,0,k-1) { sa=0; re(j,1,p-1)powreva[j]=multiply(powreva[j],reva[j],m),(sa+=powreva[j])%=m; (res+=multiply(multiply(sp,sa,m),sb[i],m))%=m; sp=multiply(sp,-p,m); } if(n%p!=0) { re(j,1,n%p)powreva[j]=1; LL sp=1,sa; re(i,0,k-1) { sa=0; re(j,1,n%p)powreva[j]=multiply(powreva[j],reva[j],m),(sa+=powreva[j])%=m; (res+=multiply(sp,sa,m))%=m; sp=multiply(sp,(-p)*(n/p),m); } } (res+=F(n/p,p,k+1,m*p)/p)%=m; return res; } int main() { freopen("math.in","r",stdin); freopen("math.out","w",stdout); LL i,p=gll(),k=gll(),n=gll(),m=1;re(i,1,k)m*=p; cout<<(F(n,p,k,m)%m+m)%m<<endl; return 0; }
