Обратим внимание, что элемент $$$a_1$$$ нельзя удалить с помощью описанных операций.

В общем случае $$$j$$$-й элемент нельзя удалить, если $$$a_j$$$ не меньше, чем все элементы $$$a_1, \ldots, a_{j - 1}$$$.

В ином случае среди элементов $$$a_1 \ldots a_{j - 1}$$$ всегда найдётся элемент $$$a_i \gt a_j$$$, а значит $$$a_j$$$ можно будет удалить.

Ограничения достаточно маленькие, поэтому для каждого элемента можно проверить данное утверждение за $$$O(n)$$$ на одну проверку — всего получается $$$O(n^2)$$$ на тест.

Идея решения за $$$O(n)$$$ целиком следующая:

• Будем обрабатывать элементы, начиная с $$$a_2$$$, до $$$a_n$$$
• Обозначим через $$$M$$$ наибольший среди уже обработанных элементов, изначально $$$M = a_1$$$
• Если $$$M \gt a_j$$$, то $$$a_j$$$ удаляется; иначе установим $$$M$$$ равным $$$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();
    }
}

Рассмотрим все группы подряд идущих нулей в зацикленной строке $$$S$$$.

Пусть количества нулей в данных группах равны $$$z_1, z_2, \ldots, z_k$$$, где $$$k$$$ — количество групп.

В таком случае ответ равен $$$MZ = \max(z_j)$$$.

Покажем, что ответ $$$MZ$$$ всегда достижим как ответ:

• $$$MZ$$$ раз выполним операцию с $$$d = 1$$$.
• После каждой операции все $$$z_j \gt 0$$$ будут уменьшаться на $$$1$$$, так как самый левый ноль в каждом блоке будет превращаться в единицу.
• Так как $$$MZ = \max(z_j)$$$, то после ровно $$$MZ$$$ операций все $$$z_j$$$ будут равны нулю.

Покажем, что нельзя получить ответ лучше, чем $$$MZ$$$:

• Предположим, что существует ответ с суммой значений циклических сдвигов $$$S = d_1 + d_2 + \dots + d_m$$$, где $$$S \lt MZ$$$.
• Рассмотрим наибольший блок из нулей (его размер равен $$$MZ$$$) и единицу слева от этого блока.
• При применении первой операции со сдвигом $$$d_1$$$ у нас не могло появиться единицы на позиции правее $$$d_1$$$ в этом блоке.
• При применении следующих циклических сдвигов данная единица также может переместиться правее только на значение сдвигов.
• Суммарно данная единица сможет переместиться только на $$$S$$$ позиций вправо от изначального положения, что строго меньше $$$MZ$$$ — таким образом последний ноль блока останется нулём, что противоречит исходному предположению.

#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();
	}
 
}

Если все имеющиеся монеты чётные, то независимо от ваших действий мешок будет пустым каждый раз.

Пусть у нас есть одна нечётная монета $$$oc$$$ и несколько чётных $$$ec_1, ec_2, \dots, ec_m$$$.

Отсортируем чётные монеты в порядке убывания: $$$ec_1 \ge ec_2 \ge \dots ec_m$$$.

В таком случае для $$$k = 1 \dots (m + 1)$$$ мы можем построить следующие ответы:

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

Но что делать, если нечётных монет несколько?

Пусть $$$oc$$$ — наибольшая из нечётных монет.

Проделаем для всех $$$k$$$ от $$$1$$$ до $$$(m + 1)$$$ описанные выше действия.

Для $$$k = (m + 2)$$$ используем следующую конструкцию:

• Вначале положим в мешок любые две нечётные монеты, не являющиеся монетой $$$oc$$$.
• Далее положим монету $$$oc$$$ и все чётные монеты, за исключением наименьшей $$$ec_m$$$.

Для $$$k = (m + 3)$$$ добавим к уже построенному ответу наименьшую чётную монету $$$ec_m$$$.

Для $$$k = (m + 4)$$$ добавим ещё две "ненужные" нечётные монеты в начало и опять уберём $$$ec_m$$$.

Для $$$k = (m + 5)$$$ опять добавим $$$ec_m$$$ и так далее.

Обратите внимание на $$$k = n$$$: если всего у вас имеется чётное количество нечётных монет, то независимо от порядка добавления мешок будет опустошён в конце.

#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();
    }
}

Давайте посчитаем $$$dp[uv]$$$ — количество путей Фибоначчи, которые начинаются с ребра ($$$u,v$$$). Переходы в такой динамике будут выглядеть как $$$dp[uv]=\sum dp[vw]$$$, по таким рёбрам ($$$v,w$$$) для которых выполняется $$$cost[w]=cost[u]+cost[v]$$$.

