题目链接

我们发现
\[f_i=2f_{i-1}+(i-1)* (i-2) \]
我们可以设
\[g_i=i * (i-1)\]
则
\[g_i=i*(i-1)=(i-1) * (i-2)+2(i-1)=g_{i-1}+2(i-1)\]
则
我们可以推出初始矩阵
\[\begin{bmatrix}& f_0&g_0&0&1 & \end{bmatrix}=\begin{bmatrix}& 0&0&0&1 & \end{bmatrix}\]
第一个数是\(f_0\)，第二个数是\(g_0\)，第三个数是i，第四个数是1。
以及转移矩阵
\[\begin{bmatrix} & 2&0&0&0 & \\\\ & 1&1&0&0 & \\\\ & 0&2&1&0 & \\\\ & 0&0&1&1 & \end{bmatrix}\]
第一行对应的是\(f_0\)，第二行对应的是\(g_0\)，第三行对应的是i，第四行对应的是1。

对于题目给出的\(10^{5000}\)，本蒟蒻因为太弱，不想打高精，想到了使用类似于快速幂的方法，因为$a^n=(a^{x})^{10}* a^y (10* x+y=n) $

可以直接避免高精

于是，我们就可以写出代码

#include<cstdio>
#include<cctype>
#include<cstring>
#define re register
#define ll long long
using namespace std;
const int N=1e7+10,mod=1e9+7;
struct matrix
{
    int n,m,g[5][5];
    inline void init(int _n,int _m)
    {
        memset(g,0,sizeof(g));
        n=_n;
        m=_m;
    }
    inline matrix operator *(const matrix &b)
    {
        matrix res;
        res.init(n,b.m);
        for(re int i=1; i<=res.n; i++)
            for(re int j=1; j<=res.m; j++)
                for(re int k=1; k<=b.n; k++)
                    res.g[i][j]=(res.g[i][j]+1ll*g[i][k]*b.g[k][j])%mod;
        return res;
    }
    template<typename T>
    inline matrix operator ^ (T b)
    {
        matrix res,a=(*this);
        res.init(4,4);
        for(re int i=1; i<=4; i++)
            res.g[i][i]=1;
        while(b)
        {
            if(b&1)
                res=res*a;
            a=a*a;
            b>>=1;
        }
        return res;
    }
    inline void print()
    {
        for(re int i=1; i<=n; i++)
        {
            for(re int j=1; j<=m; j++)
                printf("%d ",g[i][j]);
            putchar('\n');
        }
        return;
    }
} a,b,f;
int n;
char s[5010];
int main()
{
    a.init(1,4);
    a.g[1][4]=1;
    b.init(4,4);
    b.g[1][1]=2;
    b.g[2][1]=b.g[2][2]=1;
    b.g[3][2]=2,b.g[3][3]=1;
    b.g[4][3]=b.g[4][4]=1;
    f.init(4,4);
    for(re int i=1; i<=4; i++)
        f.g[i][i]=1;
    scanf("%s",s+1);
    for(re int i=1; s[i]; i++)
        f=(f^10)*(b^((int)s[i]-'0'));
    printf("%d\n",(a*f).g[1][1]);
    return 0;
}