GDKOI2016 day 2 Problem 4. 小学生数学题 - 数学题

　　题意：求

$$\sum_{i=1}^{n}i^{-1}\mod p^k,p\leq 10^5$$
　　
　　不知道是啥题。
　　我们可以注意到，当 $n<p$的时候显然是一定有解的，而当 $n\geq p$时，对于每个 $kp$项，可能单独没有逆元，但是也可能 $\sum_{k=1}^{\lfloor n/p\rfloor}1/k$是 $p$的倍数把$p$约掉，这样就存在逆元了。既然题目保证了一定存在逆元，那么我们可以把这部分分开出来算。
　　设 $$F(n)=\sum_{i=1}^{n}i^{-1},G(n)=\sum_{i=1,i\neq jp}^{n}i^{-1}$$ 那么我们有 $$F(n)=G(n)+F(\lfloor n/p\rfloor)/p$$ 因此我们每次可以将 $n$的规模缩小$p$，函数 $F$的次数是$O(\log_p n)$的。现在我们的问题是如何计算 $G(n)$。
　　为了方便分析，不妨我们先设 $p\mid n$，这样我们可以有 $$G(n)=\sum_{a=1}^{p-1}\sum_{b=0}^{\lfloor n/p \rfloor-1}\frac{1}{a+bp}$$ 这里我们还是比较难进行快速计算，可以想办法把形式比较好的 $\frac{1}{a+bp}$给拆开，比如幂级数展开： $$\frac{1}{a+bp}=a^{-1}\sum_{i=0}^{\infty}(-1)^i\frac{b^i}{a^i}p^i$$ 注意我们是在模 $p^k$的意义下做的这个玩意，因此我们可以将 $k-1$以后的项全部直接删掉。故得到 $$G(n)=\sum_{a=1}^{p-1}\sum_{b=0}^{\lfloor n/p \rfloor-1}a^{-1}\sum_{i=0}^{k-1}(-1)^i\frac{b^i}{a^i}p^i\\=\sum_{i=0}^{k-1}(-1)^i\sum_{a=1}^{p-1}\frac{p^{i}}{a^{i+1}}\sum_{b=0}^{\lfloor n/p\rfloor-1}b^i$$
　　推至此，我们可以发现前面的 $\sum_{i=0}^{k-1}\sum_{a=1}^{p-1}\frac{p^{i}}{a^{i+1}}$是可以在 $O(kp)$的时间内计算出来的，接下来的问题就是计算 $\sum_{b=0}^{\lfloor n/p\rfloor}b^i$了。
　　出题人似乎给的做法是伯努利数或者直接矩阵乘法？然而我都不会那些东西所以YY了一下 $k$次方和。
　　以下是我YY的计算$\sum_{i=0}^{n}i^k$的做法，不知道有没有错。。。
　　设 $S_k(n)$表示上面要求的这个式子，尝试用组合数来计算这串东西。
　　注意到 $$\binom{n}{k}=\frac{n^{\underline k}}{k!}=\frac{\sum_{i=0}^{k}(-1)^{i+k}s(k,i)n^{i}}{k!}$$ 其中 $s(n,k)$表示第一类斯特林数 $s(n,k)=(n-1)s(n-1,k)+s(n-1,k-1)$。那么 $$S_k(n)=\sum_{j=0}^{n}\left( k!\binom{j}{k}-\sum_{i=0}^{k-1}(-1)^{i+k}s(k,i)j^{i}\right )\\=k!\sum_{j=0}^{n}\binom{j}{k}-\sum_{i=0}^{k-1}(-1)^{i+k}s(k,i)\sum_{j=0}^{n}j^i\\=k!\binom{n+1}{k+1}-\sum_{i=0}^{k-1}(-1)^{i+k}s(k,i)S_i(n)\\=\frac{(n+1)^{\underline{k+1}}}{k+1}-\sum_{i=0}^{k-1}(-1)^{i+k}s(k,i)S_i(n)$$
　　于是我们只需要 $O(k^2)$地预处理出第一类斯特林数然后就可以按k来递推算了，当然边界是 $S_1(n)=n(n+1)/2$。
　　回到原问题中，既然我们要求的是 $\sum_{b=0}^{\lfloor n/p \rfloor-1}b^i$，那么就可以直接代进去得到： $$G(n)=\sum_{i=0}^{k-1}(-1)^i\sum_{a=1}^{p-1}\frac{p^{i}}{a^{i+1}}S_i(\lfloor n/p \rfloor-1)$$
　　到了这一步，我们就可以愉快地将 $p^i$和 $\frac{1}{a^{i+1}}$全都预处理出来，然后就可以在 $O(kp)$的时间计算 $G(n)$了。
　　补充一下当 $p\nmid n$的时候，我们可以把后面的余项划归入前面一部分，若设 $r=n\mod p$，则此时 $$G(n)=\sum_{i=0}^{k-1}(-1)^i\sum_{a=1}^{r}\frac{p^{i}}{a^{i+1}}S_i(\lfloor n/p \rfloor)+\sum_{i=0}^{k-1}(-1)^i\sum_{a=r+1}^{p-1}\frac{p^{i}}{a^{i+1}}S_i(\lfloor n/p \rfloor-1)$$
　　再回到上一步 $F(n)=G(n)+F(n/p)/p$，也就可以递归计算了。
　　至此我们只需要预处理一些东西，然后递归快速计算就可以辣。
　　时间复杂度？每一层计算 $G$的复杂度是$O(kp)$,递归层数是 $O(\log_p n)$，总复杂度是 $O(kp\log_p n)$。
　　应该有更优的做法>_<

