BM算法（Berlekamp-Massey算法）：解决线性递推问题

算法内容

Berlekamp-Massey算法，常简称为BM算法，是用来求解一个数列的最短线性递推式的算法。
BM算法可以在O(N2)的时间内求解一个长度为N的数列的最短线性递推式。

算法模板

取模，模数为质数

12124

#include<bits/stdc++.h>using namespace std;#define rep(i,a,n) for(int i=a;i<n;i++)#define SZ(x) ((int)(x).size())
typedef vector<int> VI;
typedef long long ll;const ll mod = 1e9+7;ll powmod(ll a,ll b) {ll ret = 1;a%=mod;while(b) {if(b&1) ret = ret*a%mod;a = a*a%mod;b>>=1;}return ret;
}
int _,n;namespace linear_seq {const int N = 10010;
ll res[N],base[N],_c[N],_md[N];
vector<int> Md;void mul(ll *a,ll *b,int k) {for(int i=0; i<k+k; i++)_c[i] = 0;for(int i=0; i<k; ++i)if(a[i])for(int j=0; j<k; j++)_c[i+j] = (_c[i+j]+a[i]*b[j])%mod;for(int i=k+k-1; i>=k; i--)if(_c[i])for(int j=0; j<(int) Md.size(); j++)_c[i-k+Md[j]] = (_c[i-k+Md[j]]-_c[i]*_md[Md[j]])%mod;for(int i=0; i<k; i++)a[i] = _c[i];
}int solve(ll n,VI a,VI b) {ll ans = 0,pnt = 0;int k = SZ(a);for(int i=0; i<k; i++)_md[k-1-i] = -a[i];_md[k] = 1;Md.clear();for(int i=0; i<k; i++)if(_md[i]!=0)Md.push_back(i);for(int i=0; i<k; i++)res[i] = base[i] = 0;res[0] = 1;while((1ll<<pnt)<=n)pnt++;for(int p = pnt; p>=0; p--) {mul(res,res,k);if((n>>p)&1) {for(int i=k-1; i>=0; i--) res[i+1] = res[i];res[0]  = 0;for(int j=0; j<(int) Md.size(); ++j)res[Md[j]] = (res[Md[j]]-res[k]*_md[Md[j]])%mod;}}rep(i,0,k) ans = (ans+res[i]*b[i])%mod;if(ans<0) ans+=mod;return ans;
}VI BM(VI s) {VI C(1,1),B(1,1);int L = 0,m = 1,b = 1;for(int n=0; n<(int)s.size(); ++n) {ll d = 0;for(int i=0; i<L+1; i++)d = (d+(ll)C[i]*s[n-i])%mod;if(d==0) ++m;else if(2*L<=n) {VI T=C;ll c = mod-d*powmod(b,mod-2)%mod;while(SZ(C)<SZ(B)+m)C.push_back(0);for(int i=0; i<(int) B.size(); i++)C[i+m] = (C[i+m]+c*B[i])%mod;L = n+1-L;B = T;b = d;m = 1;} else {ll c = mod-d*powmod(b,mod-2)%mod;while(SZ(C)<SZ(B)+m)C.push_back(0);for(int i=0; i<(int)B.size(); i++)C[i+m] = (C[i+m]+c*B[i])%mod;++m;}}return C;
}ll gao(VI a,ll n) {VI c = BM(a);c.erase(c.begin());for(int i=0; i<(int) c.size(); i++)c[i] = (mod-c[i])%mod;return (ll) solve(n,c,VI(a.begin(),a.begin()+SZ(c)));
}}int main() {ll t;ll nnn;VI a;a.push_back(3);a.push_back(9);a.push_back(20);a.push_back(46);a.push_back(106);a.push_back(244);a.push_back(560);a.push_back(1286);a.push_back(2956);a.push_back(6794);scanf("%lld",&t);while(t--) {scanf("%lld",&nnn);printf("%lld\n",linear_seq::gao(a,nnn-1));}
}

模数可不为质数

13704

