G - Rock (Hard)

結論から言うと、答えは $n = 1$ のとき $6$、$n ≧ 2$のときは以下の式になります。
これは行列累乗を用いて $O(logN)$ で求まるため、十分高速です。

$(答え) = 60×10^{n-2} + v_0×2 + v_1×1 + v_2×2$

この式の気持ちについて説明しましょう。
まず$k$ 桁( $k ≧ 2$ )の条件を満たす数のうち、桁の寄与を一の位とそれ以外に分けました。
すると、一の位は必ず偶数でなければならないので、$0, 2, 4, 6, 8$ のいずれかになります。
$0$ と $6$ を選ぶためには十の位以上の和が $3$ の倍数でなければなりません。同様に、$2$ と $6$ を選ぶためには $mod 3$ で$1$ 、$4$を選ぶためには $mod 3$ で $2$ である必要がありますが、十の位以上の桁の選び方のうち、$mod 3$ で $0, 1, 2$となる選び方はいずれも $3×10^{n - 2}$個になります。 よって、$(0 + 2 + 4 + 6 + 8) × 3 × 10^{n - 2} = 60 × 10^{n - 2}$となります。

次に、十の位以上の寄与を考えていきましょう。"$k$ 桁であり、十の位以上を見たときに $mod 3$ が $m(m = 1, 2, 3)$ となる数の十の位以上の部分"をすべて列挙した時の個数を $n[k][i]$ 、それらのすべての桁の和を $dp[k][i]$ とおくと、以下の漸化式が成り立ちます。 $dp[k + 1][1] = 4 * dp[k][1] + 3 * dp[k][2] + 3 * dp[k][3] + (0 + 3 + 6 + 9) * (n[k][1]) + (2 + 5 + 8) * (n[k][2]) + (1 + 4 + 7) * (n[k][3])$
$dp[k + 1][2] = 3 * dp[k][1] + 4 * dp[k][2] + 3 * dp[k][3] + (1 + 4 + 7) * (n[k][1]) + (0 + 3 + 6 + 9) * (n[k][2]) + (2 + 5 + 8) * (n[k][3])$
$dp[k + 1][3] = 3 * dp[k][1] + 3 * dp[k][2] + 4 * dp[k][3] + (2 + 5 + 8) * (n[k][1]) + (1 + 4 + 7) * (n[k][2]) + (0 + 3 + 6 + 9) * (n[k][3])$

これをそのまま行列で表現しても $6×6$ の行列累乗に帰着できるため十分高速ですが、$n[k][i]$ は $i$ によらず一定なので、その性質を用いると上に示した $4×4$ の行列に帰着することができます。

def matrixMultiplication(a,b,m): #行列の掛け算(a×b) m:mod
    I,J,K,L = len(a),len(b[0]),len(b),len(a[0])
    if(L!=K):
        return -1
    c = [[0] * J for _ in range(I)]
    for i in range(I) :
        for j in range(J) :
            for k in range(K) :
                c[i][j] += a[i][k] * b[k][j]
            c[i][j] %= m
    return c

def matrixExponentiation(x,n,m): #行列の累乗 (x^n) m:mod
    y = [[0] * len(x) for _ in range(len(x))]
    for i in range(len(x)):
        y[i][i] = 1
    while n > 0:
        if n & 1:
            y = matrixMultiplication(x,y,m)
        x = matrixMultiplication(x,x,m)
        n >>= 1
    return y

mod = 10**9 + 7
t = int(input())

P = [
    [4, 3, 3, 45],
    [3, 4, 3, 45],
    [3, 3, 4, 45],
    [0, 0, 0, 10],
    ]

A = [[12], [15], [18], [3]]

for _ in range(t):
    n = int(input())
    if(n == 1):
        print(6)
        continue

    Q = matrixExponentiation(P, n-2, mod)
    V = matrixMultiplication(Q, A, mod)
    ans = 60 * pow(10, n-2, mod) + 2* V[0][0] + V[1][0] + 2 * V[2][0]
    print(ans % mod)

#include<bits/stdc++.h>
using namespace std;
using ll = long long;

template<class T>
struct Matrix {
    int _R , _C;
    vector<vector<T>> val;

    Matrix(int r, int c, T Val): val(r,vector<T>(c,Val)), _R(r), _C(c) {}
    Matrix(int r, int c): Matrix(r, c, 0) {}

    vector<T>& operator[](int i) {
        return val[i];
    }

    Matrix operator*(const Matrix& other) {
        assert(this->_C == other._R);
        Matrix result(this->_R, other._C);

        for (int i = 0; i < this->_R; i++) {
            for (int j = 0; j < other._C; j++) {
                for (int k = 0; k < this->_C; k++) {
                    result[i][j] += this->val[i][k] * other.val[k][j];
                }
            }
        }
        return result;
    }

};

struct mint {
	static const int mod = 1'000'000'007;
	ll val;
	mint(int k = 0): val((k % mod + mod) % mod) {}
	mint operator*(const mint other) {return val * other.val % mod;}
	mint& operator+=(const mint other) {
        val = val + other.val >= mod? val + other.val - mod : val + other.val;
        return *this;
    }
};

Matrix<mint> pow_mat(Matrix<mint>a,ll n) {
    int si = a[0].size();
    Matrix<mint> res(si,si);
    for (int i = 0; i < si; i++) res[i][i] = 1;

    while (n > 0) {
        if (n & 1) res = res * a;
        a = a * a;
        n >>= 1;
    }
    return res;
}

void solve(ll n) {

    Matrix<mint> mat(4, 4);
    mat[0] = {4, 3, 3, 45};
    mat[1] = {3, 4, 3, 45};
    mat[2] = {3, 3, 4, 45};
    mat[3] = {0, 0, 0, 10};

    Matrix<mint> dp(4, 1);
    dp[0] = {18};
    dp[1] = {12};
    dp[2] = {15};
    dp[3] = {3 };

    dp = pow_mat(mat, n - 2) * dp;

    mint ans = 0;
    ans += dp[0][0] * 2;
    ans += dp[1][0] * 2;
    ans += dp[2][0] * 1;
    ans += dp[3][0] * 20;

    cout << ans.val << endl;
}

int main() {

    int q;
    cin >> q;
    for (int i = 0; i < q; i++) {
        ll n;
        cin >> n;

        if (n == 1)cout << 6 << endl;
        else solve(n);
    }
}