BZOJ3157 国王奇遇记 <倍增>

Posted on Symbols count in article: 674 Reading time ≈ 2 mins.

Problem

国王奇遇记

Time  Limit:  10  Sec\mathrm{Time\;Limit:\;10\;Sec}
Memory  Limit:  512  MB\mathrm{Memory\;Limit:\;512\;MB}

Description

Input

共一行，包括两个正整数NN和MM。

Output

共一行，为所求表达式的值对109+710^9+7取模的值。

Sample Input

1
5 3

Sample Output

1
36363

Hint

1≤N≤1091\le N\le10^9，1≤m≤2001\le m\le200

Source

Katharon\mathrm{Katharon}

标签：倍增

Solution

倍增处理幂和。
令

f(n,k)=∑i=1nikmif(n,k)=\sum_{i=1}^{n}i^km^i

则

f(n+1,k)=∑i=1n+1ikmi=m+∑i=2n+1ikmi=m+∑i=1n(i+1)kmi+1=m+m∑i=1n(i+1)kmi=m+m∑i=1nmi∑j=0k(kj)ij=m+m∑j=0k(kj)∑i=1nijmi=m+m∑j=0k(kj)f(n,j)f(2n,k)=∑i=12nikmi=∑i=1nikmi+∑i=n+12nikmi=f(n,k)+∑i=1n(i+n)kmi+n=f(n,k)+mn∑i=1n(i+n)kmi=f(n,k)+mn∑i=1nmi∑j=0k(kj)ijnk−j=f(n,k)+mn∑j=0k(kj)nk−j∑i=1nijmi=f(n,k)+mn∑j=0k(kj)nk−jf(n,j)\begin{aligned} f(n+1,k) &=\sum_{i=1}^{n+1}i^km^i\\ &=m+\sum_{i=2}^{n+1}i^km^i\\ &=m+\sum_{i=1}^{n}(i+1)^km^{i+1}\\ &=m+m\sum_{i=1}^{n}(i+1)^km^i\\ &=m+m\sum_{i=1}^{n}m^i\sum_{j=0}^{k}\binom{k}{j}i^j\\ &=m+m\sum_{j=0}^{k}\binom{k}{j}\sum_{i=1}^{n}i^jm^i\\ &=m+m\sum_{j=0}^{k}\binom{k}{j}f(n,j)\\ f(2n,k) &=\sum_{i=1}^{2n}i^km^i\\ &=\sum_{i=1}^{n}i^km^i+\sum_{i=n+1}^{2n}i^km^i\\ &=f(n,k)+\sum_{i=1}^{n}(i+n)^km^{i+n}\\ &=f(n,k)+m^n\sum_{i=1}^{n}(i+n)^km^i\\ &=f(n,k)+m^n\sum_{i=1}^{n}m^i\sum_{j=0}^{k}\binom{k}{j}i^jn^{k-j}\\ &=f(n,k)+m^n\sum_{j=0}^{k}\binom{k}{j}n^{k-j}\sum_{i=1}^{n}i^jm^i\\ &=f(n,k)+m^n\sum_{j=0}^{k}\binom{k}{j}n^{k-j}f(n,j)\\ \end{aligned}

于是可以将规模为nn的问题根据奇偶性降低到规模为n−1n-1和n2\frac{n}{2}的问题，递归处理即可。

加强版见BZOJ3516，加强版之再加强版见BZOJ4126。

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
#include <bits/stdc++.h>
#define MX 200
#define P 1000000007
using namespace std;
typedef long long lnt;
template <class T> inline void read(T &x) {
    x = 0; int c = getchar(), f = 1;
    for (; !isdigit(c); c = getchar()) if (c == 45) f = -1;
    for (; isdigit(c); c = getchar()) (x *= 10) += f*(c-'0');
}
int n, m; lnt f[60][MX+5], p[MX+5], fac[MX+5], inv[MX+5];
void init() {
    fac[0] = inv[0] = inv[1] = 1;
    for (int i = 1; i <= MX; i++) fac[i] = fac[i-1]*i%P;
    for (int i = 2; i <= MX; i++) inv[i] = (P-P/i*inv[P%i]%P)%P;
    for (int i = 2; i <= MX; i++) inv[i] = inv[i-1]*inv[i]%P;
}
lnt Pow(lnt x, int k) {
    lnt ret = 1;
    for (; k; k >>= 1, x = x*x%P)
        if (k&1) ret = ret*x%P;
    return ret;
}
lnt C(int n, int m) {return fac[n]*inv[m]%P*inv[n-m]%P;}
void calc(int n, int stp) {
    if (!n) return;
    if (n&1) {
        calc(n-1, stp+1);
        for (int i = 0; i <= m; i++) f[stp][i] = m;
        for (int i = 0; i <= m; i++) for (int j = 0; j <= i; j++)
            f[stp][i] = (f[stp][i]+m*C(i, j)%P*f[stp+1][j]%P)%P;
    } else {
        calc(n/2, stp+1); lnt pw = Pow(m, n/2); p[0] = 1;
        for (int i = 1; i <= m; i++) p[i] = p[i-1]*(n/2)%P;
        for (int i = 0; i <= m; i++) f[stp][i] = f[stp+1][i];
        for (int i = 0; i <= m; i++) for (int j = 0; j <= i; j++)
            f[stp][i] = (f[stp][i]+pw*C(i, j)%P*p[i-j]%P*f[stp+1][j]%P)%P;
    }
}
int main() {
    read(n), read(m);
    init(), calc(n, 0);
    printf("%lld\n", f[0][m]);
    return 0;
}