【USACO 2020 January Platinum】Non-Decreasing Subsequences （前缀积+矩阵的逆）

https://gmoj.net/senior/#main/show/6557

刚好前段时间学了矩阵逆的板，拿这题实践一下。

当然这个矩阵比较简单，好像是可以直接手推矩阵的逆是什么。

时间复杂度：\(O(k^4+n*k^3+Q*k^2)\)

Code

#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, _b = y; i <= _b; i ++)
#define ff(i, x, y) for(int i = x, _b = y; i <  _b; i ++)
#define fd(i, x, y) for(int i = x, _b = y; i >= _b; i --)
#define ll long long
#define pp printf
#define hh pp("\n")
using namespace std;

const int mo = 1e9 + 7;

ll ksm(ll x, ll y) {
	ll s = 1;
	for(; y; y /= 2, x = x * x % mo)
		if(y & 1) s = s * x % mo;
	return s;
}

struct P {
	int a[21][21];
	P() { memset(a, 0, sizeof a);}
} I;

int w;

P operator * (P a, P b) {
	P c = P();
	fo(k, 1, w) fo(i, 1, k) if(a.a[i][k])
		fo(j, k, w) c.a[i][j] = (c.a[i][j] + (ll) a.a[i][k] * b.a[k][j]) % mo;
	return c;
}

const int N = 50005;

int n, a[N];

P b[21], c[21];

P inv(P a) {
	P b = I;	
	fo(i, 1, w) {
		int u = -1;
		fo(j, i, w) if(a.a[i][i]) { u = j; break;}
		fo(j, 1, w) swap(a.a[i][j], a.a[u][j]), swap(b.a[i][j], b.a[u][j]);
		ll x = ksm(a.a[i][i], mo - 2);
		fo(j, 1, w) a.a[i][j] = (ll) a.a[i][j] * x % mo, b.a[i][j] = (ll) b.a[i][j] * x % mo;
		fo(j, 1, w) if(i != j && a.a[j][i]) {
			x = -a.a[j][i];
			fo(k, 1, w) {
				a.a[j][k] = (a.a[j][k] + (ll) a.a[i][k] * x) % mo;
				b.a[j][k] = (b.a[j][k] + (ll) b.a[i][k] * x) % mo;
			}
		}
	}
	return b;
}

P s1[N], s2[N];

int main() {
	scanf("%d %d", &n, &w);
	fo(i, 1, w) I.a[i][i] = 1;
	fo(i, 1, w) {
		fo(j, 1, w) b[i].a[j][j] = 1;
		fo(j, 1, i) b[i].a[j][i] ++;
		c[i] = inv(b[i]);
	}
	s1[0] = s2[0] = I;
	fo(i, 1, n) {
		scanf("%d", &a[i]);
		s1[i] = s1[i - 1] * b[a[i]];
		s2[i] = c[a[i]] * s2[i - 1];
	}
	int Q, x, y;
	scanf("%d", &Q);
	fo(ii, 1, Q) {
		scanf("%d %d", &x, &y);
		ll s = 0;
		fo(k, 1, w) fo(j, 1, w)
			s = (s + (ll) s2[x - 1].a[1][k] * s1[y].a[k][j]) % mo;
		s = (s % mo + mo) % mo;
		pp("%lld\n", s);
	}
}