常系数齐次线性递推优化矩阵快速幂

news/2024/11/18 0:17:28/

一般矩阵快速幂的形式 :
f ( n ) = ∑ i = 1 k 1 a i f ( n − i ) + ∑ i = 1 k 2 b i g ( n − i ) + c f(n)=\sum_{i=1}^{k_1} a_if(n-i)+\sum_{i=1}^{k_2}b_ig(n-i)+c f(n)=i=1k1aif(ni)+i=1k2big(ni)+c
可以做到 k 3 log ⁡ n k^3 \log n k3logn的递推,在某些特殊转移下,有更加快速的方法。

若有如下转移
f ( n ) = ∑ i = 1 k a i f ( n − i ) f(n)=\sum_{i=1}^k a_if(n-i) f(n)=i=1kaif(ni)
这就是常系数线性递推的基本形式。一般可以做到 k log ⁡ k log ⁡ n k\log k \log n klogklogn,比较简单的写法也可以做到 k 2 log ⁡ n k^2 \log n k2logn,下面介绍如何优化。

##特征多项式
前置技能:矩阵的秩

特征值,特征向量:
若有常数 λ \lambda λ,向量 v ⃗ \vec v v ,满足 λ v ⃗ = A v ⃗ \lambda\vec v=A\vec v λv =Av ,则称向量 v ⃗ \vec v v 为矩阵 A A A的一组特征向量,$\lambda 为 矩 阵 为矩阵 A$的一组特征值。
秩为k的矩阵有k组线性不相关的特征向量。

特征多项式
对关系式进行变换: ( λ I − A ) v ⃗ = 0 (\lambda I-A)\vec v=0 (λIA)v =0
这个等式有非零解的充要条件是 d e t ( λ I − A ) = 0 det(\lambda I-A)=0 det(λIA)=0,即为矩阵不满秩。
可以得到一个 k k k次多项式,这个多项式的值等于0时把这个方程称为矩阵 A A A特征方程。这个多项式称为矩阵 A A A特征多项式
k个解自然就是 k k k个特征值辣(并不需要解出来,怎么处理后面会说)!

跟据算数基本定理,最终的多项式有k个解,可以写作 f ( x ) = Π i ( λ i − x ) f(x)=\Pi_i(\lambda_i-x) f(x)=Πi(λix),而且带入 A A A矩阵本身会得到 0 0 0矩阵(Cayley-Hamilton 定理)。

证明:
(下面的证明是博主以前学OI的时候写的一种比较复杂的证明,其实可以直接用伴随矩阵来证明)
f ( A ) = Π i ( λ i I − A ) f(A)=\Pi_i(\lambda_i I -A) f(A)=Πi(λiIA)
考虑将 f ( A ) f(A) f(A)得到的矩阵分别乘特征向量,如果最后都得到了0矩阵,因为这几个特征向量线性不相关,那么可证明 f ( A ) f(A) f(A)乘以任意向量都会得到 0 0 0矩阵,从而 A A A 0 0 0矩阵。

现在问题转化为证明 f ( A ) f(A) f(A)乘任意特征向量值得到 0 0 0矩阵。首先证明 ( λ i I − A ) ( λ j I − A ) = ( λ j I − A ) ( λ i I − A ) (\lambda_i I-A)(\lambda_j I - A)=(\lambda_j I - A)(\lambda_i I-A) (λiIA)(λjIA)=(λjIA)(λiIA),这个直接展开就好了,因为矩阵满足加法交换律。证毕。

那么假设当前乘的特征向量为 v ⃗ a \vec v_a v a,有:
f ( A ) v ⃗ a = Π i ( λ i I − A ) ∗ v ⃗ a = Π i ! = a ( λ i I − A ) ∗ ( λ a I − A ) ∗ v ⃗ a \begin{aligned}f(A)\vec v_a&=\Pi_{i}(\lambda_iI-A)*\vec v_a\\&=\Pi_{i!=a}(\lambda_iI-A)*(\lambda_aI-A)*\vec v_a &&\end{aligned} f(A)v a=Πi(λiIA)v a=Πi!=a(λiIA)(λaIA)v a

