BM（Berlekamp-Massey）算法学习笔记

线性递推的题目区域赛里还是挺多的，还是有必要学一下

~ BM（Berlekamp-Massey）算法 ~

有一个$n$阶线性递推$f$，想要计算$f(m)$，有一种常用的办法是矩阵快速幂，复杂度是$O(n^3logm)$

在不少情况下这已经够用了，但是如果$n$比较大、到了$10^3$级别，这就不太适用了

而BM算法能将这个复杂度压低到$O(n^2logm)$，若加上NTT优化的话能做到$O(n^2+nlognlogm)$，十分厉害

这个算法的核心是将$f(m)$用递推的前$n$项表示

即，已知$f(0),...,f(n-1)$和递推式$f(m)=a_0f(m-1)+...+a_{n-1}f(m-n)$，该算法是求出系数$W_0,...,W_{n-1}$，使得$f(m)=W_0f(n-1)+...+W_{n-1}f(0)$

看似无从下手？实际上只要大力展开就行了

根据定义，有（只是写成$\sum$的形式而已）

\[f(m)=\sum_{i=0}^{n-1}a_i f(m-1-i)\]

而对于每一项再次展开，即

\[f(m-1-i)=\sum_{j=0}^{n-1}a_j f(m-1-i-1-j)\]

全部代入，能得到

\[f(m)=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_ia_j f(m-2-i-j)\]

把式子写的更好看一点，就是

\[f(m)=\sum_{k=0}^{2n-2}\sum_{i+j=k}a_ia_j f(m-2-k)\]

这样做之后有什么用呢？

在原本的递推式中，$f(m)$可以通过$f(m-1),...,f(m-n)$这$n$个项表示

各项展开后，就可以通过$f(m-2),...,f(m-2n)$表示

事实上，我们可以再依次对$f(m-i),2\leq i\leq n$展开，并将系数向$f(m-i-1),...,f(m-i-n)$并入，最终就能把原递推式通过$f(m-n-1),...,f(m-2n)$这$n$项表示

于是可以得到一个新的$n$阶递推式，记为$f(m)=b_0f(m-n+1),...,b_{n-1}f(m-2n)$

再用新递推式将各项展开，就可以通过$f(m-2n-2),...,f(m-4n)$表示

再用原递推式展开$f(m-2n-i),2\leq i\leq n$并向前合并系数，最终就能把原递推式通过$f(m-3n+1),...,f(m-4n)$这$n$项表示

之后都是类似的了，不再赘述

有了上面的思路，就可以用类似快速幂的方法，得到$f(m)=W_0f(m-(k-1)n+1),...,W_{n-1}f(m-kn)$这样的展开式，其中$m-kn<n$

余数$m-kn$是我们不喜欢的，但也没有必要整体再向前推，一开始计算时算出$f(0),...,f(2n-1)$就够了

按照上述思路能这样实现：

#include <cstdio>
#include <cstring>
using namespace std;typedef long long ll;
const int MOD=1000000007;
const int N=1005;int n,m;
int a[N];
int f[N<<1];int tmp[N<<1];void mul(int *y,int *x)
{memset(tmp,0,sizeof(tmp));for(int i=0;i<n;i++)for(int j=0;j<n;j++)tmp[i+j]=(tmp[i+j]+ll(y[i])*x[j])%MOD;for(int i=0;i<n-1;i++)for(int j=0;j<n;j++)tmp[i+j+1]=(tmp[i+j+1]+ll(tmp[i])*a[j])%MOD;for(int i=0;i<n;i++)y[i]=tmp[i+n-1];
}int w[N<<1],x[N<<1];int BM()
{if(m<(n<<1))return f[m];for(int i=0;i<n;i++)x[i]=a[i],w[i]=a[i];int t=(m-n)/n;int rem=m-n-t*n;while(t){if(t&1)mul(w,x);mul(x,x);t>>=1;}int res=0;for(int i=0;i<n;i++)res=(res+ll(w[i])*f[rem+n-i-1])%MOD;return res;
}int main()
{scanf("%d%d",&n,&m);for(int i=0;i<n;i++)scanf("%d",&a[i]);for(int i=0;i<n;i++)scanf("%d",&f[i]);for(int i=n;i<(n<<1);i++)for(int j=1;j<=n;j++)f[i]=(f[i]+ll(a[j-1])*f[n-j])%MOD;printf("%d\n",BM());return 0;
}

View Code

想做的更快的话，一个是要写NTT，另一个是合并系数会比较困难，待补

因为这题学的BM：牛客ACM 882B （$Eddy$ $Walker$ $2$）

$m\rightarrow \infty$时，$f(m)\rightarrow \frac{2}{k+1}$ （并不会证...）

从rls那里学了一个证明：

走$k$步，期望能走的长度是$1+2+...+k=\frac{k(k+1)}{2}$

那么在这段距离中，每个位置被走过的概率就是$\frac{k}{\frac{k(k+1)}{2}}=\frac{2}{k+1}$

在其他时候，直接套上面的板子即可

牛客的玄学评测机，同一份代码能差出500ms = =

