标签:+= can span 线性变换 就是 strong 开始 代码 inline
引入问题:给定一个对角线非零的上三角矩阵\(M\),求\(M^k\),满足\(M\)的阶\(\le 500\),\(k\le 10^9\)。
对998244353取模。
一个显而易见的算法是矩阵快速幂,然而是\(O(N^3\log k)\)的,无法通过本题。
一开始我想,既然是上三角矩阵,那么特征多项式一定不难求,那么是用CH定理+FFT多项式取模啥搞搞?
然而我naive了。
这题我们可以把\(M\)特征值分解为\(Q^{-1}AQ\)形式,其中\(A\)是一个对角矩阵。
那么\(M^k=(Q^{-1}AQ)^k=Q^{-1}A^kQ\)。
对一个对角矩阵进行幂的复杂度是\(N\log C\)的,矩阵乘法的复杂度是\(O(N^3)\)的,对一个上三角矩阵进行特征值分解可以使用高斯消元,时间复杂度也是\(O(N^3)\),具体怎么对上三角矩阵进行特征值分解??我tm怎么知道,这个得好好研究一下
upd
自己手推没推出来。观察了下std,手跑了下样例,得出来一些性质。
矩阵\(Q^{-1}\)的第\(i\)列,即为矩阵\(M\)对应第\(i\)行第\(i\)列特征值的特征向量。
这个性质通过特征值分解那套理论也不难得到--因为特征向量是\(M\)所对应“方向不变”的向量,而\(Q\)和\(Q^{-1}\)就是在这些旋转方向上的向量,通过线性变换把它们旋转过去//线代那套理论太玄学
std里在\(Q^{-1}\)上递推的,没有看得非常透彻(不过大致也观察出了一些什么),目前已经观察得比较透彻了。
求出\(Q^{-1}\)后直接上矩阵求逆板子求\(Q\),然后直接矩阵乘法就行了。
代码如下:
#include <cstdio>
using namespace std;
const int p = 998244353;
int n, k, a[500][500], l[500][500], r[500][500], v[500];
int qpow(int x, int y)
{
int res = 1;
for (x %= p; y > 0; y >>= 1, x = x * (long long)x % p)
if (y & 1) res = res * (long long)x % p;
return res;
}
int main()
{
scanf("%d%d", &n, &k);
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
scanf("%d", &a[i][j]);
//求l
for (int i = 0; i < n; i++) //枚举l矩阵的第i列,为a矩阵对应于aii的特征向量
{
l[i][i] = 1; //钦定的
for (int j = i - 1; j >= 0; j--) //求这个特征向量的第j行的值
{
int sum = 0; //这里需要满足的是\sum_{k=0}^{n-1}a_{j,k}*l{k,i}=0,此时求值为$b_{ji}$
for (int k = j + 1; k <= i; k++)
sum = (sum + a[j][k] * (long long)l[k][i]) % p;
l[j][i] = sum * (long long)qpow((a[i][i] - a[j][j] + p) % p, p - 2) % p;
//注意这里是a[i][i] - a[j][j], 相当于乘了个-1,就是我们要求的值了
}
}
//求l的逆矩阵r,注意到l是上三角矩阵
for (int i = n - 1; i >= 0; i--)
{
r[i][i] = qpow(l[i][i], p - 2);
for (int j = 0; j < i; j++)
{
int e = l[j][i] * (long long)qpow(l[j][j], p - 2) % p;
for (int k = i; k < n; k++)
r[j][k] = ((r[j][k] - r[i][k] * (long long)e % p) % p + p) % p;
}
}
//收集答案
for (int i = 0; i < n; i++) v[i] = qpow(a[i][i], k);
long long ans1 = 0, ans2 = 0;
for (int i = 0; i < n; i++)
for (int j = i; j < n; j++)
{
int sb = 0;
for (int k = i; k <= j; k++)
sb = (sb + l[i][k] * (long long)v[k] % p * r[k][j] % p) % p;
ans1 += sb, ans2 ^= sb;
}
printf("%lld %lld\n", ans1, ans2);
return 0;
}
以后再研究下一般矩阵的特征值分解,就可以弄图片压缩啥的了。
这个代码常数略大,本来可以弄小一点的
把最后收集答案时候先让对角矩阵和l乘一下再收集、以及优化一下(x%p+p)%p那部分即可。
标签:+= can span 线性变换 就是 strong 开始 代码 inline
原文地址:https://www.cnblogs.com/oier/p/10291286.html