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

原创 2016年02月23日 20:53:15

  题意:求

i=1ni1modpk,p105

  
  不知道是啥题。
  我们可以注意到,当n<p的时候显然是一定有解的,而当np时,对于每个kp项,可能单独没有逆元,但是也可能n/pk=11/kp的倍数把p约掉,这样就存在逆元了。既然题目保证了一定存在逆元,那么我们可以把这部分分开出来算。
  设
F(n)=i=1ni1,G(n)=i=1,ijpni1
那么我们有
F(n)=G(n)+F(n/p)/p
因此我们每次可以将n的规模缩小p,函数F的次数是O(logpn)的。现在我们的问题是如何计算G(n)
  为了方便分析,不妨我们先设pn,这样我们可以有
G(n)=a=1p1b=0n/p11a+bp
这里我们还是比较难进行快速计算,可以想办法把形式比较好的1a+bp给拆开,比如幂级数展开:
1a+bp=a1i=0(1)ibiaipi
注意我们是在模pk的意义下做的这个玩意,因此我们可以将k1以后的项全部直接删掉。故得到
G(n)=a=1p1b=0n/p1a1i=0k1(1)ibiaipi=i=0k1(1)ia=1p1piai+1b=0n/p1bi

  推至此,我们可以发现前面的k1i=0p1a=1piai+1是可以在O(kp)的时间内计算出来的,接下来的问题就是计算n/pb=0bi了。
  出题人似乎给的做法是伯努利数或者直接矩阵乘法?然而我都不会那些东西所以YY了一下k次方和。
  以下是我YY的计算ni=0ik的做法,不知道有没有错。。。
  设Sk(n)表示上面要求的这个式子,尝试用组合数来计算这串东西。
  注意到
(nk)=nkk!=ki=0(1)i+ks(k,i)nik!
其中s(n,k)表示第一类斯特林数s(n,k)=(n1)s(n1,k)+s(n1,k1)。那么
Sk(n)=j=0n(k!(jk)i=0k1(1)i+ks(k,i)ji)=k!j=0n(jk)i=0k1(1)i+ks(k,i)j=0nji=k!(n+1k+1)i=0k1(1)i+ks(k,i)Si(n)=(n+1)k+1k+1i=0k1(1)i+ks(k,i)Si(n)

  于是我们只需要O(k2)地预处理出第一类斯特林数然后就可以按k来递推算了,当然边界是S1(n)=n(n+1)/2
  回到原问题中,既然我们要求的是n/p1b=0bi,那么就可以直接代进去得到:
G(n)=i=0k1(1)ia=1p1piai+1Si(n/p1)

  到了这一步,我们就可以愉快地将pi1ai+1全都预处理出来,然后就可以在O(kp)的时间计算G(n)了。
  补充一下当pn的时候,我们可以把后面的余项划归入前面一部分,若设r=nmodp,则此时
G(n)=i=0k1(1)ia=1rpiai+1Si(n/p)+i=0k1(1)ia=r+1p1piai+1Si(n/p1)

  再回到上一步F(n)=G(n)+F(n/p)/p,也就可以递归计算了。
  至此我们只需要预处理一些东西,然后递归快速计算就可以辣。
  时间复杂度?每一层计算G的复杂度是O(kp),递归层数是O(logpn),总复杂度是O(kplogpn)
  应该有更优的做法>_<

  为毛一直都没有人告诉我上面的幂级数展开少了个(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;
}
版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/GEOTCBRL/article/details/50725904

【GDKOI2016】小学生数学题(附带了乘法取模黑科技)

题目描述给定n,p,kn, p, k,其中pp是质数。 求∑ni=11i(modpk)\sum_{i=1}^n \frac{1}{i} \pmod {p^k} 题目保证答案是PQ\frac{P}{...
  • Yves___
  • Yves___
  • 2016-02-24 11:20:25
  • 1655

GDKOI2016 Day2 T4 小学生数学题

T4 小学生数学题 求∑i=1n1imodpk\sum_{i=1}^{n}{1\over i}\mod p^k 首先,我们把原式分成两部分,一部分对于p有逆元,另一部分没有。 即原式= 1p∑...
  • alan_cty
  • alan_cty
  • 2016-03-05 14:47:32
  • 843

GDKOI2016 题解

day 1     Problem 1. 魔卡少女   题意:动态维护区间内所有子序列的异或和的和,单点修改。   解法:先做一个序列的异或前缀和Si=Ai⊕Si−1S_i=A...
  • GEOTCBRL
  • GEOTCBRL
  • 2016-02-22 13:24:45
  • 1395

[GDKOI2016]小学生数学题

题目大意求 ∑i=1ni−1\sum_{i=1}^ni^{-1} 对pkp^k取模 p
  • WerKeyTom_FTD
  • WerKeyTom_FTD
  • 2016-03-09 18:54:29
  • 442

[GDKOI2016] 小学生数学题

Description求∑ni=11imodpk\sum_{i=1}^n \frac{1}{i} \bmod p^k。Constraintp≤105p \le 10^5 pk⋅n≤1018p^k \...
  • AcE_DengWx
  • AcE_DengWx
  • 2016-03-03 16:19:55
  • 518

GDKOI2016 小学生数学题

题目描述给定N,K,PN,K,P,求 (∑Ni=11i)modPK(\sum_{i=1}^{N} \frac{1}{i}) \bmod P^K 若答案不存在则输出-1数据范围N∗PK≤1018N ...
  • PhilipsWeng
  • PhilipsWeng
  • 2016-02-22 21:38:48
  • 1137

[JZOJ4367]【GDKOI2016】小学生数学题(口胡)

Description Solution既然保证最后答案有逆元,那么将原式拆成有逆元和没有逆元两部分 即∑a=1⌊np⌋1a∗p+∑i=0⌊np⌋−1∑j=1p−11i∗p+j(modpk)\su...
  • hzj1054689699
  • hzj1054689699
  • 2017-12-19 17:37:07
  • 113

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

题意:求∑i=1ni−1modpk,p≤105\sum_{i=1}^{n}i^{-1}\mod p^k,p\leq 10^5
  • GEOTCBRL
  • GEOTCBRL
  • 2016-02-23 20:53:15
  • 1743

GDKOI2016 题解

#DAY1 ##P1.魔卡少女 >抽象题意:动态维护一个区间内,所有的连续子串异或值和。 方法:很显然,我们可以想到用线段树来做, ##T2.不稳定的传送门 >抽象题意:给出概率和花费,求最优花...
  • HOWARLI
  • HOWARLI
  • 2016-02-26 10:48:27
  • 1289
收藏助手
不良信息举报
您举报文章:GDKOI2016 day 2 Problem 4. 小学生数学题 - 数学题
举报原因:
原因补充:

(最多只允许输入30个字)