思路（这个解法重点在递推矩阵的推导上，可能会有些枯燥）

思路很明显，和矩阵快速幂相关，需要找到递推矩阵.

首先知道朴素求斐波那契矩阵第n项的求法，那么就用相同的方法来.

先把我们需要的$T_n$写到矩阵中，有$$\left[\begin{matrix}f_n&f_{n+1}&T_n\end{matrix}\right]$$

根据这个形式，对应的下一项就应该是$$\left[\begin{matrix}f_{n+1}&f_n+f_{n+1}&T_n+(n+1)f_{n+1}\end{matrix}\right]$$

前面两项不难，第三项缺一个$(n+1)f_{n+1}$，那么就直接在矩阵中添加这一项，变成$$\left[\begin{matrix}f_n&f_{n+1}&T_n&(n+1)f_{n+1}\end{matrix}\right]$$
因为添加了一项，所以下一项的矩阵也要对应增加一项$$\left[\begin{matrix}f_{n+1}&f_n+f_{n+1}&T_n+(n+1)f_{n+1}&(n+2)f_{n+2}\end{matrix}\right]$$
这时候递推的前三项可以比较容易直接构造矩阵了，但最后一项需要再操作一下。$(n+2)f_{n+2}$含有$n+2$的项，把它拆成$n$和$n+1$的表达形式，其实就是看着前一个矩阵中的元素，尽量往前面的矩阵靠拢
$$(n+2)f_{n+2}=(n+1)f_n+(n+1)f_{n+1}+f_n+f_{n+1}$$

可以看到只有$(n+1)f_n$在前一个矩阵是没有的，再进行一次添加，在原矩阵中添加$(n+1)f_n$，变成$$\left[\begin{matrix}f_n&f_{n+1}&T_n&(n+1)f_{n+1}&(n+1)f_n\end{matrix}\right]$$

再次求出新的下一个矩阵为$$\left[\begin{matrix}f_{n+1}&f_n+f_{n+1}&T_n+(n+1)f_{n+1}&(n+2)f_{n+2}&(n+2)f_{n+1}\end{matrix}\right]$$

这时候新添加的$(n+2)f_{n+1}=(n+1)f_{n+1}+f_{n+1}$则可以直接在原矩阵中找到！（递归推导结束！）

缓一口气

推导的结果就是，要找一个递推矩阵，使得矩阵$$\left[\begin{matrix}f_n&f_{n+1}&T_n&(n+1)f_{n+1}&(n+1)f_n\end{matrix}\right]$$能一步变成$$\left[\begin{matrix}f_{n+1}&f_n+f_{n+1}&T_n+(n+1)f_{n+1}&(n+2)f_{n+2}&(n+2)f_{n+1}\end{matrix}\right]$$

递推矩阵如下
$$ \left[\begin{matrix}f_n&f_{n+1}&T_n&(n+1)f_{n+1}&(n+1)f_n\end{matrix}\right]* \left[\begin{matrix} 0 & 1 & 0 & 1 & 0\\\\ 1 & 1 & 0 & 1 & 1\\\\ 0 & 0 & 1 & 0 & 0\\\\ 0 & 0 & 1 & 1 & 1\\\\ 0 & 0 & 0 & 1 & 0 \end{matrix}\right]= \left[\begin{matrix}f_{n+1}&f_n+f_{n+1}&T_n+(n+1)f_{n+1}&(n+2)f_{n+2}&(n+2)f_{n+1}\end{matrix}\right] $$

相同思路的另一个递推矩阵以及代码

思路主要看上面的推到就可以了
（这个构造方式已经发现可以构造三个不同的递推矩阵进行解答）
代码写的递推矩阵是7*7的，写题解的时候发现可以减少为5*5，由于本人太懒，这里只有7*7的AC代码
从矩阵
$$ \left[\begin{matrix} f_n & f_{n+1} & T_{n-1} & n+1 & 1 & nf_n & nf_{n-1} \end{matrix}\right] $$
到
$$ \left[\begin{matrix} f_{n+1} & f_n+ f_{n+1} & T_{n-1}+nf_n & n+2 & 1 & (n+1)f_{n+1} & (n+1)f_n \end{matrix}\right] $$
递推矩阵为
$$ \left[\begin{matrix} 0 & 1 & 0 & 0 & 0 & 0 & 1\\\\ 1 & 1 & 0 & 0 & 0 & 1 & 0\\\\ 0 & 0 & 1 & 0 & 0 & 0 & 0\\\\ 0 & 0 & 0 & 1 & 0 & 0 & 0\\\\ 0 & 0 & 0 & 1 & 1 & 0 & 0\\\\ 0 & 0 & 1 & 0 & 0 & 1 & 1\\\\ 0 & 0 & 0 & 0 & 0 & 1 & 0 \end{matrix}\right] $$

#include<iostream>
#include<cstring>
#include<algorithm>
#define int long long
using namespace std;

int mod;

typedef struct matrix{
    int a[7][7];

    matrix(){
        memset(a, 0, sizeof a);
    }

}matrix;

matrix multi(const matrix& m1, const matrix& m2){
    matrix ans;
    for(int i = 0; i < 7; ++i)
        for(int j = 0; j < 7; ++j)
            for(int k = 0; k < 7; ++k)
                ans.a[i][j] = (ans.a[i][j] + m1.a[i][k] * m2.a[k][j] % mod) % mod;
    return ans;
}

matrix qpow(matrix a, int b){
    if(b == 0){
        matrix ans;
        for(int i = 0; i < 7; ++i) ans.a[i][i] = 1;
        return ans;
    }
    matrix ans = a;
    b--;
    while(b){
        if(b & 1LL) ans = multi(ans, a);
        a = multi(a, a);
        b >>= 1LL;
    }
    return ans;
}

int getf(int n){
    matrix f, fk;
    f.a[0][0] = 1;
    f.a[0][1] = 2;
    f.a[0][2] = 1;
    f.a[0][3] = 3;
    f.a[0][4] = 1;
    f.a[0][5] = 2;
    f.a[0][6] = 2;

    fk.a[0][1] = fk.a[0][6] = 1;
    fk.a[1][0] = fk.a[1][1] = fk.a[1][5] = 1;
    fk.a[2][2] = 1;
    fk.a[3][3] = 1;
    fk.a[4][3] = fk.a[4][4] = 1;
    fk.a[5][2] = fk.a[5][5] = fk.a[5][6] = 1;
    fk.a[6][5] = 1;
    n--;
    fk = qpow(fk, n);
    f = multi(f, fk);
    return f.a[0][2];
}

signed main(){
    int n;
    cin >> n >> mod;
    int ans = getf(n);
    cout << ans;
    return 0;
}