【BZOJ4162shlw loves matrix II】矩阵特征多项式+拉格朗日插值解矩阵快速幂

本来矩阵快速幂应该是k^3logn的，当n足够大的时候，就有必要有一种比较高效的算法将算法时间复杂度更多的偏移到k身上了。在初步学习了线性代数和多项式之后，相信机智的你一定可以把时间复杂度变成k^4 + n^2logn。    

4162: shlw loves matrix II

Time Limit: 30 Sec  Memory Limit: 128 MB Submit: 258  Solved: 89 [Submit][Status][Discuss]

Description

给定矩阵 M，请计算 M^n，并将其中每一个元素对 1000000007 取模输出。

Input

第 1 行包含两个整数 n,k，其中 n 使用二进制表示，可能含有前导零;

余下 k 行描述了一个 k * k 的矩阵 M。

Output

输出题目描述中要求的矩阵，格式同输入。

Sample Input

010 3 5 9 5 5 4 0 8 8 8

Sample Output

110 121 65 45 61 25 144 168 104

HINT

对于 100% 数据，满足 n <= 2^10000;k <= 50; 0 <= Mij < 10^9 +7

Source

By submittersubmitter

引入-拉格朗日插值

现在有n+1个点对(xi,yi)，我们需要构造一个n次多项式F，使得F(xi) = yi。 那么有如下公式 [latex]  f(x)=\sum_{i=0}^{n-1}y_i\prod_{j=0,j\ne i}^{n-1}\frac{x-x_j}{x_i-x_j}　[/latex]

相比于高斯消元n^3,显然这是n^2的。当然，在拉格朗日插值的基础上进行一些骚操作，也可以做到nlog^2n，（即多项式多点插值）。

回到正题

我们已知对于一个矩阵A,若有非０行（列）向量y，和常数x，使得Ay = xy，那么就称y为A的特征向量，x为Ａ的特征值。我们移一个项，使得|A-Ex|y=0 (E为单位矩阵)，由于ｙ行列式不为０，那么一定|A-Ex| = 0 ，同时我们称f(x) = |A-Ex|为矩阵Ａ的特征多项式。由于某个神奇的定理（Cayley-Hamilton定理），f(A) = 0。（啥？你说为什么Ａ是一个矩阵而不是一个常数？我们其实把这个多项式f的ｘ替换成矩阵变量就可以了，感性理解的话|A-A*E| = 0) 那么，如果我们带入k个变量进去，然后k^3高斯消元得到矩阵Ａ的特征多项式之后，由于最终我们要求的是矩阵Ａ的n次幂，即A^n，那么我们对于A^n减去若干个特征多项式f最终的值应该都不会变，那么我们就可以对A^n（把A看做多项式的x)对于特征多项式取模，这样就可以将A^n用一个k-1次多项式表示出来。这样我们就得到了一个关于矩阵的k-1次多项式，把他们加起来就得到了A^n的矩阵了。 由于其时间瓶颈在于k^4带入k+1个值进去然后k^3爆算行列式子，所以显得后面多项式的操作就可以写得十分暴躁了，就可以直接n^2乘法，n^2取模。感觉为了算一个矩阵幂挺拼的orz。 code:

#include<stdio.h>
#include<algorithm>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<iostream>

using namespace std;
const int mod = 1e9+7;
int add(int x,int y) { x+=y; return x>=mod?x-mod:x; }
int mul(int x,int y) { return 1ll*x*y%mod; }
int sub(int x,int y) { x-=y; return x<0?x+mod:x; }
int k;
int ksm(int a,int b) {
    int ans = 1; a%=mod;
    for(;b;b>>=1,a=mul(a,a))
        if(b&1) ans = mul(ans,a);
    return ans;
}
struct poly{
    int o[205];
}F,kp,ANS,X;

