Note that the element $$$a_1$$$ cannot be removed using the described operations.

In general, the $$$j$$$-th element cannot be removed if $$$a_j$$$ is not less than all elements $$$a_1 \ldots a_{j - 1}$$$.

Otherwise, among the elements $$$a_1 \ldots a_{j - 1}$$$, there will always be an element $$$a_i \gt a_j$$$, which means that $$$a_j$$$ can be removed.

The constraints are small enough, so for each element, this statement can be checked in $$$O(n)$$$ for one check — resulting in a total of $$$O(n^2)$$$ for the test.

The idea for solving it in $$$O(n)$$$ is as follows:

• We will process the elements starting from $$$a_2$$$ to $$$a_n$$$.
• Let $$$M$$$ denote the largest among the already processed elements, initially $$$M = a_1$$$.
• If $$$M \gt a_j$$$, then $$$a_j$$$ is removed; otherwise, we set $$$M$$$ equal to $$$a_j$$$.

#include <iostream>
#include <vector>
#include <algorithm>
#include <map>
#include <ctime>
#include <queue>
#include <iomanip>
#include "assert.h"
#include <math.h>
#include <set>
#include <deque>
 
using namespace std;
 
 
#define vi vector<int>
#define vii vector<pair<int, int>>
#define pii pair<int, int>
#define ll long long
#define pll pair<ll, ll>
 
const ll INF = 2e18;
 
void solve() {
    int n;
    cin >> n;
    vi a(n);
    for (int i = 0; i < n; i++) cin >> a[i];
 
    int cmx = 0;
 
    int mx = 0;
 
    for (int i = 0; i < n; i++) {
        mx = max(mx, a[i]);
        if (a[i] == mx) cmx++;
    }
 
    int res = n - cmx;
    cout << res << "\n";
}
 
 
int main() {
    srand(time(0));
    ios::sync_with_stdio(0);
    cout.tie(0), cin.tie(0);
    int T = 1;
    cin >> T;
    while (T--) {
        solve();
    }
}

Consider all groups of consecutive zeros in the cyclic string $$$S$$$.

Let the number of zeros in these groups be $$$z_1, z_2, \ldots, z_k$$$, where $$$k$$$ is the number of groups.

In this case, the answer is $$$MZ = \max(z_j)$$$.

We will show that the answer $$$MZ$$$ is always achievable as a result:

• We can perform the operation with $$$d = 1$$$ exactly $$$MZ$$$ times.
• After each operation, all $$$z_j \gt 0$$$ will decrease by $$$1$$$, as the leftmost zero in each block will turn into a one.
• Since $$$MZ = \max(z_j)$$$, after exactly $$$MZ$$$ operations, all $$$z_j$$$ will be equal to zero.

We will show that it is impossible to obtain an answer better than $$$MZ$$$:

• Suppose there exists an answer with the sum of cyclic shifts $$$S = d_1 + d_2 + \dots + d_m$$$, where $$$S \lt MZ$$$.
• Consider the largest block of zeros (its size is $$$MZ$$$) and the one to the left of this block.
• When applying the first operation with shift $$$d_1$$$, we could not have created ones in positions to the right of $$$d_1$$$ in this block.
• When applying the subsequent cyclic shifts, this one can also only move to the right by the value of the shifts.
• In total, this one can only move $$$S$$$ positions to the right from its original position, which is strictly less than $$$MZ$$$ — thus the last zero of the block will remain a zero, contradicting the initial assumption.

#include <iostream>
#include <vector>
 
using namespace std;
 
#define vi vector<int>
#define ll long long
 
void solve() {
	int n;
	cin >> n;
	string s;
	cin >> s;
	s += s;
	n *= 2;
	int cur = 0;
 
	int res = 0;
 
	for (int i = 0; i < n; i++) {
		if (s[i] == '1') cur = 0;
		else cur++;
		res = max(res, cur);
	}
	cout << res << "\n";
}
 
int main() {
	int T;
	cin >> T;
	while (T--) {
		solve();
	}
 
}

If all the coins you have are even, then regardless of your actions, the bag will be empty every time.

Let's say we have one odd coin $$$oc$$$ and several even coins $$$ec_1, ec_2, \dots, ec_m$$$.

We will sort the even coins in descending order: $$$ec_1 \ge ec_2 \ge \dots ec_m$$$.

In this case, for $$$k = 1 \dots (m + 1)$$$, we can construct the following answers:

