洛谷-5175

洛谷-5175

思路

这里面的第一篇题解

根据需要啥就加啥的想法！(妙啊)

我的矩阵转换

\[\begin{pmatrix} S_i & x^2a_i^2+y^2a_{i-1}^2+2xya_ia_{i-1} & a_i^2 & xa_i^2+ya_ia_{i-1} \end{pmatrix} \\ = \begin{pmatrix} S_{i-1} & a_i^2 & a_{i-1}^2 & a_ia_{i-1} \end{pmatrix} * S \]

也即

\[ \begin{pmatrix} S_i & x^2a_i^2+y^2a_{i-1}^2+2xya_ia_{i-1} & a_i^2 & xa_i^2+ya_ia_{i-1} \end{pmatrix} \\ = \begin{pmatrix} S_1 & a_2^2 & a_1^2 & a_2a_1 \end{pmatrix} * S^{i-1} \]

其中

\[S = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 1 & x^2 & 1 & x \\ 0 & y^2 & 0 & 0 \\ 0 & 2xy & 0 & y \end{pmatrix} \]

Code

#include <bits/stdc++.h>
using namespace std;
#define _u_u_ ios::sync_with_stdio(false), cin.tie(nullptr)
#define cf int _o_o_;cin>>_o_o_;for (int Case = 1; Case <= _o_o_;Case++)
#define SZ(x) (int)(x.size())
inline void _A_A_();
signed main() {_A_A_();return 0;}

using ll = long long;
#define int long long
const int mod = 1e9 + 7;    // 如果mod不是const，会超时
const int maxn = 2e5 + 10;
const int N = 4, M = 5010;
const int inf = 0x3f3f3f3f;

struct mat
{
    int m[N][N];
};

mat operator * (const mat&a, const mat&b) {
    mat c;
    memset(c.m,0, sizeof c.m);
    for (int i = 0;i < N;i++) {
        for (int j = 0;j < N;j++) {
            for (int k = 0;k < N;k++) {
                c.m[i][j] = (c.m[i][j] + a.m[i][k] * b.m[k][j] % mod) % mod;
            }
        }
    }
    return c;
}

mat qpow(mat a, int n) {
    mat c;
    memset(c.m, 0, sizeof c.m);
    for (int i = 0;i < N;i++) c.m[i][i] = 1;
    while(n) {
        if (n & 1) {
            c = c * a;
        }
        a = a * a;
        n >>= 1;
    }
    return c;
}

ll n, a1,a2, x, y;

inline void _A_A_() {
    #ifdef LOCAL
    freopen("in.in", "r", stdin);
    #endif
    _u_u_;
    mat s;
    cf {
        cin >> n >> a1 >> a2 >> x >> y;
        memset(s.m, 0, sizeof s.m);
        s.m[0][0] = s.m[1][0] = s.m[1][2] = 1;
        s.m[1][1] = x * x % mod, s.m[1][3] = x % mod, s.m[2][1] = y * y % mod, s.m[3][1] = 2 * x * y % mod, s.m[3][3] = y % mod;
        s = qpow(s, n - 1);
        int ans = ((((a1 * a1) % mod) * s.m[0][0] % mod) + ((a2 * a2 % mod) * s.m[1][0]) % mod + ((a1 * a1 % mod) * s.m[2][0]) % mod + ((a1 * a2%mod) * s.m[3][0] ) %mod)%mod;
        cout << ans << "\n";
    }
}