　　为毛一直都没有人告诉我上面的幂级数展开少了个$(-1)^i$。。。。调了好久啊TAT
　　奇怪的Code：（懒得删注释了）

#include <bits/stdc++.h>
//using namespace std;
#define rep(i,a,b) for(int i = a , _ = b ; i <= _ ; i ++)
#define per(i,a,b) for(int i = a , _ = b ; i >= _ ; i --)
#define cr(x) memset(x , 0 , sizeof x)

#define gprintf(...) //fprintf(stderr , __VA_ARGS__)

typedef long long ll;

const int maxp = 100007;
const int maxk = 121;

ll mod , n , s[maxk][maxk] , inva[maxp];
ll powa[maxp][maxk] , S[maxk] , powp[maxk];

inline ll add(ll a , ll b , ll mod) { a = (a + b) % mod ; if (a < 0) a += mod ; if (a >= mod) a -= mod ; return a ; }
inline ll dec(ll a , ll b , ll mod) { a = (a - b) % mod ; if (a < 0) a += mod ; if (a >= mod) a -= mod ; return a ; }
inline ll qmul(ll a , ll b , ll p) {
    return (a * b - (ll)(a / (double)p * b + 1e-3) * p + p) % p;
}

#define neg(x) (((x) & 1) ? -1 : 1)

int p , k;

void input() {
    scanf("%d%d%lld" , &p , &k , &n);
    mod = 1;
    rep (i , 1 , k) mod = mod * p;
}

void init_stirling(int k , ll mod) {
//  gprintf("initalization of stirling mod %lld\n" , mod);
    s[0][0] = 1;
    s[1][0] = 0 , s[1][1] = 1;
    rep (i , 2 , k) {
        s[i][0] = 0;
        rep (j , 1 , k) {
            s[i][j] = add(qmul(i - 1 , s[i - 1][j] , mod) , s[i - 1][j - 1] , mod);
//          gprintf("%lld%c" , s[i][j] , j == k ? '\n' : ' ');
        }
    }
}

void init_inv(int k , ll mod) {
//  gprintf("calculating inverse of mod %lld\n" , mod);
    inva[1] = 1;
    powa[1][1] = 1;
    rep (i , 2 , p - 1) {
        inva[i] = dec(0 , qmul(inva[mod % i] , mod / i , mod) , mod);
//      gprintf("%lld\n" , qmul(inva[i] , i , mod));
        powa[i][1] = inva[i];
    }
    powp[0] = 1;
    rep (i , 1 , k - 1) powp[i] = powp[i - 1] * p;
    rep (a , 1 , p - 1) {
        ll t = inva[a];
        rep (i , 2 , k)
            powa[a][i] = qmul(powa[a][i - 1] , t , mod);
    }
}

namespace ASS {

    #include <assert.h>

    static ll tmp[maxk];
    static ll _pow[100 * 10][25];

