We can count the number of pairs $$$(i, j)$$$ and simply divide it with $$$n * m$$$ to get the corresponding probability. For counting the pairs we can follow the given steps:

Make a depth first search and maintain an array of hashes, which will store the hashes of the string which is in our current depth first search path.

For example, suppose our current depth first search path is $$$a_1$$$, $$$a_2$$$, $$$a_3$$$, where $$$a_1$$$ would be the root. Parent of $$$a_2$$$ is $$$a_1$$$ and parent of $$$a_3$$$ is $$$a_2$$$.For this path, our hash array length will be 3, where $$$hash_1$$$ is the hash of substring $$$a_1$$$, $$$hash_2$$$ is the hash of the substring $$$a_1 + a_2$$$ and $$$hash_3$$$ is the hash of the substring $$$a_1 + a_2 + a_3$$$. For this node we can binary search on this current hash array and a hash array of the given string to find the value of $$$upNode_{a_3}$$$. $$$upNode_x$$$ basically denotes the node up the current node till which a matching occurs with the given string.

Now, it is possible that two nodes on different branches and same depth produces the same value of $$$upNode$$$. We can remove this overcounting by using small to large merging.

#include<bits/stdc++.h>
#define ll                   long long
#define vll                  vector<ll> 
#define vi                   vector<int>
#define vvll                 vector<vector<ll>>
#define vvi                  vector<vector<int>>
#define MOD                  (ll)(1000000000 + 7)
#define all(x)               begin(x), end(x)
#define sz(x)                (int)x.size()
#define pi                   pair<int, int>

using namespace std;

template<typename T> istream& operator>>(istream &cin, pair<T, T> &v) {cin>>v.first>>v.second; return cin;}
template<typename T> istream& operator>>(istream &cin, vector<T> &v) {for(T &e: v) cin>>e; return cin;}

ll po(ll x, ll y, ll p)
{
    ll res = 1; x = x % p;   
    while (y > 0) { if (y & 1) res = (res * x) % p; y >>= 1; x = (x * x) % p; } 
    return res%p; 
}

void test_case();

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

    int test_cases = 1;
    cin>>test_cases;

    while(test_cases--) {
        test_case();
    }

    return 0;
}   

int n, s_si;
vvi g;
vi vals, s;

const vll mods = { 1027799999, 1018199999, 1000000007 };
const int P = 31;
int m = sz(mods);
vvll hashs, s_hash, pows;

vi siz, dep;
vi st, up_node;

void dfs(int par, vll cur_hash){
    
    for(int i = 0; i < m; i++){
        cur_hash[i] *= P;
        cur_hash[i] += vals[par];
        cur_hash[i]%= mods[i];
        hashs[par][i] = cur_hash[i];
    }

    st.push_back(par);
 
    if(vals[par] == s[0]){
        int l = max(0, (int)st.size() - s_si), r = (int)st.size() - 1, res = -1;
        while(l <= r){
            int mid = (l + r)/2;
            int ind_2 = (int)st.size() - 1 - mid;
            vll cur(m);
            for(int i = 0; i < m; i++){
                cur[i] = (hashs[par][i] - 
                (mid ? ((hashs[st[mid - 1]][i] * pows[st.size() - mid][i])%mods[i]) : 0) + mods[i])%mods[i]; 
            }

            assert(ind_2 >= 0 && ind_2 <= s_si - 1);
            if(cur == s_hash[ind_2]){
                res = mid;
                r = mid - 1;
            } else {
                l = mid + 1;
            }
        }
        up_node[par] = st[res];
    }

    for(int e: g[par]){
        dep[e] = dep[par] + 1;
        dfs(e, cur_hash);
        siz[par] += siz[e];
    }

    st.pop_back();
}


vector<set<pi>> v;

void dfs_get(int par, set<pi> &cur){
    if(up_node[par] != -1) cur.insert({up_node[par], dep[par]});
    for(int e: g[par]){
        dfs_get(e, cur);
    }
}

ll ans =0 ;

void dfs_sm_to_lg(int par = 0){
    int bg_ch = -1;
    for(int e: g[par]){
        if(bg_ch == -1 || siz[bg_ch] < siz[e]){
            bg_ch = e;
        }
        dfs_sm_to_lg(e);
    }
 
    if(bg_ch != -1){
        swap(v[par], v[bg_ch]);
    }
    if(up_node[par] != -1) v[par].insert({up_node[par], dep[par]});
    for(int e: g[par]){
        if(e != bg_ch){
            dfs_get(e, v[par]);
        }
    }
 
    while(!v[par].empty()){
        pi last = *v[par].rbegin();
        if(dep[last.first] <= dep[par]) break;
        v[par].erase(last);
    }

    ans += sz(v[par]);
    ans %= MOD;
}

void test_case(){
    cin>>n;

    assert(n >= 1 && n <= 5e5);

    g.assign(n, {});

    for(int i = 1; i < n; i++){
        int x;
        cin>>x;
        x--;
        g[x].push_back(i);
    }

    vals.assign(n, 0);

    for(int i = 0; i < n; i++){
        char ch;
        cin>>ch;

        vals[i] = ch - 'a' + 1;
    }

    cin>>s_si;
    
    s.assign(s_si, 0);
    for(int i = 0; i < s_si; i++){
        char ch;
        cin>>ch;

        s[i] = ch - 'a' + 1; 
    }

    hashs.assign(n, vll(m, 0));
    dep.assign(n, 0);
    siz.assign(n, 1);

    s_hash.assign(s_si, vll(m));
    pows.assign(s_si + 1, vll(m, 1));
 
    for(int i = 1; i <= s_si; i++){
        for(int j = 0; j < m; j++){
            pows[i][j] = (pows[i - 1][j] * P)%mods[j];
        }
    }
 
    vll pre(m);
    for(int i = 0; i < s_si; i++){
        for(int j = 0; j < m; j++){
            s_hash[i][j] = (pre[j] + ((s[i] * pows[i][j])%mods[j]))%mods[j];
        }
        pre = s_hash[i];
    }
 
    up_node.assign(n, -1);

    st.assign(0, {});
    dfs(0, vll(m));
 
    v.assign(n, {});

    ans = 0;

    dfs_sm_to_lg();

    ans *= po(n, MOD - 2, MOD);
    ans %= MOD;
    ans *= po(s_si, MOD - 2, MOD);
    ans %= MOD;
 
    cout<<ans<<"\n";
}

The main observation here is, For any number $$$N$$$ $$$\leq 10^{17}$$$, the smallest number that does not divide N is very small; let's assume this number is $$$X$$$. So $$$sweetness(N) = sweetness(X)+1.$$$

Now as we can observe $$$X$$$ is very small; we can easily simulate to calculate the answer for $$$X$$$.

If $$$X$$$ is the smallest number that does not divide $$$N$$$, it means every number less than $$$X$$$ will divide $$$N$$$. So $$$LCM(1,2,..., X-1)$$$ will divide $$$N$$$ and $$$LCM(1,2,...,X)$$$ will not divide $$$N$$$.

For range $$$[L, R]$$$, we have to find numbers that are divisible by $$$LCM(1,2,..., X-1)$$$ but not by $$$LCM(1,2,...,X)$$$.

To calculate $$$LCM(1,2,3,..,X)$$$ we can use property $$$LCM(LCM(1,2,3,..,X-1),X) = LCM(1,2,3..,X)$$$.

#include<bits/stdc++.h>
using namespace std;
#pragma GCC optimize("O3")
#pragma GCC target("avx2")
#define M 1000000007
#define ll long long int
#define pb push_back
#define pf push_front
#define all(x) x.begin(), x.end()
#define pob pop_back
#define pof pop_front
#define debug1(x) cout<<#x<<" "<<x<<endl;
#define debug2(x,y) cout<<#x<<" "<<x<<"   "<<#y<<" "<<y<<endl;
#define debug3(x,y,z) cout<<#x<<" "<<x<<"   "<<#y<<" "<<y<<"   "<<#z<<" "<<z<<endl;
#define debug4(x,y,z,p) cout<<#x<<" "<<x<<"   "<<#y<<" "<<y<<"   "<<#z<<" "<<z<<" "<<#p<<" "<<p<<endl;
#define vll      vector<ll> 
#define vvll     vector<vector<ll>>
#define fi first
#define se second
#define boost ios_base::sync_with_stdio(false);cin.tie(NULL);cout.tie(NULL);
#define inf 1e18
#define endl '\n'
typedef pair<int, int>		pii;
typedef pair<ll, ll>		pll;
#define pi (3.14159265358979323846264338327950288)
template<typename T>
void printv(const T& t) {std::copy(t.cbegin(), t.cend(), std::ostream_iterator<typename T::value_type>(std::cout, ", "));cout<<endl;}
ll modpower(ll a, ll b, ll mod) { ll res = 1; for (; b; b >>= 1) { if (b & 1)res = res * a; res %= mod; a = a * a; a %= mod; }return res; }
struct custom_hash {
    static uint64_t splitmix64(uint64_t x) {
        // http://xorshift.di.unimi.it/splitmix64.c
        x += 0x9e3779b97f4a7c15;
        x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9;
        x = (x ^ (x >> 27)) * 0x94d049bb133111eb;
        return x ^ (x >> 31);
    }