而矩阵乘法满足结合律,对于后半部分,展开得:
( λ a I − A ) ∗ v ⃗ a = λ a I v ⃗ a − A v ⃗ a (\lambda_aI -A)*\vec v_a=\lambda_aI\vec v_a-A\vec v_a (λaIA)v a=λaIv aAv a

根据定义, λ a I v ⃗ a − A v ⃗ a = 0 \lambda_aI\vec v_a-A\vec v_a=0 λaIv aAv a=0,那么代回原式可得 0 0 0矩阵。

它同时是一个关于 A A A k k k次多项式。那么假设我们求 A w A^w Aw,可以做多项式取模,模掉多余的0矩阵得到一个关于A的 k − 1 k-1 k1次多项式: g ( A ) = A w % f ( A ) = ∑ i = 0 k − 1 a i A i g(A)=A^w\% f(A) =\sum_{i=0}^{k-1}a_iA^i g(A)=Aw%f(A)=i=0k1aiAi

设初始递推矩阵为 h = h k , h k − 1 , h k − 2 , . . . , h 1 h={h_k,h_{k-1},h_{k-2},...,h_1} h=hk,hk1,hk2,...,h1,那么答案为:
∑ a i ∗ [ ( h ∗ A i ) ] 11 \sum a_i*[(h*A^i)]_{11} ai[(hAi)]11,而 h ∗ A i h*A^i hAi的第一项就是 h i + 1 h_{i+1} hi+1,那么可以直接 O ( n ) O(n) O(n)计算了。因为多项式取模可以用快速傅里叶变换做到 k log ⁡ k k\log k klogk,加上快速幂的时间复杂度,总复杂度为 O ( k log ⁡ k log ⁡ n ) O(k\log k\log n) O(klogklogn)。为了方便,可以将多项式取模写 k 2 k^2 k2算法,那么总复杂度是 O ( k 2 log ⁡ n ) O(k^2\log n) O(k2logn)的。

##求解特征多项式
要做多项式取模的前提是要知道多项式。显然解出 k k k个特征值是不现实的,考虑到常系数齐次递推的特点,可以快速求出矩阵的特征多项式:
首先,列出转移方程:

( a 1 a 2 a 3 ⋯ a k − 2 a k − 1 a k 1 0 0 ⋯ 0 0 0 0 1 0 ⋯ 0 0 0 0 0 1 ⋯ 0 0 0 ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ 0 0 0 ⋯ 1 0 0 0 0 0 ⋯ 0 1 0 ) ( h k − 1 h k − 2 h k − 3 ⋮ h 2 h 1 h 0 ) = ( h k h k − 1 h k − 2 ⋮ h 3 h 2 h 1 ) \begin{pmatrix} a_1 & a_2 & a_3 & \cdots & a_{k - 2} & a_{k - 1} & a_k \\ 1 & 0 & 0 & \cdots & 0 & 0 & 0 \\ 0 & 1 & 0 & \cdots & 0 & 0 & 0 \\ 0 & 0 & 1 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 1 & 0 & 0 \\ 0 & 0& 0 & \cdots & 0 & 1 & 0 \end{pmatrix} \begin{pmatrix}h_{k - 1} \\ h_{k - 2} \\ h_{k - 3} \\ \vdots \\ h_2 \\ h_1 \\ h_0 \end{pmatrix} = \begin{pmatrix}h_{k} \\ h_{k - 1} \\ h_{k - 2} \\ \vdots \\ h_3 \\ h_2 \\ h_1 \end{pmatrix} a110000a201000a300100ak200010ak100001ak00000hk1hk2hk3h2h1h0=hkhk1hk2h3h2h1

