STQ - EDITORIAL

PROBLEM LINK:

Practice
Contest: Division 3
Contest: Division 2
Contest: Division 1

Author: Rahim Mammadli
Tester: Anshu Garg
Editorialist: Istvan Nagy

DIFFICULTY:

HARD

PREREQUISITES:

PROBLEM:

You are given an array A of N positive integers. You are allowed to perform any of the following operations any number of times:

  • Pick an index i (1 \le i \le n), and a prime number p. Assign A_i \times p to A_i.
  • Pick an index i (1 \le i \le n), and a prime number p such that A_i is divisible by p. Assign A_i \div p to a_i.

The goal is to make all the numbers in array equal using the minimum number of operations. Among all sequences of operations using the minimum number of operations, denote S to be the smallest final value of all elements in A. You are given Q queries of the form (L_{i}, R_{i}). For each query, you need to calculate the S value of the contiguous subarray A_{L_i \ldots R_i}.

EXPLANATION

Lets start with a simpler version of the problem when all A_i = p^{b_i}. We have to do the minimum number of operations if we choose the median of the {b_i} array. When the array has odd number of elements the median is unambiguous, but when the array has even number of elements, then all the power between the two middle element requires the same number of operations, so to get the smallest final value we have to pick the smaller one of the two.

Observation: A_i has at most 7 prime factor.
Proof: The first 8 prime product = 2\times 3\times 5\times 7\times 11\times 13\times 17\times 19\times = 8314020

Observation 2: The smallest final value has at most 13 prime factor.
Proof: When a prime is a factor of the final value, it means this prime is a divisor more than a half of the range values. E.g the Pigeonhole principle gives us the 13 upper bound.