    size_t operator()(uint64_t x) const {
        static const uint64_t FIXED_RANDOM = chrono::steady_clock::now().time_since_epoch().count();
        return splitmix64(x + FIXED_RANDOM);
    }
};
//unordered_map<ll,ll, custom_hash>m1;
ll expo(ll a, ll b, ll mod) {ll res = 1; while (b > 0) {if (b & 1)res = (res * a) % mod; a = (a * a) % mod; b = b >> 1;} return res;}
ll mod_inv(ll a, ll b) {return expo(a, b - 2, b);}
ll mod_add(ll a, ll b, ll m) {a = a % m; b = b % m; return (((a + b) % m) + m) % m;}
ll mod_mul(ll a, ll b, ll m) {a = a % m; b = b % m; return (((a * b) % m) + m) % m;}
ll mod_sub(ll a, ll b, ll m) {a = a % m; b = b % m; return (((a - b) % m) + m) % m;}
ll mod_div(ll a, ll b, ll m) {a = a % m; b = b % m; return (mod_mul(a, mod_inv(b, m), m) + m) % m;}
//-------------------------------Template--Above-----------------------------------------------
vll lc;
vll val;
ll gcd(ll a, ll b)
{
  if (b == 0)
    return a;
  return gcd(b, a % b);
}
 
ll lcm(ll a, ll b)
{
    return (a / gcd(a, b)) * b;
}

void solve(){
	lc.assign(60,1);
	
	val.assign(60,1);
	
	val[2]=1;
	val[1]=1;
	for(ll i=3;i<=59;i++){
		ll x=i;
		ll curr=0;
		while(x>2){
			ll j=2;
			while(j<x and (x%j==0)){
				j++;
			}
			curr++;
			x=j;
		}
		val[i]=curr+1;
	}
	
	for(ll i=2;i<=59;i++){
		lc[i]=lcm(lc[i-1],i);
	}
	
	ll l,r;cin>>l>>r;
	
	assert(3<=l<=r and r<=100000000000000000);
	
	ll ans=0;
	for(ll i=2;i<=59;i++){
		ll prev=(r/lc[i-1])-((l-1)/lc[i-1]);
		ll curr=(r/lc[i])-((l-1)/lc[i]);
		
		ll x=((prev-curr)*val[i])+(prev-curr);
		
		if(prev>curr)
		ans+=((prev-curr)*val[i])+(prev-curr);
		
	}
	
	cout<<ans<<endl;
	
    return;
}
int main()
{
    boost
    /*
    freopen("moocast.in", "r", stdin);
	freopen("moocast.out", "w", stdout);
	*/
    ll t;
    t=1;
   // cin>>t;
    for(ll i=1;i<=t;i++){
    	//cout << "Case #" << i << ": "; solve();
    	solve();
	}
    return 0;
}

Probabilty that Alice gets to a pile of $$$1$$$ stone strictly before Bob can be written as

#include "bits/stdc++.h"
using namespace std;
    
#pragma GCC target ("avx2")
#pragma GCC optimize ("O3")
#pragma GCC optimize ("unroll-loops")
    
#define ll                    long long
#define ull                   unsigned long long
#define loop(i,a,b)           for(int i=(int)a;i<(int)b;++i)
#define loopr(i,a,b)          for(int i=(int)a;i>=(int)b;--i)
#define count_1(n)            __builtin_popcountll(n)
#define pb                    push_back
#define vi                    vector<int>
#define vvi                   vector<vector<int>>
const ll MOD = 998244353; //1000000007 998244353
const ll INF = 1e17  + 10; //careful
    
inline ll po(ll a, ll b) { ll res = 1; for (; b; b >>= 1) { if (b & 1)res = res * a; a = a * a; }return res; }
inline ll gcd (ll a, ll b) { ll r; while(b) { r = a%b; a=b; b=r; } return a; }
    
template<typename T> inline void debug(T a) { cout << a << endl; }
template<typename T, typename... Args> inline void debug(T a, Args... args) { cout << a << " "; debug(args...); }
template<typename T> inline void printlist(T& a) { for(auto i : a) cout << i << " "; cout << endl; }
 
inline ll modpow(ll a, ll b, ll mod) { ll res = 1; for (; b; b >>= 1) { if (b & 1)res = (res * a)%mod; a = (a * a)%mod; }return res; }
 
template <ll mod = MOD>
struct modint {
    ll val;
    modint() : val(0) {}
    modint(ll a) : val(a) { val += mod; if(val >= mod) val -= mod; }
    void operator = (modint a) { val = a.val; }
    bool operator == (modint a) { return val == a.val; } 
    bool operator != (modint a) { return val != a.val; }
    inline modint& operator += (modint a) { val = val + a.val; if(val >= mod) val-=mod; return *this; }
    inline modint& operator -= (modint a) { val = val - a.val + mod; if(val >= mod) val -= mod; return *this; }
    inline modint& operator *= (modint a) { val = (val*a.val)%mod; return *this; }
    friend modint operator + (modint a, modint b) { return (a += b); }
    friend modint operator - (modint a, modint b) { return (a -= b); } 
    friend modint operator * (modint a, modint b) { return (a *= b); }
    friend modint operator ^ (modint a, ll b) { modint res = 1; for (; b; b >>= 1) { if (b & 1)res*=a; a*=a; } return res; }
    friend modint operator / (modint num, modint den) { den = den ^ (mod-2); return (num*den); }
    friend ostream & operator << (ostream &out, modint a) { out << a.val; return out; }
    friend istream & operator >> (istream &in, modint& a) { in >> a.val; return in; }
};
 
using mint = modint<MOD>;
 
int lim = 500000;
vector<vector<mint>> dp(lim+1);
vector<vector<int>> primes(lim+1);
void preprocess() {
    loop(i,2,lim+1) {
        if(primes[i].size()) continue;
        for(int j = i; j <= lim; j+=i) primes[j].pb(i);
    }
    dp[1].pb(1);
    loop(i,1,lim+1) {
        mint den = 1;
        if(i > 1) {
            for(int j : primes[i]) { int res = 0; int ii = i; while(ii%j==0) { res++; ii/=j; } den *= res+1; }
            den = mint(1) / (den - 1);
            loop(j,0,dp[i].size()) dp[i][j] *= den;
        }
        for(int j = 2*i; j <= lim; j+=i) { if(dp[j].size() < dp[i].size() + 1) dp[j].resize(dp[i].size() + 1); loop(k,0,dp[i].size()) dp[j][k+1] += dp[i][k]; }
    }
}
 
void testcases(int test) {
    int n; cin>>n;
    mint ans = 0; loop(i,0,dp[n].size()) ans += dp[n][i] * dp[n][i];
    ans = (1 - ans) / mint(2);
    debug(ans);
}   
    
int32_t main()
{
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);
    
    
    auto start = chrono::high_resolution_clock::now();
    
    int t = 1;
    cin >> t;
    preprocess();
    loop(i,0,t) { testcases(i); }
    
    auto stop = chrono::high_resolution_clock::now();
    auto duration = chrono::duration_cast<chrono::milliseconds>(stop - start);
    //cout << duration.count() << " ms\n";
}

Here we can make pairs of drawing only if their colours are different. So if we have $$$K$$$ colours available, we will try to use all the colours So that we have as many different colour drawings available as possible.

It's also obvious (and easy to proof) divide the paintings equally among K colours to maximise the answer. If we can't divide drawings equally into K colours (if N is not divisible by k), We should divide the available drawing equally and then colour the remaining drawings with different colours.

