PURIFYIT - Editorial

PROBLEM LINK:

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

Author: mamaleshrh
Tester: wasd2401
Editorialist: iceknight1093

DIFFICULTY:

TBD

PREREQUISITES:

Offline queries, segment tree/fenwick tree, sieve of Eratosthenes

PROBLEM:

You’re given L and R. To purify an integer x between L and R, you can do the following:

  • Choose some sequence (a_1, a_2, \ldots, a_k) such that each a_i lies between L and R, and the LCM of the sequence equals x.
  • The cost of this is the maximum a_i.

Find the minimum cost to purify everything between L and R.

EXPLANATION:

Firstly, observe that since \text{lcm}(a_1, a_2, \ldots, a_k) = x, each a_i should be a divisor of x.

Let’s first attempt to solve subtask 1: L = 1 always, meaning we have every divisor of x available to us when trying to purify it.
As it turns out, in this case, the cost for a fixed x can be found greedily: keep taking factors of x from smallest to largest, till their LCM reaches x.

Proof

This is fairly easy to prove once you observe it.
Since we’re only taking factors of x, their LCM will never exceed x (after all, x is a common multiple of all its factors).

So, if c_x is the actual answer for x, notice that simply taking every divisor that’s \leq c_x still leaves us with an LCM of x, meaning we might as well do it.

Let \text{ans}[x] denote the minimum cost necessary to purify x using some of its divisors.
As noted by the greedy argument, this can be found by just listing out all the divisors of x in sorted order, then iterating through them till the LCM becomes equal to x.
Computing the list of divisors of every number from 1 to N can be done with the sieve of Eratosthenes in \mathcal{O}(N\log N) time; and of course computing the answer once you have these lists takes \mathcal{O}(N\log N) LCM computations at worst, so it’s pretty fast for N = 3\cdot 10^5.
The constraints are low enough (and operations simple enough) that even brute-force square-root factorization is fast enough here, really.

Note that once you have \text{ans}[x] for every x, the answer for the range [1, R] is simply

\sum_{i=1}^R \text{ans}[i]

This can be found in \mathcal{O}(1) time, being just a prefix sum of the \text{ans} array.


Now, let’s move on to the full problem.
Recall the greedy idea from the first subtask? Well, it still works here!
More specifically, for any range [L, R] and integer x in that range, the minimum cost to purify x is obtained by just iterating over the divisors of x from smallest to largest (of course, only those that are \geq L in the first place), till the LCM reaches x.
The proof is exactly the same.

Unfortunately, this doesn’t let us immediately solve the problem; since L and R can vary wildly (and the sum of R-L+1 isn’t bounded).

However, observe that the cost of x doesn’t really depend on R at all: it only depends on L.
So, for a fixed L, x\geq L will contribute the same cost to any [L, R] such that R\geq x.

Using this observation, the problem can be solved offline, i.e, we’ll find the answers to all queries not necessarily in the order they’re given in the input (which is fine as long as they’re printed in the right order).
In our case, we’ll answer queries in descending order of L.

Let \text{cost}[x] denote, well, the cost of x.
For a fixed L, suppose we’ve answered all queries with left endpoints \gt L.
What happens when L is allowed to be a left endpoint?
It can be observed that:

  • If x is not a multiple of L, the cost is unchanged.
  • If x is a multiple of L, the cost may change.
    • Specifically, you might want to consider some range of factors of x starting from L whose LCM is x, and see if it ends somewhere less than the current \text{cost}[x].
  • Once all the \text{cost}[x] values of multiples of L are updated appropriately, the answer to a query [L, R] is just
    \sum_{x=L}^R \text{cost}[x]

To compute these range sums quickly, notice that you need a data structure that supports point updates (changing some specific \text{cost}[x] value), and range sums.
This is easily done with a segment tree/fenwick tree in \mathcal{O}(\log N) time.
You require \mathcal{O}(N\log N) updates in total (\text{cost}[x] gets updated once for each of its divisors), and one range query for each [L, R], for a total of \mathcal{O}(N\log^2 N) (N = 3\cdot 10^5 here).

The only piece remaining is the actual recomputation of \text{cost}[x] values for multiples of L, i.e, finding the shortest range of divisors of x starting from L whose LCM is x.
This can be done with, for example, two pointers (remembering the count of each maximum prime power in x, and updating as the endpoints move) for \mathcal{O}(N\log^2 N) in total again; but there’s a much simpler method that happens to run fast here: brute force!

