Pitchers (Editorial)

Idea

For each query (l, r, k), we need to compute \displaystyle\sum_{i=l}^{r} \frac{k}{a_i^{p_i}}, where p_i — the largest integer such that $k$ is divisible by a_i^{p_i}.

The naive approach (for each query, traverse the segment and compute the exponent separately for each jug) works in O(Q \cdot n \cdot \log k), which is far too slow for n, Q \le 3 \cdot 10^6.

Key observation: reformulation via divisors

Fix $k$. For each possible value a_i = d (where $d$ — a divisor of $k$), the value \frac{k}{d^p} is the same and depends only on $k$ and $d$. Denote this value by f(k, d).

If a_i does not divide $k$, then p_i = 0 and the contribution is simply $k$.

Then the answer to the query can be written as:

where \text{cnt}(d, l, r) — the number of elements a_i = d on the segment [l, r].

Precomputing f(k, d)

For each $k$ from $2$ to $mx$ and each of its divisors $d$, we compute f(k, d) = \frac{k}{d^p}, where $p$ — the maximum power of $d$ dividing $k$. This can be done recursively:

• If \frac{k}{d} is divisible by $d$, then f(k, d) = f\!\left(\frac{k}{d},\; d\right).

• Otherwise, f(k, d) = \frac{k}{d}.

We find the divisors of each number $k$ using the classic "sieve": for each $i$ from $2$ to $mx$, we add $i$ to the divisor list of all its multiples. The total number of pairs $(k,d)$ is O(\mathit{mx} \ln \mathit{mx}).

Offline query processing (sweep-line)

To efficiently find \text{cnt}(d, l, r) for each query, we use prefix sums and an offline sweep from left to right over the array $a$:

1. For each query $i$ with segment [l_i, r_i], create two "events":

   • $+$ at position r_i (add the contribution);

   • $-$ at position l_i - 1 (subtract the contribution).

2. Traverse the array from left to right, maintaining a counter $\text{cnt}[v]$ — how many times the value $v$ has appeared among a_0, a_1, \ldots, a_i.

3. When we encounter a query event $\text{num}$ at position $i$, iterate over all divisors $d$ of k_{\text{num}} and add (or subtract) the contribution f(k_{\text{num}}, d) \cdot \text{cnt}[d] to the answer.

Thus, for each query we obtain \text{cnt}(d, l, r) = \text{cnt}(d, 0, r) - \text{cnt}(d, 0, l-1) without explicit computation.

Final formula

After the sweep-line, for each query $i$ we have a partial sum \mathit{res}_i (contributions from those jugs whose $a$-value divides k_i) and a counter c_i (the number of such jugs). We complete the answer with the contribution from the "neutral" jugs:

Complexity

• Building divisor lists (sieve): O(\mathit{mx} \ln \mathit{mx}).

• Precomputing $f$: O(\mathit{mx} \ln \mathit{mx}).

• Sweep-line: for each query we iterate over the divisors of k_i, whose number is O(\ln \mathit{mx}) on average, so in total O(Q \cdot d(k)), where $d(k)$ — the number of divisors of $k$.

Overall: O\!\left((n + \mathit{mx}) \ln \mathit{mx} + Q \cdot d(k)\right), where \mathit{mx} = \max(a_i, k_i) \le 5 \cdot 10^5.

#include <bits/stdc++.h>

#define endl '\n'

using namespace std;

using ll = long long;

const ll DIM = 500007;

struct query {
    int l, r, k;
};

vector < int > d[DIM];
vector < int > f[DIM];

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

    int n, q;
    cin >> n >> q;

    vector < int > a(n);
    for(int i = 0; i < n; i++) cin >> a[i];

    vector < query > Q(q);
    for(int i = 0; i < q; i++) {
        cin >> Q[i].l >> Q[i].r >> Q[i].k;
    }

    int mx = 0; 
    for(int i = 0; i < n; i++) mx = max(mx, a[i]);
    for(int i = 0; i < q; i++) mx = max(mx, Q[i].k);

    // sieve: build divisor lists
    for(int i = 2; i <= mx; i++) {
        for(int j = i; j <= mx; j += i) {
            d[j].push_back(i);
        }
    }
    vector < int > ptr(mx + 1, 0);

    // precompute f[k][j] = k / d[k][j]^p
    for(int i = 2; i <= mx; i++) {
        f[i].resize(d[i].size());
        for(int j = 0; j < (int)d[i].size(); j++) {
            if((i / d[i][j]) % d[i][j] == 0) {
                int to = i / d[i][j];
                while(d[to][ptr[to]] < d[i][j]) ptr[to]++;
                f[i][j] = f[to][ptr[to]];
            }
            else f[i][j] = i / d[i][j];
        }
    }

    // offline sweep
    vector < vector < int > > E(n);
    
    for(int i = 0; i < q; i++) {
        if(Q[i].l != 0) E[Q[i].l - 1].push_back(-(i + 1));
        E[Q[i].r].push_back((i + 1));
    }

    vector < ll > res(q);
    vector < ll > c(q);
    vector < int > cnt(mx + 1, 0);

    for(int i = 0; i < n; i++) {
        cnt[a[i]]++;

        for(auto e: E[i]) {
            int num = abs(e) - 1;
            for(int j = 0; j < (int)d[Q[num].k].size(); j++) {
                if(e > 0) {
                    c[num] += (ll)cnt[d[Q[num].k][j]];
                    res[num] += (ll)(f[Q[num].k][j])
                              * (ll)cnt[d[Q[num].k][j]];
                }
                if(e < 0) {
                    c[num] -= (ll)cnt[d[Q[num].k][j]];
                    res[num] -= (ll)(f[Q[num].k][j])
                              * (ll)cnt[d[Q[num].k][j]];
                }
            }
        }
    }

    for(int i = 0; i < q; i++) {
        res[i] += ll((Q[i].r - Q[i].l + 1) - c[i])
                * (ll)Q[i].k;
        cout << res[i] << endl;
    }

    return 0;
}