#include<bits/stdc++.h>
using namespace std;
#pragma GCC optimize("O3")
#pragma GCC target("avx2")
#define M 1000000007
#define ll long long int
#define pb push_back
#define pf push_front
#define all(x) x.begin(), x.end()
#define pob pop_back
#define pof pop_front
#define debug1(x) cout<<#x<<" "<<x<<endl;
#define debug2(x,y) cout<<#x<<" "<<x<<"   "<<#y<<" "<<y<<endl;
#define debug3(x,y,z) cout<<#x<<" "<<x<<"   "<<#y<<" "<<y<<"   "<<#z<<" "<<z<<endl;
#define debug4(x,y,z,p) cout<<#x<<" "<<x<<"   "<<#y<<" "<<y<<"   "<<#z<<" "<<z<<" "<<#p<<" "<<p<<endl;
#define vll      vector<ll> 
#define vvll     vector<vector<ll>>
#define fi first
#define se second
#define boost ios_base::sync_with_stdio(false);cin.tie(NULL);cout.tie(NULL);
#define inf 1e18
#define endl '\n'
typedef pair<int, int>		pii;
typedef pair<ll, ll>		pll;
#define pi (3.14159265358979323846264338327950288)
template<typename T>
void printv(const T& t) {std::copy(t.cbegin(), t.cend(), std::ostream_iterator<typename T::value_type>(std::cout, ", "));cout<<endl;}
ll modpower(ll a, ll b, ll mod) { ll res = 1; for (; b; b >>= 1) { if (b & 1)res = res * a; res %= mod; a = a * a; a %= mod; }return res; }
struct custom_hash {
    static uint64_t splitmix64(uint64_t x) {
        // http://xorshift.di.unimi.it/splitmix64.c
        x += 0x9e3779b97f4a7c15;
        x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9;
        x = (x ^ (x >> 27)) * 0x94d049bb133111eb;
        return x ^ (x >> 31);
    }

    size_t operator()(uint64_t x) const {
        static const uint64_t FIXED_RANDOM = chrono::steady_clock::now().time_since_epoch().count();
        return splitmix64(x + FIXED_RANDOM);
    }
};
//unordered_map<ll,ll, custom_hash>m1;
ll expo(ll a, ll b, ll mod) {ll res = 1; while (b > 0) {if (b & 1)res = (res * a) % mod; a = (a * a) % mod; b = b >> 1;} return res;}
ll mod_inv(ll a, ll b) {return expo(a, b - 2, b);}
ll mod_add(ll a, ll b, ll m) {a = a % m; b = b % m; return (((a + b) % m) + m) % m;}
ll mod_mul(ll a, ll b, ll m) {a = a % m; b = b % m; return (((a * b) % m) + m) % m;}
ll mod_sub(ll a, ll b, ll m) {a = a % m; b = b % m; return (((a - b) % m) + m) % m;}
ll mod_div(ll a, ll b, ll m) {a = a % m; b = b % m; return (mod_mul(a, mod_inv(b, m), m) + m) % m;}
//-------------------------------Template--Above-----------------------------------------------
void solve(){
	ll n,k;cin>>n>>k;
	assert(2<=k and k<=n);
	assert(n<=100000);
	vll v(k,0);
	ll curr=n;
	ll i=0;
	while(curr>0){
		v[i]++;
		curr--;
		i++;
		i%=k;
	}
	
	ll ans=0;
	for(ll i=0;i<k;i++){
		ll rem=n-v[i];
		ans+=(rem*v[i]);
	}
	
	ans/=2;
	cout<<ans<<endl;
	
    return;
}
int main()
{
    boost
    /*
    freopen("moocast.in", "r", stdin);
	freopen("moocast.out", "w", stdout);
	*/
    ll t;
    t=1;
   // cin>>t;
    for(ll i=1;i<=t;i++){
    	//cout << "Case #" << i << ": "; solve();
    	solve();
	}
    return 0;
}

We want to build an initial sequence from which we can add or subtract terms to get any value in $$$[1,W]$$$. We can also think of this as multiplying each item in the list with $$$1, 0$$$ or $$$-1$$$ to get values $$$[1,W]$$$. This is simply ternary counting, so the powers of 3 ($$${log}_3(2W+1)$$$ of them) are an optimal choice for our initial sequence.

#include "bits/stdc++.h"
using namespace std;
    
#pragma GCC target ("avx2")
#pragma GCC optimize ("O3")
#pragma GCC optimize ("unroll-loops")
    
#define ll                    long long
#define ull                   unsigned long long
#define loop(i,a,b)           for(int i=(int)a;i<(int)b;++i)
#define loopr(i,a,b)          for(int i=(int)a;i>=(int)b;--i)
#define count_1(n)            __builtin_popcountll(n)
#define pb                    push_back
#define vi                    vector<int>
#define vvi                   vector<vector<int>>
const ll MOD = 1000000007; //1000000007 998244353
const ll INF = 1e17  + 10; //careful
    
inline ll po(ll a, ll b) { ll res = 1; for (; b; b >>= 1) { if (b & 1)res = res * a; a = a * a; }return res; }
inline ll gcd (ll a, ll b) { ll r; while(b) { r = a%b; a=b; b=r; } return a; }
    
template<typename T> inline void debug(T a) { cout << a << endl; }
template<typename T, typename... Args> inline void debug(T a, Args... args) { cout << a << " "; debug(args...); }
template<typename T> inline void printlist(T& a) { for(auto i : a) cout << i << " "; cout << endl; }
 
vector<int> pot(18);
void preprocess() {
    pot[0] = 1;
    loop(i,1,18) pot[i] = 3*pot[i-1]; 
}
 
void testcases(int test) {
    int w; cin>>w; int W = 0;
    int s = 0; while(pot[s] < 2*w+1) { W += pot[s]; s++; }
    cout << s << "\n";
    loop(i,0,s) cout << pot[i] << " "; cout << "\n";
    loop(i,1,w+1) {
        int key = W + i;
        vector<int> vals(s);
        loop(j,0,s) { vals[j] = key%3; key/=3; }
        vector<int> s1, s2; loop(j,0,s) { if(vals[j] == 2) s2.pb(pot[j]); else if(vals[j] == 0) s1.pb(pot[j]); }
        cout << s1.size() << " ";
        printlist(s1);
        cout << s2.size() << " ";
        printlist(s2);
    }                
}   
    
int32_t main()
{
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);
    
    
    auto start = chrono::high_resolution_clock::now();
    
    int t = 1;
    cin >> t;
    preprocess();
    loop(i,0,t) { testcases(i); }
    
    auto stop = chrono::high_resolution_clock::now();
    auto duration = chrono::duration_cast<chrono::milliseconds>(stop - start);
    //cout << duration.count() << " ms\n";
}

For the subarray $$$a_l, a_{l + 1}, ... a_{r}$$$, the number of possible games will be:

$$$(r - l + 1)!/(A! * B!)$$$

where $$$A$$$ is the number of elements less than $$$x$$$ and $$$B$$$ is the number of elements greater than $$$x$$$.

For each query, we just need to find the value of $$$A$$$ and $$$B$$$. We can solve these queries offline and use segment tree to find there respective values.

#include<bits/stdc++.h>
#define ll                   long long
#define vll                  vector<ll> 
#define vi                   vector<int>
#define MOD                  (ll)(1000000000 + 7)
#define all(x)               begin(x), end(x)
#define sz(x)                (int)x.size()
#define pi                   pair<int, int>

using namespace std;

ll inline mo(ll a){ return a%MOD;}

ll po(ll x, ll y, ll p)
{
    ll res = 1; x = x % p;   
    while (y > 0) { if (y & 1) res = (res * x) % p; y >>= 1; x = (x * x) % p; } 
    return res%p; 
}

void test_case();

const int maxn = 500000 + 5;
vll fact;

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

    fact.assign(maxn, 0);
    fact[0] = 1;

    for(int i = 1; i < maxn; i++){
        fact[i] = (fact[i - 1] * i)%MOD;
    }

    int test_cases = 1;
    cin>>test_cases;

    while(test_cases--) {
        test_case();
    }

    return 0;
}   


struct segtree
{
    int size;
    vi sum;

    void init(int n){
        size = 1;
        while(size<=n) size *= 2;
        sum.assign(2*size, neutral());
    }

    /* Change this function for updating nodes */

    int update(int f, int s){
        return (f + s);
    }

    /* Return neutral element */
    int neutral(){
        return 0;
    }

    int get(int l, int r){
        return get(l, r, 0, 0, size);
    }

    int get(int l, int r, int x, int lx, int rx){
        if(lx>=r||rx<=l) return neutral();
        if(lx>=l&&rx<=r){
            return sum[x];
        }

        ll m = (lx + rx)/2;

        ll a = get(l, r, 2*x + 1, lx, m);
        ll b = get(l, r, 2*x + 2, m, rx);

        return update(a, b);
    }

    void set(int l, int r, int v){
        set(l, r, v, 0, 0, size);
    }

    void set(int l, int r, int v, int x, int lx, int rx){
        if( lx >= r || rx <= l) return;
        if( lx >= l && rx <= r){
            sum[x] += v;
            return;
        }

        ll m = (lx + rx)/2;

        set(l, r, v, 2*x + 1, lx, m);
        set(l, r, v, 2*x + 2, m, rx);

        sum[x] = update(sum[2*x + 1], sum[2*x + 2]);
    }
    
};