• $$$oc$$$
• $$$oc + ec_1$$$
• $$$oc + ec_1 + ec_2$$$
• $$$\dots$$$
• $$$oc + ec_1 + ec_2 + \dots + ec_m$$$.

But what to do if there are several odd coins?

Let $$$oc$$$ be the largest of the odd coins.

We will perform the above actions for all $$$k$$$ from $$$1$$$ to $$$(m + 1)$$$.

For $$$k = (m + 2)$$$, we will use the following construction:

• First, we will put any two odd coins that are not the coin $$$oc$$$ into the bag.
• Next, we will put the coin $$$oc$$$ and all even coins except the smallest one $$$ec_m$$$.

For $$$k = (m + 3)$$$, we will add the smallest even coin $$$ec_m$$$ to the already constructed answer.

For $$$k = (m + 4)$$$, we will add two more "unnecessary" odd coins at the beginning and again remove $$$ec_m$$$.

For $$$k = (m + 5)$$$, we will add $$$ec_m$$$ again, and so on.

Note for $$$k = n$$$: if you have an even number of odd coins in total, then regardless of the order of addition, the bag will be emptied in the end.

#include <iostream>
#include <vector>
#include <algorithm>
#include <map>
#include <ctime>
#include <queue>
#include <iomanip>
#include "assert.h"
#include <math.h>
#include <set>
#include <deque>
 
using namespace std;
 
 
#define vi vector<int>
#define vii vector<pair<int, int>>
#define pii pair<int, int>
#define ll long long
#define pll pair<ll, ll>
 
const ll INF = 2e18;
 
void solve() {
    int n;
    cin >> n;
    vi a(n);
    for (int i = 0; i < n; i++) cin >> a[i];
 
 
    vi odd, even;
 
    for (int i = 0; i < n; i++) {
        if (a[i] & 1) odd.push_back(a[i]);
        else even.push_back(a[i]);
    }
 
    sort(odd.begin(), odd.end());
    sort(even.begin(), even.end());
    reverse(even.begin(), even.end());
 
    vector<ll> prefeven((int)even.size() + 1);
 
    if (even.size() > 0) {
        prefeven[0] = 0;
        for (int i = 0; i < even.size(); i++) prefeven[i + 1] = prefeven[i] + even[i];
    }
 
    ll ans = -INF;
 
    int cntodd = (int)odd.size();
    int cnteven = 0;
    for (int i = 0; i < n; i++) {
        if (a[i] % 2 == 0) cnteven++;
    }
 
    int odd1 = 1, even1 = 0;
 
    if (odd.size() == 0) odd1 = 0, even1 = 1;
 
    for (int k = 1; k <= n; k++) {
        if (k > 1) {
            if (even1 < even.size()) even1++;
            else {
                if (odd1 + 2 <= odd.size() && even1 > 0) {
                    odd1 += 2;
                    even1--;
                }
                else {
                    odd1++;
                }
            }
        }
        if (odd1 & 1) {
            cout << odd.back() + prefeven[even1] << " ";
        }
        else {
            cout << 0 << " ";
        }
    }
    cout << '\n';
}
 
 
int main() {
    srand(time(0));
    ios::sync_with_stdio(0);
    cout.tie(0), cin.tie(0);
    int T = 1;
    cin >> T;
    while (T--) {
        solve();
    }
}

Let's count $$$dp[uv]$$$ — the number of Fibonacci paths that start with the edge ($$$u,v$$$). The transitions in such dynamics will look like $$$dp[uv]=\sum dp[vw]$$$, for such edges ($$$v,w$$$) for which $$$cost[w]=cost[u]+cost[v]$$$.

How to compute such dynamics if it depends on the order of processing edges? We could introduce a parameter $$$k$$$ — fix the length of the Fibonacci path, and then the dynamics would look like $$$dp[uv][k]=\sum dp[vw][k-1]$$$. But there is a simpler way.

Notice that in a Fibonacci path, for each subsequent edge, the sum of the numbers at its ends is strictly greater than the sum of the numbers at the ends of the previous edge. This means that edges with a greater sum of numbers should be processed strictly before edges with a smaller sum of numbers, and the order of processing edges with the same sum can be arbitrary, as they cannot be part of the same Fibonacci path.

Then we will sort the edges, and it will remain only to quickly recalculate the sum $$$dp[uv]=\sum dp[vw]$$$ from the dynamics. Naively searching for edges ($$$v, w$$$) where $$$cost[w]=cost[u]+cost[v]$$$ will not work, as the asymptotic complexity of such a solution will easily be $$$O(E^2)$$$ or worse.

