洛谷P9510 「STA - R3」高维立方体题解

题意：

如下定义斐波那契数列：
$\operatorname{fib}(n)=\begin{cases}1&n\le 2\\\operatorname{fib}(n-1)+\operatorname{fib}(n-2)&n>2\end{cases}$
现在我们定义一个函数（注意在 $n<1$ 时这个函数的值是 $0$ ）：
$f(n)=\sum_{i=1}^n\operatorname{fib}^2(i)$
由于求斐波那契数列的前缀和太简单了，你需要求出：
$\sum_{i=1}^n\operatorname{fib}(i)\cdot(f(i-2)+\operatorname{fib}^2(i)+\operatorname{fib}(i))$
的值，答案对输入的 $p$ 取模。

注：$\operatorname{fib}^2(x)$ 表示 $\operatorname{fib}(x)$ 的平方。

输入格式：

本题有多组数据。

第一行一个整数 $T$，表示数据的组数。

对于每组数据，一行两个整数 $n,p$。

输出格式：

对于每组数据，输出一行一个整数，表示答案对 $p$ 取模后的结果。

数据范围：

$1\le T\le 2\times 10^5,~1\le n\le 10^{18},~2\le p\le 10^9+7$。

考虑把 $\operatorname{fib}^2(x)$ 看作边长为 $\operatorname{fib}(x)$ 的正方形，如图所示

可以发现这个长方形的长是 $\operatorname{fib}(n + 1)$ ，宽是 $\operatorname{fib}(n)$ ，则

$\sum_{i=1}^n \operatorname{fib}^2(i)=\operatorname{fib}(n) \operatorname{fib}(n+1)$

同理，$\operatorname{fib}^3(x)$ 可以看作棱长为 $\operatorname{fib}(x)$ 的立方体，如图所示

其中红色部分的体积为 $\operatorname{fib}(n) \operatorname{fib}(n-1) \operatorname{fib}(n-2)$ ，那么

$\begin{aligned} \mathrm{Ans} &=\sum_{i=1}^n \mathrm{fib}(i) \cdot\left(f(i-2)+\operatorname{fib}^2(i)+\operatorname{fib}(i)\right) \\[12pt]& =\sum_{i=1}^n \left(\operatorname{fib}(i) f(i-2)+\operatorname{fib}^3(i)\right)+\sum_{i=1}^n \mathrm{fib}^2(i) \\[12pt]& =\operatorname{fib}^2(n) \operatorname{fib}(n+1)+\operatorname{fib}(n) \operatorname{fib}(n+1) \end{aligned}$

用矩阵快速幂求出 $\operatorname{fib}(n)$ 和 $\operatorname{fib}(n+1)$ 即可。

时间复杂度 $\mathcal{O}(8T \log n)$ ，本题卡常。

代码：

#include <bits/stdc++.h>
using namespace std;
// #define int long long
// #define INF 0x3f3f3f3f3f3f3f3f
typedef long long ll;
// typedef vector< vector<ll> > mat;
void up(int &x, int y) { x < y ? x = y : 0; }
void down(int &x, int y) { x > y ? x = y : 0; }
#define rep(i, a, b) for(int i = (a), i##END = (b); i <= i##END; i++)
#define Rep(i, a, b) for(int i = (a), i##END = (b); i >= i##END; i--)
#define N ((int)())

int mod; struct mat { ll a[2][2]; };
mat mul(const mat &A, const mat &B)
{
    static mat C; C.a[0][0] = C.a[0][1] = C.a[1][0] = C.a[1][1] = 0;
    rep(i, 0, 1) rep(j, 0, 1) {
        ll t = 0; rep(k, 0, 1) t += A.a[i][k] * B.a[k][j];
        C.a[i][j] = (C.a[i][j] + t) % mod;
    }
    return C;
}
mat qpow(mat &A, ll k)
{
    static mat B; B.a[0][0] = B.a[1][1] = 1; B.a[0][1] = B.a[1][0] = 0;
    for(; k; A = mul(A, A), k >>= 1) if(k & 1) B = mul(B, A);
    return B;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    static mat A,B; B.a[0][0] = B.a[1][0] = 1; int qwq; cin >> qwq;
    while(qwq--)
    {
        static ll n, a, b; cin >> n >> mod;
        A.a[0][0] = A.a[0][1] = A.a[1][0] = 1; A.a[1][1] = 0;
        A = qpow(A, n - 1); A = mul(A, B); a = A.a[1][0], b = A.a[0][0];
        cout << (a * a % mod + a) % mod * b % mod << '\n';
    }
    return 0;
}