ll C(int n, int r){
    if(n < r){
        return 0;
    }

    return (fact[n] * ((po(fact[n - r], MOD - 2, MOD) * po(fact[r], MOD - 2, MOD))%MOD))%MOD;
}

void test_case(){
    int n;
    cin>>n;

    vll a(n);
    for(ll &e: a){
        cin>>e;
    }

    ll q;
    cin>>q;

    vector<pair<int, pair<pi, int>>> qu(q);

    for(int query = 0; query < q; query++){
        ll l, r, x;
        cin>>l>>r>>x;
        l--; r--;

        qu[query] = {x, {{l, r}, query}};
    }

    sort(all(qu));

    vector<pair<int, int>> vals(n);
    for(int i = 0; i < n; i++){
        vals[i] = {a[i], i};
    }

    sort(all(vals));
    int ind = 0;

    segtree sg;
    sg.init(n);

    vi ans(q);

    for(int i = 0; i < q;){
        int j = i;
        while(j < q && qu[j].first == qu[i].first){
            j++;
        }

        int x = qu[i].first;

        while(ind < n && vals[ind].first < x){
            sg.set(vals[ind].second, vals[ind].second + 1, 1);
            ind++;
        }

        while(i < j){
            int index = qu[i].second.second;
            int l = qu[i].second.first.first;
            int r = qu[i].second.first.second;

            int less = sg.get(l, r + 1);

            ans[index] = C(r - l + 1, less);
            i++;
        }
    }

    for(int e: ans){
        cout<<e<<"\n";
    }
}

We can first convert the problem into $$$1$$$-dimension, imagine the player moving left or right on the non-negative number line, which will be reflected by his $$$y$$$-coordinate in the actual board, and the passage of time will be maintained by his $$$x$$$-coordinate (they move $$$1$$$ row down each second). Now, the problem becomes: find the expected number of moves to get to position $$$z$$$ in a number line from position $$$y$$$, where each move you go back one step with probability $$$b$$$, or go forward one step with probability $$$a$$$; and at the end we just add this expected number of moves to the starting row number to get the answer.

We can use linearity of expectation here, as we can instead look at the expected number of steps it takes to get to position $$$x$$$ from position $$$x-1$$$. If we define this as $$$p_x$$$, our required expected value will be $$$p_{y+1} + p_{y+2} + .. + p_{z}$$$. We can also quiet easily find a recurrence for $$$p_x$$$:

Given that $$$p_0 = 0$$$, we can calculate the sum,

If you expand this out, and group terms favourably, it will break down into a sum of a gemotric sequence and an arithmetico-geomteric sequcence, thus this sum can be calculated in $$$O(log(n))$$$ time.

#include "bits/stdc++.h"
using namespace std;
    
#pragma GCC target ("avx2")
#pragma GCC optimize ("O3")
#pragma GCC optimize ("unroll-loops")
    
#define ll                    long long
#define ull                   unsigned long long
#define loop(i,a,b)           for(int i=(int)a;i<(int)b;++i)
#define loopr(i,a,b)          for(int i=(int)a;i>=(int)b;--i)
#define count_1(n)            __builtin_popcountll(n)
#define pb                    push_back
#define vi                    vector<int>
#define vvi                   vector<vector<int>>
const ll MOD = 998244353; //1000000007 998244353
const ll INF = 1e17  + 10; //careful
    
inline ll po(ll a, ll b) { ll res = 1; for (; b; b >>= 1) { if (b & 1)res = res * a; a = a * a; }return res; }
inline ll gcd (ll a, ll b) { ll r; while(b) { r = a%b; a=b; b=r; } return a; }
    
template<typename T> inline void debug(T a) { cout << a << endl; }
template<typename T, typename... Args> inline void debug(T a, Args... args) { cout << a << " "; debug(args...); }
template<typename T> inline void printlist(T& a) { for(auto i : a) cout << i << " "; cout << endl; }
 
inline ll modpow(ll a, ll b, ll mod) { ll res = 1; for (; b; b >>= 1) { if (b & 1)res = (res * a)%mod; a = (a * a)%mod; }return res; }
 
template <ll mod = MOD>
struct modint {
    ll val;
    modint() : val(0) {}
    modint(ll a) : val(a) { val += mod; if(val >= mod) val -= mod; }
    void operator = (modint a) { val = a.val; }
    bool operator == (modint a) { return val == a.val; } 
    bool operator != (modint a) { return val != a.val; }
    inline modint& operator += (modint a) { val = val + a.val; if(val >= mod) val-=mod; return *this; }
    inline modint& operator -= (modint a) { val = val - a.val + mod; if(val >= mod) val -= mod; return *this; }
    inline modint& operator *= (modint a) { val = (val*a.val)%mod; return *this; }
    friend modint operator + (modint a, modint b) { return (a += b); }
    friend modint operator - (modint a, modint b) { return (a -= b); } 
    friend modint operator * (modint a, modint b) { return (a *= b); }
    friend modint operator ^ (modint a, ll b) { if(b < 0) return ((modint(1) / a) ^ (-b)); modint res = 1; for (; b; b >>= 1) { if (b & 1)res*=a; a*=a; } return res; }
    friend modint operator / (modint num, modint den) { den = den ^ (mod-2); return (num*den); }
    friend ostream & operator << (ostream &out, modint a) { out << a.val; return out; }
    friend istream & operator >> (istream &in, modint& a) { in >> a.val; return in; }
};
 
using mint = modint<MOD>;
 
mint calc(ll n, mint a, mint b) {
    if(n == 0) return 0;
 
    if(a == b) {
        mint t1 = mint(n) * 2;
        mint t2 = (mint(n) * (n-1));
        return t1 + t2;
    }
 
    mint t1 = mint(n) / a;
    mint t2 = ((mint(n) * b) / ((a ^ n) * (b-a))) * ((b^(n-1)) - (a^(n-1)));
    mint t3 = b / (a * (a - b));
    mint t4 = ((b^2)) / (a * (a-b) * (a-b));
    mint t5 = (b^n) / ((a^(n-1)) * (a-b) * (a-b));
    mint t6 = ((b^n) * (n-1)) / ((a^n) * (a-b));
 
    mint ans = t1 + t2 - t3 - t4 + t5 + t6;
    return ans;
}
 
void testcases(int test) {
    ll x,y,z; cin>>x>>y>>z;
    mint a,b; cin>>a>>b; mint c = a+b; a = a/c; b = b/c;
    mint ans = calc(z,a,b) - calc(y,a,b);
    ans = ans + x;
    debug(ans);
}   
    
int32_t main()
{
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);
    
    
    auto start = chrono::high_resolution_clock::now();
    
    int t = 1;
    cin >> t;
    loop(i,0,t) { testcases(i); }
    
    auto stop = chrono::high_resolution_clock::now();
    auto duration = chrono::duration_cast<chrono::milliseconds>(stop - start);
    //cout << duration.count() << " ms\n";
}

Observation 1 : A diameter of a tree through a node passes through one of the ancestors of  . Observation 2 : when we join $$$2$$$ trees , the diameter of the combined tree can be any one of the individual diameters of any of the trees or the one passing through the newly added edge . If we know for each node , the maximum depth in its left and right subtree . then the diameter passing through this node can be easily calculated . But keeping $$$left[i]$$$ and $$$right[i]$$$ for each node would be too much memory, because there are order of $$$2^30$$$ nodes in a single tree .

Observation 3 : we do not require to keep the $$$right[i]$$$ and $$$left[i]$$$ for each and every node of the tree, but only for those whose $$$left[i]$$$ and $$$right[i]$$$ are changed from the initial default values .  After removing a subtree at node $$$u$$$, update all $$$left[i]$$$ and $$$right[i]$$$ for the nodes who are ancestors of $$$u$$$. Thus we did at most $$$Q*30$$$ operations. 

Now the main task is to know the diameter of the tree at any point of time . For each tree keep a multiset  for each level  denoting the $$$left[i] + right[i]$$$ values in that level whose $$$left[i]$$$ and $$$right[i]$$$ have been changed from initial default values. This gives the diameter passing through the node $$$i$$$ in its subtree. While finding diameter, go through each level and see if the multiset is full or not (having all enteries from this level). If its full, then take the max value from this multiset, else take the default value of diameter from this level. 

See the code for more details

