zoukankan      html  css  js  c++  java
  • bzoj#4161-Shlw loves matrixI【常系数线性齐次递推】

    正题

    题目链接:https://darkbzoj.tk/problem/4161


    题目大意

    给出序列(a),和(h)(0sim k-1)项,满足

    [h_n=sum_{i=1}^na_ih_{n-i} ]

    (h_n)

    (1leq nleq 10^9,1leq kleq 2000)


    解题思路

    显然(kleq 2000)我们就不能矩阵乘法了,然后就去研究了一下常系数线性齐次递推,而且这题(k^2)能过所以不用(O(nlog n))的多项式取模。

    首先对于一个(k imes k)的矩阵(A)的特征多项式(f(x))的定义是

    [f(x)=sum_{i=0}^kdet(A-klambda)x^i ]

    而我们知道对于这种递推式的特征多项式就是

    [f(x)=x^k-sum_{i=0}^{k-1}a_ix^i ]

    然后有个Cayley-Hamilton定理就是说对于(A)的特征多项式(f(x))满足(f(A)=0)

    如果按照这种特殊的特征多项式去理解的话挺好懂的,但是实际上的就不会证了。

    那么我们如果要求一个转移矩阵(A^n=A^nmod f(A)),我们可以求出一个多项式(r(x)=x^nmod f(x)),然后直接把(A)带入这个多项式就好了。

    这个东西要用到快速幂+多项式取模,(O(k^2))的话实现起来很简单,考虑一个(x^{k+i}(igeq 0)mod f(x))会发生什么,因为(f(x)[k]=1)所以(lfloor frac{x^{k+i}}{f(x)} floor=x^i)相当于把(x^i(f(x)-x^k))再丢回去就好了。

    然后我们得到了一个多项式(r(x)),考虑怎么计算答案,只需要计算(A^{n-k}vec{h})的最后一项,所以我们实际上求的(r(x)=x^{n-k}mod f(x)),记(r(x))(i)次系数为(b_i),那么我们有

    [r(A)vec{h}=sum_{i=0}^{k-1}b_iA^{i}vec{h} ]

    所以我们算出(h)(ksim 2k-1)项然后就有

    [h_n=sum_{i=0}^{k-1}r(x)[i]h_{k+i} ]

    时间复杂度:(O(k^2log n))

    值得一提的是如果转移矩阵(A)不是一个递推式的矩阵,我们也可以用拉格朗日插值插出它的特征多项式从而降低平均的复杂度。


    code

    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #define ll long long
    using namespace std;
    const ll N=4100,P=1e9+7;
    ll n,k,a[N],p[N],h[N],b[N],f[N],tmp[N];
    void Mul(ll *a,ll *b,ll *c){
    	memset(tmp,0,sizeof(tmp));
    	for(ll i=0;i<k;i++)
    		for(ll j=0;j<k;j++)
    			(tmp[i+j]+=a[i]*b[j]%P)%=P;
    	for(ll i=2*k-2;i>=k;i--)
    		for(ll j=0;j<k;j++)
    			(tmp[i+j-k]+=P-tmp[i]*p[j]%P)%=P;
    	for(ll i=0;i<k;i++)c[i]=tmp[i];
    	return;
    }
    signed main()
    {
    	scanf("%lld%lld",&n,&k);
    	for(ll i=1;i<=k;i++)
    		scanf("%lld",&a[i]),a[i]=(a[i]+P)%P,p[k-i]=P-a[i];
    	for(ll i=0;i<k;i++)
    		scanf("%lld",&h[i]),h[i]=(h[i]+P)%P;
    	for(ll i=k;i<(k<<1);i++)
    		for(ll j=1;j<=k;j++)
    			(h[i]+=h[i-j]*a[j]%P)%=P;
    	if(n<(k<<1))return printf("%lld
    ",h[n])&0;
    	b[1]=p[k]=f[0]=1;
    	ll B=n-k,ans=0;
    	while(B){
    		if(B&1)Mul(f,b,f);
    		Mul(b,b,b);B>>=1;
    	}
    	for(ll i=0;i<k;i++)
    		(ans+=f[i]*h[i+k]%P)%=P;
    	printf("%lld
    ",ans);
    	return 0;
    }
    
  • 相关阅读:
    页面布局
    Vue学习指南
    《前端JavaScript重点》学习笔记 6-12
    复习3----作用域和闭包
    复习1-变量类型和计算
    复习2--js原型与原型链2
    慕课网《前端JavaScript面试技巧》学习笔记(2)-原型和原型链
    旋转图片
    UITextView添加行距
    YYKit之YYText
  • 原文地址:https://www.cnblogs.com/QuantAsk/p/15346227.html
Copyright © 2011-2022 走看看