That is, simply start from L and keep taking the LCM till you attain x.
The complexity of this is bounded above by the sum of the squares of the divisor function, which is \mathcal{O}(N\log^3 N) (see OEIS), but in practice it happens to run much faster because you likely won’t have to iterate long before you reach the required LCM (you can verify by running locally that the total number of LCM calls is about 1.6\cdot 10^7).

TIME COMPLEXITY:

\mathcal{O}(N\log^2 N + T\log N), where N = 3\cdot 10^5 here.

CODE:

Author's code (C++)
#include <bits/stdc++.h>
using namespace std;
// #include<ext/pb_ds/assoc_container.hpp>
// #include<ext/pb_ds/tree_policy.hpp>
// using namespace __gnu_pbds;
// template<class T> using oset =tree<T, null_type, less<T>, rb_tree_tag,tree_order_statistics_node_update> ;//oset name name.order_of_key(#ele<k) *name.find_by_order(index) less_equal greater greater_equal

#define vi vector<int>
#define pii pair<int, int>
#define mii map<int, int>
#define int long long
#define ld long double
#define pb push_back
#define all(v) v.begin(), v.end()
#define yes cout << "YES" \
                 << "\n";
#define no cout << "NO" \
                << "\n";
#define nl "\n"
#define FastIO                    \
    ios_base::sync_with_stdio(0); \
    cin.tie(0);                   \
    cout.tie(0)
#define mod 1000000007
const int oo = 1e18;
vi cost(1000001, 0);
vector<vi> factor(1000001);
vector<vi> store;

//-----------------------------------SEGMENT TREE----------------------------------------//
class SegmentTree
{
private:
    vector<int> tree;
    int n;

    // Function to build the segment tree
    void build(const vector<int> &arr, int v, int tl, int tr)
    {
        if (tl == tr)
        {
            tree[v] = arr[tl];
        }
        else
        {
            int tm = (tl + tr) / 2;
            build(arr, v * 2, tl, tm);
            build(arr, v * 2 + 1, tm + 1, tr);
            tree[v] = tree[v * 2] + tree[v * 2 + 1]; // Modify this line based on the query
        }
    }

public:
    SegmentTree(const vector<int> &arr)
    {
        n = arr.size();
        tree.resize(4 * n); // Adjust the size based on the maximum size of your input array
        build(arr, 1, 0, n - 1);
    }

    // Function to update a value at index idx to val
    void update(int idx, int val)
    {
        update(1, 0, n - 1, idx, val);
    }

    // Function to query the range [l, r]
    int query(int l, int r)
    {
        return query(1, 0, n - 1, l, r);
    }

private:
    // Function to update a value at index idx to val in the segment tree
    void update(int v, int tl, int tr, int idx, int val)
    {
        if (tl == tr)
        {
            tree[v] = val;
        }
        else
        {
            int tm = (tl + tr) / 2;
            if (idx <= tm)
            {
                update(v * 2, tl, tm, idx, val);
            }
            else
            {
                update(v * 2 + 1, tm + 1, tr, idx, val);
            }
            tree[v] = tree[v * 2] + tree[v * 2 + 1]; // Modify this line based on the query
        }
    }

    // Function to query the range [l, r] in the segment tree
    int query(int v, int tl, int tr, int l, int r)
    {
        if (l > r)
        {
            return 0; // Modify this line based on the query
        }
        if (l == tl && r == tr)
        {
            return tree[v];
        }
        int tm = (tl + tr) / 2;
        return query(v * 2, tl, tm, l, min(r, tm)) + query(v * 2 + 1, tm + 1, tr, max(l, tm + 1), r);
        // Modify the line above based on the query
    }
};

void solve()
{
    vi ans(store.size());
    int j = 0;
    SegmentTree sgt(cost);

    for (int i = 3e5; i >= 1; i--)
    {
        cost[i] = i;
        sgt.update(i, cost[i]);
        factor[i].pb(i);
        for (int j = 2 * i; j <= 3e5; j += i)
        {
            factor[j].pb(i);
            int ans = j;
            int lc = 1;
            for(int ii=factor[j].size()-1;ii>=0;ii--){
                int it=factor[j][ii];
                lc=lcm(lc,it);
                if (lc == j)
                {
                    cost[j] = it;
                    sgt.update(j, cost[j]);
                    break;
                }

            }
        }
        while (j < store.size() && store[j][0] == i)
        {
            ans[store[j][2]] = sgt.query(store[j][0], store[j][1]);
            // cerr<<i<<nl;
            // for(int jjj=store[j][0];jjj<=store[j][1];jjj++){
            //     cout<<cost[jjj]<<nl;
            // }
            j++;
        }
    }
    for (auto it : ans)
    {
        cout << it << nl;
    }
}
signed main()
{
    FastIO;
    // freopen("P10_7.in", "r", stdin);
    // freopen("P10_7.out", "w", stdout);
    int test = 1;
    cin >> test;
    for (int i = 0; i < test; i++)
    {
        int l, r;
        cin >> l >> r;
        store.pb({l, r, i});
    }
    sort(all(store));
    reverse(all(store));
    solve();
    return 0;
}


