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