【TCO 2013 3A】TrichyInequality

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);}};