    void _S(ll n , ll mod) {
        rep (i , 1 , n) {
            _pow[i][1] = i;
            rep (j , 2 , k)
                _pow[i][j] = qmul(_pow[i][j - 1] , i , mod);
        }
        rep (i , 0 , k) {
            tmp[i] = 0;
            rep (j , 1 , n)
                tmp[i] = add(tmp[i] , _pow[j][i] , mod);
            if (!i) tmp[i] = (n + 1) % mod;
            gprintf("%lld%c" , tmp[i] , i == k ? '\n' : ' ');
            assert(tmp[i] == S[i]);
        }
    }

}

void get_S(ll n , ll mod , int k) {
//  gprintf("calculating sum %lld of power %d\n" , n , k);
    S[0] = (n + 1) % mod;
//  if (!n) S[0] = 1;
    if (n < k) {
        static ll pown[maxk][maxk];
        rep (i , 1 , n) {
            pown[i][1] = i;
            rep (j , 2 , k)
                pown[i][j] = qmul(pown[i][j - 1] , i , mod);
        }
        rep (i , 1 , k) {
            S[i] = 0;
            rep (j , 1 , n)
                S[i] = add(S[i] , pown[j][i] , mod);
        }
    } else {
        int cur_k = -1;
        rep (K , 1 , k) {
            S[K] = 1;
            rep (j , 0 , K) {
                if ((n + 1 - j) % (K + 1) == 0)
                    S[K] = qmul((n + 1 - j) / (K + 1) , S[K] , mod);
                else
                    S[K] = qmul(n + 1 - j , S[K] , mod);
            }
            int cur_i = cur_k;
            rep (i , 0 , K - 1) {
                S[K] = dec(S[K] , cur_i * qmul(s[K][i] , S[i] , mod) , mod);
                cur_i = - cur_i;
            }
            cur_k = - cur_k;
        }
    }
    #ifdef DEBUG
        rep (i , 0 , k) {
            gprintf("%lld%c" , S[i] , i == k ? '\n' : ' ');
        }
        ASS::_S(n , mod);
    #endif
}

ll G(ll n , ll mod , int k) {
    init_inv(k , mod);
    init_stirling(k , mod);
//  gprintf("calculating G(%lld , %lld)\n" , n , mod);
    ll ret = 0;
    if (n % p == 0) {
        get_S(n / p - 1 , mod , k);
        rep (a , 1 , p - 1)
            rep (i , 0 , k - 1) {
//              gprintf("p^%d : %lld , 1/%d^(%d+1) %lld , S[%d] : %lld\n" , i , powp[i] , a , i , powa[a][i + 1] , i , S[i]);               
                ret = add(ret , neg(i) * qmul(qmul(powp[i] , powa[a][i + 1] , mod) , S[i] , mod) , mod);
            }
    } else {
        if (n < p) {
            rep (i , 1 , n)
                ret = add(ret , inva[i] , mod);
            gprintf("result of G(%lld , %lld) = %lld\n" , n , mod , ret);
            return ret;
        }
        ll r = n % p;
        get_S(n / p , mod , k);
        rep (a , 1 , r)
            rep (i , 0 , k - 1) {
//              gprintf("p^%d : %lld , 1/%d^(%d+1) %lld , S[%d] : %lld\n" , i , powp[i] , a , i , powa[a][i + 1] , i , S[i]);
                ret = add(ret , neg(i) * qmul(qmul(powp[i] , powa[a][i + 1] , mod) , S[i] , mod) , mod);
            }
        get_S(n / p - 1 , mod , k);
        rep (a , r + 1 , p - 1)
            rep (i , 0 , k - 1)
                ret = add(ret , neg(i) * qmul(qmul(powp[i] , powa[a][i + 1] , mod) , S[i] , mod) , mod);
    }
    gprintf("result of G(%lld , %lld) = %lld\n" , n , mod , ret);
    return ret;
}

ll F(ll n , ll mod , int k) {
    if (!n) return 0;
    return (G(n , mod , k) + F(n / p , mod * p , k + 1) / p) % mod;
}

void solve() {
//  printf("%lld\n" , G(102728 , 823543));
    ll ans = F(n , mod , k);
    printf("%lld\n" , ans);
}

int main() {
    #ifndef ONLINE_JUDGE
        freopen("math.in" , "r" , stdin);
    #endif
    input();
    solve();
    return 0;
}