#include             <bits/stdc++.h>
#include             <ext/pb_ds/assoc_container.hpp>
#include             <ext/pb_ds/tree_policy.hpp>
#define PRE(x,p)     cout<<setprecision(x)<<p; 
#define pb           push_back
#define mp           make_pair
#define f            first
#define s            second
#define pi           3.14159265358979
#define mod          (ll)(1e9 + 7)
#define endl         "\n"
#define high         1e18
#define low          -1e18
#define ll           long long int
#define ld           long double
#define mem(x,val)   memset(x,0,sizeof(x));
#define rep(i,l,r)   for(ll i=l;i<=r;i++)
#define p(a)         for(auto i:a) cout<<i<<' '; cout<<endl;
#define vll          vector<ll> 
#define vb           vector<bool>
#define vpll         vector<pair<ll,ll>>
#define vi           vector<int>
#define vpi          vector<pair<int, int>>
#define vvll         vector<vector<ll>>
#define vvi          vector<vector<int>>
#define vvvll        vector<vector<vector<ll>>>
#define pll          pair<ll,ll>
#define vs           vector<string>
#define vvpll        vector<vector<pair<ll, ll>>>
#define vvpi         vector<vector<pair<int, int>>>
#define vpii         vector<pair<int, int>>
#define sz(a)        (ll)a.size()
#define po(x)        (ll)(1ll<<x)
#define all(x)       begin(x), end(x)
#define speed        ios_base::sync_with_stdio(false);cin.tie(NULL);cout.tie(NULL);
#define yes          {cout<<"YES"<<endl;return;}
#define no           {cout<<"NO"<<endl; return;}
#define ok           cout<<"ok"<<endl;
#define ordered_set  tree<int, null_type,less<int>, rb_tree_tag,tree_order_statistics_node_update>


using namespace std;
using namespace __gnu_pbds;

void showa(ll a[],ll n)  { for(ll i=1;i<=n;i++) cout<<a[i]<<' ';  cout<<endl;  }
ll ison(ll w ,ll i)    {return  w&(1ll<<i);}
void amax(ll &a, ll b){ a=max(a,b); }
void amin(ll &a, ll b){ a=min(a,b);}
void modadd(ll &a , ll b) {a=((a%mod)+(b%mod))%mod;}
void modsub(ll &a , ll b) {a=((a%mod)-(b%mod)+mod)%mod;}
void modmul(ll &a , ll b) {a=((a%mod)*(b%mod))%mod;}

#ifndef ONLINE_JUDGE
#define debug(x) cerr << #x <<" "; _print(x); cerr << endl;
#else
#define debug(x)
#endif

void _print(ll t)     {cerr << t<<' ';}
void _print(int t)    {cerr << t<<' ';}
void _print(string t) {cerr << t<<' ';}
void _print(char t)   {cerr << t<<' ';}
void _print(ld t)     {cerr << t<<' ';}
void _print(double t) {cerr << t<<' ';}
template<class T,class V> void _print(pair <T, V> p);
template<class T>void _print(vector <T> v);
template<class T>void _print(vector <T> v);
template<class T>void _print(set <T> v);
template<class T,class V> void _print(map <T, V> v);
template<class T>void _print(multiset <T> v);
template<class T,class V> void _print(pair <T, V> p) {cerr << "{"; _print(p.f); cerr << ","; _print(p.s); cerr << "}";}
template<class T>void _print(vector <T> v) {cerr << "[ "; for (T i : v) {_print(i); cerr << " ";} cerr << "]";}
template<class T>void _print(set <T> v) {cerr << "[ "; for (T i : v) {_print(i); cerr << " ";} cerr << "]";}
template<class T>void _print(multiset <T> v) {cerr << "[ "; for (T i : v) {_print(i); cerr << " ";} cerr << "]";}
template<class T,class V> void _print(map <T, V> v) {cerr << "[ "; for (auto i : v) {_print(i); cerr << " ";} cerr << "]";}
void check(ll x , ll l , ll r){
  assert(x>=l && x<=r);
}
//const ll l=30;   //log2(n)
const ll N=200005;
ll h[N];  

unordered_map<ll,ll>removed[N];

unordered_map<ll,ll>lefty[N];
unordered_map<ll,ll>righty[N];

set<ll>right_left[N][31];  //for each tree and at each level
ll def_depth(ll x, ll y)
{
  // default depth from tree x's node y to a leaf node down
  ll node=y;
  ll ans=0;
  while(node>=1)
  {
    ans++; node/=2;
  }
  return h[x]-ans;
}

ll get_right(ll x, ll node)
{
  if(righty[x].count(node)) return righty[x][node];
  else return def_depth(x, node);
}
ll get_left(ll x, ll node)
{
  if(lefty[x].count(node)) return lefty[x][node];
  else return def_depth(x, node);
}
void update_right(ll x, ll node , ll par)
{
  
  ll val=max(get_right(x,node), get_left(x,node)) +1;
  righty[x][par]=val;
}
void update_left(ll x, ll node , ll par)
{
  ll val=max(get_right(x,node), get_left(x,node)) +1 ; 
  lefty[x][par]=val;
}

ll get_maxi(ll x, ll node)
{
  ll ans=0;
  ll base=1;
  ans=max(get_right(x,node),get_left(x,node));
  // debug(def_depth(2,5))
  // debug(node) debug(ans)
  while(node>1)
  { 
    ll par=node/2;
    if(par*2 == node) ans=max(ans, base + get_right(x,par));
    else  ans=max(ans, base + get_left(x,par));
    base++;
    node/=2;
  }
  return ans;
}

ll getdiam(ll x)
{
  ll height=h[x];
  ll ans=0;
  ll totnodes=1;
  ll default_right_left=2*(h[x]-1);
  rep(i,1,height)
  {
    ll here;
    if(right_left[x][i].size()==totnodes) here=*(--right_left[x][i].end());
    else here=default_right_left;  
    
    ans=max(ans,here);
      totnodes*=2;
      default_right_left-=2;
  }

  return ans;
}

void remove_right_left(ll x, ll node , ll level)
{
  ll val=get_left(x,node) + get_right(x,node);
  if(val== 2 * def_depth(x,node)) return; 
  right_left[x][level].erase(right_left[x][level].lower_bound(val));
}
void add_right_left(ll x, ll node , ll level )
{
  ll val=get_left(x,node) + get_right(x,node);
  right_left[x][level].insert(val);
}

void solve()
{
   ll n,q;
   assert(cin>>n>>q);
   check(n,1,2e5);
   check(q,1,2e5);


   rep(i,1,n) { assert(cin>>h[i]); check(h[i],1,30); }
     rep(i,1,n)
   {
     lefty[i].clear();
     righty[i].clear();
     rep(j,1,h[i]) right_left[i][j].clear();
   }
   rep(jj,1,q)
   {
     ll type;
     assert(cin>>type); check(type,1,2);
     if(type==1)
     {
      ll x,node;
      assert(cin>>x>>node);
      check(x,1,n);
      check(node,1, (1ll<<h[x]) -1 );
      
           ll k=node;
           while(k>=1)
           {
             assert(removed[x].count(k)==0);
             k/=2;
           }
          ll level= h[x]- def_depth(x,node);
      
         
      lefty[x][node]=-1;
      righty[x][node]=-1;
      removed[x][node]=1;
      

      while(node>1)
      {
        ll par=node/2;
        level--;
        remove_right_left(x,par,level);
        if(par*2 +1 == node) update_right(x,node,par);
        else update_left(x,node,par);
        add_right_left(x,par,level);
        node/=2;
      }

     }
     else
     {
       ll x,y,u,v;
       assert(cin>>x>>y>>u>>v);
       check(x,1,n); check(y,1,n); assert(x!=y);
       check(u,1, (1ll<<h[x]) -1);
       check(v,1, (1ll<<h[y]) -1);
     
       ll maxi1= get_maxi(x,u);        

       ll maxi2= get_maxi(y,v);
       ll ans1=maxi1+1+maxi2;
       cout<<max(ans1, max(getdiam(x) , getdiam(y)))<<endl;  

     }
   }
   
}
 
signed main()
{
  speed
  solve();
  return 0;
}  

It is clear from the problem statement that we need a mechanism to do the following stuff:

1. Keep track of the obstacles on the number line.

2. To check whether we can cover k contiguous points in the range, we need to find the maximum contiguous segment in the range given.

We can use set to track the obstacles on the number line. To find the maximum contiguous segment length, we can use a Data Structure like segment tree. For every obstacle placed on the number line, we can update the point with the distance between the current and the next obstacle and update the previous obstacle point with the distance between the current and previous obstacle. Similarly, we can update the segment lengths during the removal query.

Now we can simply query the maximum contiguous segment in the range l to r using the segment tree. Let ob1 and ob2 be the first and last obstacle in the range l to r then max contiguous segment = $$$max(ob1-l, r-ob2, query(ob1,ob2)).$$$

