常系数齐次线性递推 & 拉格朗日插值

常系数齐次线性递推

具体记在笔记本上了，以后可能补照片，这里稍微写一下，主要贴代码。

---

概述

形式：
\[ h_n = a_1 h_{n-1}+a_2h_{n-2}+...+a_kh_{n-k} \]
矩阵乘法是\(O(k^3 \log n)\)

利用特征多项式可以做到\(O(k^2\log n)\)

特征多项式

特征值和特征向量

特征多项式
\[ f(\lambda) = \mid M - \lambda I\mid \]
是关于\(\lambda\)的\(n\)次多项式

根据\(Cayley-hamilton\)定理得到 \(f(M)=0\)(zero matrix)，所以就能得到\(M^k\)与\(M^0...M^{k-1}\)的线性递推关系，\(M^n\)也可以用\(M^0...M^{k-1}\)线性表示，多项式乘法快速幂求这个递推关系的系数。最后代入就行了。

其实就是：
\[ M^n = M^n \mod f(M) \]
所以还需要多项式求余。

这玩意卡我好长时间，然后发现就是手算的原理。当然用毒瘤算法可以做到nlogn

两种形式

对于常系数齐次线性递推关系，我们可以一眼看出它的特征多项式
\[ f(\lambda) = \lambda^k - a_1 \lambda^{k-1} -...-a_k \]
并且最高次系数为1，求余很简单

对于一般的矩阵，需要求行列式得到n+1个点的值，然后拉格朗日插值，这一步复杂度\(O(n^4)\)

拉格朗日插值公式：
\[ \sum_{j=0}^n y_j \prod_{i=0 , i\neq j}^n \frac{x-x_i}{x_j - x_i} \]

---

代码

例题为bzoj 4161 & bzoj 4162

