To solve this problem optimally, we need to ensure that we maximize the size of the set $$$S$$$ while satisfying the given constraints on XOR and OR operations.

A necessary condition for a valid set is that $$$ y $$$ must be a subset of $$$ v $$$. In terms of bits, this means every bit that is set in $$$ y $$$ must also be set in $$$ v $$$. If $$$ y \gt x $$$, then at least one bit in $$$ y $$$ is missing from $$$ x $$$, making it impossible to construct a valid set where the XOR of elements results in $$$ y $$$. Thus, in such cases, the answer is immediately $$$ -1 $$$.

To understand how we construct valid sets, consider forming all possible subsets of the bits set in $$$ v $$$. The total number of such subsets is $$$ 2^{pc} $$$, where $$$ pc $$$ is the number of set bits in $$$ v $$$. Each individual bit contributes exactly $$$ 2^{(pc - 1)} $$$ times across all subsets, meaning the overall XOR of the entire set cancels out to $$$ 0 $$$ when $$$ pc \geq 2 $$$.

We begin by ensuring that our OR value is exactly $$$ y $$$, so that all required bits are present in our final set. Then, we greedily attempt to maximize the number of set bits by starting from the least significant bit (LSB) and setting as many bits as possible while ensuring of $$$ v \leq x $$$.

Now, consider the number of set bits, denoted as $$$ pc $$$:

1. Case $$$ pc = 1 $$$: The set will only contain $$$ { 0, v } $$$, meaning the XOR remains nonzero. Thus, we return $$$ 1 $$$ if $$$ y = 0 $$$, otherwise, $$$ 2 $$$.

1. Case $$$ pc \geq 2 $$$: The full set of $$$ 2^{pc} $$$ elements XORs to $$$ 0 $$$. Since we need to ensure the XOR of the subset is $$$ y $$$, we must exclude exactly one element when $$$ y \neq 0 $$$ to achieve the correct XOR sum.

This ensures that the largest possible set size is achieved while maintaining the required XOR condition.

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

void solve() 
{
    ll x, y;
    cin >> x >> y;
    
    if (y > x) 
    {
        cout << -1 << endl;
        return;
    }
    
    ll v = y;

    for(ll bit = 0; bit < 61; bit++)
    {
        if ((v | (1LL << bit)) <= x)   
            v |= (1LL << bit);
    }

    ll answer = (1LL << __builtin_popcountll(v));

    if ((__builtin_popcountll(v) == 1 && y == 0)) answer--;
    if ((__builtin_popcountll(v) > 1 && y != 0))  answer--;

    cout << answer << endl;
}

int main()
{
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);
    cout.tie(NULL);

    ll t;
    cin >> t;
    while (t--)
    {
        solve();
    }

    return 0;
}

Main Observation: A solution always exists.

We need to construct a permutation $$$a$$$ of length $$$n$$$ such that:

We can model the problem as a bipartite matching problem. Consider a graph with $$$n$$$ nodes where:
1. The odd numbers from $$$1$$$ to $$$n$$$ form one partition.
2. The even numbers from $$$1$$$ to $$$n$$$ form the other partition.
3. An edge exists between an odd node $$$u$$$ and an even node $$$v$$$ if $$$u + v$$$ is prime.

Using Hall’s theorem, we can prove that a perfect matching always exists. Specifically:
1. The largest prime gap up to $$$n = 10^6$$$ is $$$114$$$.
2. More generally, for any $$$n$$$, the largest prime gap is much smaller than $$$n$$$.
3. Thus, for any subset of odd vertices $$$X$$$, the number of adjacent even vertices is always at least $$$|X|$$$, ensuring a perfect matching.

Construction

There are a bunch of constructions for this problem. Here are some:

Construction 1

1. Find the smallest $$$i$$$ such that $$$i + n$$$ is prime.
2. Assign the numbers in pairs:
   • Iteratively increase $$$k$$$ and set $$$a_k = i+k$$$ and $$$a_{n-k} = n-k$$$.
3. Perform bipartite matching on the remaining numbers from $$$0$$$ to $$$i-1$$$.

Time Complexity: $$$O({P_n}^3+n)$$$ where $$$P_n$$$ is highest prime gap till $$$n$$$.

Construction 2 (IceKnight1093’s Approach)

1. Find the smallest $$$i$$$ such that $$$i + n$$$ is prime.
2. Assign numbers from $$$i$$$ to $$$n$$$ using same approach as before.
3. Recursively apply the same strategy for numbers from $$$1$$$ to $$$i-1$$$.
4. This approach can be justified using Bertrand’s Postulate, which ensures the existence of a prime in necessary ranges.

Time Complexity: $$$O(n)$$$

#include<bits/stdc++.h>
using namespace std;
 
typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef set<int> si;
typedef pair<int, int> pi;
typedef pair<ll, ll> pl;
#define ordered_set tree<ll, null_type,less<ll>, rb_tree_tag,tree_order_statistics_node_update> 
#define fo(i, a, b) for (ll i = a; i < b; i++)
#define F first
#define S second
#define pb push_back
#define mp make_pair
#define no cout << "No" << "\n"
#define yes cout << "Yes" << "\n"
#define cr(a,b) a.clear(); a.resize(b);
ll p = 998244353;
#define mod p
#define MAX 100001
#define NIL 0
#define INF (1<<28)
 
vector<bool> prime;
void sieve(int n)
{
    prime.clear();
    prime.resize(n + 1);
    for (int p = 0; p <= n; p++) {
        prime[p] = 1;
    }
    for (int p = 2; p * p <= n; p++) {
        if (prime[p] == true) {
            for (int i = p * p; i <= n; i += p)
                prime[i] = false;
        }
    }
}

vector< int > G[MAX];
int n, m, match[MAX], dist[MAX];
// n: number of nodes on left side, nodes are numbered 1 to n
// m: number of nodes on right side, nodes are numbered n+1 to n+m
// G = NIL[0] ∪ G1[G[1---n]] ∪ G2[G[n+1---n+m]]
 
bool bfs() {
    int i, u, v, len;
    queue< int > Q;
    for(i=1; i<=n; i++) {
        if(match[i]==NIL) {
            dist[i] = 0;
            Q.push(i);
        }
        else dist[i] = INF;
    }
    dist[NIL] = INF;
    while(!Q.empty()) {
        u = Q.front(); Q.pop();
        if(u!=NIL) {
            len = G[u].size();
            for(i=0; i<len; i++) {
                v = G[u][i];
                if(dist[match[v]]==INF) {
                    dist[match[v]] = dist[u] + 1;
                    Q.push(match[v]);
                }
            }
        }
    }
    return (dist[NIL]!=INF);
}
 
bool dfs(int u) {
    int i, v, len;
    if(u!=NIL) {
        len = G[u].size();
        for(i=0; i<len; i++) {
            v = G[u][i];
            if(dist[match[v]]==dist[u]+1) {
                if(dfs(match[v])) {
                    match[v] = u;
                    match[u] = v;
                    return true;
                }
            }
        }
        dist[u] = INF;
        return false;
    }
    return true;
}
 
int hopcroft_karp() {
    int matching = 0, i;
    // match[] is assumed NIL for all vertex in G
    while(bfs())
        for(i=1; i<=n; i++)
            if(match[i]==NIL && dfs(i))
                matching++;
    return matching;
}
 
int main() {
    cin.tie(0)->sync_with_stdio(0);
    cin.exceptions(cin.failbit);
    ll k=2;
    int m=3e6;
    sieve(m);
    memset(match,0,sizeof(match));
    fo(i,0,MAX)G[i].clear();
    ll T, n1;
    cin >> T;
    while (T--) {
        fo(i,0,n1+1){
            G[i].clear();
        }
        fo(i,0,n1+1)match[i]=0;
        
        cin >> n1;
        int k2=k;
        while(!prime[k2+n1]){
            k2++;
        }
        if(n1<=k2){
            vi ans(n1);
            vi odd,even;
            if(n1 %2==1){
                ans[n1/2]=1;
                for(int i=3;i<=n1;i+=2){
                    odd.pb(i);
                }
                for(int i=2;i<=n1;i+=2){
                    even.pb(i); 
                }
            }
            else{
                for(int i=1;i<=n1;i+=2){
                    odd.pb(i);
                }
                for(int i=2;i<=n1;i+=2){
                    even.pb(i); 
                }
            }
            int j=0;
            int i1=0,j1=0;
            n=n1/2;
            m=n1/2;
            for(auto x:odd){
                i1++;
                j1=0;
                for(auto y:even){
                    j1++;
                    if(prime[x+y]){
                        G[i1].pb(j1+n1/2);
                    }
                }
            }
            
            int matching=hopcroft_karp();
            fo(i,1,n1/2+1){
                ans[j]=odd[i-1];
                ans[n1-1-j]=even[match[i]-n1/2-1];
                j++;
            }
            fo(i,0,n1){
                cout <<ans[i]<< " ";
            }
            cout << "\n";
        }
        else{
            vi ans(n1);
            vi odd,even;
            k2-=1;
            if(n1%2==1){
                ans[n1/2]=1;
                for(int i=3;i<=k2;i+=2){
                    odd.pb(i);
                }
                for(int i=2;i<=k2;i+=2){
                    even.pb(i); 
                }
            }
            else{
                for(int i=1;i<=k2;i+=2){
                    odd.pb(i);
                }
                for(int i=2;i<=k2;i+=2){
                    even.pb(i); 
                }
            }
            int j=0;
            int i1=0,j1=0;
            n=k2/2;
            m=k2/2;
            for(auto x:odd){
                i1++;
                j1=0;
                for(auto y:even){
                    j1++;
                    if(prime[x+y]){
                        G[i1].pb(j1+k2/2);
                    }
                }
            }
            
            int matching=hopcroft_karp();
            fo(i,1,k2/2+1){
                ans[j]=odd[i-1];
                ans[n1-1-j]=even[match[i]-k2/2-1];
                j++;
            }
            k2++;
            int k3=n1;
            while(k2<k3){
                ans[j]=k2;
                ans[n1-1-j]=k3;
                k2++;
                k3--;
                j++;
            }
            fo(i,0,n1){
                cout <<ans[i]<< " ";
            }
            cout << "\n";
 
        }
    }
}

