原来是矩阵快速幂，我已经忘光了....

设R， G， B各为一个三维向量，那么初始值为：

i=0, R（1， 0， 0）， G（0， 1， 0）， B（0， 0， 1） 那么对每个i，都进行一遍操作，得到R， G， B的三维向量

如例子：

i = 1， （1， 1 / 3， 2 / 3），（0， 2 / 3， 0）， （0， 0， 1 / 3）

i= 2， （1/3， 1/9， 2/9）， （2/3， 8/9， 5/9）， (0, 0, 2/9)

那么对于每个询问l, r，我们只需要知道当i = l - 1时的三维向量和i=r时的三维向量，然后把当i=l-1的三维向量看成基底来表示当i = r时的三维向量，则可以得到l-r区间的操作和

即当i=l-1, R(a, b, c), G(d, e, f), B(g, h, i)

当i = r， R1（A， B， C), G1(D, E, F), B1(G, H, I)

设R1 = xR + yG + z * B 那么（x, y, z）为l-r区间的R的向量，其值就是x+y+z

对于求x, y, z可以用线性代数学到的自己推一下，G1，B1同理

然后就AC了...

#include <bits/stdc++.h>
#define endl '\n'
#define pb push_back

using namespace std;
using ll = long long;
using db = double;
using pii = pair<int , int>;

const int mod = 1e9 + 7;
const int maxn  = 2e5 + 10;

struct node
{
	ll x, y, z;
	node operator * (const ll &d)
	{
		return {x * d % mod, y * d % mod, z * d % mod};
	}
	node operator + (const node &p)
	{
		return {x + p.x % mod, y + p.y % mod, z + p.z % mod};
	}
	node operator % (const ll &d)
	{
		return {x % mod, y % mod, z % mod};
	}
}R[maxn], G[maxn], B[maxn];

ll qpow(ll base, ll num)
{
	ll ans = 1;
	while (num)
	{
		if (num & 1)
			ans = ans * base % mod;
		base = base * base % mod;
		num >>= 1;
	}
	return ans;
}

ll gg(ll a, ll b, ll c, ll d)
{
	return (((a * b % mod) - (c * d % mod) + mod) % mod);
}

ll f(ll a, ll b, ll c, ll d, ll e, ll f, ll g, ll h, ll i, ll aa, ll bb, ll cc)
{
	ll t = a * (e * i % mod - h * f % mod + mod) % mod;
	t = (t + b * (g * f % mod - d * i % mod + mod) % mod) % mod;
	t = (t + c * (d * h % mod - g * e % mod + mod) % mod) % mod;
	ll up = aa * gg(e, i, h, f) % mod + bb * gg(h, c, b, i) % mod + cc * gg(b, f, e, c);
	up %= mod;
	up = (up + aa * gg(g, f, d, i) % mod + bb * gg(a, i, g, c) + cc * gg(d, c, a, f) % mod) % mod;
	up = (up + aa * gg(d, h, g, e) % mod + bb * gg(g, b, a, h) + cc * gg(a, e, b, d) % mod) % mod;
	return up * qpow(t, mod - 2) % mod;
}

void solve()
{
	ll n, m, i, j;
	cin >> n >> m;
	R[0] = {1, 0, 0};
	G[0] = {0, 1, 0};
	B[0] = {0, 0, 1};
	ll t = qpow(3, mod - 2);
	for (i = 1; i <= n; i++)
	{
		cin >> j;
		if (j == 1)
		{
			R[i] = R[i - 1] + G[i - 1] * t % mod + B[i - 1] * 2 * t % mod;
			G[i] = G[i - 1] * 2 * t % mod;
			B[i] = B[i - 1] * t % mod;
		}
		else if (j == 2)
		{
			G[i] = G[i - 1] + B[i - 1] * t % mod + R[i - 1] * 2 * t % mod;
			R[i] = R[i - 1] * t % mod;
			B[i] = B[i - 1] * 2 * t % mod; 
		}
		else
		{
			B[i] = B[i - 1] + R[i - 1] * t % mod + G[i - 1] * 2 * t % mod;
			G[i] = G[i - 1] * t % mod;
			R[i] = R[i - 1] * 2 * t % mod; 
		}
		R[i] = R[i] % mod;
		G[i] = G[i] % mod;
		B[i] = B[i] % mod;
	}
	int l, r;
	while (m--)
	{
		cin >> l >> r;
		cout << f(R[l - 1].x, G[l - 1].x, B[l - 1].x, 
				  R[l - 1].y, G[l - 1].y, B[l - 1].y,
				  R[l - 1].z, G[l - 1].z, B[l - 1].z,
				  R[r].x, R[r].y, R[r].z) << ' ';
		cout << f(R[l - 1].x, G[l - 1].x, B[l - 1].x, 
				  R[l - 1].y, G[l - 1].y, B[l - 1].y,
				  R[l - 1].z, G[l - 1].z, B[l - 1].z,
				  G[r].x, G[r].y, G[r].z) << ' ';
		cout << f(R[l - 1].x, G[l - 1].x, B[l - 1].x, 
				  R[l - 1].y, G[l - 1].y, B[l - 1].y,
				  R[l - 1].z, G[l - 1].z, B[l - 1].z,
				  B[r].x, B[r].y, B[r].z) << '\n';
	}
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    int T = 1;
   // cin >> T;
    while (T--)
        solve();
    
    return 0;
}

%%% //已超链接添加到题解页面(没想到还有 O(n+m) 的做法)

wawawa orz