常系数齐次线性递推

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <queue>
using namespace std;
typedef long long ll;
const int N = 4005, mo = 1e9+7;
inline int read(){
    char c=getchar(); int x=0,f=1;
    while(c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9') {x=x*10+c-'0';c=getchar();}
    return x*f;
}

int n, k, a[N], f[N], b[N], c[N], ans;
void mul_mod(int *x, int *y) {
    static int c[N];
    for(int i=0; i<k<<1; i++) c[i] = 0;
    for(int i=0; i<k; i++)
        for(int j=0; j<k; j++) c[i+j] = (c[i+j] + (ll) x[i] * y[j]) %mo;
    // mod f(M)
    for(int i=2*k-2; i>=k; i--)
        for(int j=1; j<=k; j++) c[i-j] = (c[i-j] + (ll) c[i] * a[j]) %mo;
    for(int i=0; i<k; i++) x[i] = c[i];
}

void Pow(int *a, int b, int *ans) {
    for(; b; b>>=1, mul_mod(a, a)) if(b&1) mul_mod(ans, a);
}
int main() {
    freopen("in", "r", stdin);
    n=read(); k=read();
    for(int i=1; i<=k; i++) a[i] = read(), a[i] += a[i] < 0 ? mo : 0;
    for(int i=0; i<k; i++)  f[i] = read(), f[i] += f[i] < 0 ? mo : 0;
    c[1] = 1; b[0] = 1;
    Pow(c, n, b);
    for(int i=0; i<k; i++) ans = (ans + (ll) b[i] * f[i]) %mo;
    printf("%d\n", ans);
}

一般矩阵快速幂

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
#define pll pair<ll, ll>
#define fir first
#define sec second
const int N = 52, mo = 1e9+7;
inline int read(){
    char c=getchar(); int x=0,f=1;
    while(c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9') {x=x*10+c-'0';c=getchar();}
    return x*f;
}

ll Pow(ll a, int b) {
    ll ans = 1;
    for(; b; b>>=1, a=a*a%mo)
        if(b&1) ans=ans*a%mo;
    return ans;
}

int n, k, a[N][N], t[N][N], t2[N][N], ans[N][N];
char s[10005];

void mul(int a[N][N], int b[N][N], int c[N][N]) {
    for(int i=1; i<=n; i++) 
        for(int k=1; k<=n; k++) if(a[i][k])
            for(int j=1; j<=n; j++) if(b[k][j])
                c[i][j] = (c[i][j] + (ll) a[i][k] * b[k][j]) %mo;
}

int det(int a[N][N]) {
    int s = 0;
    for(int i=1; i<=n; i++) {
        int r;
        for(r=i; r<=n; r++) if(a[r][i]) break;
        if(r == n+1) return 0;
        if(r != i) {
            s ^= 1;
            for(int j=1; j<=n; j++) swap(a[r][j], a[i][j]);
        }
        ll inv = Pow(a[i][i], mo-2);
        for(int k=i+1; k<=n; k++) {
            ll t = (ll) (mo - a[k][i]) * inv %mo;
            for(int j=i; j<=n; j++) a[k][j] = (a[k][j] + t * a[i][j]) %mo;
        }
    }
    ll ans = 1;
    for(int i=1; i<=n; i++) ans = ans * a[i][i] %mo;
    return (s&1) ? mo - ans : ans;
}

struct poly {
    int a[N<<1];
    poly(int x=0) {memset(a, 0, sizeof(a)); a[0] = x;}
    int& operator [](int x) {return a[x];}

    poly operator + (poly b) {
        poly c;
        for(int i=0; i<=n; i++) c[i] = (a[i] + b[i]) %mo;
        return c;
    }
    poly operator * (ll b) {
        poly c;
        for(int i=0; i<=n; i++) c[i] = (ll) a[i] * b %mo;
        return c;
    }
    poly operator * (poly b) {
        poly c;
        for(int i=0; i<=n; i++) 
            for(int j=0; j<=n; j++) c[i+j] = (c[i+j] + (ll) a[i] * b[j]) %mo;
        return c;
    }
    poly operator * (pll t) {
        ll k = t.fir, b = t.sec;
        poly c(a[0] * b %mo);
        for(int i=1; i<=n; i++) c[i] = (a[i-1] * k + a[i] * b) %mo;
        return c;
    }
    friend poly operator % (poly a, poly b) {
        for(int i=n; i>=0; i--) {
            ll t = (ll) (mo - a[i+n]) * Pow(b[n], mo-2) %mo;
            for(int j=0; j<=n; j++) a[i+j] = (a[i+j] + b[j] * t) %mo;
        }
        return a;
    }
} ;

int y[N];
int main() {
    freopen("in", "r", stdin);
    scanf("%s",s);
    n=read();
    for(int i=1; i<=n; i++) 
        for(int j=1; j<=n; j++) a[i][j] = read();
    for(int x=0; x<=n; x++) {
        memcpy(t, a, sizeof(a));
        for(int i=1; i<=n; i++) t[i][i] = (t[i][i] + mo - x) %mo;
        y[x] = det(t);
    }
    poly f;
    for(int j=0; j<=n; j++) {
        poly t(1);
        for(int i=0; i<=n; i++) if(i != j) 
            t = (t * make_pair(1, mo-i) ) * Pow(j-i+mo, mo-2);
        t = t * y[j];
        f = f + t;
    }

    poly p, b(1); p[1] = 1;
    for(int i=strlen(s)-1; i>=0; i--, p = p * p % f)
        if(s[i] == '1') b = b * p % f;

    memset(t, 0, sizeof(t));
    for(int i=1; i<=n; i++) t[i][i] = 1;
    for(int p=0; p<n; p++) {
        for(int i=1; i<=n; i++) 
            for(int j=1; j<=n; j++) 
                ans[i][j] = (ans[i][j] + (ll) t[i][j] * b[p]) %mo;
        memset(t2, 0, sizeof(t2));
        mul(t, a, t2);  
        memcpy(t, t2, sizeof(t2));
    }
    for(int i=1; i<=n; i++) 
        for(int j=1; j<=n; j++) printf("%d%c", ans[i][j], " \n"[j==n]);
    return 0;
}