// #include <bits/allocator.h>
// #pragma GCC optimize("O3,unroll-loops")
// #pragma GCC target("avx2,bmi,bmi2,lzcnt,popcnt")
#include "bits/stdc++.h"
using namespace std;
using ll = long long int;
mt19937_64 RNG(chrono::high_resolution_clock::now().time_since_epoch().count());

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

    int t; cin >> t;
    while (t--) {
        int n; cin >> n;

        vector<int> mark(2*n + 1), primes, ans(n, 1);
        for (int i = 2; i <= 2*n; ++i) {
            if (mark[i]) continue;
            primes.push_back(i);
            for (int j = 2*i; j <= 2*n; j += i)
                mark[j] = 1;
        }
        
        int L = 0, R = n-1, x = n;
        while (x > 1) {
            int p = *ranges::upper_bound(primes, x);
            for (int i = 0; x-i > p-x+i; ++i) {
                ans[L++] = x-i;
                ans[R--] = p-x+i;
            }
            x = p-x-1;
        }
        for (int y : ans) cout << y << ' ';
        cout << '\n';
    }
}

Due to the condition of sum of diagonals being equal, the difference between adjacent elements is equal for each row and each column. Hence rectangle can be uniquely identified by any starting number(here we are using the top left as starting number) and the $$$l-1$$$ operations along rows and $$$b-1$$$ operations along columns. For the other condition ($$$max-min=(l-1)+(b-1)$$$) to satisfy, a bad rectangle must have a monotonically increasing/decreasing sequence along both it's rows and columns.

We need to count over all patterns of length $$$n$$$, here we assume starting from $$$0$$$ and extend digits from $$$-9$$$ to $$$+9$$$ along with maintaining the highest and lowest number that was achieved in the sequence. For any starting digit $$$s$$$ and following pattern $$$p1$$$ along length and $$$p2$$$ along breadth of grid, whose maximum were $$$mx1$$$ and $$$mx2$$$ respectively will achieve a highest value of $$$s+mx1+mx2$$$ which should be $$$ \lt 10$$$ and whose lowest were $$$mn1$$$ and $$$mn2$$$ respectively will achieve the lowest value of $$$x+mn1+mn2( \gt =0)$$$. Along with keeping fix the range of digits $$$0-9$$$, we also need to maintain the maximum length streak (max consecutive increasing or decreasing) in the pattern. $$$2$$$ patterns having streaks $$$s1$$$ and $$$s2$$$ arranged along two sides of this grid will have $$$s1*s2( \lt =c)$$$ as the area of the largest bad rectangle.

The final dp state is dp[ith element from the start][current diff from starting][max value upto i][min value][cur streak][max streak] representing the count of patterns for first $$$i$$$ digits starting from $$$0$$$ extending from min to max element which is currently in streak of cur streak and has achieved a maximum streak of max streak upto $$$i$$$.

Using this we can find the sum over all starting digits $$$s=0...9$$$ and over all pairs of patterns $$$p1,p2$$$ satisfying $$$s+mx1+mx2 \lt 10$$$ and $$$s+mn1+mn2 \gt =0$$$ and $$$s1*s2 \lt =c$$$.