#include <bits/stdc++.h>
using namespace std;
#define ll int
const int req = 500005;
ll A[req + 5], tree[4 * req + 5];
void build(ll node, ll start, ll end)
{
    if (start == end)
    {
        tree[node] = A[start];
    }
    else
    {
        int mid = (start + end) / 2;
        build(2 * node, start, mid);
        build(2 * node + 1, mid + 1, end);
        tree[node] = max(tree[2 * node], tree[2 * node + 1]);
    }
}
void update(ll node, ll start, ll end, ll idx, ll val)
{
    if (start == end)
    {
        A[idx] = val;
        tree[node] = val;
    }
    else
    {
        ll mid = (start + end) / 2;
        if (start <= idx and idx <= mid)
        {
            update(2 * node, start, mid, idx, val);
        }
        else
        {
            update(2 * node + 1, mid + 1, end, idx, val);
        }
        tree[node] = max(tree[2 * node], tree[2 * node + 1]);
    }
}
int query(ll node, ll start, ll end, ll l, ll r)
{
    if (r < start or end < l)
    {
        return 0;
    }
    if (l <= start and end <= r)
    {
        return tree[node];
    }
    int mid = (start + end) / 2;
    int p1 = query(2 * node, start, mid, l, r);
    int p2 = query(2 * node + 1, mid + 1, end, l, r);
    return max(p1, p2);
}
int main()
{
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);
    cout.tie(NULL);
    ll t = 1;
    // cin >> t;
    for (ll o = 1; o <= t; o++)
    {
        ll n;
        ll i, j, k;
        ll q;
        cin  >> n >> q;
        map<ll, ll> obs;
        obs.insert({1, n});
        obs.insert({n + 3, 0});
        A[1] = n;
        build(1, 0, n);
        for (i = 0; i < q; i++)
        {
            cin >> j;
            if (j == 1)
            {
                cin >> k;
                auto itr = obs.lower_bound(k + 1);
                itr--;
                update(1, 0, n, (*itr).first, k - (*itr).first);
                int prev = (*itr).second;
                obs[(*itr).first] = k - (*itr).first;
                update(1, 0, n, k + 1, prev - (*itr).second - 1);
                obs.insert({k + 1, prev - (*itr).second - 1});
            }
            else if (j == 2)
            {
                cin >> k;
                auto itr = obs.find(k + 1);
                int val = query(1, 0, n, (*itr).first, (*itr).first);
                update(1, 0, n, (*itr).first, 0);
                itr--;
                update(1, 0, n, (*itr).first, (*itr).second + 1 + val);
                val += (*itr).second + 1;
                int ind = (*itr).first;
                obs.erase(k + 1);
                obs[ind] = val;
            }
            else
            {
                int l, r;
                cin >> l >> r >> k;
                int st, en;
                auto itr = obs.lower_bound(l + 1);
                int pos = (*itr).first - 1;
                st = pos - l;
                itr = obs.lower_bound(r + 2);
                itr--;
                int pos2 = (*itr).first - 1;
                en = r - pos2;
                int val;
                if (pos > r)
                {
                    val = r - l + 1;
                }
                else if (pos == pos2)
                {
                    val = max(pos - l, r - pos);
                }
                else
                {
                    val = query(1, 0, n, pos, pos2);
                    val = max({st, en, val});
                }
                if (val >= k)
                {
                    cout << "Possible"
                         << "\n";
                }
                else
                {
                    cout << "Not Possible"
                         << "\n";
                }
            }
        }
    }
}

Let us denote the parity of each node's depth in the array $$$parity$$$. Now, by solving the equation we can find that the solution for any subtree is equals to the xor of all the values whose parity is equal to the parity of the root node of that particular subtree, i.e.,

$$$ans$$$ = $$$val_{x_1}$$$ ^ $$$val_{x_2}$$$, ... ^ $$$val_{x_k}$$$, such that $$$parity_{x_i}$$$ = $$$parity_Y$$$, where $$$Y$$$ is the root of the corresponding subtree ans $$$x_i$$$ is present in the corresponding subtree.

Now, to solve the above problem we can first construct the entire tree and make an euler tour. Initially nodes which are not attached will have a value equals to 0 and whenever we add a node we can simply change its corresponding value. This is how we handle query of type 1.

To handle query of type 2, we can simply use segment tree on the array obtained by euler tour to answer range queries.

#include<bits/stdc++.h>
#define ll long long
#define vi                   vector<int>
#define vpi                  vector<pair<int, int>>
#define vvi                  vector<vector<int>>
#define vvpi                 vector<vector<pair<int, int>>>
#define MOD                  (ll)(1000000000 + 7)
#define all(x)               begin(x), end(x)
#define sz(x)                (int)x.size()
#define pi                   pair<int, int>

using namespace std;

template<typename T> istream& operator>>(istream &cin, pair<T, T> &v) {cin>>v.first>>v.second; return cin;}
template<typename T> istream& operator>>(istream &cin, vector<T> &v) {for(T &e: v) cin>>e; return cin;}

ll inline mo(ll a){ return a%MOD;}

ll po(ll x, ll y, ll p)
{
    ll res = 1; x = x % p;   
    while (y > 0) { if (y & 1) res = (res * x) % p; y >>= 1; x = (x * x) % p; } 
    return res%p; 
}

void test_case();

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

    int test_cases = 1;
    cin>>test_cases;

    while(test_cases--) {
        test_case();
    }

    return 0;
}   

struct segtree
{
    int size;
    // {even_depth_xor_val, odd_depth_xor_val}
    vpi xor_val;

    void init(ll n){
        size = 1;
        while(size<=n) size *= 2;
        xor_val.assign(2*size, neutral());
    }

    /* Change this function for updating nodes */
    pi update(pi f, pi s){
        return {(f.first ^ s.first), (f.second ^ s.second)};
    }

    /* Return neutral element */
    pi neutral(){
        return {0, 0};
    }

    pi get(int l, int r){
        return get(l, r, 0, 0, size);
    }

    pi get(int l, int r, int x, int lx, int rx)
    {
        if(lx>=r||rx<=l) return neutral();
        if(lx>=l&&rx<=r){
            return xor_val[x];
        }

        int m = (lx + rx)/2;

        pi a = get(l, r, 2*x + 1, lx, m);
        pi b = get(l, r, 2*x + 2, m, rx);

        return update(a, b);
    }

    void set(int l, int r, pi v){
        set(l, r, v, 0, 0, size);
    }

    void set(int l, int r, pi v, int x, int lx, int rx){
        if( lx >= r || rx <= l) return;
        if( lx >= l && rx <= r){
            xor_val[x] = v;
            return;
        }

        ll m = (lx + rx)/2;

        set(l, r, v, 2*x + 1, lx, m);
        set(l, r, v, 2*x + 2, m, rx);

        xor_val[x] = update(xor_val[2*x + 1], xor_val[2*x + 2]);
    }
};

vi vals, parent, dep;
vvi g;

int ti;
vpi range;

void euler(int par = 0){
    range[par].first = ti++;
    for(int e: g[par]){
        dep[e] = dep[par] + 1;
        euler(e);
    }
    range[par].second = ti;
}

void test_case(){
    vals.assign(1, -1);

    int q;
    cin>>q>>vals[0];

    assert(q >= 1 && q <= 5e5);
    assert(vals[0] >= 0 && vals[0] <= 1e9);

    vvi queries(q);
    int n = 1;

    parent.assign(1, -1);

    for(int que = 0; que < q; que++){
        int ty;
        cin>>ty;

        if(ty == 1){
            int par, val;
            cin>>par>>val;
            assert(par >= 1 && par <= n);
            assert(val >= 0 && val <= 1e9);
            par--;
            n++;
            parent.push_back(par);
            vals.push_back(val);
            queries[que] = {ty, par, val};
        } else {
            int node;
            cin>>node;
            assert(node >= 1 && node <= n);
            node--;
            queries[que] = {ty, node};
        }
    }

    assert(n == sz(vals) && n == sz(parent));

    g.assign(n, {});
    for(int i = 1; i < n; i++){
        g[parent[i]].push_back(i);
    }

    ti = 0;
    dep.assign(n, 0);
    range.assign(n, {-1, -1});
    euler();

    int cur_node = 1;
    segtree sg;
    sg.init(n);

    sg.set(0, 1, {vals[0], 0});

    vi ans;

    for(int i = 0; i < q; i++){
        vi &cur = queries[i];
        int ty = cur[0];
        if(ty == 1){
            pi update = {0, 0};
            int val = cur[2];
            if(dep[cur_node]&1){
                update.second = val;
            } else {
                update.first = val;
            }
            sg.set(range[cur_node].first, range[cur_node].first + 1, update);
            cur_node++;
        } else {
            int node = cur[1];
            pi res = sg.get(range[node].first, range[node].second);
            if(dep[node]&1){
                ans.push_back(res.second);
            } else {
                ans.push_back(res.first);
            }
        }
    }

    for(int e: ans){
        cout<<e<<"\n";
    }
}