poly operator^(const poly&a,int x){//插值(x+xi)
    poly tp = kp;
    for(int i=(k<<1);i;i--)
        tp.o[i]=add(a.o[i-1],mul(x,a.o[i]));
    tp.o[0] = mul(a.o[0],x);
    return tp;
}
poly operator*(const poly&a,int x) {
    poly tp = kp;
    for(int i=0;i<=(k<<1);i++) 
        tp.o[i] = mul(a.o[i],x);
    return tp;
}
poly operator*(const poly&aa,const poly&bb) {
    poly tp = kp;
    for(int i=0;i<=k;i++) {
        for(int j=0;j<=k;j++) {
            tp.o[i+j] = add(tp.o[i+j],mul(aa.o[i],bb.o[j]));
        }
    }
    return tp;
}
poly operator+(poly aa,poly bb) {
    poly tp = kp;
    for(int i=0;i<=(k<<1);i++) tp.o[i] = add(aa.o[i],bb.o[i]);
    return tp;
}
poly operator%(const poly &aa,const poly&bb) {
    poly tp = aa;
    for(int i=k;i>=0;i--) {
        int t = mul(tp.o[i+k],ksm(bb.o[k],mod-2));
        for(int j=0;j<=k;j++) {
            tp.o[i+j] = sub(tp.o[i+j],mul(bb.o[j],t));
        }
    }
    return tp;
}
struct mat{
    int o[55][55];
}a,b,kong,ans;
mat operator*(mat aa,mat bb) {
    mat tmp = kong;
    for(int x=1;x<=k;x++) {
        for(int y=1;y<=k;y++) {
            for(int z=1;z<=k;z++) {
                tmp.o[x][z] = add(tmp.o[x][z],mul(aa.o[x][y],bb.o[y][z]));
            }
        }
    }
    return tmp;
}
int calc(mat oo) {
    int mk = 1; int ans = 1;
    for(int i=1;i<=k;i++) {
        int id = -1;
        for(int j=i;j<=k;j++) {
            if(oo.o[j][i]) { id = j; break; }
        }
        if(id==-1) return 0;
        if(id!=i) {
            mk^=1;
            for(int j=i;j<=k;j++) swap(oo.o[i][j],oo.o[id][j]);
        }
        int iv = ksm(oo.o[i][i],mod-2);
        for(int j=i+1;j<=k;j++) {
            if(oo.o[j][i]) {
                int tt = mul(oo.o[j][i],iv);
                for(int z=1;z<=k;z++) {
                    oo.o[j][z] = sub(oo.o[j][z],mul(oo.o[i][z],tt));
                }
            }
        }
        ans = mul(ans,oo.o[i][i]);
    }
    return mk?ans:(mod-ans)%mod;
}
char ss[100005];
int main() {
    scanf("%s%d",&ss[1],&k);
    for(int i=1;i<=k;i++) {
        for(int j=1;j<=k;j++) {
            int x; scanf("%d",&x);
            a.o[i][j] = x;
        }
    }
    for(int i=0;i<=k;i++) {
        b = a; poly tmp = kp;
        tmp.o[0] = 1;
        for(int j=1;j<=k;j++) b.o[j][j] = add(b.o[j][j],mod-i);
        int y = calc(b);
        for(int j=0;j<=k;j++) {
            if(i==j) continue;
            tmp = tmp^( (mod-j)%mod );
            tmp = tmp*ksm((i-j+mod)%mod,mod-2);     
        }
        tmp = tmp*y;
        F = F + tmp;
    }
    ANS.o[0]=1; X.o[1] = 1;
    int slen = strlen(ss+1);
    for(int i=slen;i;i--) {
        if(ss[i]=='1') ANS = ANS*X%F;
        X = X*X%F;
    }

    b = kong;
    for(int i=1;i<=k;i++) b.o[i][i] = 1;
    for(int t=0;t<k;t++,b=b*a) {
        for(int i=1;i<=k;i++) {
            for(int j=1;j<=k;j++) {
                ans.o[i][j] = add(ans.o[i][j],mul(ANS.o[t],b.o[i][j]));
            }
        }
    }
    for(int i=1;i<=k;i++) {
        for(int j=1;j<=k;j++) {
            printf("%d ",ans.o[i][j]);
        }
        puts("");
    }
}