3A TrichyInequality
Description
求出满足 ∑mi=1xi≤s,∀i≤m,xi>0,∀i≤n,xi≤t . 的向量 X 的解数。
Difficulty
★★★
MainAlgorithm
矩乘加速
Complexity
O((m−n)3logn)
Solution
标解给的
O((m−n)2)
太科幻了……不管了,写了个矩乘。
首先我们枚举前
n
个
x
的取值。
则根据简单的组合数学原理我们知道答案是
∑∑ni=1xi=k,∀i,1≤xi≤t(s−km−n)
.
注意到那个组合数对于
k
是一个
m−n
次多项式,我们可以比较方便地展开并求出每一项的系数。
将系数提前,就变成了
∑m−ni=1ai∗∑∑ni=1xi=k,∀i,1≤xi≤tki
.
如何求
ki
?
不妨设
fn,m=∑∑ni=1xi=k,∀i,1≤xi≤tkm
.
可以很容易地将
k
写成
k′+xn
,然后二项式展开得到递推式。
矩乘加速将
fn,0∼fn,m−n
求出来即可。
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#define Rep(i, x, y) for (int i = x; i <= y; i ++)
#define Dwn(i, x, y) for (int i = x; i >= y; i --)
#define RepE(i, x) for(int i = pos[x]; i; i = g[i].nex)
using namespace std;
typedef long long LL;
const int N = 105, mod = 1000000007, M = 100005;
LL ans, C[N][N], h[N][N]; int n, pw[M][N];
struct Mt {
LL a[N][N];
Mt() { memset(a, 0, sizeof(a)); }
}p, q;
Mt operator* (Mt x, Mt y) {
Mt z;
Rep(i, 0, n)
Rep(k, 0, n)
Rep(j, 0, n) (z.a[i][j] += x.a[i][k] * y.a[k][j]) %= mod;
return z;
}
LL Mult(LL x, LL y) {
LL z = 0;
while (y) {
if (y&1) (z += x) %= mod;
(x += x) %= mod, y >>= 1;
}
return z;
}
class TrickyInequality {
public:
void Pow(int y) {
while (y) {
if (y & 1) q = q * p;
p = p * p, y >>= 1;
}
}
LL pow2(LL x, int y) {
LL z = 1;
while(y) {
if (y&1) z = z * x % mod;
x = x * x % mod, y >>= 1;
}
return z;
}
int countSolutions(long long s, int t, int n0, int m) {
n = m - n0;
h[0][0] = 1;
Rep(i, 1, n) {
h[i][0] = Mult(h[i - 1][0], (s - i + 1));
Rep(j, 1, n) {
h[i][j] = ((Mult(h[i - 1][j], (s - i + 1)) + h[i - 1][j - 1] * (-1)) % mod + mod) % mod;
}
}
Rep(i, 1, t) {
pw[i][0] = 1;
Rep(j, 1, n) pw[i][j] = (LL)pw[i][j - 1] * i % mod;
}
Rep(j, 0, n) {
Rep(i, 1, t) (pw[i][j] += pw[i - 1][j]) %= mod;
} // we will use pw[t][ .. ]
C[0][0] = 1;
Rep(i, 1, n) {
C[i][0] = 1;
Rep(j, 1, i) C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % mod;
}
Rep(i, 0, n) {
Rep(j, 0, i) p.a[i - j][i] = C[i][j] * (LL)pw[t][j] % mod;
q.a[0][i] = pw[t][i];
}
Pow(n0 - 1);
Rep(i, 0, n) {
(ans += h[n][i] * q.a[0][i]) %= mod;
}
LL p0 = 1;
Rep(i, 1, n) (p0 *= i) %= mod;
return int(ans * pow2(p0, mod - 2) % mod);
}
};