Therefore, we will compress edges with the same value at the end and coming from the same vertex together. For each vertex $$$u$$$, we will maintain a dictionary $$$dp[u][cost]$$$ — the number of Fibonacci paths from vertex $$$u$$$ that go to some vertex $$$v$$$ with value $$$cost$$$.

The recalculation of this dynamics is such that $$$dp[u][cost[v]] = dp[u][cost[v]] + (1 + dp[v][cost[u] + cost[v]])$$$ — the value of the dynamics $$$dp[u][cost[v]]$$$ is incremented by the number of paths from vertex $$$v$$$ with cost $$$cost[u]+cost[v]$$$, plus one for the path from a single edge ($$$u,v$$$).

The asymptotic complexity of the solution is $$$O(E\cdot\operatorname{log}(E))$$$.

#include <bits/stdc++.h>
 
using namespace std;
 
#define vi vector<int>
#define vii vector<pair<int, int>>
#define pii pair<int, int>
#define ll long long
#define pll pair<ll, ll>
 
const ll INF = 2e18;
const int MOD = 998244353;
 
struct Edge {
    int u, v;
    ll s;
    Edge(){}
    Edge(int u, int v, ll s) : u(u), v(v), s(s){}
    bool operator < (const Edge &other) const {
        return s < other.s;
    }
};
 
void add(int &a, int b) {
    a += b;
    if (a >= MOD) a -= MOD;
}
 
 
void solve() {
    int n, m;
    cin >> n >> m;
    vector <ll> a(n);
    for (int i = 0; i < n; i++) cin >> a[i];
    vector <vi> g(n);
    vector <Edge> edges;
 
    for (int i = 0; i < m; i++) {
        int u, v;
        cin >> u >> v;
        u--;
        v--;
        g[u].push_back(v);
        edges.push_back(Edge(u, v, a[u] + a[v]));
    }
 
    sort(edges.begin(), edges.end());
    reverse(edges.begin(), edges.end());
 
    vector<map <ll, int>> sumdp(n);
 
    int ans = 0;
 
    for (auto e : edges) {
        int u = e.u;
        int v = e.v;
        ll s = e.s;
        //cout << u << " " << v << " " << s << endl;
        int curdp = sumdp[v][s];
        add(curdp, 1);
        add(sumdp[u][a[v]], curdp);
        add(ans, curdp);
    }
 
    cout << ans << "\n";
}
 
 
int main() {
    srand(time(0));
    ios::sync_with_stdio(0);
    cout.tie(0), cin.tie(0);
    int T = 1;
    cin >> T;
    while (T--) {
        solve();
    }
}

Let's first provide a theoretical estimate for the answer, followed by an algorithm that achieves this estimate.

Let's introduce $$$\operatorname{f}(l, r, x)$$$, which means the following: suppose we only consider the elements $$$a_l, a_{l + 1}, \ldots, a_{r - 1}$$$, and we can eliminate each of them using an element with a removal cost of $$$x$$$. What is the minimum removal cost for the entire range? Let's assume that we pass the best possible value to $$$x$$$.

Then let $$$p_1, p_2, \ldots, p_m$$$ be the positions of all the maxima on the interval. Note the following facts:

• For any $$$i, j$$$ for which there exists $$$k$$$ such that $$$l \leq i \lt p_k \lt j \lt r$$$, there is no sequence of operations in which one of these elements removes the other

• At least one of the elements $$$p$$$ must be removed by an element outside this interval

• The elements $$$p_1, p_2, \ldots, p_m$$$ can be removed either for the cost of $$$x$$$ or for the cost of another maximum element