#include <bits/stdc++.h>
using namespace std;
#define FOR(i,a,b) for(int i(a);i<=(b);++i)
#define FOL(i,a,b) for(int i(a);i>=(b);--i)
#define SZ(x) ((long long)(x).size())
#define REW(a,b) memset(a,b,sizeof(a))const int64_t Mod (1000000000) ;int64_t powEx(int64_t base,int64_t n,int64_t Mod = ::Mod )
{int64_t ret(1);while(n){if(n&1) ret = ret*base%Mod ;base = base * base %Mod ;n >>= 1 ;}return ret % Mod ;
}class Linear_Seq
{using VI = vector<int64_t> ;
public:static const int N=10010;int64_t res[N],base[N],c[N],md[N];vector<int> Md;inline void mulEx(int64_t *a,int64_t *b,int k) {for(int i(0);i<k+k;++i) c[i]=0;for(int i(0);i<k;++i)if(a[i])for(int j(0);j<k;++j)c[i+j]=(c[i+j]+a[i]*b[j])%Mod;for (int i(k+k-1);i>=k;--i) if (c[i])for(int j(0);j<Md.size();++j)c[i-k+Md[j]]=(c[i-k+Md[j]]-c[i]*md[Md[j]])%Mod;copy(c,c+k,a) ;}int solve(int64_t n,VI a,VI b) { // a 系数 b 初值 b[n+1]=a[0]*b[n]+...int64_t ans(0),cnt(0);int k(a.size());for(int i(0);i<k;++i) md[k-1-i]=-a[i];md[k]=1 ;  Md.clear() ;for(int i(0);i<k;++i) if (md[i]) Md.push_back(i);for(int i(0);i<k;++i) res[i] = base[i] = 0;res[0]=1;while ((1LL<<cnt)<=n) ++ cnt;for (int p(cnt);~p;-- p) {mulEx(res,res,k);if ((n>>p)&1) {for (int i(k-1);~i;--i) res[i+1]=res[i];res[0]=0;for(int j(0);j<Md.size();++j)res[Md[j]]=(res[Md[j]]-res[k]*md[Md[j]])%Mod;}}for(int i(0);i<k;++i) ans=(ans+res[i]*b[i])%Mod;return ans+(ans<0?Mod:0);}VI BM(VI s) {VI ret(1,1),B(1,1);int L(0),m(1),b(1);for(int n(0);n<s.size();++n) {int64_t d(0);for(int i(0);i<=L;++i)d=(d+(int64_t)ret[i]*s[n-i])%Mod;if (!d) ++m;else if (2*L<=n) {VI T=ret;int64_t c(Mod-d*powEx(b,Mod-2)%Mod);while (ret.size()<B.size()+m) ret.push_back(0);for (int i(0);i<B.size();++i)ret[i+m]=(ret[i+m]+c*B[i])%Mod;L=n+1-L; B=T; b=d; m=1;} else {int64_t c(Mod-d*powEx(b,Mod-2)%Mod);while (ret.size()<B.size()+m) ret.push_back(0);for(int i(0);i<B.size();++i)ret[i+m]=(ret[i+m]+c*B[i])%Mod;++m;}}return ret;}inline static void extand(VI &a, size_t d, int64_t value = 0) {if (d <= a.size()) return;a.resize(d, value);}inline static int64_t gcdEx(int64_t a, int64_t b, int64_t&x, int64_t& y){if(!b) {x=1;y=0;return a;}int64_t d = gcdEx(b,a%b,y,x);y -= (a/b)*x;return d;}static int64_t CRT(const VI &c, const VI &m) {int n(c.size());int64_t M(1), ans(0);for (int i = 0; i < n; ++i) M *= m[i];for (int i = 0; i < n; ++i) {int64_t x,y,tM(M / m[i]);gcdEx(tM, m[i], x, y);ans = (ans + tM * x * c[i] % M) % M;}return (ans + M) % M;}static VI ReedsSloane(const VI &s, int64_t Mod) {auto Inv = [](int64_t a, int64_t Mod) {int64_t x, y;return gcdEx(a, Mod, x, y)==1?(x%Mod+Mod)%Mod:-1;};auto L = [](const VI &a, const VI &b) {int da = (a.size()>1||(a.size()== 1&&a[0]))?a.size()-1:-1000;int db = (b.size()>1||(b.size()== 1&&b[0]))?b.size()-1:-1000;return max(da, db + 1);};auto prime_power = [&](const VI &s, int64_t Mod, int64_t p, int64_t e) {// linear feedback shift register Mod p^e, p is primevector<VI> a(e), b(e), an(e), bn(e), ao(e), bo(e);VI t(e), u(e), r(e), to(e, 1), uo(e), pw(e + 1);;pw[0] = 1;for (int i(pw[0] = 1); i <= e; ++i) pw[i] = pw[i - 1] * p;for (int64_t i(0); i < e; ++i) {a[i] = {pw[i]}; an[i] = {pw[i]};b[i] = {0}; bn[i] = {s[0] * pw[i] % Mod};t[i] = s[0] * pw[i] % Mod;if (!t[i]) {t[i] = 1; u[i] = e;}else for (u[i] = 0; t[i] % p == 0; t[i] /= p, ++u[i]);}for (size_t k(1);k < s.size(); ++k) {for (int g(0); g < e; ++g) {if (L(an[g], bn[g]) > L(a[g], b[g])) {ao[g] = a[e-1-u[g]];bo[g] = b[e-1-u[g]];to[g] = t[e-1-u[g]];uo[g] = u[e-1-u[g]];r[g] = k - 1;}}a = an; b = bn;for (int o(0); o < e; ++o) {int64_t d(0);for (size_t i(0); i < a[o].size() && i <= k; ++i)d = (d + a[o][i] * s[k - i]) % Mod;if (d == 0) {t[o] = 1;u[o] = e;}else {for (u[o]=0, t[o]=d;!(t[o]%p);t[o]/=p ,++u[o]);int g (e-1-u[o]);if (!L(a[g], b[g])) {extand(bn[o], k + 1);bn[o][k] = (bn[o][k] + d) % Mod;} else {int64_t coef = t[o]*Inv(to[g],Mod)%Mod*pw[u[o]-uo[g]]%Mod;int m(k-r[g]);extand(an[o],ao[g].size()+m);extand(bn[o],bo[g].size()+m);for (size_t i(0);i < ao[g].size(); ++i) {an[o][i+m] -= coef*ao[g][i]%Mod;if (an[o][i + m]<0) an[o][i+m] += Mod;}while (an[o].size() && !an[o].back()) an[o].pop_back();for (size_t i(0); i < bo[g].size(); ++i) {bn[o][i+m] -= coef*bo[g][i]%Mod;if (bn[o][i + m] < 0) bn[o][i + m] -= Mod;}while (bn[o].size()&& !bn[o].back()) bn[o].pop_back();}}}}return make_pair(an[0], bn[0]);};vector<tuple<int64_t, int64_t, int> > fac;for (int64_t i(2); i*i <= Mod; ++i)if (!(Mod % i)) {int64_t cnt(0),pw(1);while (!(Mod % i)) {Mod /= i; ++cnt; pw *= i;}fac.emplace_back(pw, i, cnt);}if (Mod > 1) fac.emplace_back(Mod, Mod, 1);vector<VI> as;size_t n = 0;for (auto &&x: fac) {int64_t Mod, p, e;VI a, b;std::tie(Mod, p, e) = x;auto ss = s;for (auto &&x: ss) x %= Mod;std::tie(a, b) = prime_power(ss, Mod, p, e);as.emplace_back(a);n = max(n, a.size());}VI a(n),c(as.size()),m(as.size());for (size_t i(0); i < n; ++i) {for (size_t j(0); j < as.size(); ++j) {m[j] = std::get<0>(fac[j]);c[j] = i < as[j].size() ? as[j][i] : 0;}a[i] = CRT(c, m);}return a;}int64_t solve(VI a,int64_t n,int64_t Mod,bool prime=true) {VI c;if(prime) c = BM(a);else c = ReedsSloane(a,Mod);c.erase(c.begin());for(int i(0);i<c.size();++i) c[i] = (Mod-c[i])%Mod;return solve(n,c,VI(a.begin(),a.begin()+c.size()));}
}BM;
typedef long long ll;
ll sum;
ll a[1005];
int main()
{cin.tie(0);cout.tie(0);int64_t n, m;cin >> n >> m;a[0]=a[1]=1;vector<int64_t> f;//({0, 1});//f.push_back(0);f.push_back(1);f.push_back(2);sum=2;for (int i = 2; i < 1000; i++){a[i]=(a[i-1]+a[i-2])%Mod;sum =(sum+powEx(a[i],m,Mod))%Mod;f.push_back(sum);                         //这里是因为题意要求出和的值所以导入数列和}printf("%lld\n",BM.solve(f,n-1,Mod,false));   //false-模数为合数return 0;
}

实数域

代码来源

#include<bits/stdc++.h>
using namespace std;
const int MAXN = 2005;
const double eps = 1e-8;
template <typename T> void read(T &x) {x = 0; int f = 1;char c = getchar();for (; !isdigit(c); c = getchar()) if (c == '-') f = -f;for (; isdigit(c); c = getchar()) x = x * 10 + c - '0';x *= f;
}
int cnt, fail[MAXN];
double val[MAXN], delta[MAXN];
vector <double> ans[MAXN];
int main() {int n; read(n);for (int i = 1; i <= n; i++)scanf("%lf", &val[i]);for (int i = 1; i <= n; i++) {double tmp = val[i];for (unsigned j = 0; j < ans[cnt].size(); j++)tmp -= ans[cnt][j] * val[i - j - 1];delta[i] = tmp;if (fabs(tmp) <= eps) continue;fail[cnt] = i;if (cnt == 0) {ans[++cnt].resize(i);continue;}double mul = delta[i] / delta[fail[cnt - 1]];cnt++; ans[cnt].resize(i - fail[cnt - 2] - 1);ans[cnt].push_back(mul);for (unsigned j = 0; j < ans[cnt - 2].size(); j++)ans[cnt].push_back(ans[cnt - 2][j] * -mul);if (ans[cnt].size() < ans[cnt - 1].size()) ans[cnt].resize(ans[cnt - 1].size());for (unsigned j = 0; j < ans[cnt - 1].size(); j++)ans[cnt][j] += ans[cnt - 1][j];}for (unsigned i = 0; i < ans[cnt].size(); i++)cout << ans[cnt][i] << ' ';return 0;
}

BM算法（Berlekamp-Massey算法）：解决线性递推问题相关推荐

（其他算法与技巧）【线性递推 Berlekamp-Massey算法】
( 其他算法与技巧 )[ 线性递推 Berlekamp-Massey算法 ] 原理请看:https://www.cnblogs.com/p-b-p-b/p/10844127.html 和 https ...
最短线性递推式求解与有理函数重建
这一算法来自于我们对"线性递推式拟合"的视角转换,其后得到的算法是自然的. 引理 1. 如果两个有理分式 p1/q1,p2/q2p_1/q_1, p_2/q_2p1/q1,p2 ...
Berlekamp–Massey算法简要介绍
这是一篇翻译向的文章,笔者整理了一些有关Berlekamp–Massey算法的笔记,还增加了一些自己的理解. 下面列出了笔者写此文时所参考的一些资料: wikipedia fjzzq2002 别人的博 ...
[线性代数学习笔记] 线性递推数列及 Berlekamp-Massey 算法的详细推导过程
线性递推数列线性递推对于无限数列 {a0,a1,...}\{a_0,a_1,...\}{a0,a1,...} 和有限非空数列 {r0,r1,...,rm−1}\{r_{0},r_1,...,r ...
L. Poor God Water(ACM-ICPC 2018 焦作赛区网络预赛,ac自动机+矩阵快速幂或 BM线性递推)
描述 God Water likes to eat meat, fish and chocolate very much, but unfortunately, the doctor tells hi ...
杜教板子（BM) 线性递推式（板子）
杜教板子(BM) 线性递推式解决传统线性递推式神器 1 #include <cstdio> 2 #include <cstring> 3 #include <cmath ...
【THUSC 2017】如果奇迹有颜色【polya引理】【矩阵】【计数dp】【BM打表+线性递推】
题意:长度为 nnn 的环染 mmm 种颜色,要求任意相邻 mmm 个元素不能包含全部的颜色.求方案数模 109+710^9+7109+7,循环同构. n≤109,m≤7n\leq 10^9,m\l ...
【算法设计与分析】09 递推方程与算法分析
关于什么是递推方程,这里就不再多说了.本文主要讲讲简单的递推方程来求解算法的时间复杂度文章目录 1. 递推方程的引入 1.1 插入排序时间复杂度求解 1.2 二分归并排序时间复杂度求解 2 总结 1 ...
【模板】BM + CH（线性递推式的求解，常系数齐次线性递推）
这里所有的内容都将有关于一个线性递推: $f_{n} = \sum\limits_{i = 1}^{k} a_{i} * f_{n - i}$,其中$f_{0}, f_{1}, ... , f_{k ...

BM算法（Berlekamp-Massey算法）：解决线性递推问题

算法内容

算法模板

取模，模数为质数

模数可不为质数

实数域

BM算法（Berlekamp-Massey算法）：解决线性递推问题相关推荐

最新文章

热门文章