Yukicoder No.426 往復漸化式

,

$b_{i-1}$を$b_{i+1}$だと思い込み$b_n$には気付いたが$b_0$のtypoだと無視した結果、誤読で時間を潰し、editorialを見て気付いた。 ついでに $ \left( \begin{matrix} 6i & 6i+1 & 6i+2 \
6i+3 & 6i+4 & 6i+5 \end{matrix} \right) $ の$i$は虚数だと思ってた。

solution

segment木。$O(N + Q \log N)$。

漸化式を整理すると以下のようになる。

$$ \begin{array}{ccl} a_i & = & A_{i-1} A_{i-2} \dots A_0 a_0 \
b_i & = & B_{i+1} A_{i+2} \dots B_n b_n \
& & + \; S_{i+1} a_{i+1} \
& & + \; B_{i+1} S_{i+2} A_{i+1} a_{i+1} \
& & + \; B_{i+1} B_{i+2} S_{i+3} A_{i+2} A_{i+1} a_{i+1} \
& & + \; \vdots \
& & + \; B_{i+1} B_{i+2} \dots B_{n-1} S_n A_{n-1} \dots A_{i+2} A_{i+1} a_{i+1} \end{array} $$

ただし

$$ S_i = \left( \begin{matrix} 6i & 6i+1 & 6i+2 \
6i+3 & 6i+4 & 6i+5 \end{matrix} \right) $$

以下をsegment木で求めるのは容易である。

  • $A_{i-1} A_{i-2} \dots A_0$
  • $B_{i+1} A_{i+2} \dots B_n$

残る部分について、$a_{i+1}$をくくり出せば、次を求めればよい。

$$ \begin{array}{cc} & S_{i+1} \
+ & B_{i+1} S_{i+2} A_{i+1} \
+ & B_{i+1} B_{i+2} S_{i+3} A_{i+2} A_{i+1} \
+ & \vdots \
+ & B_{i+1} B_{i+2} \dots B_{n-2} S_n A_{n-2} \dots A_{i+2} A_{i+1} \
+ & B_{i+1} B_{i+2} \dots B_{n-2} B_{n-1} S_n A_{n-1} A_{n-2} \dots A_{i+2} A_{i+1} \end{array} $$

これを変形すると例えば以下のようになる。

$$ \begin{array}{crcl} & & S_{i+1} & \
+ & & B_{i+1} S_{i+2} A_{i+1} & \
+ & & B_{i+1} B_{i+2} S_{i+3} A_{i+2} A_{i+1} & \
+ & B_{i+1} B_{i+2} B_{i+3} & S_{i+4} & A_{i+3} A_{i+2} A_{i+1} \
+ & B_{i+1} B_{i+2} B_{i+3} & B_{i+4} S_{i+5} A_{i+4} & A_{i+3} A_{i+2} A_{i+1} \
+ & \vdots & \vdots & \vdots \
+ & B_{i+1} B_{i+2} B_{i+3} & B_{i+4} \dots B_{n-2} S_n A_{n-2} \dots A_{i+4} & A_{i+3} A_{i+2} A_{i+1} \
+ & B_{i+1} B_{i+2} B_{i+3} & B_{i+4} \dots B_{n-2} B_{n-1} S_n A_{n-1} A_{n-2} \dots A_{i+4} & A_{i+3} A_{i+2} A_{i+1} \end{array} $$

こう分割することにより、$A_i, B_i$単体の区間を用いて、三角形のそれの区間を合成できる。

implementation

segment木を再帰なしで実装するやつをやってみたらほぼ$2$倍速ぐらいになった: recursion $\to$ loop

#include <iostream>
#include <vector>
#include <array>
#include <functional>
#define repeat(i,n) for (int i = 0; (i) < int(n); ++(i))
using ll = long long;
using namespace std;

const int mod = 1e9+7;

