[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}^{\frac{n}{p}-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));
}