B-M (特征方程教程)模板

1 篇文章 0 订阅

 出处 https://jhcloud.top/blog/?p=1271

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define rep(i,a,n)for(int i=a;i<n;i++)
#define per(i,a,n)for(int i=n-1;i>=a;i--)
#define pb push_back
#define mp make_pair
#define fi first
#define se second
#define SZ(x) ((ll)(x).size())
#define all(x) (x).begin(),(x).end()
typedef vector<ll> VI;
typedef pair<ll,ll>PII;
const ll mod = 1e9+7;

ll powmod(ll A,ll b)
{
    ll res = 1;
    A%=mod;
    for(; b; b>>=1)
    {
        if(b&1)
            res  =res*A%mod;
        A = A*A%mod;
    }
    return res;
}

namespace liner_seq
{
    const int N = 100010;
    ll res[N],base[N],_c[N],_md[N];
    vector<ll>Md;
    void mul(ll *a,ll *b,ll k)
    {
        rep(i,0,k+k)_c[i] = 0;
        rep(i,0,k)if(a[i])
            rep(j,0,k)_c[i+j] = (_c[i+j]+a[i]*b[j])%mod;
        for(ll i=k+k-1; i>=k; i--)
        {
            if(_c[i])
                rep(j,0,SZ(Md))_c[i-k+Md[j]] = (_c[i-k+Md[j]]-_c[i]*_md[Md[j]])%mod;
        }
        rep(i,0,k)a[i]=_c[i];
    }

    ll solve(ll n,VI a,VI b)
    {
        ll ans =0,pnt =0;
        ll k = SZ(a);
        rep(i,0,k) _md[k-1-i] = -a[i];
        _md[k] = 1;
        Md.clear();
        rep(i,0,k)if(_md[i]!=0)
            Md.push_back(i);
        rep(i,0,k)res[i ] =base[i] =0;
        res[0] = 1;
        while((1ll<<pnt)<=n)
            pnt++;
        for(ll p=pnt; p>=0; p--)
        {
            mul(res,res,k);
            if((n>>p)&1)
            {
                for(ll i=k-1; i>=0; i--)
                    res[i+1] = res[i];
                res[0 ] =0;
                rep(j,0,SZ(Md))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);
        ll L  =0,m=1,b = 1;
        rep(n,0,SZ(s))
        {
            ll d =0;
            rep(i,0,L+1)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.pb(0);
                rep(i,0,SZ(B))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.pb(0);
                rep(i,0,SZ(B))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());
        rep(i,0,SZ(c))c[i] = (mod-c[i])%mod;
        return solve(n,c,VI(a.begin(),a.begin()+SZ(c)));
    }
}

VI a;

int main()
{
    for(int i=0; i<18; i++)
    {
        a.push_back(i);
    }
    int T;
    cin>>T;
    while(T--)
    {
        long long n;
        cin>>n;
        printf("%lld\n",liner_seq::gao(a,n)%mod);
    }
    return 0;
}

 

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
可以处理无限数据,求其线性复杂度 //by史瑞 //LFSR线性移位寄存器以及求异或运算OK unsigned char LFSRB_M(unsigned char *fun,unsigned char *seq,unsigned long Cont){ unsigned long x; unsigned char ch=0x00,t=0x00,*array; array=(unsigned char *)malloc(Cont*sizeof(unsigned char)); for(x=0;x<Cont;x++) *(array+x)=(*(seq+x))&(*(fun+x)); for(x=0;x<Cont*8;x++){ ch=(*array)&0x80; LeftShift(array,Cont); t^=ch; } free(array); return t; } #define word(ln) ((ln-1)/8+1) #define place(ln) ((ln-1)%8) //Berlek_Massey求生成任意序列的联接多项式OK unsigned char *Berlek_Massey(unsigned char *seq,unsigned long *Rank,unsigned long Cont,FILE *fmm){ unsigned long n=0,m,lm,ln=0; int d; unsigned long lfm; unsigned char *fun,*fm,*array,ch,t=0x00; unsigned char *func,*fmc; unsigned long x,y,k; fun=(unsigned char *)malloc(sizeof(unsigned char)); *(fun)=0x00; for(x=0;x<Cont*8;x++){ if(ln!=0){ array=(unsigned char *)malloc(word(ln)); InitialDSR(array,word(ln)); for(y=0;y<ln;y++){ ch=((*(seq+(n-(y+1))/8))<<((n-(y+1))%8))&0x80; if(ch) *(array+y/8)^=(ch>>(y%8)); } t=LFSRB_M(fun,array,word(ln)); d=((((*(seq+n/8))<<(n%8))&0x80;)^t)?1:0; free(array); } else d=(((*(seq+n/8))<<(n%8))&0x80;)?1:0; if(d){ if(ln!=0){ lm=ln; func=(unsigned char *)malloc(word(ln)*sizeof(unsigned char)); memcpy(func,fun,word(ln)); if(ln<(n+1-ln)){ ln=n+1-ln; } fmc=(unsigned char *)malloc(word(ln)*sizeof(unsigned char)); InitialDSR(fmc,word(ln)); memcpy(fmc,fm,word(lfm)); for(k=0;k<n-m;k++){ RightShift(fmc,word(ln)); } ch=0x80; *(fmc+(n-m-1)/8)^=(ch>>((n-m-1)%8)); fun=(unsigned char *)realloc(fun,word(ln)*sizeof(unsigned char)); for(k=word(lm);k<word(ln);k++) *(fun+k)=0x00; XorDSR(fun,fmc,word(ln)); free(fmc); if(lm<(n+1-lm)){ m=n; lfm=lm; fm=(unsigned char *)realloc(fm,word(lm)); memcpy(fm,func,word(lfm)); } free(func); } else{ ln=n+1; m=n; if(m!=0){ fm=(unsigned char *)malloc(word(m)*sizeof(unsigned char)); InitialDSR(fm,word(m)); lfm=m; } else{ fm=(unsigned char *)malloc(sizeof(unsigned char)); *fm=0x00; lfm=1; } fun=(unsigned char *)realloc(fun,word(ln)); InitialDSR(fun,word(ln)); ch=0x80; ch>>=(place(ln)); *(fun+(word(ln)-1))^=ch; } } n++; printf("\t<%d,%d>",n,ln); fprintf(fmm,"\t<n,d,ln>=<%d,%d,%d>",n,d,ln); if(n%3==0){ printf("\n"); fprintf(fmm,"\n"); } } printf("\nFn="); // for(k=0;k<word(ln);k++) // FromBytetoBit(*(fun+k)); *Rank=ln; return fun; }

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值