Как посчитать такую динамику, если она зависит от порядка обработки рёбер? Мы могли бы ввести параметр $$$k$$$ — зафиксировать длину пути Фибоначчи, и тогда динамика выглядела бы как $$$dp[uv][k]=\sum dp[vw][k-1]$$$. Но есть более простой способ.

Заметим, что в пути Фибоначчи у каждого следующего ребра сумма чисел на его концах строго больше, чем сумма чисел на концах предыдущего ребра. Это означает, что рёбра с большей суммой чисел выгодно обработать строго раньше, чем рёбра с меньшей суммой чисел, а порядок обработки рёбер с одинаковой суммой может быть произвольным, так как они не могут входить в один путь Фибоначчи.

Тогда мы отсортируем рёбра, и останется только быстро пересчитывать сумму $$$dp[uv]=\sum dp[vw]$$$ из динамики. Наивно искать рёбра ($$$v, w$$$), у которых $$$cost[w]=cost[u]+cost[v]$$$ не получится, асимптотика такого решения легко будет $$$O(E^2)$$$ или хуже.

Поэтому сожмём рёбра с одинаковым значением на конце и идущие из одной вершины вместе. Для каждой вершины $$$u$$$ заведём словарь $$$dp[u][cost]$$$ — число путей Фибоначчи из вершины $$$u$$$, которые идут в какую-то вершину $$$v$$$ со значением $$$cost$$$.

Пересчет этой динамики такой $$$dp[u][cost[v]] = dp[u][cost[v]] + (1 + dp[v][cost[u] + cost[v]])$$$ — к значению динамики $$$dp[u][cost[v]]$$$ прибавляется число путей из вершины $$$v$$$ со стоимостью $$$cost[u]+cost[v]$$$, и ещё плюс единица за путь из одного ребра ($$$u,v$$$).

Асимптотика решения $$$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();
    }
}

Давайте вначале дадим теоретическую оценку ответа, затем алгоритм, который достигает этой оценки.

Введем $$$\operatorname{f}(l, r, x)$$$, которая означает следующее: пускай мы рассматриваем только элементы $$$a_l, a_{l + 1}, \ldots, a_{r - 1}$$$ и мы можем каждый из них вычеркнуть при помощи элемента, который имеет стоимость удаления $$$x$$$. Тогда какая минимальная стоимость удаления всего данного отрезка? Будем считать, что мы передаем в $$$x$$$ самое лучшее подходящее значение.

Тогда пусть $$$p_1, p_2, \ldots, p_m$$$ — это позиции всех максимумов на отрезке. Заметим следующие факты:

• Для любых $$$i, j$$$, для которых существует $$$k$$$ такой, что $$$l \leq i \lt p_k \lt j \lt r$$$, не существует последовательности операций, при которой один из этих элементов удаляет другой

• Хотя бы один из элементов $$$p$$$ обязан быть удален элементом вне этого отрезка

• Элементы $$$p_1, p_2, \ldots, p_m$$$ могут быть удалены либо за стоимость $$$x$$$, либо за стоимость другого максимального элемента

Определим $$$x' = min(x, c_{p_1}, c_{p_2}, \ldots, c_{_m})$$$. Тогда из этих трех замечаний можно сделать вывод, что при оптимальном выборе $$$x$$$ ранее, $$$\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')$$$.

Это так, ведь по пункту 2 должен быть хотя бы один элемент, удаляемый элементом вне отрезка. По пункту 3 мы оставшиеся $$$m - 1$$$ максимумов можем удалить только при помощи $$$x$$$ и других элементов $$$p_i$$$, а оставшийся максимум — за $$$x'$$$(по условию). А по пункту 1 все отрезки между собой независимые, если не учитывать максимумы.

И в этом случае ответом на задачу будет $$$\operatorname{f}(0, n, inf) - inf$$$.

Теперь покажем, что такая последовательность удалений действительно существует.

Будем также рекурсивно удалять отрезки, но с условием, что этот минимальный элемент вне отрезка будет лежать строго слева либо справа от отрезка $$$l, r$$$. Тогда, когда мы выделили $$$x'$$$, мы можем делать следующий процесс:

У нас есть множество отрезков $$$(l_1, r_1), \ldots, (l_{m + 1}, r_{m + 1})$$$. У нас найдется отрезок, который является соседом элемента $$$x'$$$, рекурсивно удалим его. После этого у этого оптимального соседа будет соседом один из максимальных элементов (либо же мы удалили весь отрезок). Мы его удаляем за $$$x'$$$.

Этот процесс закончится, когда останется ровно один элемент, который можно удалить за $$$x'$$$. И несложно заметить, что мы достигли нашей оценки, а значит, рекурсивный алгоритм выше действительно возвращает ответ на задачу.