特征多项式:
f ( λ ) = ∣ λ I − A ∣ = ( λ − a 1 − a 2 − a 3 ⋯ − a k − 2 − a k − 1 − a k − 1 λ 0 ⋯ 0 0 0 0 − 1 λ ⋯ 0 0 0 0 0 − 1 ⋯ 0 0 0 ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ 0 0 0 ⋯ − 1 λ 0 0 0 0 ⋯ 0 − 1 λ ) f(\lambda) = |\lambda I - A| = \begin{pmatrix} \lambda - a_1 & -a_2 & -a_3 & \cdots & -a_{k - 2} & -a_{k - 1} & -a_k \\ -1 & \lambda & 0 & \cdots & 0 & 0 & 0 \\ 0 & -1 & \lambda & \cdots & 0 & 0 & 0 \\ 0 & 0 & -1 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -1 & \lambda & 0 \\ 0 & 0& 0 & \cdots & 0 & -1 & \lambda \end{pmatrix} f(λ)=λIA=λa110000a2λ1000a30λ100ak200010ak1000λ1ak0000λ

考虑对第一行展开分别求和,有:

f ( λ ) = ( λ − a 1 ) A 11 + ( − a 2 ) A 12 + ⋯ + ( − a k ) A 1 n = λ k − a 1 λ k − 1 − a 2 λ k − 2 − ⋯ − a k f(\lambda) = (\lambda - a_1)A_{11} + (-a_2)A_{12} + \cdots + (-a_k)A_{1n} = \lambda ^ k - a_1 \lambda ^ {k - 1} - a_2 \lambda ^ {k - 2} - \cdots - a_k f(λ)=(λa1)A11+(a2)A12++(ak)A1n=λka1λk1a2λk2ak

其中 A 1 i A_{1i} A1i A A A的代数余子式。

至此,常系数齐次线性递推也就可以优化矩阵快速幂了。

##小优化
1.是不是觉得计算出结果后很麻烦还要一项一项相乘??不用!!
直接计算 A n A^{n} An而不是 A n − k A^{n-k} Ank那么最后结果会是向量的最后一位。此时 A i ( i ≤ k ) A^{i}(i\le k) Ai(ik)的最后一位是原递推式的一项,可以直接求取。

2.注意 k 2 k^2 k2计算并取模时枚举的顺序,依次枚举 i , j i,j i,j加到 i + j i+j i+j上而不是依次枚举 i , j i,j i,j, a i = a j ∗ a i − j a_i=a_j*a_{i-j} ai=ajaij

附上BZOJ4161的代码,这是一道常系数齐次线性递推的裸题。可以用 k 2 log ⁡ n k^2\log n k2logn的方法过。

#include<bits/stdc++.h>
using namespace std;
inline int read(){char ch=getchar();int i=0,f=1;while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}while(isdigit(ch)){i=(i<<1)+(i<<3)+ch-'0';ch=getchar();}return i*f;
}
const int Maxn=2e3+20,Mod=1000000007;
#define rg register 
int n,k,a[Maxn<<1],b[Maxn<<1],c[Maxn<<1],h[Maxn],tp[Maxn<<1],ans,lim;
inline void add(int &x,int y){x+=y;(x>Mod?x-=Mod:(x<0?x+=Mod:0));
}
inline void mul(int *A,int *B){for(rg int i=0;i<=lim;i++)tp[i]=0;for(rg int i=0;i<=k;i++)for(rg int j=0;j<=k;j++)tp[i+j]=(tp[i+j]+1ll*A[i]*B[j])%Mod;for(rg int i=lim;i>=k;i--){for(rg int j=0;j<k;j++)tp[i-k+j]=(tp[i-k+j]+1ll*tp[i]*a[k-j])%Mod;tp[i]=0;}for(int i=0;i<=lim;i++)A[i]=tp[i];
}
inline void power(int *A,int B,int *C){for(;B;B>>=1,mul(A,A))if(B&1)mul(C,A);
}
int main(){n=read(),k=read();lim=k<<1;for(int i=1;i<=k;i++)a[i]=read();for(int i=1;i<=k;i++)h[i]=(read()%Mod+Mod)%Mod;if(n<k){printf("%d\n",h[n+1]);return 0;}b[0]=1;c[1]=1;power(c,n,b);for(int i=0;i<k;i++)add(ans,1ll*b[i]*h[i+1]%Mod);printf("%d\n",ans);
}