#include <cstdio>
#include <cstring>
using namespace std;typedef long long ll;
const int MOD=1000000007;
const int N=1100;inline int quickpow(int x,int t)
{int res=1;while(t){if(t&1)res=ll(res)*x%MOD;x=ll(x)*x%MOD;t>>=1;}return res;
}inline int rev(int x)
{return quickpow(x,MOD-2);
}int n,rn;
ll m;
int a[N];
int f[N<<1];int tmp[N<<1];void mul(int *y,int *x)
{memset(tmp,0,sizeof(tmp));for(int i=0;i<n;i++)for(int j=0;j<n;j++)tmp[i+j]=(tmp[i+j]+ll(y[i])*x[j])%MOD;for(int i=0;i<n-1;i++)for(int j=0;j<n;j++)tmp[i+j+1]=(tmp[i+j+1]+ll(tmp[i])*a[j])%MOD;for(int i=0;i<n;i++)y[i]=tmp[i+n-1];
}int w[N<<1],x[N<<1];int BM()
{if(m<(n<<1))return f[m];for(int i=0;i<n;i++)x[i]=a[i],w[i]=a[i];ll t=(m-n)/n;int rem=m-n-t*n;while(t){if(t&1)mul(w,x);mul(x,x);t>>=1;}int res=0;for(int i=0;i<n;i++)res=(res+ll(w[i])*f[rem+n-i-1])%MOD;return res;
}int main()
{int T;scanf("%d",&T);while(T--){scanf("%d%lld",&n,&m);if(m==-1){printf("%d\n",2LL*rev(n+1)%MOD);continue;}rn=rev(n);for(int i=0;i<n;i++)a[i]=rn;memset(f,0,sizeof(f));f[0]=1;for(int i=1;i<(n<<1);i++)for(int j=1;j<=n && j<=i;j++)f[i]=(f[i]+ll(rn)*f[i-j])%MOD;printf("%d\n",BM());}return 0;
}

View Code

比较特定的知识点吧，以后遇到就是赚到（然后发现强制NTT，直接白给= =）

（完）

转载于:https://www.cnblogs.com/LiuRunky/p/Berlekamp_Massey.html

BM（Berlekamp-Massey）算法学习笔记相关推荐

Berlekamp–Massey算法简要介绍
这是一篇翻译向的文章,笔者整理了一些有关Berlekamp–Massey算法的笔记,还增加了一些自己的理解. 下面列出了笔者写此文时所参考的一些资料: wikipedia fjzzq2002 别人的博 ...
基于MVS的三维重建算法学习笔记（四）— 立体匹配经典算法Semi-Global Matching（SGM）论文翻译及要点解读
基于MVS的三维重建算法学习笔记(四)- 立体匹配经典算法Semi-Global Matching(SGM)论文翻译及要点解读声明 SGM概述 Cost Calculation(像素代价计算)--M ...
大顶堆删除最大值_算法学习笔记(47): 二叉堆
堆(Heap)是一类数据结构,它们拥有树状结构,且能够保证父节点比子节点大(或小).当根节点保存堆中最大值时,称为大根堆:反之,则称为小根堆. 二叉堆(Binary Heap)是最简单.常用的堆,是一 ...
Manacher算法学习笔记 | LeetCode#5
Manacher算法学习笔记 DECLARATION 引用来源:https://www.cnblogs.com/grandyang/p/4475985.html CONTENT 用途:寻找一个字符串的 ...
数据结构与算法学习笔记之从0编号的数组
数据结构与算法学习笔记之从0编号的数组前言数组看似简单,但掌握精髓的却没有多少:他既是编程语言中的数据类型,又是最基础的数据结构: 一个小问题: 为什么数据要从0开始编号,而不是从1开始呢? ...
输出dag的所有拓扑排序序列_算法学习笔记(53): 拓扑排序
拓扑排序是对DAG(有向无环图)上的节点进行排序,使得对于每一条有向边 , 都在之前出现.简单地说,是在不破坏节点先后顺序的前提下,把DAG拉成一条链.如果以游戏中的科技树(虽然名字带树,其实常常 ...
算法学习笔记：对指定金额计算最少钞票数
算法学习笔记:对指定金额计算最少钞票数一.引出问题财务人员给员工发工资时经常遇到这样一个问题,即根据每个人的工资额(以元作为单位)计算出各种面值的钞票的张数,且要求总张数最少.例如,某职工工资为3 ...
matlab中x从0到5不含0,关于MATLAB的数学建模算法学习笔记
关于MATLAB的数学建模算法学习笔记目录线性规划中应用: (3) 非线性规划: (3) 指派问题;投资问题:(0-1问题) (3) 1)应用fmincon命令语句 (3) 2)应用指令函数:bi ...
机器学习篇01：在线学习的支持向量机算法学习笔记
在线学习的支持向量机算法学习笔记 oisvm算法实现说明 oisvm算法实现说明 % 本程序是用于实现基于在线学习的调制信号识别的程序 % % % 第一步:调制信号的生成 % 首先是7个信号:2ASK ...
数据结构与算法学习笔记之提高读取性能的链表（上）
数据结构与算法学习笔记之提高读取性能的链表(上) 前言链表(Linked list)比数组稍微复杂一点,在我们生活中用到最常见的应该是缓存,它是一种提高数据读取性能的技术,常见的如cpu缓存,浏览 ...

BM（Berlekamp-Massey）算法学习笔记

BM（Berlekamp-Massey）算法学习笔记相关推荐

最新文章

热门文章