Observation 3: Lets consider an (L,R) query. Lets pick uniform randomly k integer from the query range the probability of the final value (of the query) has a prime factor what is not represented any of the k picked item is less then 13\cdot \frac{1}{2^k}.
Proof: For one picked item the probability is less then 1/2 because the prime factor appears more than a half of the items. Because of [Boole’s inequality] (Boole's inequality - Wikipedia) we can get the upper bound of the failure.

The setter choose k=40 what guarantees the probability of the fail is less then 13/2^{40}. The probability of success for all query is :

\left(1-Q\cdot \frac{13}{2^{40}}\right)\geq 0.999997635

We have to check at most 105 primes per query because at 25 primes are less than 100. Every query point can check at most 2 primes which is greater than 100, so we get the number of primes per query cannot exceed 25+2*40.

Create an auxiliary array for the queries, it has N item (all item is a vector) and push the the query index to the l_i-1 and also r_i.
Scan through the A array use another auxiliary for countine prime power p^j, how many times occurs increment p_j counter only if p^j \mid A_i and p^{j+1} \nmid A_i. After we processed the A_i item and updated the auxiliary array counters, we can check is there a query for the actual position. In case there is a query starting point then we subtract the belonging counters for the 105 primes and its powers, if its an ending point then we add the counters. So at the end we know how much times a prime power occurs in the query range, so we can calculate the median easily. The number of 105 prime powers of the 40 primes cannot exceed (18 + 12+ ...+ 80 \cdot 1)= 93+80=173, this number of steps also enough to calculate the median.

TIME COMPLEXITY:

O(N\cdot log(N)+ Q \cdot (105 + 173))

SOLUTIONS:

Setter's Solution
#pragma GCC optimize("O3")
#pragma GCC target("sse4,avx2,abm,fma,tune=native")
#pragma GCC optimize("unroll-loops")
 
#include <bits/stdc++.h>
using namespace std;
#define mp(a,b) make_pair(a,b)
#define ff first
#define setp(a) setprecision(a)<<fixed
#define ss second
#define fori(v) for(ll i=0; i<v; i++)
#define forj(v) for(ll j=0; j<v; j++)
#define fork(v) for(ll k=0; k<v; k++)
#define forl(v) for(ll l=0; l<v; l++)
#define fort(v) for(ll t=0; t<v; t++)
#define forz(v) for(ll z=0; z<v; z++)
#define forx(v) for(ll x=0; x<v; x++)
#define fory(v) for(ll y=0; y<v; y++)
#define ll long long
#define lll __int128
#define pb(a) push_back(a)
typedef long double ld;
const ll INF = 0x3f3f3f3f;
const ll inf =  pow(10,9);
const ll modulo = pow(10,9) + 7;
 
#define MAX 1'000'002 

bool isp[MAX];
int dvsz[MAX];
int dv[MAX][7];
int dvpv[MAX][7];
int maxpov[MAX];

void pre(){
    fori(MAX){
        isp[i] = 1;
    }
    isp[0] = isp[1] = 0;
    for(ll i = 2; i<MAX; i++){
        if(isp[i]){
            for(ll j = i; j<MAX; j+=i){
                dv[j][dvsz[j]] = i;
                dvsz[j]++;
                isp[j] = 0;
            }
            int cur= i;
            maxpov[i] = 1;
            while(cur < MAX){
                cur*=i;
                maxpov[i]++;
            }
        }
    }
    fori(MAX){
        forj(dvsz[i]){
            int pr = dv[i][j];
            int cur = pr;
            int cnt = 0;
            while(i%cur == 0){
                ++cnt;
                cur *= pr;
            }
            dvpv[i][j] = cnt;
        }
    }
}

#define MAXN 1'000'010
#define MAXQ 200'010

int arr[MAXN];

int quepr[MAXQ][108];
int quesz[MAXQ];
int cntpr[MAXQ][108];
int seenpr[MAX];

vector<int> qs[MAXN];    // queries
int ls[MAXQ];
int rs[MAXQ];

int prpv[MAX][20];


int satpr[MAXQ][14];    // primes appearing more than half
int satsz[MAXQ];        
vector<int> pvs[MAXQ][14];


// #define cin in
// #define cout out


void deal(){    
    mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
    pre();
    
    ifstream in;
    ofstream out;
    int n , q;
    cin>>n>>q;
    assert(n < MAXN);
    assert(q < MAXQ);
    for(ll i = 1; i<=n; i++){
        cin>>arr[i];
        assert(arr[i] < MAX);
    }
    for(ll i = 1; i<=q; i++){
        int li, ri;
        cin>>li>>ri;
        assert(li <= ri);
        assert(li >= 1);
        assert(ri <= n);

        ls[i] = li, rs[i] = ri;
        qs[li-1].pb(i);
        qs[ri].pb(i);
        forj(40){
            int cur = li + rng()%(ri-li+1);
            int num = arr[cur];
            int sz = dvsz[num];
            fork(sz){
                int pr = dv[num][k];
                if(seenpr[pr] != i){
                    seenpr[pr] = i;
                    quepr[i][quesz[i]] = pr;
                    ++quesz[i];
                }
            }
        }
    }
    
    fori(n+1){
        int cur = arr[i];
        {
            int sz = dvsz[cur];
            forj(sz){
                int pr = dv[cur][j];
                prpv[pr][0]++;
            }
        }
        for(auto& el : qs[i]){
            int mlt = 1;
            if(i == ls[el] - 1){
                mlt = -1;
            }
            int sz = quesz[el];
            forj(sz){
                int pr = quepr[el][j];
                cntpr[el][j]+=mlt*prpv[pr][0];
            }
        }
    }
    
    
    for(ll i = 1; i<= q; i++){
        int sz = quesz[i];
        int range = (rs[i] - ls[i] + 1);
        forj(sz){
            if(cntpr[i][j] >= range/2+1){
                int id = satsz[i];
                int pr = quepr[i][j];
                satpr[i][id] = pr;
                pvs[i][id] = vector<int>(maxpov[pr], 0);
                pvs[i][id][0] = range  - cntpr[i][j];
                ++satsz[i];
            }
        }
    }
    
    
    fori(MAX){
        prpv[i][0] = 0;
    }
    
    fori(n+1){
        int cur = arr[i];
        {
            int sz = dvsz[cur];
            forj(sz){
                int pr = dv[cur][j];
                int num = dvpv[cur][j];
                prpv[pr][num]++;
            }
        }
        for(auto& el : qs[i]){
            int mlt = 1;
            if(i == ls[el] - 1){
                mlt = -1;
            }
            int sz = satsz[el];
            forj(sz){
                int pr = satpr[el][j];
                int pv = maxpov[pr];
                fork(pv){
                    pvs[el][j][k]+=mlt*prpv[pr][k];
                }
            }
        }
    }
    
    for(ll i = 1; i<=q; i++){
        ll ans = 1;
        int sz = satsz[i];
        int range = (rs[i] - ls[i] + 1);
        forj(sz){
            int pr = satpr[i][j];
            int s = pvs[i][j][0];
            for(int k = 1; s<range-s; k++){
                s+=pvs[i][j][k];
                ans*=pr;
            }
        }
        cout<<ans<<'\n';
    } 
}
 
int main()  {
    cin.tie(0);
    ios_base::sync_with_stdio(0);
    deal();
}
Tester's Solution
#include<bits/stdc++.h>
using namespace std ;

#define ll              long long 
#define pb              push_back
#define all(v)          v.begin(),v.end()
#define sz(a)           (ll)a.size()
#define F               first
#define S               second
#define INF             2000000000000000000
#define popcount(x)     __builtin_popcountll(x)
#define pll             pair<ll,ll>
#define pii             pair<int,int>
#define ld              long double

const int M = 1000000007;
const int MM = 998244353;

template<typename T, typename U> static inline void amin(T &x, U y){ if(y<x) x=y; }
template<typename T, typename U> static inline void amax(T &x, U y){ if(x<y) x=y; }

#ifdef LOCAL
#define debug(...) debug_out(#__VA_ARGS__, __VA_ARGS__)
#else
#define debug(...) 2351
#endif

const int N = 1e6 + 2;

#define runSieve
bool prime[N];
int pf[N];
void sieve()
{
    int k=sqrt(N-2);
    for(int i=1;i<=N-2;++i)
        prime[i]=true;
    for(int i=2;i<=k;++i)
    {
        if(prime[i])
        {
            for(int k=i*i;k<=N-2;k+=i)
            {
                prime[k]=false;
                pf[k]=i;
            }
        }
    }
}

mt19937 rng(chrono::high_resolution_clock::now().time_since_epoch().count());

vector<int> q[N];

int pt2[N];
int divs[N][8], fq[N][8], divcnt[N];
int cnt[N];

int f[N][20];


const int QR = 4e5 + 2;

int ans[QR];
int a[N];

int qry[QR][90];

int pt[QR];
int prs[QR][90];

int aux[100];

bool mp[N];
int l[QR], r[QR];

int mxp[N],pw[N][13];

int fnl[QR][13][20];


int _runtimeTerror_()
{
    int n, Q;
    cin >> n >> Q;
    for(int i=1;i<=n;++i)
        cin >> a[i];
    for(int i=2;i<=1e6;++i) {
        ll x = i;
        if(!prime[i])
            continue;
        mxp[i] = 0;
        while(x <= N) {
            ++mxp[i];
            x *= i;
        }
        int f = 1;
        int cnt = 0;
        x = 1;
        while(x <= N) {
            pw[i][cnt] = x;
            ++cnt;
            x *= i;
        }
    }
    for(int i=2;i<=1e6;++i) {
        int tmp = i;
        while(tmp > 1 && !prime[tmp]) {
            int z = pf[tmp];
            divs[i][divcnt[i]] = z;
            int &cnt = fq[i][divcnt[i]];
            while(tmp % z == 0) {
                tmp /= z;
                ++cnt;
            }
            ++divcnt[i];
        }
        if(tmp > 1) {
            divs[i][divcnt[i]] = tmp;
            fq[i][divcnt[i]] = 1;
            ++divcnt[i];
        }
    }
    auto get = [&](int x,int i) {
        for(int j=0;j<divcnt[x];++j) {
            if(mp[divs[x][j]] == 0) {
                mp[divs[x][j]] = 1;
                prs[i][pt[i]] = divs[x][j];
                ++pt[i];
            }
        }
    };
    for(int i=1;i<=Q;++i) {
        ans[i] = 1;
        cin >> l[i] >> r[i];
        q[r[i]].push_back(i);
        q[l[i]-1].push_back(-i);
        for(int j=0;j<25;++j) {
            int id = l[i] + rng() % (r[i] - l[i] + 1);
            get(a[id],i);
        }
        for(int j=0;j<pt[i];++j) {
            mp[prs[i][j]] = 0;
        }
    }
    for(int i=1;i<=n;++i) {
        for(int j=0;j<divcnt[a[i]];++j) {
            ++cnt[divs[a[i]][j]];
        }
        for(auto &j:q[i]) {
            int p = j > 0 ? j : -j;
            int fac = j / p;
            for(int k=0;k<pt[p];++k) {
                qry[p][k] += cnt[prs[p][k]] * (fac);
            }
        }
    }
    for(int i=1;i<=Q;++i) {
        for(int j=0;j<pt[i];++j) {
            aux[j] = prs[i][j];
        }
        int sz = pt[i];
        pt[i] = 0;
        for(int j=0;j<sz;++j) {
            if(qry[i][j] > (r[i] - l[i] + 1)/2) {
                prs[i][pt[i]] = aux[j];
                fnl[i][pt[i]][0] += r[i] - l[i] + 1 - qry[i][j];
                ++pt[i];
            }
        }
        assert(pt[i] < 12);
    }
    for(int i=1;i<=n;++i) {
        for(int j=0;j<divcnt[a[i]];++j) {
            ++f[divs[a[i]][j]][fq[a[i]][j]];
        }
        for(auto &j:q[i]) {
            int p = j > 0 ? j : -j;
            int fac = j / p;
            for(int k=0;k<pt[p];++k) {
                for(int l=0;l<=mxp[prs[p][k]];++l) {
                    fnl[p][k][l] += f[prs[p][k]][l] * fac;
                }   
            }
            if(j > 0) {
                for(int k=0;k<pt[p];++k) {
                    int g = 0;
                    int len = (r[p] - l[p] + 2)/2;
                    for(int l=0;l<=mxp[prs[p][k]];++l) {
                        g += fnl[p][k][l];
                        if(g >= len) {
                            if(prs[p][k] == 2) {
                                ans[p] *= 1 << l;
                            }
                            else {
                                ans[p] *= pw[prs[p][k]][l];
                            }
                            break;
                        }
                    }
                }
            }
        }
    }
    for(int i=1;i<=Q;++i) {
        cout << ans[i] << "\n";
    }
    return 0;
}

int main()
{
    ios_base::sync_with_stdio(0);cin.tie(0);cout.tie(0);
    #ifdef runSieve
        sieve();
    #endif
    #ifdef NCR
        initialize();
    #endif
    int TESTS = 1;
    //cin >> TESTS;
    while(TESTS--)
        _runtimeTerror_();
    return 0;
}

If you have other approaches or solutions, let’s discuss in comments.

1 Like

I understood all the observations but didn’t get the part “We have to check at most 13*7=91 primes per query”. Can anyone help me understand it better?
Also If anyone has any deterministic solution, Can you please share it.

That part should be 105, so there are 25 primes < 100 , assume we hit all those primes as a result of checking 40 random points. Now if a prime is > 100 then we can have at most 2 of those per number, so we hit at most 40 * 2 + 25 = 105 primes as a result of taking all distinct primes from the 40 random points we check.

Thank you very much. Got it.

Hey!

The problem is very nice. Enjoyed upsolving it. The constraints seem very tight though, had to squeeze my solution to get it accepted.

The bound of 173 on the number of distinct prime powers seems wrong. Could you clarify a bit more where did you get that number from? I asserted it in my code and it got RE. Maybe the problem is in my code though…

Also, the editorial has a lot of typos here and there. The text overall was also very hard for me to understand. I hope you will rewrite it in a more cleaner way.

Anyway, thanks very much for the amazing problem!

Hey, glad you enjoyed solving the problem. I think the author was going to update the editorial. I am guessing your code might be a bit suboptimal if it has a hard time fitting in TL. Going of from your description, and the editorial missing the required information, I ll put some explanation here for the last part of the solution.

So we know that max number of primes we hit per query is 105. We push the index of each query to li -1 and ri. And first scan the array to get the required primes that appear more than half for each query. During this first scan, we do not consider the actual power of each prime, we just look at whether it is zero or nonzero, e.g. we answer the question how many array indices in given range divide a non-zero power of the prime without any remainders. Afterwards we go through all queries and find the actual primes that appear more than half, there s at most 13 (actually i think i can prove an upper bound of 11) of those per query. Let us call these 13 primes “interesting”. Now we allocate the needed space for maximum powers of these interesting primes (the sum of these “maximum powers” is at most 82), and scan through array one more time doing almost the same thing, but here we actually consider the actual powers. And afterwards, we find the median power for each interesting prime of each query. Hope this helps. So the final thing should be Q * ( 105 + 82 )

1 Like

Ah, now I get it completely. Indeed, my solution was far from optimal. Guess I will rewrite it now. Thank you so much for devoting your time and making it clear!

You are welcome :+1:

This problem can be solved deterministically using the solution of the following problem(a well-known problem but I don’t know the name in English).

There is an array a with length n. Answer q queries with the following format:
given l,r find any value of the subarray a[l…r] which appears strictly more than (r-l+1)/2 times or report that it doesn’t exist.
O(n+q log n)

I think you meant “all values” instead of “any value” .
Feel free to share your solution/the solution you know which has better time complexity to this “well known” problem.