杜教BM递推板子

Posted 2022-07-29 forgenvueory

Berlekamp-Massey 算法用于求解常系数线性递推式

#include<bits/stdc++.h>

typedef std::vector<int> VI;
typedef long long ll;
typedef std::pair<int, int> PII;

const ll mod = 1000000007;

ll powmod(ll a, ll b) 
  ll res = 1;
  a %= mod;
  assert(b >= 0);
  for(; b; b >>= 1) 
    if(b & 1)
      res = res * a % mod;
    a = a * a % mod;
  
  return res;


ll _, n;

namespace linear_seq 
const int N = 10010;
ll res[N], base[N], _c[N], _md[N];
std::vector<ll> 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 = ((int)(a).size());
  assert(((int)(a).size()) == ((int)(b).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] != 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;
    
  
  for (int i = 0; i < k; i++)
    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 (((int)(C).size()) < ((int)(B).size()) + 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 (((int)(C).size()) < ((int)(B).size()) + 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;

int 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 solve(n, c, VI(a.begin(), a.begin() + ((int)(c).size())));

;
int main() 
  int t;
  scanf("%d", &t);
  while(t--) 
    scanf("%lld", &n);
    std::vector<int>v = 
      1, 1, 0, 3, 0, 3,
      5, 4, 1, 9, 1, 6,
      9, 7, 2, 15, 2, 9,
      13, 10, 3, 21, 3, 12
    ;
    //至少8项，越多越好。
    printf("%lld\n", linear_seq::gao(v, n - 1) % mod);
  

数据大时都改为 long long
若mod不为质数，则需替换 powmod ,并将BM中的 d*powmod(b, mod-2) 改为 powmod(d, b)

void exgcd(ll a, ll b, ll &x, ll &y) 
  if (!b)
    x = 1, y = 0;
  else
    exgcd(b, a % b, y, x), y -= x * (a / b);

ll inv(ll a, ll p) 
  ll x, y;
  exgcd(a, p, x, y);
  return(x + p) % p;

VI prime, g;
void getPrime() 
  ll qw = mod;
  for(ll i = 2; i * i <= qw; i++) 
    if(qw % i == 0)
      prime.pb(i);
    while(qw % i == 0)
      qw /= i;
  
  if(qw > 1)
    prime.push_back(qw);

ll powmod(ll fz, ll fm) 
  ll ret = 1ll;
  ll cnt[5] = 0;
  for(int k = 0; k < prime.size(); k++) 
    if(fz % prime[k] == 0) 
      while(fz % prime[k] == 0) 
        fz /= prime[k];
        cnt[k]++;
      
    
  
  ret = fz % mod;
  for(int k = 0; k < prime.size(); k++) 
    if(fm % prime[k] == 0) 
      while(fm % prime[k] == 0) 
        fm /= prime[k];
        cnt[k]--;
      
    
  
  if(fm > 1)
    ret = (ret * inv(fm, mod)) % mod;
  for(int k = 0; k < prime.size(); k++)
    for(int kk = 1; kk <= cnt[k]; kk++)
      ret = (ret * prime[k]) % mod;
  return ret;