template <typename T, size_t H, size_t W>
using matrix = array<array<T,W>,H>;
template <typename T, size_t A, size_t B, size_t C>
matrix<T,A,C> operator * (matrix<T,A,B> const & a, matrix<T,B,C> const & b) {
    matrix<T,A,C> c = {};
    repeat (y,A) {
        repeat (z,B) {
            repeat (x,C) {
                c[y][x] = (c[y][x] + a[y][z] *(ll) b[z][x]) % mod;
            }
        }
    }
    return c;
}
template <typename T, size_t H, size_t W>
array<T,H> operator * (matrix<T,H,W> const & a, array<T,W> const & b) {
    array<T,H> c = {};
    repeat (y,H) {
        repeat (z,W) {
            c[y] = (c[y] + a[y][z] *(ll) b[z]) % mod;
        }
    }
    return c;
}
template <typename T, size_t H, size_t W>
matrix<T,H,W> operator + (matrix<T,H,W> const & a, matrix<T,H,W> const & b) {
    matrix<T,H,W> c;
    repeat (y,H) {
        repeat (x,W) {
            c[y][x] = (a[y][x] + b[y][x]) % mod;
        }
    }
    return c;
}
template <typename T, size_t N>
array<T,N> operator + (array<T,N> const & a, array<T,N> const & b) {
    array<T,N> c;
    repeat (i,N) {
        c[i] = (a[i] + b[i]) % mod;
    }
    return c;
}
template <typename T, size_t H, size_t W>
matrix<T,H,W> matrix_zero() { return {}; }
template <typename T, size_t N>
matrix<T,N,N> matrix_unit() { matrix<T,N,N> a = {}; repeat (i,N) a[i][i] = 1; return a; }

template <typename T>
struct segment_tree { // on monoid
    int n;
    vector<T> a;
    function<T (T,T)> append; // associative
    T unit; // unit
    segment_tree() = default;
    segment_tree(int a_n, T a_unit, function<T (T,T)> a_append) {
        n = 1; while (n < a_n) n *= 2;
        a.resize(2*n-1, a_unit);
        unit = a_unit;
        append = a_append;
    }
    template <typename F>
    void point_update(int i, F func) { // 0-based
        a[i+n-1] = func(a[i+n-1]);
        for (i = (i+n)/2; i > 0; i /= 2) {
            a[i-1] = append(a[2*i-1], a[2*i]);
        }
    }
    T point_get(int i) {
        return range_concat(i, i+1);
    }
    T range_concat(int l, int r) { // 0-based, [l, r)
        T lacc = unit, racc = unit;
        for (l += n, r += n; l < r; l /= 2, r /= 2) {
            if (l % 2 == 1) lacc = append(lacc, a[(l ++) - 1]);
            if (r % 2 == 1) racc = append(a[(-- r) - 1], racc);
        }
        return append(lacc, racc);
    }
};

struct range_t {
    matrix<int,3,3> A;
    matrix<int,2,2> B;
    matrix<int,2,3> S;
};

int main() {
    int n; cin >> n;
    array<int,3> a0; repeat (i,3) cin >> a0[i];
    array<int,2> bn; repeat (i,2) cin >> bn[i];
    segment_tree<range_t> segtree(n+1, (range_t) { matrix_unit<int,3>(), matrix_unit<int,2>(), matrix_zero<int,2,3>() }, [&](range_t const & l, range_t const & r) {
        range_t m;
        m.A = r.A * l.A;
        m.B = l.B * r.B;
        m.S = l.S + l.B * r.S * l.A;
        return m;
    });
    repeat (i,n+1) {
        segtree.point_update(i, [&](range_t m) {
            repeat (y,2) repeat (x,3) m.S[y][x] = 6*i + 3*y+x;
            return m;
        });
    }
    int query; cin >> query;
    while (query --) {
        string type; int i; cin >> type >> i;
        if (type == "a") {
            segtree.point_update(i, [](range_t m) {
                repeat (y,3) repeat (x,3) cin >> m.A[y][x];
                return m;
            });
        } else if (type == "b") {
            segtree.point_update(i, [](range_t m) {
                repeat (y,2) repeat (x,2) cin >> m.B[y][x];
                return m;
            });
        } else if (type == "ga") {
            range_t m = segtree.range_concat(0, i);
            array<int,3> ai = m.A * a0;
            cout << ai[0] << ' ' << ai[1] << ' ' << ai[2] << endl;
        } else if (type == "gb") {
            range_t l = segtree.range_concat(0, i+1);
            range_t r = segtree.range_concat(i+1, n+1);
            array<int,2> bi = r.B * bn + r.S * l.A * a0;
            cout << bi[0] << ' ' << bi[1] << endl;
        }
    }
    return 0;
}