http://www.ppmy.cn/news/244047.html

相关文章

我的电流测试板(一)

最近做了块电流板&#xff0c;这是本人第一次独立的负责软件和硬件&#xff0c;本人以前做的产品都是比较大的系统&#xff0c;双层板&#xff0c;类似于万用表&#xff0c;可以测试电流&#xff0c;不过跟万用表不同的是我这块板子可以累加电流&#xff0c;可以测试其他电源。…

jtag和swd有什么不同_jtag和swd的区别

本文为大家介绍jtag和swd的区别。 jtag和swd有什么不同----引脚对比 对于仿真ARM&#xff0c;TKScope仿真器家庭的AK100/AK100Pro、K8/K9等仿真器提供标准的20PIN调试接口。接口管脚定义如下。其中ARM芯片有两种调试模式&#xff0c;一种是JTAG&#xff0c;一种是SWD&#xf…

5v继电器模块实物接线_高手教你玩传感器系列之继电器的使用

原标题:高手教你玩传感器系列之继电器的使用 单片机、嵌入式系统等是一个弱电器件,一般情况下它们大都工作在5V甚至更低,驱动电流在mA级以下。而要把它用于一些大功率场合,比如控制电动机,显然是不行的。所以,就要有一个环节来衔接,这个环节就是所谓的”功率驱动”。继电…

单片机的四种烧写方式

参考&#xff1a;单片机的四种烧写方式 作者&#xff1a;爱学习的小王呀 发布时间&#xff1a;2020-11-27 20:05:12 网址&#xff1a;https://blog.csdn.net/hongliwong/article/details/110245095?spm1001.2014.3001.5501 参考&#xff1a;单片机3种烧录方式解析 作者&#x…

单片机3种烧录方式解析

内容包括ISP、IAP、ICP三种烧录方式的详细介绍&#xff0c;STM32单片机与宏晶STC单片机烧录方法&#xff0c;STM32单片机自动ISP的详细介绍(附电路原理图)。紫色文字是超链接&#xff0c;点击自动跳转至相关博文。持续更新&#xff0c;原创不易&#xff01; 目录&#xff1a; 一…

ISP、IAP、ICP、JTAG、SWD的编程特点

电子工程师都知道&#xff0c;半导体技术发展迅猛&#xff0c;带动了各种芯片技术的不断升级。在数据存储方面&#xff0c;从最初的掩膜ROM&#xff0c;发展到现在的Flash技术&#xff0c;存储技术的不断改进&#xff0c;相对应的编程技术也在不断发展。 记得老一辈工程师在烧…

5.11 gethostname() ---- 我是谁?

原文&#xff1a;https://beej.us/guide/bgnet/html/#gethostnamewho-am-i 5.11 gethostname() ---- 我是谁&#xff1f; 甚至比getpeername() 还简单的函数是 gethostname()。它返回本机名。 然后可以使用 gethostbyname() 以获得本机 IP 地址。 下面是定义&#xff1a; #in…

JTAG、JLink、ULINK、ST-LINK仿真器区别

首先要了解一下JTAG。 JTAG协议 JTAG&#xff08;Joint Test Action Group&#xff0c;联合测试行动小组&#xff09;是一种国际标准测试协议&#xff08;IEEE 1149.1兼容&#xff09;&#xff0c;主要用于芯片内部测试。现在多数的高级器件都支持JTAG协议&#xff0c;如ARM、D…