Что представляет из себя структура рекурсивных вызовов функции $$$\operatorname{f}$$$? Можем ли мы запомнить её, а затем по ней эффективно обрабатывать запросы? Для каких элементов изменится стоимость, за которую мы их удаляем, при обнулении некоторого $$$c_i$$$?

Заметим, что структура вызовов функции $$$\operatorname{f}$$$ представляет собой дерево. Давайте сохраним его и запомним, за какую цену мы удаляли каждый элемент, пусть это $$$r_i$$$.

Затем, когда приходит запрос обнуления $$$c_{p_i}$$$, возьмем узел, для которого $$$p_i$$$ является одним из максимумов, и запустим поиск в глубину по поддереву этой вершины, не заходя в вершины, для которых уже было выполнено обнуление. Таким образом, поменяем соответствующие значения $$$r_i$$$ на $$$0$$$ и пересчитаем глобальный ответ.

Сложность решения: $$$O(n \log{n})$$$, так как для реализации функции $$$\operatorname{f}$$$ нам нужно вычислять максимум на отрезке, для чего может понадобиться дерево отрезков или sparce table.

Возможны и другие подходы по задаче.

#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;
}

Обозначим за $$$K$$$ наибольшее количество уникальных простых в разложении любого числа в заданных ограничениях. В этой задаче $$$K \leq 6$$$, и асимптотически $$$K=O(\operatorname{log}(maxA))$$$.

Давайте попробуем посчитать $$$dp[g][sum_{len}]$$$ — количество упорядоченных пар ($$$i,j$$$) таких, что $$$\operatorname{gcd}(a_i,a_j)=g$$$, а $$$\omega(a_i)+\omega(a_j)=sum_{len}$$$. Если мы умеем считать эту динамику, то ответ на задачу — это $$$\sum_{g=1}^{g=n}\sum_{len=1}^{len=2 \cdot K}dp[g][sum_{len}] \cdot (sum_{len}-\omega(g))^k$$$.

Довольно известным способом можно посчитать количество пар чисел массива ($$$i,j$$$), у которых $$$\operatorname{gcd}(a_i,a_j)=g$$$ для любого числа $$$g$$$ от $$$1$$$ до $$$maxA$$$. Давайте модифицируем этот способ для нашей задачи.

Пусть $$$cnt[x][len]$$$ — количество чисел исходного массива, которые делятся на $$$x$$$, и $$$\omega(x)$$$ для которых равен $$$len$$$.

Чтобы посчитать эту вспомогательную динамику, мы можем проитерироваться по каждому $$$x$$$ от $$$1$$$ до $$$maxA$$$, для каждого $$$x$$$ проитерироваться по его кратным и прибавить к $$$cnt[x][\omega(i \cdot x)]$$$ количество кратных $$$i \cdot x$$$ из массива. Асимптотика этого подсчета $$$O(n \cdot \operatorname{log}(n))$$$ (частичные суммы гармонического ряда).

Теперь мы готовы посчитать динамику $$$dp[g][sum_{len}]$$$. Давайте проитерируемся по наибольшему общему делителю $$$g$$$ от $$$maxA$$$ до $$$1$$$. Если бы мы просто хотели посчитать число пар с $$$\operatorname{gcd}=g$$$, то мы бы вычли уже посчитанные пары с большим кратным $$$\operatorname{gcd}$$$.

Но теперь у нас еще фиксирована суммарная длина $$$\omega(x)$$$ обоих чисел пары. Поэтому при фиксированном $$$g$$$ переберем еще 2 параметра — $$$len_a$$$ и $$$len_b$$$, то есть $$$\omega(a)$$$ и $$$\omega(b)$$$ для обоих чисел.

Тогда $$$dp[g][sum_{len}]=dp[g][\omega(a) + \omega(b)]=dp[g][i+j]=cnt[g][i] \cdot cnt[g][j]$$$, где $$$i,j$$$ перебираются от $$$1$$$ до $$$K$$$. Здесь также надо быть аккуратным и учесть случай, когда $$$i=j$$$, в таком случае к $$$dp[g][i+j]$$$ прибавляется $$$\frac{cnt[g][i] \cdot (cnt[g][i] - 1)}{2}$$$.

Далее мы просто вложенным циклом перебираем кратные $$$g$$$: $$$2 \cdot g, 3 \cdot g, \ldots$$$, чтобы вычесть пересечения $$$dp[g][sum_{len}] = dp[g][sum_{len}] - dp[i \cdot g][sum_{len}]$$$.

Асимптотика решения $$$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;
}

...You are given an array of up to 500000 numbers which can be up to 10^30. Define the cost of the subset of the array as sum of elements of this subset plus gcd of elements of this subset. For every k from 1 to n print the maximum possible cost of some subset of size k. TL:8 sec. (By the way on the municipal stage there was only k=2).