Tester's code (C++)
/*

*       *  *  ***       *       *****
 *     *   *  *  *     * *        *
  *   *    *  ***     *****       *
   * *     *  * *    *     *   *  *
    *      *  *  *  *       *   **

                                 *
                                * *
                               *****
                              *     *
        *****                *       *
      _*     *_
     | * * * * |                ***
     |_*  _  *_|               *   *
       *     *                 *  
        *****                  *  **
       *     *                  ***
  {===*       *===}
      *  IS   *                 ***
      *  IT   *                *   *
      * RATED?*                *  
      *       *                *  **
      *       *                 ***
       *     *
        *****                  *   *
                               *   *
                               *   *
                               *   *
                                ***   

*/

// 2 Years Tribute to Competitive Programming

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace __gnu_pbds;
using namespace std;

#define osl tree<ll, null_type, less<ll>, rb_tree_tag, tree_order_statistics_node_update>
#define ll long long
#define ld long double
#define forl(i, a, b) for(ll i = a; i < b; i++)
#define rofl(i, a, b) for(ll i = a; i > b; i--)
#define fors(i, a, b, c) for(ll i = a; i < b; i += c)
#define fora(x, v) for(auto x : v)
#define vl vector<ll>
#define vb vector<bool>
#define pub push_back
#define pob pop_back
#define fbo find_by_order
#define ook order_of_key
#define yesno(x) cout << ((x) ? "YES" : "NO")
#define all(v) v.begin(), v.end()

const ll N = 3e5 + 4;
const ll mod = 1e9 + 7;
// const ll mod = 998244353;

vl factor[N];
struct FenwickTree {
    vector<ll> bit;  // binary indexed tree
    ll n;

    FenwickTree(ll n) {
        this->n = n;
        bit.assign(n, 0);
    }

    FenwickTree(vector<ll> const &a) : FenwickTree(a.size()) {
        for (size_t i = 0; i < a.size(); i++)
            add(i, a[i]);
    }

    ll sum(ll r) {
        ll ret = 0;
        for (; r >= 0; r = (r & (r + 1)) - 1)
            ret += bit[r];
        return ret;
    }

    ll sum(ll l, ll r) {
        return sum(r) - sum(l - 1);
    }

    void add(ll idx, ll delta) {
        for (; idx < n; idx = idx | (idx + 1))
            bit[idx] += delta;
    }
};

ll modinverse(ll a) {
	ll m = mod, y = 0, x = 1;
	while (a > 1) {
		ll q = a / m;
		ll t = m;
		m = a % m;
		a = t;
		t = y;
		y = x - q * y;
		x = t;
	}
	if (x < 0) x += mod;
	return x;
}
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;
}
bool poweroftwo(ll n) {
	return !(n & (n - 1));
}
ll power(ll a, ll b, ll md = mod) {
	ll product = 1;
	a %= md;
	while (b) {
		if (b & 1) product = (product * a) % md;
		a = (a * a) % md;
		b /= 2;
	}
	return product % md;
}
void panipuri() {
	ll n, m = 0, k = -1, c = 0, sum = 0, q = 0, ans = 0, p = 1;
	string s;
	bool ch = true;
	cin >> n;
	forl(i,1,N){
		fors(j,i,N,i) factor[j].pub(i);
	}
	vl curr(N),curr1(N);
	vector<vl> v;
	forl(i,0,n){
		ll l,r;
		cin>>l>>r;
		v.pub({l,r,i});
	}
	sort(all(v));
	FenwickTree ft(N);
	forl(i,1,N){
		p=1;
		forl(j,0,factor[i].size()){
			p=lcm(p,factor[i][j]);
			if(p==i){
				curr1[i]=factor[i][j];
				ft.add(i,curr1[i]);
				break;
			}
		}
	}
	vl v1(n);
	c=1;
	forl(i,0,n){
		while(c<v[i][0]){
			ft.add(c,-c);
			fors(j,c*2,N,c){
				ft.add(j,-curr1[j]);
				curr[j]++;
				p=1;
				forl(i1,curr[j],factor[j].size()){
					p=lcm(p,factor[j][i1]);
					if(p==j){
						curr1[j]=factor[j][i1];
						ft.add(j,curr1[j]);
						break;
					}
				}
			}
			c++;
		}
		v1[v[i][2]]=ft.sum(v[i][1]);
	}
	fora(x,v1) cout<<x<<'\n';
	return;
}
int main() {
	ios::sync_with_stdio(false);
	cin.tie(NULL);
	#ifndef ONLINE_JUDGE
	freopen("input.txt", "r", stdin);
	freopen("output.txt", "w", stdout);
	#endif
	int laddu = 1;
	// cin >> laddu;
	forl(i, 1, laddu + 1) {
		// cout << "Case #" << i << ": ";
		panipuri();
		cout << '\n';
	}
}
Editorialist's code (C++, subtask 1)
#include <bits/stdc++.h>
using namespace std;
using ll = long long int;