Final Time complexity: $$$n*20*10*10*20*10$$$

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
# define FAST ios_base::sync_with_stdio(false); cin.tie(NULL); cout.tie(NULL);
const ll MOD = 1e9+7;
int dp[20][10][10][20][10]={0};
int dpnew[20][10][10][20][10]={0};
int main(){
    FAST
    ll n, c;
    cin >> n >> c;
 
    dp[9][0][0][9][0] = 1;
 
    for(int a = 1; a < n; a++){
        memset(dpnew,0,sizeof(dpnew));
 
        for(int lastd = -9; lastd <= 9; lastd++){
            for(int mx = 0; mx < 10; mx++){
                for(int mn = 0; mn < 10; mn++){
                    for(int streak = -9; streak <= 9; streak++){
                        for(int mxstrk = 0; mxstrk < 10; mxstrk++){
                            ll val = dp[lastd + 9][mx][mn][streak + 9][mxstrk];
                            if((lastd < 9) && (streak < 9)){
                                dpnew[lastd + 1 + 9][max(mx, lastd + 1)][mn][max(1, streak + 1) + 9][max(mxstrk, max(1, streak + 1))] = (dpnew[lastd + 1 + 9][max(mx, lastd + 1)][mn][max(1, streak + 1) + 9][max(mxstrk, max(1, streak + 1))] + val) % MOD;
                            }
                            dpnew[lastd + 9][mx][mn][0 + 9][mxstrk] = (dpnew[lastd + 9][mx][mn][0 + 9][mxstrk] + val) % MOD;
                            if(lastd > -9 && streak > -9){
                                dpnew[lastd - 1 + 9][mx][max(mn, -lastd + 1)][min(-1, streak - 1) + 9][max(mxstrk, abs(min(-1, streak - 1)))] = (dpnew[lastd - 1 + 9][mx][max(mn, -lastd + 1)][min(-1, streak - 1) + 9][max(mxstrk, abs(min(-1, streak - 1)))] + val) % MOD;
                            }
                        }
                    }
                }
            }
        }
        swap(dp, dpnew);
    }
    
    ll val[10][10][10] = {};
    for(int i = 0; i < 19; i++){
        for(int mx = 0; mx < 10; mx++){
            for(int mn = 0; mn < 10; mn++){
                for(int cstrk = 0; cstrk < 19; cstrk++){
                    for(int mxstk = 0; mxstk < 10; mxstk++){
                        val[mx][mn][mxstk] = (val[mx][mn][mxstk] + dp[i][mx][mn][cstrk][mxstk]) % MOD;
                    }
                }
            }
        }
    }
 
    ll ans = 0;
    for(int i = 0; i < 10; i++){
        for(int mx1 = 0; mx1 < 10; mx1++){
            for(int mx2 = 0; mx2 < 10; mx2++){
                for(int mn1 = 0; mn1 < 10; mn1++){
                    for(int mn2 = 0; mn2 < 10; mn2++){
                        if(((i + mx1 + mx2) < 10) && ((i - mn1 - mn2) >= 0)){
                            for(int mxstr = 1; mxstr <= 10; mxstr++){
                                for(int mxstr2 = 1; mxstr2 <= min(c / mxstr, 10LL); mxstr2++){
                                    ans = (ans + val[mx1][mn1][mxstr - 1] * val[mx2][mn2][mxstr2 - 1] % MOD) % MOD;
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    cout << ans;
}

This task can be solved in many ways. We will discuss a simple approach here.

The key idea is to find the first $$$4$$$ entries of the permutation by making $$$ {4}\choose{3} $$$ queries as follows :

• query 1 — $$$2$$$ $$$3$$$ $$$4$$$, let the result of this query is $$$q_1$$$.

• query 2 — $$$1$$$ $$$3$$$ $$$4$$$, let the result of this query is $$$q_2$$$.

• query 3 — $$$1$$$ $$$2$$$ $$$4$$$, let the result of this query is $$$q_3$$$.

• query 4 — $$$1$$$ $$$2$$$ $$$3$$$, let the result of this query is $$$q_4$$$.

Now, it is easy to see that $$$p_i$$$ $$$=$$$ $$$ ( \frac{q_1 + q_2 + q_3 + q_4}{3} ) - q_i $$$ for $$$ 1 \le i \le 4 $$$.

After this, we can find $$$p_i$$$ $$$(5 \le i \le n)$$$ by querying $$$1$$$ $$$2$$$ $$$i$$$ as we already know $$$p_1$$$ and $$$p_2$$$.

#include<bits/stdc++.h>
using namespace std;
int query(int i,int j,int k){
    int x;
    cout << "? " << i << ' ' << j << ' ' << k << flush << endl;
    cin >> x;
    return x;
}
int main(){
    int t;
    cin >> t;

    while(t--){
        int n;
        cin >> n;
        vector<int> p(n+1);

        int q[5];
        q[1] = query(2,3,4);
        q[2] = query(1,3,4);
        q[3] = query(1,2,4);
        q[4] = query(1,2,3);

        for(int i=1;i<=4;i++) p[i] = (q[1] + q[2] + q[3] + q[4])/3 - q[i];

        for(int i=5;i<=n;i++) p[i] = query(1,2,i) - p[1] - p[2];

        cout << "! ";
        for(int i=1;i<=n;i++) cout << p[i] << ' ';
        cout << flush << endl;

    }
    return 0;
}

A natural way to look at the string is to decompose it into contiguous groups (or blocks) of identical characters. Suppose the string can be represented as:

where each $$$c_i$$$ is either $$$0$$$ or $$$1$$$, and $$$c_i \ne c_{i+1}$$$. The score of the original string is:

since each block contributes $$$(\text{length} - 1)$$$ to the score.

The idea to improve the score is to create longer consecutive sequences by “merging” blocks of the same character by removing some characters. The main observations are as follows:

1. We will never remove a block partially because it never increases the beauty. So, we only remove a block completely.
2. We only remove a block if $$$n_i=1$$$. In this case, the answer increases by $$$1$$$ by merging left and right blocks. If $$$n_i \gt 1$$$, the answer either stays the same or decreases upon removal.

Time Complexity: $$$O(n)$$$

#include<bits/stdc++.h>
using namespace std;
 
typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef set<int> si;
typedef pair<int, int> pi;
typedef pair<ll, ll> pl;
 
#define ordered_set tree<ll, null_type,less<ll>, rb_tree_tag,tree_order_statistics_node_update> 
#define fo(i, a, b) for (ll i = a; i < b; i++)
#define F first
#define S second
#define pb push_back
#define mp make_pair
#define no cout << "No" << "\n"
#define yes cout << "Yes" << "\n"
#define cr(a,b) a.clear(); a.resize(b);
#define rep(i, a, b) for(int i = a; i < (b); ++i)
#define all(x) begin(x), end(x)
#define sz(x) (int)(x).size()
typedef long long ll;
typedef pair<int, int> pii;
typedef vector<int> vi;
typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef set<int> si;
typedef pair<int, int> pi;
typedef pair<ll, ll> pl;
ll p = 1e9+7;
#define mod p

int main() {
    cin.tie(0)->sync_with_stdio(0);
    cin.exceptions(cin.failbit);
    ll T,n;
    cin >> T;
    while (T--) {
        cin >>n;
        string s;
        cin >> s;
        ll ans=0;
        fo(i,1,n){
            if(s[i]!=s[i-1]){
                if(i<n-1){
                    if(s[i]!=s[i+1]){
                        s[i]=s[i-1];
                    }
                }
            }
            else{
                ans++;
            }
        }
        cout <<ans << "\n";
    }
}

Main Observation: It is optimal to convert the entire array to its minimum element $$$m$$$.

Let the minimum only occur once, then we can imagine that for any operation $$$l$$$ to $$$r$$$, each element except the minimum, contributes a cost of $$$1$$$. Thus, answer always stays same or increases if that the element is greater than or equal to min plus one. If there are multiple elements equal to min, then we can break down the operation into multiple operations.

For example, given array $$$[2,3,3,2,3]$$$, instead of performing $$$l=1$$$ and $$$r=5$$$. We break it down into $$$l=1$$$ to $$$r=3$$$ and $$$l=4$$$ to $$$r=5$$$. The final approach is to iterate on array and add min to answer if element is equal to min, else add min plus one.

Time Complexity: $$$O(n)$$$

#include<bits/stdc++.h>
using namespace std;
 
typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef set<int> si;
typedef pair<int, int> pi;
typedef pair<ll, ll> pl;
 
#define ordered_set tree<ll, null_type,less<ll>, rb_tree_tag,tree_order_statistics_node_update> 
#define fo(i, a, b) for (ll i = a; i < b; i++)
#define F first
#define S second
#define pb push_back
#define mp make_pair
#define no cout << "No" << "\n"
#define yes cout << "Yes" << "\n"
#define cr(a,b) a.clear(); a.resize(b);
#define rep(i, a, b) for(int i = a; i < (b); ++i)
#define all(x) begin(x), end(x)
#define sz(x) (int)(x).size()
typedef long long ll;
typedef pair<int, int> pii;
typedef vector<int> vi;
typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef set<int> si;
typedef pair<int, int> pi;
typedef pair<ll, ll> pl;
ll p = 1e9+7;
#define mod p
 
int main() {
    cin.tie(0)->sync_with_stdio(0);
    cin.exceptions(cin.failbit);
    ll T, n,m;
    cin >> T;
    while (T--) {
        cin >> n;
        vl a(n);
        ll e=1e9;
        fo(i,0,n){
            cin >> a[i];
            e=min(e,a[i]);
        }
        ll ans=0;
        fo(i,0,n){
            if(a[i]==e){
                ans+=e;
            }
            else{
                ans+=e+1;
            }
        }
        cout << ans << "\n";
    }
}

Let product of a subsequence $$${a_{i_1},a_{i_2},...,a_{i_k}}$$$ be defined as $$$\prod_{j=1}^{k} \frac{a_{i_j}}{\phi(a_{i_j})}$$$. $$$\sum_{g=1}^{\max a} \frac{1}{g} \cdot (\text{sum of products of all subsequences whose } \gcd \text{ is } g)$$$

We can write in terms of divisors $$$d \mid g$$$ , $$$h(g)=\frac{1}{g} = \sum_{d \mid g} f(d)$$$

$$$f$$$ can be calculated using the mobius inversion theorem $$$\newline$$$ $$$f(n)=\sum_{d \mid n} \mu(d) \cdot h(n/d)=\sum_{d \mid n} \frac{\mu(d)}{(n/d)}$$$

Finally, the original sum can be transformed into $$$\sum_{d=1}^{\max a} f(d) \cdot (\text{sum of products of all subsequences whose elements are multiples of } d)$$$

$$$= \gt \sum_{d=1}^{\max a} f(d) \cdot \left[\left(\sum_{\text{m:d|m }} \left[ \left(1+\frac{m}{\phi(m)}\right)^{\mathrm{freq}(m)} \right] \right)-1\right]$$$

Alternatively, a standard method involves iterating over the possible gcd values directly. For each divisor $$$d$$$, we compute the product

which represents the sum of contributions of all subsequences whose elements are multiples of $$$d$$$ (i.e., those with a gcd that is a multiple of $$$d$$$). By processing these values in descending order of $$$d$$$, we can subtract the contributions from all higher multiples of $$$d$$$ to isolate the contribution of subsequences whose gcd is exactly $$$d$$$.

Final Time Complexity: $$$O(nlogn)$$$

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define in(a, b) for (ll i = (a); i <= (b); i++)                // in using i
#define inj(a, b) for (ll j = (a); j <= (b); j++)               // in using j
#define ink(a, b) for (ll k = (a); k <= (b); k++)               // in using k
#define inl(a, b) for (ll l = (a); l <= (b); l++)               // in using l
#define inr(a, b) for(ll i = (a); i >= (b); i--)                // in reverse
#define inrj(a, b) for(ll j = (a); j >= (b); j--)                // in reverse
#define inrk(a, b) for(ll k = (a); k >= (b); k--)                // in reverse
#define inrl(a, b) for(ll l = (a); l >= (b); l--)                // in reverse
#define tt ll tcs; cin>>tcs; while(tcs--)                       // include test cases
#define ina(arr,n) ll arr[(n+1)]={0}; in(1,n) cin>>arr[i]       // input arr of n elements
#define inv(vec,n) vector<ll> vec(n+1); vec[0]=0; in(1,n) cin>>vec[i]; // input vector of n elements
#define pb push_back
#define lb lower_bound
#define ub upper_bound
#define pll pair <ll,ll>
#define vpll vector <pll>
#define sll set <ll>
#define vll vector<ll>
#define mll map <ll,ll>
#define all(x) x.begin(), x.end()
const ll mod=1e9+7;
#define vvll vector<vll>
#define pref(p,a,n) vll p(n+1); in(1,n) p[i]=p[i-1]+a[i];
#define vec2(a,n,m) vvll a(n+1,vll(m+1))
// #define vec2(a,n,m,val) vvll a(n+1,vll(m+1,val))

#define vec3(a,l,m,n) vector<vvll>a(l+1,vvll(m+1,vll(n+1)));
// #define vec3(a,l,m,n,val) vector<vvll>a(l+1,vvll(m+1,vll(n+1,val)));
# define FAST ios_base::sync_with_stdio(false); cin.tie(NULL); cout.tie(NULL);
//Dont forget to reset MOD if some other MOD is needed
const ll MOD = mod;
 
//Comment the line above and uncomment the line below if problem requires more than 1 MOD
//After uncommenting the below line, declaration of mint becomes [ mint<mod> M; ]
 
//template<ll MOD>
class mint
{
    //WARNING:
    //Be very careful not to use two mints with different mods for any operation
    //No guarantee of behavior in this case
    public:
    ll val;
    static ll mod_exp(ll a, ll b){ ll res=1;   a=a%MOD; while(b>0){ if(b%2==1) res=(res*a)%MOD; b/=2; a=(a*a)%MOD; } return res; }
    static ll gcdExtended(ll a, ll b, ll *x, ll *y) { if (a == 0) { *x = 0, *y = 1; return b; } ll x1, y1; ll gcd = gcdExtended(b%a, a, &x1, &y1);*x = y1 - (b/a) * x1; *y = x1; return gcd; }
    static ll modInverse(ll a) { ll x, y; 	ll g = gcdExtended(a, MOD, &x, &y); g++; ll res = (x%MOD);	if(res < 0) res += MOD;	return res;} 
    mint(){	val = 0;} 
    mint(ll x){	val = x%MOD;	if(val < 0) val += MOD;}
    mint& operator +=(const mint &other){	val += other.val;	if(val >= MOD) val -= MOD; return (*this); }
    mint& operator -=(const mint &other){   val -= other.val;if(val < 0) val += MOD;  return (*this); }
    mint& operator *=(const mint &other){	val = (val * other.val)%MOD; return (*this); }
    mint& operator /=(const mint &other){	val = (val * modInverse(other.val)) % MOD; return (*this); }
    mint& operator =(const mint &other) { 	val = other.val; return (*this); }
    mint operator +(const mint &other) const {	return mint(*this) += other; }
    mint operator -(const mint &other) const {	return mint(*this) -= other; }
    mint operator *(const mint &other) const {	return mint(*this) *= other; }
    mint operator /(const mint &other) const {	return mint(*this) /= other; }
    bool operator ==(const mint &other) const {   return val == other.val; }
 
    mint operator ++() { ++val; if(val == MOD) val = 0; return (*this); }
    mint operator ++(int) { val++; if(val == MOD) val = 0; return mint(val-1); }
    mint operator --() { --val; if(val == -1) val = MOD-1; return (*this); }
    mint operator --(int) { val--; if(val == -1) val = MOD-1; return mint(val+1); }
 
    // ^ has very low precedence, careful!!
 
    template<typename T>
    mint& operator ^=(const T &other){   val = mod_exp(val, other); return (*this); }
    template<typename T>
    mint operator ^(const T &other) const {  return mint(*this) ^= other; }
 
    mint& operator ^=(const mint &other){   val = mod_exp(val, other.val); return (*this); }
    mint operator ^(const mint &other) const {  return mint(*this) ^= other; }
 
    template<typename T>
    explicit operator T() {	return (T)val; }
    template<typename T>
    friend mint operator +(T other, const mint &M){	return mint(other) + M; }
    template<typename T>
    friend mint operator -(T other, const mint &M){	return mint(other) - M; }
    template<typename T>
    friend mint operator *(T other, const mint &M){	return mint(other) * M; }
    template<typename T>
    friend mint operator /(T other, const mint &M){	return mint(other) / M; }
    template<typename T>
    friend mint operator ^(T other, const mint &M){	return mint(other) ^ M; }
 
 
    friend std::ostream &operator << (std::ostream &output, const mint &M){  return output << M.val; }
    friend std::istream &operator >> (std::istream &input, mint &M) { input >> M.val;	M.val %= MOD;	return input;}
};
vll phi; //totient function of n
void totient_sieve(ll n){
    // totient function = no. of integers upto n co-prime to n
    phi.clear();
    phi.resize(n+1);
    for (ll i=0;i<=n;i++) phi[i]=i;
    for (ll i=2;i<=n;i++) {
        if (phi[i]==i){
            for (ll j=i;j<=n;j+=i) phi[j]=phi[j]-phi[j]/i;
        }
    }
}
 
vector<mint>cnt(1e6+1);
vector<mint>dp(1e6+1);
int main(){
    FAST
    totient_sieve(1e6);
    ll n;
    cin>>n;
    inv(a,n);
    in(1,n){
        cnt[a[i]]++;
    }
    for(int i=1;i<=1e6;i++){
        cnt[i]=(1+(mint(i)/mint(phi[i])))^cnt[i];
    }
    for(int d=1;d<=1e6;d++){
        for(int n=2*d;n<=1e6;n+=d){
            cnt[d]*=cnt[n];
        }
    }
    mint ans=0;
    for(int d=1e6;d>=1;d--){
        for(int n=2*d;n<=1e6;n+=d){
            cnt[d]-=cnt[n];
        }
        cnt[d]--;
        ans+=cnt[d]/d;
    }
    cout<<ans;
    
}

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define in(a, b) for (ll i = (a); i <= (b); i++)       
#define inj(a, b) for (ll j = (a); j <= (b); j++)              
#define ink(a, b) for (ll k = (a); k <= (b); k++)
#define vll vector<ll>
#define all(x) x.begin(), x.end()
const ll mod=1e9+7; 
# define FAST ios_base::sync_with_stdio(false); cin.tie(NULL); cout.tie(NULL);
 
//Dont forget to reset MOD if some other MOD is needed
const ll MOD = mod;
 
//Comment the line above and uncomment the line below if problem requires more than 1 MOD
//After uncommenting the below line, declaration of mint becomes [ mint<mod> M; ]
 
//template<ll MOD>
class mint
{
    //WARNING:
    //Be very careful not to use two mints with different mods for any operation
    //No guarantee of behavior in this case
    public:
    ll val;
    static ll mod_exp(ll a, ll b){ ll res=1;   a=a%MOD; while(b>0){ if(b%2==1) res=(res*a)%MOD; b/=2; a=(a*a)%MOD; } return res; }
    static ll gcdExtended(ll a, ll b, ll *x, ll *y) { if (a == 0) { *x = 0, *y = 1; return b; } ll x1, y1; ll gcd = gcdExtended(b%a, a, &x1, &y1);*x = y1 - (b/a) * x1; *y = x1; return gcd; }
    static ll modInverse(ll a) { ll x, y; 	ll g = gcdExtended(a, MOD, &x, &y); g++; ll res = (x%MOD);	if(res < 0) res += MOD;	return res;} 
    mint(){	val = 0;} 
    mint(ll x){	val = x%MOD;	if(val < 0) val += MOD;}
    mint& operator +=(const mint &other){	val += other.val;	if(val >= MOD) val -= MOD; return (*this); }
    mint& operator -=(const mint &other){   val -= other.val;if(val < 0) val += MOD;  return (*this); }
    mint& operator *=(const mint &other){	val = (val * other.val)%MOD; return (*this); }
    mint& operator /=(const mint &other){	val = (val * modInverse(other.val)) % MOD; return (*this); }
    mint& operator =(const mint &other) { 	val = other.val; return (*this); }
    mint operator +(const mint &other) const {	return mint(*this) += other; }
    mint operator -(const mint &other) const {	return mint(*this) -= other; }
    mint operator *(const mint &other) const {	return mint(*this) *= other; }
    mint operator /(const mint &other) const {	return mint(*this) /= other; }
    bool operator ==(const mint &other) const {   return val == other.val; }
 
    mint operator ++() { ++val; if(val == MOD) val = 0; return (*this); }
    mint operator ++(int) { val++; if(val == MOD) val = 0; return mint(val-1); }
    mint operator --() { --val; if(val == -1) val = MOD-1; return (*this); }
    mint operator --(int) { val--; if(val == -1) val = MOD-1; return mint(val+1); }
 
    // ^ has very low precedence, careful!!
 
    template<typename T>
    mint& operator ^=(const T &other){   val = mod_exp(val, other); return (*this); }
    template<typename T>
    mint operator ^(const T &other) const {  return mint(*this) ^= other; }
 
    mint& operator ^=(const mint &other){   val = mod_exp(val, other.val); return (*this); }
    mint operator ^(const mint &other) const {  return mint(*this) ^= other; }
 
    template<typename T>
    explicit operator T() {	return (T)val; }
    template<typename T>
    friend mint operator +(T other, const mint &M){	return mint(other) + M; }
    template<typename T>
    friend mint operator -(T other, const mint &M){	return mint(other) - M; }
    template<typename T>
    friend mint operator *(T other, const mint &M){	return mint(other) * M; }
    template<typename T>
    friend mint operator /(T other, const mint &M){	return mint(other) / M; }
    template<typename T>
    friend mint operator ^(T other, const mint &M){	return mint(other) ^ M; }
 
 
    friend std::ostream &operator << (std::ostream &output, const mint &M){  return output << M.val; }
    friend std::istream &operator >> (std::istream &input, mint &M) { input >> M.val;	M.val %= MOD;	return input;}
};
 
 
 
ll modInverse(ll A, ll M)
{
    ll m0 = M;
    ll y = 0, x = 1;
    if (M == 1)
        return 0;
    while (A > 1) {
        int q = A / M;
        int t = M;
        M = A % M, A = t;
        t = y;
        y = x - q * y;
        x = t;
    }
    if (x < 0)
        x += m0;
    return x;
}
 
 
ll binpow(ll a,ll b,ll m = LLONG_MAX){
    ll res = 1;
    while (b){
        if (b&1) res = (res*a)%m;
        a = (a*a)%m;
        b>>=1;
    }
    return res;
}
 
 
 
vll mob;
vll pr;
vll phi; //totient function of n
 
const int mxn=1e6;
vll cnt(mxn+1);
vll inv(mxn+1);
vector<mint> f(mxn+1);
vector<mint> val(mxn+1);
void pre_compute(ll n){
    mob.clear();
    mob.assign(n+1,0);
    pr.assign(n+1,0);
    for(ll p=2;p<=n;p++){
        if(!pr[p]){// p is prime
            pr[p]=p;
            for(ll j=(ll)p*p;j<=n;j+=p){
                if(!pr[j]){//no smaller factor
                    pr[j]=p;
                }
            }
        }
    }
    for(ll d=1;d<=n;d++){
        if(d==1){
            mob[d]=1;
        }
        else{
            if(pr[d/pr[d]]==pr[d]){
                mob[d]=0;
            }
            else{
                mob[d]=(-1)*mob[d/pr[d]];//agar 0 hoga(pichhle prime pe square) toh bhi direct 0 aajayega
            }
        }
    }
    // totient function = no. of integers upto n co-prime to n
    phi.clear();
    phi.resize(n+1);
    for (ll i=0;i<=n;i++) phi[i]=i;
    for (ll i=2;i<=n;i++) {
        if (phi[i]==i){
            for (ll j=i;j<=n;j+=i) phi[j]=phi[j]-phi[j]/i;
        }
    }
    //h(n)=1/n === sigma(d divides n) f(n) using mobius inversion f[n]=sigma mod[d]*(1/(n/d))
    for(int d=1;d<=mxn;d++){
        for(int n=d;n<=mxn;n+=d){
            f[n]+=(mob[d] * d);
        }
    }
    for(int n=1;n<=mxn;n++){
        f[n]/=n;
    }
 

}
 
 
int main(){
    FAST
    pre_compute(mxn+10);
 
    int n;cin>>n;
    vll a(n+1);
    in(1,n) {
        cin>>a[i];
        cnt[a[i]]++;
    }
    for(int i=1;i<=mxn;i++){
        val[i]=(1+i/mint(phi[i]))^cnt[i];
    }
    mint sum=0;
    for(int d=1;d<=mxn;d++){
        mint prod=1;
        for(int j=d;j<=mxn;j+=d){
            prod*=val[j];
        }
        sum+=(prod-1)*f[d];
    }
    cout<<sum;
}

For any element $$$x$$$, Let $$$s$$$ be the set of elements whose bits are a subset of bits of $$$x$$$, any elements outside of $$$s$$$ will have atleast one extra bit not in $$$x$$$ and hence cant be used in the set making $$$x$$$. Among all subsets $$$y$$$ of $$$s$$$, $$$OR(y)$$$ will be a subset of $$$x$$$. Hence subtracting the answers for all subsets of $$$x$$$ will give the answer for $$$x$$$. Here we only consider subsets of size $$$k$$$. After calculating answer for number of sets whose $$$OR$$$ is a subset $$$x$$$. The answer is just running sum over subsets DP in reverse order to get the answer for every $$$x$$$.$$$\newline$$$ The final steps are:-

1.Construct the $$$dp[x]$$$ array equal to count of elements $$$=x$$$.

2.Run Sum over Subsets to get $$$dp[x]=$$$ count of all subsets of $$$x$$$.

3.for every element $$$x$$$ get $$$\binom{dp[x]}{k}$$$ to be the number of ways to pick $$$k$$$ numbers whose or is a subset of $$$x$$$.

4.in the current state $$$dp[x] =$$$ sum of number of ways over all subsets of $$$x$$$. Running SOS Dp in reverse order cancel out the number of ways that are subsets of $$$x$$$ but not exactly $$$x$$$.

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define in(a, b) for (ll i = (a); i <= (b); i++)                // in using i
#define inj(a, b) for (ll j = (a); j <= (b); j++)               // in using j
#define ink(a, b) for (ll k = (a); k <= (b); k++)               // in using k
#define inl(a, b) for (ll l = (a); l <= (b); l++)               // in using l
#define inr(a, b) for(ll i = (a); i >= (b); i--)                // in reverse
#define inrj(a, b) for(ll j = (a); j >= (b); j--)                // in reverse
#define inrk(a, b) for(ll k = (a); k >= (b); k--)                // in reverse
#define inrl(a, b) for(ll l = (a); l >= (b); l--)                // in reverse
#define tt ll tcs; cin>>tcs; while(tcs--)                       // include test cases
#define ina(arr,n) ll arr[(n+1)]={0}; in(1,n) cin>>arr[i]       // input arr of n elements
#define inv(vec,n) vector<ll> vec(n+1); vec[0]=0; in(1,n) cin>>vec[i]; // input vector of n elements
#define pb push_back
#define lb lower_bound
#define ub upper_bound
#define pll pair <ll,ll>
#define vpll vector <pll>
#define sll set <ll>
const ll mod=1e9+7;
 
const ll MOD = mod;
 
//Comment the line above and uncomment the line below if problem requires more than 1 MOD
//After uncommenting the below line, declaration of mint becomes [ mint<mod> M; ]
 
//template<ll MOD>
class mint
{
    //WARNING:
    //Be very careful not to use two mints with different mods for any operation
    //No guarantee of behavior in this case
    public:
    ll val;
    static ll mod_exp(ll a, ll b){ ll res=1;   a=a%MOD; while(b>0){ if(b%2==1) res=(res*a)%MOD; b/=2; a=(a*a)%MOD; } return res; }
    static ll gcdExtended(ll a, ll b, ll *x, ll *y) { if (a == 0) { *x = 0, *y = 1; return b; } ll x1, y1; ll gcd = gcdExtended(b%a, a, &x1, &y1);*x = y1 - (b/a) * x1; *y = x1; return gcd; }
    static ll modInverse(ll a) { ll x, y; 	ll g = gcdExtended(a, MOD, &x, &y); g++; ll res = (x%MOD);	if(res < 0) res += MOD;	return res;} 
    mint(){	val = 0;} 
    mint(ll x){	val = x%MOD;	if(val < 0) val += MOD;}
    mint& operator +=(const mint &other){	val += other.val;	if(val >= MOD) val -= MOD; return (*this); }
    mint& operator -=(const mint &other){   val -= other.val;if(val < 0) val += MOD;  return (*this); }
    mint& operator *=(const mint &other){	val = (val * other.val)%MOD; return (*this); }
    mint& operator /=(const mint &other){	val = (val * modInverse(other.val)) % MOD; return (*this); }
    mint& operator =(const mint &other) { 	val = other.val; return (*this); }
    mint operator +(const mint &other) const {	return mint(*this) += other; }
    mint operator -(const mint &other) const {	return mint(*this) -= other; }
    mint operator *(const mint &other) const {	return mint(*this) *= other; }
    mint operator /(const mint &other) const {	return mint(*this) /= other; }
    bool operator ==(const mint &other) const {   return val == other.val; }
 
    mint operator ++() { ++val; if(val == MOD) val = 0; return (*this); }
    mint operator ++(int) { val++; if(val == MOD) val = 0; return mint(val-1); }
    mint operator --() { --val; if(val == -1) val = MOD-1; return (*this); }
    mint operator --(int) { val--; if(val == -1) val = MOD-1; return mint(val+1); }
 
    // ^ has very low precedence, careful!!
 
    template<typename T>
    mint& operator ^=(const T &other){   val = mod_exp(val, other); return (*this); }
    template<typename T>
    mint operator ^(const T &other) const {  return mint(*this) ^= other; }
 
    mint& operator ^=(const mint &other){   val = mod_exp(val, other.val); return (*this); }
    mint operator ^(const mint &other) const {  return mint(*this) ^= other; }
 
    template<typename T>
    explicit operator T() {	return (T)val; }
    template<typename T>
    friend mint operator +(T other, const mint &M){	return mint(other) + M; }
    template<typename T>
    friend mint operator -(T other, const mint &M){	return mint(other) - M; }
    template<typename T>
    friend mint operator *(T other, const mint &M){	return mint(other) * M; }
    template<typename T>
    friend mint operator /(T other, const mint &M){	return mint(other) / M; }
    template<typename T>
    friend mint operator ^(T other, const mint &M){	return mint(other) ^ M; }
 
 
    friend std::ostream &operator << (std::ostream &output, const mint &M){  return output << M.val; }
    friend std::istream &operator >> (std::istream &input, mint &M) { input >> M.val;	M.val %= MOD;	return input;}
};
 
#define vll vector<ll>
#define mll map <ll,ll>
#define all(x) x.begin(), x.end()
#define vvll vector<vll>
#define pref(p,a,n) vll p(n+1); in(1,n) p[i]=p[i-1]+a[i];
#define vec2(a,n,m) vvll a(n+1,vll(m+1))
// #define vec2(a,n,m,val) vvll a(n+1,vll(m+1,val))
 
#define vec3(a,l,m,n) vector<vvll>a(l+1,vvll(m+1,vll(n+1)));
// #define vec3(a,l,m,n,val) vector<vvll>a(l+1,vvll(m+1,vll(n+1,val)));
# define FAST ios_base::sync_with_stdio(false); cin.tie(NULL); cout.tie(NULL);
//Dont forget to reset MOD if some other MOD is needed
vector<mint>f(1e6+1);
vector<mint>inf(1e6+1);
void pre(){
    f[0]=1;
    for(int i=1;i<=1e6;i++){
        f[i]=f[i-1]*i;
    }
    inf[1e6]=1/f[1e6];
    for(int i=1e6-1;i>=0;i--){
        inf[i]=(inf[i+1]*(i+1));
    }
}
mint ncr(int n, int r){
    if(r>n) return 0;
    return f[n]*inf[r]*inf[n-r];
}

int main(){
    FAST
    pre();
    ll n,k;
    cin>>n>>k;
    vector<int>cnt(1<<20);
    vector<ll> a(n+1);
    in(1,n){
        cin>>a[i];
        cnt[a[i]]++;
    }
 
    for(int lg=0;lg<20;lg++){
        for(int i=1;i<(1<<20);i++){
            if((i&(1<<lg))){
                cnt[i]+=cnt[i^(1<<lg)];
            }
        }
    }
 
    vector<mint>dp(1<<20);
    for(int j=1;j<(1<<20);j++){
        dp[j]=ncr(cnt[j],k);
    }
    for(int lg=0;lg<20;lg++){
        for(int j=1;j<(1<<20);j++){
            if((j&(1<<lg))){
                dp[j]-=dp[j^(1<<lg)];
            }
        }
    }
    

    int m;
    cin>>m;
    while(m--){
        int x;
        cin>>x;
        cout<<dp[x]<<endl;
    }
}

Task 1

This construction is built on three observations:

1. The XUM value of a permutation $$$p$$$ can be expressed as $$$\sum_{i=0}^{n} \sum_{j=i+1}^{n} a_i \oplus a_j$$$, where each $$$a_i$$$ is the prefix XOR of the first $$$i$$$ elements of the permutation.

2. For any fixed bit (say the $$$i$$$th bit) that appears $$$c_i$$$ times among the numbers $$$0$$$ to $$$n-1$$$, it must appear at least $$$\lceil \tfrac{c_i}{2} \rceil$$$ times in the prefix XOR array. If we can make a construction that achieves this for each bit, then it will be optimal.

3. The crucial observation for the construction is that for any four consecutive numbers of the form $$$4k, 4k+1, 4k+2, 4k+3,$$$ their XOR is zero. We can group numbers into blocks of $$$4$$$, such that the xor of this block cancels out for future numbers.

Using these observations and some case work on how to arrange an individual 4 size block, the final construction is as follows:

• Case 1: $$$n \bmod 4 \neq 3$$$
  Arrange the blocks so that within each block the order is:

This ordering minimizes the occurrence of each bit.

• Case 2: $$$n \bmod 4 = 3$$$
  Here the final block consists of only three terms. In this situation, arranging the block as:

(with the appropriate adjustment for the three-term block) helps to avoid one excess occurrence of the 1st bit in the prefix XOR of the last $$$3$$$ terms, keeping it within the bound of $$$\lceil \tfrac{c_1}{2} \rceil$$$.

Task 2

Using the first and second observation in task $$$1$$$, the answer for the $$$i$$$-th bit is the number of elements in prefix xor having the $$$i$$$-th bit set multiplied by the number of bits having it unset, you can write out this sum as follows:

Here, $$$c_i= \frac{n}{2^{i+1}} \times 2^i + max(0,n \mod 2^{i+1} -2^i)$$$

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const ll mod = 1e9 + 7;

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);
    int t;
    cin >> t;
    while (t--) {
        ll c, n;
        cin >> c >> n;
        if (c == 2) 
        {
            ll sum = 0;
            ll d=1,f;
            for(ll i=0;d<=n-1;i++){
                f=(n/(2*d))*d+max(ll(0),n%(2*d)-d);
                f=(f+1)/2;
                sum+=(((f*(n+1-f))%mod)*d)%mod;
                sum%=mod;
                d*=2;
            }
        cout << sum << '\n';
        } else {
            vector<ll> v;
            if(n%4!=3)
                v={0,1,3,2};
            else
                v={1,0,2,3};
            for(int i=0;i<n;i++){
                cout<<(i/4)*4+v[i%4]<<" ";
            }
            cout << "\n";
        }
    }
    return 0;
}

Main observation: The operation never increases the gcd of the whole array because $$$gcd(a,b)=gcd(a,b-a)$$$.

Thus, we only need to bruteforce, removing each element. To do this, save the prefix and suffix gcd and take max over $$$gcd(pref_{i-1},suf_{i+1})$$$.

#include<bits/stdc++.h>
using namespace std;
 
typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef set<int> si;
typedef pair<int, int> pi;
typedef pair<ll, ll> pl;
 
#define ordered_set tree<ll, null_type,less<ll>, rb_tree_tag,tree_order_statistics_node_update> 
#define fo(i, a, b) for (ll i = a; i < b; i++)
#define F first
#define S second
#define pb push_back
#define mp make_pair
#define no cout << "No" << "\n"
#define yes cout << "Yes" << "\n"
#define cr(a,b) a.clear(); a.resize(b);
#define rep(i, a, b) for(int i = a; i < (b); ++i)
#define all(x) begin(x), end(x)
#define sz(x) (int)(x).size()
typedef long long ll;
typedef pair<int, int> pii;
typedef vector<int> vi;
typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef set<int> si;
typedef pair<int, int> pi;
typedef pair<ll, ll> pl;
ll p = 1e9+7;
#define mod p

int main() {
    cin.tie(0)->sync_with_stdio(0);
    cin.exceptions(cin.failbit);
    ll T, n,m;
    cin >> T;
    while (T--) {
        cin >> n;
        vl a(n);
        fo(i,0,n){
            cin >> a[i];
        }
        vl b(n);
        b[n-1]=a[n-1];
        for(int i=n-2;i>=0;i--){
            b[i]=gcd(b[i+1],a[i]);
        }
        ll ans=0;
        ll f=0;
        fo(i,0,n-1){
            ans=max(ans,gcd(f,b[i+1]));
            f=gcd(f,a[i]);
        }
        ans=max(ans,f);
        cout << ans << "\n";
    }
}

Let $$$val[v]$$$ be the amount of land $$$v$$$ received from it's parent, if all nodes only distributed their lands equally among their children. Let $$$child[v]$$$ be the number of children of node $$$v$$$, if $$$v$$$ receives exactly $$$val[v]$$$ land, it needs to pass on $$$\lfloor \frac{val[v]}{child[v]} \rfloor$$$ land to all it's children and $$$k=val[v] \bmod child[v]$$$ children will receive an additional acre of land.

Note that even if $$$v$$$ receives an additional acre of land, it's children doesn't receive more than $$$\lfloor \frac{val[v]}{child[v]} \rfloor+1$$$ land, since, $$$k+1 \leq child[v]$$$.

Finally, each child can receive either $$$val[v]$$$ or $$$val[v]+1$$$ land, giving us two dp states $$$dp[v][0]$$$ if $$$v$$$ receives exactly $$$val[v]$$$ and $$$dp[v][1]$$$ if it receives $$$val[v]+1$$$. Finally, $$$dp[v][0]$$$ can be trivially calculated as $$$\sum_{i} dp[c_i][0]$$$ plus the $$$k$$$ maximum values of $$$dp[c_i][1]-dp[c_i][0]$$$, where $$$c_i$$$ is the $$$i$$$th children of $$$v$$$. For $$$dp[v][1]$$$, we only need to add one more additional value to this.

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define in(a, b) for (ll i = (a); i <= (b); i++)                // in using i
#define inj(a, b) for (ll j = (a); j <= (b); j++)               // in using j
#define ink(a, b) for (ll k = (a); k <= (b); k++)               // in using k
#define inl(a, b) for (ll l = (a); l <= (b); l++)               // in using l
#define inr(a, b) for(ll i = (a); i >= (b); i--)                // in reverse
#define inrj(a, b) for(ll j = (a); j >= (b); j--)                // in reverse
#define inrk(a, b) for(ll k = (a); k >= (b); k--)                // in reverse
#define inrl(a, b) for(ll l = (a); l >= (b); l--)                // in reverse
#define tt ll tcs; cin>>tcs; while(tcs--)                       // include test cases
#define ina(arr,n) ll arr[(n+1)]={0}; in(1,n) cin>>arr[i]       // input arr of n elements
#define inv(vec,n) vector<ll> vec(n+1); vec[0]=0; in(1,n) cin>>vec[i]; // input vector of n elements
#define pb push_back
#define lb lower_bound
#define ub upper_bound
#define pll pair <ll,ll>
#define vpll vector <pll>
#define sll set <ll>
#define vll vector<ll>
#define mll map <ll,ll>
#define all(x) x.begin(), x.end()
const ll mod=1e9+7;
#define vvll vector<vll>
#define pref(p,a,n) vll p(n+1); in(1,n) p[i]=p[i-1]+a[i];
#define vec2(a,n,m) vvll a(n+1,vll(m+1))
// #define vec2(a,n,m,val) vvll a(n+1,vll(m+1,val))

#define vec3(a,l,m,n) vector<vvll>a(l+1,vvll(m+1,vll(n+1)));
// #define vec3(a,l,m,n,val) vector<vvll>a(l+1,vvll(m+1,vll(n+1,val)));
# define FAST ios_base::sync_with_stdio(false); cin.tie(NULL); cout.tie(NULL);
int main(){
    FAST
    tt{
        ll n;
        cin>>n;
        ll k;
        cin>>k;
        vvll adj(n+1);
        in(2,n){
            ll x;cin>>x;
            adj[x].push_back(i);
        }
        vll val(n+1);
        val[1]=k;
        vec2(dp,n,2);
        function<void(int)> dfs=[&](int v){
            if(adj[v].size()==0){
                dp[v][0]=val[v]^v;
                dp[v][1]=(val[v]+1)^v;
                return;
            }
            ll sbk=val[v]/adj[v].size();
            ll xtr=val[v]%adj[v].size();
            vll mor;
            for(auto x:adj[v]){
                val[x]=sbk;
                dfs(x);
                dp[v][0]+=dp[x][0]; 
                mor.push_back(dp[x][1]-dp[x][0]);
            }
            sort(all(mor),greater<ll>());
            for(int i=0;i<xtr;i++){
                dp[v][0]+=mor[i];
            }
            dp[v][1]=dp[v][0]+mor[xtr];
            dp[v][0]+=v^(val[v]);
            dp[v][1]+=v^(val[v]+1);
        };
        dfs(1);
        cout<<dp[1][0]<<endl;

    }
}

Probalistic Solution

To determine whether nodes $$$u$$$ and $$$v$$$ are compatible, we check the following conditions:

• First Condition → $$$v$$$ must be in the subtree of $$$u$$$.
• Second Condition → The difference between the subtree size of $$$u$$$ and $$$v$$$ must be even.
• Third Condition → The maximum frequency of any value in the remaining nodes must be at most half of the remaining nodes.

How to Check the Third Condition Efficiently:

• Perform DFS to compute the Euler tour, entry, and exit time of each node.
• Maintain a 2D array, where for each unique value, we store the indices of nodes having that value in the Euler tour.
• Randomly pick nodes from the remaining set and check their frequency. To find the frequency of a value in the remaining nodes, we count its occurrences in $$$u$$$ ’s subtree using its entry and exit times, then exclude occurrences in $$$v$$$ ’s subtree. The difference gives the frequency. If the frequency is more than half of the remaining nodes, pairing is impossible.
• If we repeat the previous step approximately $$$30$$$ times, then if there is a node value with a frequency greater than half of the remaining nodes, the probability that it is not picked is at most $$$2^{-30}$$$, which is extremely small.

Time Complexity: $$$O(n + q*30*log(n))$$$

Deterministic Solution

Flatten out the tree using a Euler tour, such that each subarray forms a continuous subarray. For pair of nodes $$$(u,v)$$$, the nodes form two disjoint subarrays $$$(tin[u],tin[v]-1)$$$ and $$$(tout[v]+1,tout[u])$$$ in the flattened array.

Main Observation: For the most frequent element to constitute at atleast half of the remaining nodes, it must be one of the maximum frequency elements from either the first or second subarray. This is because if its frequency is less than half in both segments, it cannot exceed half overall.

You can find the max frequent element and its frequency using a segment tree and query for other segment using an ordered set.

Time Complexity: $$$O((n + q)logn)$$$

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

mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
int rand(int l, int r) {
    uniform_int_distribution<int> uid(l, r);
    return uid(rng);
}

int sz = 1e5 + 10;
vector<vector<int>> adj;
vector<int> val;
vector<int> euler;
vector<int> tin, tout, cnt;
vector<vector<int>> v1(sz);

void dfs(int v, int p = -1) {
    cnt[v] = 1;
    tin[v] = euler.size();
    euler.push_back(v);
    for (int u : adj[v]) {
        if (u != p) {
            dfs(u, v);
            cnt[v] += cnt[u];
        }
    }
    tout[v] = euler.size() - 1;
}

int getCnt(int x, int l, int r) {
    int lb = lower_bound(begin(v1[x]), end(v1[x]), l) - begin(v1[x]);
    int rb = upper_bound(begin(v1[x]), end(v1[x]), r) - begin(v1[x]);
    return rb - lb;
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);

    int n;
    cin >> n;
    adj.assign(n + 1, {});
    val.assign(n + 1, 0);
    tin.assign(n + 1, 0);
    tout.assign(n + 1, 0);
    cnt.assign(n + 1, 0);

    for (int i = 1; i <= n; i++) {
        cin >> val[i];
    }

    for (int i = 1; i <= n - 1; i++) {
        int u, v;
        cin >> u >> v;
        adj[u].push_back(v);
        adj[v].push_back(u);
    }

    dfs(1);

    for (int i = 0; i < euler.size(); i++) {
        int x = euler[i];
        v1[val[x]].push_back(i);
    }

    int q;
    cin >> q;
    while (q--) {
        int u, v;
        cin >> u >> v;

        int L1 = tin[u];
        int R1 = tout[u];

        int L2 = tin[v];
        int R2 = tout[v];

        if (!((L1 <= L2) && (R2 <= R1))) {
            cout << "NO" << '\n';
            continue;
        }
        int sz = cnt[u] - cnt[v];
        if (sz % 2 != 0) {
            cout << "NO" << '\n';
            continue;
        }

        int trial = 30;
        bool flag = 1;
        while (trial--) {
            int i = rand(L1, L2);
            int j = rand(R2, R1);
            int x = val[euler[i]];
            if (rand(0, 1) == 1) {
                x = val[euler[j]];
            }

            int f = getCnt(x, L1, R1) - getCnt(x, L2, R2);
            if (f > (sz / 2)) {
                flag = 0;
                break;
            }
        }
        if (flag) {
            cout << "YES" << '\n';
        } else {
            cout << "NO" << '\n';
        }
    }
    return 0;
}

// In the name of Allah.
// We are nothing and you are everything.
// Ya Ali!

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

const int MAX_N = 1e5 + 14, LG = 17;

int n, a[MAX_N], st[MAX_N], ft[MAX_N], sp[LG][MAX_N];
vector<int> g[MAX_N], appear[MAX_N];

void dfs(int v = 0, int p = -1) {
    static int time = 0;
    sp[0][time] = a[v];
    appear[a[v]].push_back(time);
    st[v] = time++;
    for (int u : g[v])
        if (u != p)
            dfs(u, v);
    ft[v] = time;
}

int count(int x, int l, int r) {
    return ranges::lower_bound(appear[x], r) - ranges::lower_bound(appear[x], l);
}

int get(int l, int r, int c1, int c2) {
    return get<1>(max(tuple{count(c1, l, r), c1}, tuple{count(c2, l, r), c2}));
}

int king(int l, int r) {
    int ans = 0;
    int p = l;
    for (int i = 0; i < LG; ++i)
        if (r - p >> i & 1) {
            ans = get(l, r, ans, sp[i][p]);
            p += 1 << i;
        }
    return ans;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cin >> n;
    for (int i = 0; i < n; ++i)
        cin >> a[i];
    for (int i = 0; i < n - 1; ++i) {
        int v, u;
        cin >> v >> u;
        --v;
        --u;
        g[v].push_back(u);
        g[u].push_back(v);
    }
    dfs();
    for (int k = 1; k < LG; ++k)
        for (int i = 0; i + (1 << k) <= n; ++i)
            sp[k][i] = get(i, i + (1 << k), sp[k - 1][i], sp[k - 1][i + (1 << k - 1)]);
    int q;
    cin >> q;
    while (q--) {
        int u, v;
        cin >> u >> v;
        --u;
        --v;
        int sz = st[v] - st[u] + ft[u] - ft[v];
        if (!(st[u] <= st[v] && ft[v] <= ft[u]) || sz % 2) {
            cout << "NO\n";
            continue;
        }
        int c1 = king(st[u], st[v]), c2 = king(ft[v], ft[u]);
        int max_count = max(count(c1, st[u], st[v]) + count(c1, ft[v], ft[u]),
                            count(c2, st[u], st[v]) + count(c2, ft[v], ft[u]));
        // cerr << max_count << endl;
        cout << (max_count * 2 > sz ? "NO\n" : "YES\n");
    }
}

The first step is to realize the solution to a well known dp problem for finding $$$f(A)$$$ for a given array. We define $$$dp[i][0]$$$ to be the maximum sum for elements $$$A_1 ... A_i$$$ where we do not pick element $$$i$$$ and can pick element $$$i+1$$$. In $$$dp[i][1]$$$, we may or may not pick $$$i^{\text{th}}$$$ element. This ensures that $$$dp[i][1] \gt =dp[i][0]$$$ $$$\newline$$$ The transition will be such $$$dp[i][0]=dp[i-1][1]$$$ and $$$dp[i][1]=max(dp[i][0],dp[i-1][0]+a[i])$$$ where $$$dp[n][1]$$$ will be the answer for entire array. $$$\newline$$$ We can maintain $$$dp[i][1]-dp[i][0]$$$ as one of the states of DP of the original problem where only on basis of this difference, future elements would be picked. Out of the excess $$$dp[i][1]-dp[i][0]$$$, we take into account the future dp transitions, which is guaranteed to be within the range of just $$$1$$$ additional element $$$A_i$$$ and hence the difference is within the range of $$$A_i \lt =2*10^5$$$ $$$\newline$$$ So, for the original problem we only maintain the $$$dp[i][1]-dp[i][0]$$$ value and add the sum $$$dp[i][0]$$$ for all possible remaining sets($$$2^{n-i}$$$) to the answer and forgetting $$$dp[i][0]$$$ completely $$$\newline$$$ The final DP states are $$$ansdp[i][j]$$$ representing the number of subsets of the first i numbers having a difference of $$$j$$$ between $$$dp[i][1]$$$ and $$$dp[i][0]$$$, note that it is not necessary that the subsequences have the $$$i^{th}$$$ element in them. $$$\newline$$$ To transition from $$$i$$$ to $$$i+1$$$ for all subsequences which is not taking $$$A_{i+1}$$$ we can directly update for them $$$ansdp[i+1][j]+=ansdp[i][j]$$$ later as they will have the same difference now considering all subsets having $$$A_{i+1}$$$ in them

Since, we have counted the contribution of $$$dp[i][0]$$$ already, the above expressions for $$$dp[i][1]-dp[i][0]=j$$$ remain

Here we add $$$dp[i+1][0]=j$$$ to the answer for all $$$2^{n-(i+1)}$$$ remaining ways to complete the set and the difference for the next state is

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

const ll MOD = 1000000007;
class mint
{
    public:
    ll val;
    static ll mod_exp(ll a, ll b){ ll res=1;   a=a%MOD; while(b>0){ if(b%2==1) res=(res*a)%MOD; b/=2; a=(a*a)%MOD; } return res; }
    static ll gcdExtended(ll a, ll b, ll *x, ll *y) { if (a == 0) { *x = 0, *y = 1; return b; } ll x1, y1; ll gcd = gcdExtended(b%a, a, &x1, &y1);*x = y1 - (b/a) * x1; *y = x1; return gcd; }
    static ll modInverse(ll a) { ll x, y; 	ll g = gcdExtended(a, MOD, &x, &y); g++; ll res = (x%MOD);	if(res < 0) res += MOD;	return res;} 
    mint(){	val = 0;} 
    mint(ll x){	val = x%MOD;	if(val < 0) val += MOD;}
    mint& operator +=(const mint &other){	val += other.val;	if(val >= MOD) val -= MOD; return (*this); }
    mint& operator -=(const mint &other){   val -= other.val;if(val < 0) val += MOD;  return (*this); }
    mint& operator *=(const mint &other){	val = (val * other.val)%MOD; return (*this); }
    mint& operator /=(const mint &other){	val = (val * modInverse(other.val)) % MOD; return (*this); }
    mint& operator =(const mint &other) { 	val = other.val; return (*this); }
    mint operator +(const mint &other) const {	return mint(*this) += other; }
    mint operator -(const mint &other) const {	return mint(*this) -= other; }
    mint operator *(const mint &other) const {	return mint(*this) *= other; }
    mint operator /(const mint &other) const {	return mint(*this) /= other; }
    bool operator ==(const mint &other) const {   return val == other.val; }
    mint operator ++() { ++val; if(val == MOD) val = 0; return (*this); }
    mint operator ++(int) { val++; if(val == MOD) val = 0; return mint(val-1); }
    mint operator --() { --val; if(val == -1) val = MOD-1; return (*this); }
    mint operator --(int) { val--; if(val == -1) val = MOD-1; return mint(val+1); }
    template<typename T>
    mint& operator ^=(const T &other){   val = mod_exp(val, other); return (*this); }
    template<typename T>
    mint operator ^(const T &other) const {  return mint(*this) ^= other; }
    mint& operator ^=(const mint &other){   val = mod_exp(val, other.val); return (*this); }
    mint operator ^(const mint &other) const {  return mint(*this) ^= other; }
    template<typename T>
    explicit operator T() {	return (T)val; }
    template<typename T>
    friend mint operator +(T other, const mint &M){	return mint(other) + M; }
    template<typename T>
    friend mint operator -(T other, const mint &M){	return mint(other) - M; }
    template<typename T>
    friend mint operator *(T other, const mint &M){	return mint(other) * M; }
    template<typename T>
    friend mint operator /(T other, const mint &M){	return mint(other) / M; }
    template<typename T>
    friend mint operator ^(T other, const mint &M){	return mint(other) ^ M; }
    friend std::ostream &operator << (std::ostream &output, const mint &M){  return output << M.val; }
    friend std::istream &operator >> (std::istream &input, mint &M) { input >> M.val;	M.val %= MOD;	return input;}
};
int mxa=2e5;
int main(){
    ios_base::sync_with_stdio(false); cin.tie(NULL); cout.tie(NULL);
    int t;
    cin>>t;
    vector<mint>pow2(1001);
    pow2[0]=1;
    for(int i=1;i<=1000;i++){
        pow2[i]=pow2[i-1]*2;
    }
    while(t--){
        int n;cin>>n;
        vector<int> a(n+1);
        for(int i=1;i<=n;i++){
            cin>>a[i];
        }
        vector<mint>ansdp(mxa+1);
        vector<mint>nwdp(mxa+1);
        ansdp[0]=1;
        mint ans=0;
        for(int i=1;i<=n;i++){
            for(int j=0;j<=mxa;j++){
                nwdp[j]=0;
            }
            for(int j=0;j<=mxa;j++){
                nwdp[max(0,a[i]-j)]+=ansdp[j];
                ans+=ansdp[j]*(pow2[n-i])*j;
                nwdp[j]+=ansdp[j];
            }
            swap(ansdp,nwdp);
        }
        for(int j=0;j<=mxa;j++){
            ans+=ansdp[j]*j;
        }
        cout<<ans<<endl;
    }
}

Case 1: $$$k \geq 3$$$

If we choose three or more vertices, notice that every pair of chosen vertices will yield a unique simple path. For the MEX on any path to be greater than $$$1$$$, every such path must contain both $$$0$$$ and $$$1$$$. However, in a tree the intersection of all these paths (considering every pair among $$$k$$$ vertices) can only be a single vertex.

Thus, for $$$k \geq 3$$$, the answer will be $$$1$$$ till $$$k$$$ is equal to the degree of node $$$0$$$ plus one, else it is $$$0$$$.

Case 2: $$$k = 2$$$

The problem becomes simpler when $$$k = 2$$$. In this case, we are just looking for a pair of vertices such that the unique path between them has the maximum possible MEX. A simple approach for this is to perform a binary search on answer i.e. final mex value. When querying for the first $$$k$$$ nodes, find the node with highest dfs in-time. If a path exists, this node is guaranteed to be an endpoint. We can run a dfs from the node with the highest dfs in-time and check for path.

Time Complexity: $$$O(nlogn)$$$

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

typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef set<int> si;
typedef pair<int, int> pi;
typedef pair<ll, ll> pl;

#define ordered_set tree<ll, null_type,less<ll>, rb_tree_tag,tree_order_statistics_node_update> 
#define fo(i, a, b) for (ll i = a; i < b; i++)
#define F first
#define S second
#define pb push_back
#define mp make_pair
#define no cout << "No" << "\n"
#define yes cout << "Yes" << "\n"
#define cr(a,b) a.clear(); a.resize(b);
#define rep(i, a, b) for(int i = a; i < (b); ++i)
#define all(x) begin(x), end(x)
#define sz(x) (int)(x).size()
typedef long long ll;
typedef pair<int, int> pii;
typedef vector<int> vi;
typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef set<int> si;
typedef pair<int, int> pi;
typedef pair<ll, ll> pl;
ll p = 1e9+7;
#define mod p

vector<vi> edges;
vector<bool> vis;
int timer=0;
vi intime;
vi node;
ll n;
int ma=0;
void dfs(int k){
    vis[k]=1;
    intime[k]=timer++;
    for(auto x:edges[k]){
        if(!vis[x]){
            dfs(x);
        }
    }
}
void dfs1(int k,int f,int t){
    vis[k]=1;
    if(k<=t){
        f++;
    }
    ma=max(ma,f);
    for(auto x:edges[k]){
        if(!vis[x]){
            dfs1(x,f,t);
        }
    }
}
bool f(int k){
    vis.clear();
    vis.resize(n);
    ma=0;
    dfs1(node[k],0,k);
    if(ma==k+1)return 1;
    else return 0;
}
int main() {
    cin.tie(0)->sync_with_stdio(0);
    cin.exceptions(cin.failbit);
    int T,a,b,c;
    cin >> T;
    while (T--) {
        cin >> n;
        c=1;
        edges.clear();
        edges.resize(n);
        fo(i,0,n-1){
            cin >>a >> b;
            edges[a].pb(b);
            edges[b].pb(a);
            if(min(a,b)==0){
                c++;
            }
        }
        vis.clear();
        vis.resize(n);
        intime.clear();
        intime.resize(n);
        timer=0;
        dfs(0);
        node.clear();
        node.resize(n);
        int r=intime[0];
        fo(i,1,n){
            node[i]=node[i-1];
            if(r<intime[i]){
                r=intime[i];
                node[i]=i;
            }
        }
        int start=0,end=n-1,mid;
        while(start<end){
            mid=(start+end+1)/2;
            if(f(mid)){
                start=mid;
            }
            else end=mid-1;
        }
        cout << start+1 << " ";
        fo(i,3,c+1){
            cout << 1 << " ";
        }
        fo(i,c+1,n+1){
            cout << 0 << " ";
        }
        cout << "\n";
    }
}