We will first detail the $$$O(k^2log(n))$$$ solution, after which it should be easy to optimize to $$$O(klog(k)log(n))$$$ using NTT. We will calculate $$$dp[i][x]$$$ = probability of a number which is $$$2^i$$$ digits long to give a remainder of $$$x$$$ modulo $$$k$$$. The transition for this dp is

Thus, we can calculate this dp for all necessary $$$i$$$ in $$$O(k^2log(n))$$$ time. We can now observe that this is very similar to multiplying two polynomials:

This is exactly like multiplying two polynomials, and we can now calculate dp in $$$O(klog(k)log(n))$$$ time, and the answer can be constructed using this dp in a similar manner.

#include "bits/stdc++.h"
using namespace std;
    
#pragma GCC target ("avx2")
#pragma GCC optimize ("O3")
#pragma GCC optimize ("unroll-loops")
    
#define ll                    long long
#define ull                   unsigned long long
#define loop(i,a,b)           for(int i=(int)a;i<(int)b;++i)
#define loopr(i,a,b)          for(int i=(int)a;i>=(int)b;--i)
#define count_1(n)            __builtin_popcountll(n)
#define pb                    push_back
#define vi                    vector<int>
#define vvi                   vector<vector<int>>
const ll MOD = 998244353; //1000000007 998244353
const ll INF = 1e17  + 10; //careful
    
inline ll po(ll a, ll b) { ll res = 1; for (; b; b >>= 1) { if (b & 1)res = res * a; a = a * a; }return res; }
inline ll gcd (ll a, ll b) { ll r; while(b) { r = a%b; a=b; b=r; } return a; }
    
template<typename T> inline void debug(T a) { cout << a << endl; }
template<typename T, typename... Args> inline void debug(T a, Args... args) { cout << a << " "; debug(args...); }
template<typename T> inline void printlist(T& a) { for(auto i : a) cout << i << " "; cout << endl; }
 
inline ll modpow(ll a, ll b, ll mod) { ll res = 1; for (; b; b >>= 1) { if (b & 1)res = (res * a)%mod; a = (a * a)%mod; }return res; }
 
template <ll mod = MOD>
struct modint {
    ll val;
    modint() : val(0) {}
    modint(ll a) : val(a) { val += mod; if(val >= mod) val -= mod; }
    void operator = (modint a) { val = a.val; }
    bool operator == (modint a) { return val == a.val; } 
    bool operator != (modint a) { return val != a.val; }
    inline modint& operator += (modint a) { val = val + a.val; if(val >= mod) val-=mod; return *this; }
    inline modint& operator -= (modint a) { val = val - a.val + mod; if(val >= mod) val -= mod; return *this; }
    inline modint& operator *= (modint a) { val = (val*a.val)%mod; return *this; }
    friend modint operator + (modint a, modint b) { return (a += b); }
    friend modint operator - (modint a, modint b) { return (a -= b); } 
    friend modint operator * (modint a, modint b) { return (a *= b); }
    friend modint operator ^ (modint a, ll b) { modint res = 1; for (; b; b >>= 1) { if (b & 1)res*=a; a*=a; } return res; }
    friend modint operator / (modint num, modint den) { den = den ^ (mod-2); return (num*den); }
    friend ostream & operator << (ostream &out, modint a) { out << a.val; return out; }
    friend istream & operator >> (istream &in, modint& a) { in >> a.val; return in; }
};
 
using mint = modint<MOD>;
 
 
const int root = 565042129;
const int root_1 = 950391366;
const int root_pw = 1 << 20;
 
//MOD 998244353
void fft(vector<mint> & a, bool invert) {
    int n = a.size();
 
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1)
            j ^= bit;
        j ^= bit;
 
        if (i < j)
            swap(a[i], a[j]);
    }
 
    for (int len = 2; len <= n; len <<= 1) {
        mint wlen = invert ? root_1 : root;
        for (int i = len; i < root_pw; i <<= 1)
            wlen = wlen * wlen;
 
        for (int i = 0; i < n; i += len) {
            mint w = 1;
            for (int j = 0; j < len / 2; j++) {
                mint u = a[i+j], v = a[i+j+len/2] * w;
                a[i+j] = u + v;
                a[i+j+len/2] = u - v;
                w = w * wlen;
            }
        }
    }
 
    if (invert) {
        mint n_1 = mint(1) / n;
        for (mint & x : a)
            x = x * n_1;
    }
}
 
//result stored in a
void multiply(vector<mint> & a, vector<mint> const& b) {
    int n = 1; vector<mint> fb = b;
    while (n < a.size() + b.size()) n <<= 1;
    a.resize(n);
    fb.resize(n);
 
    fft(a, false);
    fft(fb, false);
    for (int i = 0; i < n; i++)
        a[i] *= fb[i];
    fft(a, true);
    
}
 
 
void testcases(int test) {
    int n,k; cin>>n>>k; ll sum = 0;
    vector<ll> po102(30);  po102[0] = 10%k; loop(i,1,30) { po102[i] = (po102[i-1] * po102[i-1])%k; }
    vector<mint> p(10); loop(i,0,10) { ll x; cin>>x; p[i] = x; sum += x; }
    loop(i,0,10) p[i] = p[i] / sum;
    vector<vector<mint>> dp(30, vector<mint>(k));
    loop(i,0,10) { dp[0][i%k] += p[i]; }
    loop(i,1,30) {
        ll run = 0;
        vector<mint> temp(k); loop(j,0,k) { temp[run] += dp[i-1][j]; run += po102[i-1]; if(run >= k) run -= k; }
        multiply(temp, dp[i-1]);
        run = 0;
        loop(j,0,temp.size()) { dp[i][run] += temp[j]; run++; if(run >= k) run -= k; }
    } 
    vector<mint> probs(k); int flag = 0;
    int pot = 0;
    while(n) {
        if(n&1) {
            if(!flag) {
                probs = dp[pot];
                flag = 1;
            }
            else {
                vector<mint> nprobs(k);
                ll run = 0;
                vector<mint> temp(k); loop(j,0,k) { temp[run] += probs[j]; run += po102[pot]; if(run >= k) run -= k; }
                multiply(temp, dp[pot]);
                run = 0;
                loop(j,0,temp.size()) { nprobs[run] += temp[j]; run++; if(run >= k) run -= k; }
                swap(probs, nprobs);
            }
        }
        n>>=1; pot++;
    }
    cout << probs[0] << "\n";
}   
    
int32_t main()
{
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);
    
    
    auto start = chrono::high_resolution_clock::now();
    
    int t = 1;
    cin >> t;
    loop(i,0,t) { testcases(i); }
    
    auto stop = chrono::high_resolution_clock::now();
    auto duration = chrono::duration_cast<chrono::milliseconds>(stop - start);
    //cout << duration.count() << " ms\n";
}

Here let's say the sum of time of all students is $$$S$$$, So we can clearly observe that the lower bound for the answer is $$$S$$$.

Here we will have two cases:

1. When largest element of array $$$\leq S/2$$$ — In this case our answer will always be $$$S$$$. There are many ways to achieve this answer, one is to take maths viva from smallest to largest and Science from largest to smallest.

2. When the largest element of array $$$>$$$ $$$S/2$$$ — In this case, let's say the largest element is $$$X$$$. So our answer will be $$$2*X$$$ always. We can achieve this answer Let's say Alice is taking the viva of $$$X$$$; in that time, Bob will take the viva of all other students, and then Bob will take the viva of $$$X$$$, and Alice will take the viva of all other students.