Define $$$x' = min(x, c_{p_1}, c_{p_2}, \ldots, c_{_m})$$$. Then from these three remarks, we can conclude that, with the optimal choice of $$$x$$$ earlier, $$$\operatorname{f}(l, r, x) = x' \cdot m + \operatorname{f}(l, p_1, x') + \operatorname{f}(p_1 + 1, p_2, x') + \ldots + \operatorname{f}(p_m + 1, r, x')$$$.

This is true because, according to point 2, there must be at least one element that is removed by an element outside the interval. According to point 3, we can remove the remaining $$$m - 1$$$ maxima only using $$$x$$$ and other elements $$$p_i$$$, and the remaining maximum is $$$x'$$$(according to the condition). And according to point 1, all the segments are independent of each other, if we do not consider the maxima.

In this case, the answer to the problem is $$$\operatorname{f}(0, n, inf) - inf$$$.

Now, we will show that this sequence of deletions is valid.

We will also recursively remove segments, but with the condition that this minimum element outside the segment will lie strictly to the left or to the right of the segment $$$l, r$$$. Then, when we have selected $$$x'$$$, we can perform the following process:

We have a set of segments $$$(l_1, r_1), \ldots, (l_{m + 1}, r_{m + 1})$$$. We will find a segment that is a neighbor of the element $$$x'$$$, and we will recursively remove it. After that, this optimal neighbor will have one of the maximum elements as a neighbor (or we will have removed the entire segment). We remove it for $$$x'$$$.

This process will end when there is exactly one element that can be removed for $$$x'$$$. It is easy to see that we have reached our estimate, which means that the recursive algorithm above actually returns the answer to the problem.

What is the structure of the recursive calls to the function $$$\operatorname{f}$$$? Can we remember it and then use it to efficiently process requests? For which elements will the cost of deleting them change when we set it to zero?

Note that the structure of the function calls $$$\operatorname{f}$$$ is a tree. Let's save it and remember the price we paid for removing each element, let's call it $$$r_i$$$.

Then, when the request to zero $$$c_{p_i}$$$ comes, we take the node for which $$$p_i$$$ is one of the maximums and perform a depth-first search on the subtree of this node, avoiding the nodes that have already been zeroed. Thus, we change the corresponding values of $$$r_i$$$ to $$$0$$$ and recalculate the global answer.

The complexity of the solution is $$$O(n \log{n})$$$, as we need to calculate the maximum on a segment to implement the function $$$\operatorname{f}$$$, which may require a segment tree or a sparse table.

There may be other approaches to the problem.

#include <iostream>
#include <vector>
#include <set>
#include <map>
#include <algorithm>
#include <deque>
#include <queue>
#include <iomanip>
 
using namespace std;
 
#define ll long long
#define vi vector<int>
#define pii pair<int, int>
#define vii vector<pii>
 
 
const int N = 505010;
const int MOD = 1e9 + 7;
 
vi g[3 * N];
int cost1[3 * N], cost2[3 * N];
int tree[4 * N], pos[4 * N], a[N], c[N];
int nxtL[N], nxtR[N];
int num[N];
int n, q;
 
void build(int v, int tl, int tr) {
    if (tl == tr) {
        tree[v] = a[tl];
        pos[v] = tl;
        return;
    }
    int tm = (tl + tr) / 2;
    build(v * 2, tl, tm);
    build(v * 2 + 1, tm + 1, tr);
    if (tree[v * 2] > tree[v * 2 + 1]) {
        tree[v] = tree[v * 2];
        pos[v] = pos[v * 2];
    }
    else {
        tree[v] = tree[v * 2 + 1];
        pos[v] = pos[v * 2 + 1];
    }
}
 
pii getmax(int v, int tl, int tr, int l, int r) {
    if (l > r) return { -1, -1 };
    if (l == tl && r == tr) {
        return { tree[v], pos[v] };
    }
    int tm = (tl + tr) / 2; 
    pii tmp = max(getmax(v * 2, tl, tm, l, min(r, tm)), getmax(v * 2 + 1, tm + 1, tr, max(l, tm + 1), r));
    return tmp;
}
 
 
int state = -1;
ll ans;
 
void addEdge(int v, int u) {
    if (v == -1 || u == -1) return;
    g[v].push_back(u);
}
 
int it = 0;
 
int precalc(int l, int r, int minc) {
    int xx = minc;
    state++;
    int vres = state;
    if (l > r) return vres;
    pii p = getmax(1, 0, n - 1, l, r);

    vi t;
 
    int pos = p.second;
    while (pos >= l) {
        t.push_back(pos);
        pos = nxtL[pos];
    }
 
    reverse(t.begin(), t.end());
 
    pos = nxtR[p.second];
    while (pos <= r) {
        t.push_back(pos);
        pos = nxtR[pos];
    }
 
    for (int x : t) minc = min(minc, c[x]);
    for (int x : t) num[x] = vres;
    if (xx <= (int)1e9) ans += minc;
 
    cost1[vres] = minc;
 
    ans += 1ll * minc * ((int)t.size() - 1);
    addEdge(vres, precalc(l, t[0] - 1, minc));
 
    for (int i = 1; i < t.size(); i++) {
        addEdge(vres, precalc(t[i - 1] + 1, t[i] - 1, minc));
    }
 
    addEdge(vres, precalc(t.back() + 1, r, minc));
 
    for (int v : t) num[v] = vres;
    return vres;
}
 
void go(int v) {
    if (!cost1[v]) return;
    if (!g[v].size()) return;
    if (v == 1) ans -= 1ll * ((int)g[v].size() - 2) * cost1[v];
    else ans -= 1ll * ((int)g[v].size() - 1) * cost1[v];
    cost1[v] = 0;
    for (int u : g[v]) {
        go(u);
    }
 
 
}
 
void solve() {
    cin >> n;
    for (int i = 0; i < n; i++) cin >> a[i];
    for (int i = 0; i < n; i++) cin >> c[i];
    vi p(n);
    for (int i = 0; i < n; i++) cin >> p[i], p[i]--;
    map <int, int> lst;
    for (int i = 0; i < n; i++) {
        if (lst.find(a[i]) == lst.end()) nxtL[i] = -1;
        else nxtL[i] = lst[a[i]];
        lst[a[i]] = i;
    }
    lst.clear();
    for (int i = n - 1; i >= 0; i--) {
        if (lst.find(a[i]) == lst.end()) nxtR[i] = n;
        else nxtR[i] = lst[a[i]];
        lst[a[i]] = i;
    }
 
    build(1, 0, n - 1);
    for (int i = 0; i <= 3 * n; i++) g[i].clear();
    ans = 0;
    state = 0;
    int root = precalc(0, n - 1, (int)1e9 + 1);
    cout << ans << " ";
 
    for (int i = 0; i < n; i++) {
        int j = num[p[i]];
        go(j);
        cout << ans << " ";
    }
    
    cout << "\n";
}
 
int main()
{
    int T = 1;
    cin >> T;
    while (T--) {
        solve();
    }
    return 0;
}

Let $$$K$$$ denote the maximum number of unique primes in the factorization of any number under the given constraints. In this problem, $$$K \leq 6$$$, and asymptotically $$$K=O(\operatorname{log}(maxA))$$$.

Let's try to compute $$$dp[g][sum_{len}]$$$ — the number of ordered pairs ($$$i,j$$$) such that $$$\operatorname{gcd}(a_i,a_j)=g$$$, and $$$\omega(a_i)+\omega(a_j)=sum_{len}$$$. If we can compute this dynamic programming, then the answer to the problem is: $$$\sum_{g=1}^{g=n}\sum_{len=1}^{len=2 \cdot K}dp[g][sum_{len}] \cdot (sum_{len}-\omega(g))^k$$$.

A well-known method can be used to count the number of pairs of numbers in the array ($$$i,j$$$) for which $$$\operatorname{gcd}(a_i,a_j)=g$$$ for any number $$$g$$$ from $$$1$$$ to $$$maxA$$$. Let's modify this method for our problem.

Let $$$cnt[x][len]$$$ be the number of numbers in the original array that are divisible by $$$x$$$, and for which $$$\omega(x)$$$ equals $$$len$$$.

To compute this auxiliary dynamic programming, we can iterate over each $$$x$$$ from $$$1$$$ to $$$maxA$$$, and for each $$$x$$$, iterate over its multiples and add to $$$cnt[x][\omega(i \cdot x)]$$$ the count of multiples $$$i \cdot x$$$ from the array. The asymptotic complexity of this counting is $$$O(n \cdot \operatorname{log}(n))$$$ (partial sums of the harmonic series).

Now we are ready to compute the dynamic programming $$$dp[g][sum_{len}]$$$. Let's iterate over the greatest common divisor $$$g$$$ from $$$maxA$$$ down to $$$1$$$. If we simply wanted to count the number of pairs with $$$\operatorname{gcd}=g$$$, we would subtract the already counted pairs with larger multiples of $$$\operatorname{gcd}$$$.

But now we also have a fixed total length $$$\omega(x)$$$ of both numbers in the pair. Therefore, for a fixed $$$g$$$, we will also iterate over 2 parameters — $$$len_a$$$ and $$$len_b$$$, that is, $$$\omega(a)$$$ and $$$\omega(b)$$$ for both numbers.

Then $$$dp[g][sum_{len}]=dp[g][\omega(a) + \omega(b)]=dp[g][i+j]=cnt[g][i] \cdot cnt[g][j]$$$, where $$$i,j$$$ are iterated from $$$1$$$ to $$$K$$$. Here we also need to be careful and account for the case when $$$i=j$$$, in which case we add $$$\frac{cnt[g][i] \cdot (cnt[g][i] - 1)}{2}$$$ to $$$dp[g][i+j]$$$.

Next, we simply iterate over the multiples of $$$g$$$: $$$2 \cdot g, 3 \cdot g, \ldots$$$, to subtract the intersections $$$dp[g][sum_{len}] = dp[g][sum_{len}] - dp[i \cdot g][sum_{len}]$$$.

The asymptotic complexity of the solution is $$$O(maxA \cdot (K^2 + K \cdot \operatorname{log}(maxA)))$$$.

#include <bits/stdc++.h>
 
using std::cin;
using std::cout;
using vi = std::vector<int>;
using vvi = std::vector<vi>;
using ll = long long;
 
const auto ready = []()
    {
        cin.tie(0);
        std::ios_base::sync_with_stdio(false);
        return true;
    }();
 
ll binpow(ll a, ll p, ll mod)
{
    ll res = 1;
    ll mult = a;
    while (p) {
        if (p & 1) res = res * mult % mod;
        mult = mult * mult % mod;
        p >>= 1;
    }
    return res;
}
 
const ll mod = 998244353;
const int max_a = 2e5 + 10;
const int max_len = 6;
vi fi_div(max_a);
vi num_len(max_a);
 
void precalc() {
    // Here we precompute the number of unique primes in each number, i.e. w(n)
    for (int i = 2; i < max_a; ++i) {
        if (fi_div[i] == 0) {
            for (int j = i; j < max_a; j += i) {
                if (fi_div[j] == 0) fi_div[j] = i;
            }
        }
    }
 
    for (int i = 2; i < max_a; ++i) {
        int x = i;
        int len = 0;
        while (x > 1) {
            ++len;
            int tmp = fi_div[x];
            while (fi_div[x] == tmp) x /= tmp;
        }
        num_len[i] = len;
    }
}
 
void solve() {
    int n, k;
    cin >> n >> k;
 
    vi vec(n);
    for(int i = 0; i < n; ++i) cin >> vec[i];
 
    vvi cnt(n + 1, vi(max_len + 1));
 
    for(int i = 0; i < n; ++i) {
        int x = vec[i];
        cnt[vec[i]][num_len[x]] += 1;
    }
    
    // Here we calcualte auxilary dp - cnt[x][len] - how many numbers with w(n) = len are divisible by x
    for (int i = 1; i < n + 1; ++i) {
        vi dp(max_len + 1);
 
        for (int j = i; j < n + 1; j += i) {
            for(int s =0 ; s < max_len + 1; ++s) dp[s] += cnt[j][s];
        }
        cnt[i] = dp;
    }
    
    vvi dp(n + 1, vi(2 * max_len + 1));
    ll ans = 0LL;
    
    // Here we calcualte dp[g][len]
    for (int g = n; g >= 1; --g) {
        // Account for cases len_a != len_b
        for (int len1 = 0; len1 <= max_len; ++len1) {
            for (int len2 = len1 + 1; len2 <= max_len; ++len2) {
                dp[g][len1 + len2] += 1LL * cnt[g][len1] * cnt[g][len2] % mod;
                dp[g][len1 + len2] %= mod;
            }
        }
        // Account for cases len_a = len_b
        for (int len = 0; len <= max_len; ++len) {
            dp[g][len + len] += 1LL * cnt[g][len] * (cnt[g][len] - 1) / 2 % mod;
            dp[g][len + len] %= mod;
        }
        
        // Subtruct states with where gcd is multiple of current g
        for (int j = 2 * g; j < n + 1; j += g) {
            for(int s = 0; s < 2 * max_len + 1; ++s) {
                dp[g][s] -= dp[j][s];
                dp[g][s] %= mod;
            }
        }
        
        for(int s = 0; s <  2 * max_len + 1; ++s) {
            dp[g][s] %= mod;
            dp[g][s] += mod;
            dp[g][s] %= mod;
        }
    
        // Add k-th powers to the answer
        int gcd_len = num_len[g];
        for (int len = 0; len < 2 * max_len + 1; ++len) {
            int rad = len - gcd_len;
            ans += dp[g][len] % mod * binpow(rad, k, mod) % mod;
            ans %= mod;
        }
 
    }
    cout << ans << "\n";
}
 
int main()
{
    precalc();
    int t;
    cin >> t;
    while (t--) solve();
    return 0;
}