(Analysis by Sujay Konda)

First, let's analyze what the operation does to the bit string in the compressed format (the lengths of maximal consecutive blocks of equal characters). The operation is equivalent to taking two adjacent runs of the same size and reversing that part of the bitstring. Once you reverse the bitstring, the runs will merge to the runs outside (unless there are no runs outside). Let's look at example of this: 

11(0011)00
2,(2,2),2
->
11(1100)00
2+2,2+2

For the edge case, then: 

(1100)111
(2,2),3
->
(0011)111
2,2+3

We can actually make the edge case the same as the normal case, by padding the runs with $0$s at the start and at the end, so now 

(1100)111
0,(2,2),3,0
->
(0011)111
0+2,2+3,0

This motivates a range DP, where $dp[l][r] = \text{# of ways to remove all the runs in strictly between } l \text{ and } r$. Note that we don't need to store the minimum number of operations, since every operation erases two runs, so the min ops is just $(r - l - 1) / 2$.

For the transitions, since right before reaching the state where everything between $l$ and $r$ is removed, we must be in a state where everything between $l$ and $r$ except $a$ and $b$ has been removed. To remove $v_a$ and $v_b$, we need $v_a = v_b$. To find the value of $v_a$ and $v_b$, we can return to the normal bitstring format and notice that will look like 

[run of 1s][run of 0s][run of 1s][run of 0s] 
   v_l,       v_a,       v_b,        v_r
(or inverted).

Note that all the $0$s between $l$ and $a$ will connect with the run of $0$s, the $1$s between $b$ and $r$ will connect with the run of 1s and both the $0$s and $1$s in between $a$ and $b$ will connect with the run of 0s and run of 1s respectively. Therefore, $v_a = \text{# of } 0\text{s in } [l,b]$ and $v_b = \text{# of }1\text{s in }[a,r]$ (to check if it is the inverted case, it depends on the character for the run at $a$).

Now define $F(a,b,c) = \frac{(a+b+c)!}{a!b!c!}$ (i.e. the number of ways to interleave 3 sequences of operations of length $a, b, c$), then the transitions are

$$dp[l][r] = \sum_{l < a < b < r, v_a = v_b} dp[l][a] \cdot dp[a][b] \cdot dp[b][r] \cdot F((a-l-1)/2, (b-a-1)/2, (r-b-1)/2)$$

To extract the answer, we loop through all $l,r$ such that everything outside $l$ and $r$ is padding, i.e. $l$ and $r$ are either in the padding or are the leftmost/rightmost original runs. We can update min ops by maintaining a second boolean dp to see if removing $l, r$ is possible while also keeping the original # of ways dp to count the number of ways.

This almost gives us an $O(R^4)$ algorithm, but not entirely. We still haven't shown how much padding we need to add for the algorithm to work. For subtask 3, since the min ops equals $R/2 - 1$, we don't need any padding because we remove exactly $R-2$ runs leaving exactly $2$ runs left. For the other subtasks, you can easily bound the padding by $R$, but you can actually do better. Suppose you do some edge operation 

0,0,a,a,b,...
->
0,a,b+a,....
Now note that to ever do another edge operation, we first need to at least delete the b+a element, which will at least turn $a$ into $2a + b$. Therefore, the value at the edge at least doubles every time we need one extra padding. Therefore, we can bound the padding by $\log N \approx 30$. This gives us a $O((R + \log N)^4)$ algorithm.

To speed it up further, notice that $(a,b)$ works with $(l,r)$ if

$$\text{# of ones before }a + \text{# of zeros before }b = \text{# of ones before }r + \text{# of zeros before }l.$$
This means as we increase $a$, we increase the $\text{# of ones before }a$, so we must decrease $b$ to decrease $\text{# of zeros before }b$ to balance. Therefore, every $a$ pairs with one $b$ (other than padding cases). Note that the formula and logic is the exact same in the inverted case, except the $0$s are replaced with $1$s and vice versa. This leads to an $O((R + \log N)^3 + \log^4 N)$ algorithm (the padding case contributes the $\log^4 N$). Note that you only need to test even-sized ranges because you remove 2 at a time, which improves the constant factor dramatically, which is why $R=800$ here.

My code: 

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

const int LGN = 30;
const int INF = 1e9;
const int MOD = 1e9 + 7;
using ll = long long;

const int MXM = 1000;

int f[MXM + 1], invf[MXM + 1];

int mulm(int x, int y) {
    return (ll)x * y % MOD;
}

int bpow(int x, int y) {
    return (y == 0 ? 1 : mulm(bpow(mulm(x, x), y / 2), (y % 2 ? x : 1)));
}

int choose(int n, int k) {
    if(k > n) return 0;
    if(k < 0) return 0;
    assert(n >= 0);
    return mulm(mulm(f[n], invf[k]), invf[n - k]);
}

void tc() {
    int M; cin >> M; char c; cin >> c;
    vector<int> a;
    for(int i = 0; i < LGN; i++)
        a.push_back(0);
    for(int i = 0; i < M; i++) {
        int ai; cin >> ai;
        a.push_back(ai);
    }
    for(int i = 0; i < LGN; i++)
        a.push_back(0);

    vector<int> p(M + 2 * LGN);
    for(int i = 2; i < M + 2 * LGN; i++) {
        p[i] = a[i] + p[i - 2];
    }

    vector<vector<bool>> dp(M + 2 * LGN, vector<bool>(M + 2 * LGN));
    vector<vector<int>> dp2(M + 2 * LGN, vector<int>(M + 2 * LGN));

    map<int, vector<pair<int, int>>> trans;

    // add possible transitions from l, r
    auto add_trans = [&] (int l, int r) {
        if(r < LGN || l >= M + LGN) return;

        trans[p[r - 1] + p[l - 1]].push_back({l, r});
    };
    
    for(int i = 0; i < M + 2 * LGN - 1; i++) {
        dp[i][i + 1] = true;
        dp2[i][i + 1] = 1;
        add_trans(i, i + 1);
    }

    for(int sz = 4; sz <= M + 2 * LGN; sz += 2) {
        for(int l = 1; l + sz - 1 < M + 2 * LGN; l++) {
            int r = l + sz - 1;
            for(auto [u, v] : trans[p[r - 1] + p[l - 1]]) {
               if(l < u && v < r && dp[l][u] && dp[v][r]) {
                    dp[l][r] = true;
                    int ways = mulm(f[(r - l) / 2 - 1], mulm(invf[(r - v) / 2], mulm(invf[(v - u) / 2], invf[(u - l) / 2])));
                    dp2[l][r] += mulm(ways, mulm(dp2[u][v], mulm(dp2[l][u], dp2[v][r])));
                    dp2[l][r] %= MOD;
                } 
            }
            if(dp[l][r])
                add_trans(l, r);
        }
    }
    int ans = INF;
    int ans2 = 0;
    for(int l = 0; l <= LGN; l++) {
        for(int r = M + LGN - 1; r < M + 2 * LGN; r++) {
            if(dp[l][r] && (LGN - l) % 2 == (c == '0')) {
                if ((r - l) / 2 < ans) {
                    ans = (r - l) / 2;
                    ans2 = 0;
                }
                if((r - l) / 2 == ans) {
                    ans2 += dp2[l][r];
                    ans2 %= MOD;
                }
            }
        }
    }
    cout << ((ans == INF) ? -1 : ans) << " ";
    cout << ans2 << endl;
}

int main() {

    f[0] = 1;
    for(int i = 1; i <= MXM; i++) {
        f[i] = mulm(f[i - 1], i);
    }
    invf[MXM] = bpow(f[MXM], MOD - 2);
    for(int i = MXM; i >= 1; i--) {
        invf[i - 1] = mulm(invf[i], i);
    }

    ios::sync_with_stdio(false), cin.tie(nullptr);
    int T; cin >> T;
    while(T--) tc();
}

Bonus: Solve the problem in $O(R^3)$ time.