#include<bits/stdc++.h>
using namespace std;
#pragma GCC optimize("O3")
#pragma GCC target("avx2")
#define M 1000000007
#define ll long long int
#define pb push_back
#define pf push_front
#define all(x) x.begin(), x.end()
#define pob pop_back
#define pof pop_front
#define debug1(x) cout<<#x<<" "<<x<<endl;
#define debug2(x,y) cout<<#x<<" "<<x<<"   "<<#y<<" "<<y<<endl;
#define debug3(x,y,z) cout<<#x<<" "<<x<<"   "<<#y<<" "<<y<<"   "<<#z<<" "<<z<<endl;
#define debug4(x,y,z,p) cout<<#x<<" "<<x<<"   "<<#y<<" "<<y<<"   "<<#z<<" "<<z<<" "<<#p<<" "<<p<<endl;
#define vll      vector<ll> 
#define vvll     vector<vector<ll>>
#define fi first
#define se second
#define boost ios_base::sync_with_stdio(false);cin.tie(NULL);cout.tie(NULL);
#define inf 1e18
#define endl '\n'
typedef pair<int, int>		pii;
typedef pair<ll, ll>		pll;
#define pi (3.14159265358979323846264338327950288)
template<typename T>
void printv(const T& t) {std::copy(t.cbegin(), t.cend(), std::ostream_iterator<typename T::value_type>(std::cout, ", "));cout<<endl;}
ll modpower(ll a, ll b, ll mod) { ll res = 1; for (; b; b >>= 1) { if (b & 1)res = res * a; res %= mod; a = a * a; a %= mod; }return res; }
struct custom_hash {
    static uint64_t splitmix64(uint64_t x) {
        // http://xorshift.di.unimi.it/splitmix64.c
        x += 0x9e3779b97f4a7c15;
        x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9;
        x = (x ^ (x >> 27)) * 0x94d049bb133111eb;
        return x ^ (x >> 31);
    }

    size_t operator()(uint64_t x) const {
        static const uint64_t FIXED_RANDOM = chrono::steady_clock::now().time_since_epoch().count();
        return splitmix64(x + FIXED_RANDOM);
    }
};
//unordered_map<ll,ll, custom_hash>m1;
ll expo(ll a, ll b, ll mod) {ll res = 1; while (b > 0) {if (b & 1)res = (res * a) % mod; a = (a * a) % mod; b = b >> 1;} return res;}
ll mod_inv(ll a, ll b) {return expo(a, b - 2, b);}
ll mod_add(ll a, ll b, ll m) {a = a % m; b = b % m; return (((a + b) % m) + m) % m;}
ll mod_mul(ll a, ll b, ll m) {a = a % m; b = b % m; return (((a * b) % m) + m) % m;}
ll mod_sub(ll a, ll b, ll m) {a = a % m; b = b % m; return (((a - b) % m) + m) % m;}
ll mod_div(ll a, ll b, ll m) {a = a % m; b = b % m; return (mod_mul(a, mod_inv(b, m), m) + m) % m;}
//-------------------------------Template--Above-----------------------------------------------
void solve(){
	ll n;cin>>n;
	vll v(n);
	ll sum=0;
	for(ll i=0;i<n;i++){
		cin>>v[i];
		sum+=v[i];
	}
	sort(all(v));
	
	ll x=v[n-1];
	if(x<=sum/2){
		cout<<sum<<endl;
	}
	else{
		cout<<2*x<<endl;
	}
	
    return;
}
int main()
{
    boost
    /*
    freopen("moocast.in", "r", stdin);
	freopen("moocast.out", "w", stdout);
	*/
    ll t;
    t=1;
   // cin>>t;
    for(ll i=1;i<=t;i++){
    	//cout << "Case #" << i << ": "; solve();
    	solve();
	}
    return 0;
}

The Boruvka algorithm for finding an MST can pretty much be directly applied in this problem to find the MST. The main task in this algorithm is to: for each connected component (in the MST being built), find the lowest weight edge that connects this component to another; and add all these edges to the MST. We can do this by maintaining a set that initially contains all vertices $$$[1,n]$$$ during each iteration of Boruvka, and:

For each component,

• Remove all members of this component from the set.
• For each member vertex $$$i$$$, find the lowest weight edge connecting it to a vertex outside its component by simply scanning the set from left to right and stopping at the first vertex $$$j$$$ such that edge $$$(i,j)$$$ is not deleted.
• Take the cheapest of all such edges found, and that s the required edge for this component.
• Add all the member vertices back to the set, and repeat this for all other connected components.

This solution will run in $$$((n+m)log^2)$$$ amortized time, as each 'removed' edge is inspected during an iteration of Boruvka atmost twice, and Boruvka itself will run for atmost $$$O(log(n))$$$ iterations. This was intended to TLE, as there is a simple optimization here, we can maintain for each vertex $$$i$$$, the last vertex that it found in the above process, and simply start its search from this vertex in the next iteration, as the vertex pairs before are already found to be not in the graph. This allows the solution to run in $$$O(nlog^2(n) + mlog(n))$$$ time.

#include "bits/stdc++.h"
using namespace std;
    
#pragma GCC target ("avx2")
#pragma GCC optimize ("O3")
#pragma GCC optimize ("unroll-loops")
    
#define ll                    long long
#define ull                   unsigned long long
#define loop(i,a,b)           for(int i=(int)a;i<(int)b;++i)
#define loopr(i,a,b)          for(int i=(int)a;i>=(int)b;--i)
#define count_1(n)            __builtin_popcountll(n)
#define pb                    push_back
#define vi                    vector<int>
#define vvi                   vector<vector<int>>
const ll MOD = 1000000007; //1000000007 998244353
const ll INF = 1e17  + 10; //careful
    
inline ll po(ll a, ll b) { ll res = 1; for (; b; b >>= 1) { if (b & 1)res = res * a; a = a * a; }return res; }
inline ll gcd (ll a, ll b) { ll r; while(b) { r = a%b; a=b; b=r; } return a; }
    
template<typename T> inline void debug(T a) { cout << a << endl; }
template<typename T, typename... Args> inline void debug(T a, Args... args) { cout << a << " "; debug(args...); }
template<typename T> inline void printlist(T& a) { for(auto i : a) cout << i << " "; cout << endl; }
 
inline ll modpow(ll a, ll b, ll mod) { ll res = 1; for (; b; b >>= 1) { if (b & 1)res = (res * a)%mod; a = (a * a)%mod; }return res; }
 
struct dsu {
    vector<int> p,r;
    dsu(int n) {
        loop(i,0,n+1) {
            p.pb(i); r.pb(0);
        }
    }
    int get(int a) {
        return (p[a] = (p[a]==a)?a:get(p[a]));
    }
    void uni(int a, int b) {
        a = get(a); b = get(b);
        if(a==b) return;
        if(r[a]==r[b]) r[a]++;
        if(r[a] < r[b]) swap(a,b);
        p[b] = a;
    }
};
 
void testcases(int test) {
    int n,m; cin>>n>>m;
    set<pair<int,int>> edges;
    vector<ll> a(n); loop(i,0,n) cin>>a[i];
    loop(i,0,m) {
        int x,y; cin>>x>>y; x--; y--;
        edges.insert({x,y}); edges.insert({y,x});
    }   
    dsu d(n+1); int cc = n;
    set<pair<ll,ll>> hold; loop(i,0,n) hold.insert({a[i],i}); int flag = 1;
    vector<pair<int,int>> mst; vector<pair<int,int>> last(n, {-1,-1});
    while(cc > 1) {
        vvi components(n);
        loop(i,0,n) components[d.get(i)].pb(i);
        vector<pair<int,int>> minedge(n, {-1,-1});
        loop(i,0,n) {
            if(components[i].empty()) continue;
            for(int j : components[i]) hold.erase({a[j],j});
            for(int j : components[i]) {
                auto itr = hold.lower_bound(last[j]);
                while(itr != hold.end() && edges.find({j, itr->second}) != edges.end()) itr++;
                if(itr == hold.end()) { last[j] = {100000000,100000000}; continue; }
                last[j] = *itr;
                if(minedge[i].first == -1) minedge[i] = {j, itr->second};
                else if(a[j] * itr->first < a[minedge[i].first] * a[minedge[i].second]) minedge[i] = {j, itr->second}; 
            }
            for(int j : components[i]) hold.insert({a[j],j});
        }
        loop(i,0,n) {
            if(minedge[i].first == -1 && !components[i].empty()) { flag = 0; break; }
            if(components[i].empty()) continue;
            if(d.get(minedge[i].first) != d.get(minedge[i].second)) { 
                cc--; 
                mst.pb({minedge[i].first, minedge[i].second});
                d.uni(minedge[i].first, minedge[i].second);
            }
        }
        if(!flag) break;
    }
    if(!flag) { cout << -1 << "\n"; return; }
    ll ans = 0; for(pair<int,int> i : mst) { ans += a[i.first] * a[i.second]; }
    debug(ans);
    for(pair<int,int> i : mst) debug(i.first + 1, i.second + 1);
}   
    
int32_t main()
{
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);
    
    
    auto start = chrono::high_resolution_clock::now();
    
    int t = 1;
    cin >> t;
    loop(i,0,t) { testcases(i); }
    
    auto stop = chrono::high_resolution_clock::now();
    auto duration = chrono::duration_cast<chrono::milliseconds>(stop - start);
    //cout << duration.count() << " ms\n";
}