int main() {
    const int N = 3e5 + 10;
    vector<long long> ans(N);
    for (int i = 1; i < N; ++i) {
        int L = 1;
        ans[i] = ans[i-1];
        vector<int> facs;
        for (int j = 1; j*j <= i; ++j) {
            if (i%j) continue;
            facs.push_back(j);
            if (j*j != i) facs.push_back(i/j);
        }
        sort(begin(facs), end(facs));
        for (int d : facs) {
            L = lcm(L, d);
            if (L == i) {
                ans[i] += d;
                break;
            }
        }
    }
    
    ios::sync_with_stdio(false);
    cin.tie(0);
	int t; cin >> t;
	while (t--) {
	    int L, R; cin >> L >> R;
	    cout << ans[R] << '\n';
	}
}

Editorialist's code (C++, full)
// #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());

template<class T, T unit = T()>
struct SegTree {
    T f(T a, T b) { return a+b; }
    vector<T> s; int n;
    SegTree(int _n = 0, T def = unit) : s(2*_n, def), n(_n) {}
    void update(int pos, T val) {
        for (s[pos += n] = val; pos /= 2;)
            s[pos] = f(s[pos * 2], s[pos * 2 + 1]);
    }
    T query(int b, int e) {
        T ra = unit, rb = unit;
        for (b += n, e += n; b < e; b /= 2, e /= 2) {
            if (b % 2) ra = f(ra, s[b++]);
            if (e % 2) rb = f(s[--e], rb);
        }
        return f(ra, rb);
    }
};

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

    const int N = 3e5 + 10;
    int t; cin >> t;
    vector<vector<array<int, 2>>> queries(N);
    for (int i = 0; i < t; ++i) {
        int L, R; cin >> L >> R;
        queries[L].push_back({R, i});
    }

    SegTree<ll> seg(N);
    vector<ll> ans(t);

    vector<vector<int>> facs(N);
    for (int i = 1; i < N; ++i)
        for (int j = i; j < N; j += i)
            facs[j].push_back(i);
    vector<int> pos(N);
    for (int i = 1; i < N; ++i) pos[i] = facs[i].size();
    
    for (int L = N-1; L >= 1; --L) {
        for (int x = L; x < N; x += L) {
            int l = 1;
            --pos[x];
            for (int i = pos[x]; ; ++i) {
                l = lcm(l, facs[x][i]);
                if (l == x) {
                    seg.update(x, facs[x][i]);
                    break;
                }
            }
        }
        for (auto [R, id] : queries[L]) {
            ans[id] = seg.query(0, R+1);
        }
    }
    for (auto x : ans) cout << x << '\n';
}
1 Like

Let’s say the last greater L in query which we calculated to be 9. If x = 24, the cost[x] is also equal to 24, (for L = 9). Let the next decreasing L in the query be 7. Now 24 is not divisible by 7 so according to this statement

the cost should remain same, but that is not the case, the cost[24] now becomes 8.

What am I missing?

You don’t iterate only across only those L that appear in queries, you consider all of them - if there happen to be no queries with that L, doesn’t matter, it doesn’t worsen the complexity anyway (if your implementation really depends on every L being queried, just add fake queries of the form (L, L) for every L and don’t print out their answers at the end).

In the example you gave, you’d also pass through 8 when moving from 9 to 7, so cost[24] will be its correct value (which is 12 and not 8 for L = 7, btw).