常系数齐次线性递推

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

概述

形式：

\[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;

}

巴特西

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

常系数齐次线性递推

概述

特征多项式

两种形式

代码

最新文章

热门文章