数列 题解【矩阵加速递推】【数学】

作者: wjyyy 分类: 数列,数学,矩阵乘法,矩阵加速递推,解题报告,递推 发布时间: 2018-10-28 15:22

点击量:36

 

    带关于\(n\)的高次项也可以矩阵加速。

 

题目描述

这是一道很简单的题。

给你两个数\(a_1,a_2\),求\(a_n\)在模\(2147439647\)下的值。

满足递推关系:\(a_n=a_{n-1}+2a_{n-2}+n^4\)

保证\(1\le a_1,a_2\le 2^{31}\)

输入描述

三个数字,\(n,a_1,a_2\)。

输出描述

一个数\(a_n\)。

输入输出样例

Sample Input

3 1 2

Sample Output

85

数据范围

对于\(30\%\)的数据,\(n\le 1000\),
对于\(50\%\) 的数据,\(n\le 100000\),
对于\(100\%\)的数据,\(n\le 10^{15}\)。

题解:

    看到递推式,第一个想到的就是矩阵加速了,但是发现每次有一个会改变的高次项,感觉非常不好处理。

    不过\(n^4\)既然不是常数项,而且变化有规律,那么就也是可以递推的。可以发现\(n^4\)和\((n-1)^4\)的关系,但是这样要配方,不如研究\(n^4\)和\((n+1)^4\)的关系。

    那么根据二项式定理化简式子,\((n+1)^4=n^4+4n^3+6n^2+4n+1\),则每次转移需要调用\(n^3,n^2,\)以及\(n\)。同时推出\(n^3\)的关系为\((n+1)^3=n^3+3n^2+3n+1\),\(n^2\)为\((n+1)^2=n^2+2n+1\),而\(n\)又是可以直接从\(n-1\)推过来的,所以\(n^4\)就是能递推的了。

    我们把这些状态放在矩阵里,就是\(a:\begin{bmatrix} a_n\\ a_{n-1}\\ n^4\\ n^3\\ n^2\\ n\\ 1 \end{bmatrix}\),然后递推矩阵为\(b:\begin{bmatrix}1 & 2 & 1 & 4 & 6 & 4 & 1 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 4 & 6 & 4 & 1 \\ 0 & 0 & 0 & 1 & 3 & 3 & 1 \\ 0 & 0 & 0 & 0 & 1 & 2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1\end{bmatrix}\)。第\(i\)行表示那一行\(a_i\)的状态由上一个状态转移过来的式子(式子中的\(n\)指上一个状态的\(n\)):\(b_{(i,1)}a_n+b_{(i,2)}a_{n-1}+b_{(i,3)}n^4+b_{(i,4)}n^3+b_{(i,5)}n^2+b_{(i,6)}n+b_{(i,7)}\)。

    最后递推得到的矩阵就是\(c=b^{n-2}\times a\),答案在\(c_{(1,1)}\)处取,时间复杂度为\(O(7^3\log n)\)。模数为231+9999什么鬼啊。

Code:

#include<cstdio>
#include<cstring>
long long p=2147493647;
struct matrix
{
    long long a[7][7];
    int x,y;
    matrix(int x,int y)
    {
        this->x=x;
        this->y=y;
        memset(a,0,sizeof(a));
    }
    matrix(){}
    friend matrix operator *(matrix u,matrix v)//矩乘
    {
        matrix ans(u.x,v.y);
        for(int i=0;i<u.x;++i)
            for(int j=0;j<u.y;++j)
                for(int k=0;k<u.y;++k)
                    ans.a[i][j]=(ans.a[i][j]+u.a[i][k]*v.a[k][j])%p;
        return ans;
    }
    friend matrix operator ^(matrix u,long long v)//快速幂
    {
        matrix ans(u.x,u.y),m=u;
        for(int i=0;i<u.x;++i)
           ans.a[i][i]=1;
        for(;v;v>>=1)
        {
            if(v&1)
                ans=ans*m;
            m=m*m;
        }
        return ans;
    }
};
int st[7][7]={//转移矩阵
    {1,2,1,4,6,4,1},
    {1},
    {0,0,1,4,6,4,1},
    {0,0,0,1,3,3,1},
    {0,0,0,0,1,2,1},
    {0,0,0,0,0,1,1},
    {0,0,0,0,0,0,1}};
int main()
{
    long long n,f1,f2;
    scanf("%lld%lld%lld",&n,&f1,&f2);
    matrix t(7,7),ans(7,1);
    ans.a[0][0]=f2;//初始化状态矩阵
    ans.a[1][0]=f1;
    ans.a[2][0]=16;
    ans.a[3][0]=8;
    ans.a[4][0]=4;
    ans.a[5][0]=2;
    ans.a[6][0]=1;
    for(int i=0;i<7;++i)
        for(int j=0;j<7;++j)
            t.a[i][j]=st[i][j];
    t=t^(n-2);
    ans=t*ans;
    printf("%lld\n",ans.a[0][0]);
    return 0;
}

 

说点什么

avatar
  Subscribe  
提醒
/* */