665F - Four Divisors

The editorial for this problem is a little modification of the materials from the lecture of Mikhail Tikhomirov Endagorion of the autumn of 2015 in Moscow Institute of Physics and Technology. Thanks a lot to Endagorion for that materials.

Easy to see that only the numbers of the form p·q and p³ (for different prime p, q) have exactly four positive divisors.

We can easily count the numbers of the form p³ in , where n is the number from the problem statement.

Now let p < q and π(k) be the number of primes from 1 to k. Let's iterate over all the values p. Easy to see that . So for fixed p we should increase the answer by the value .

So the task is ot to find ~--- the number of primes not exceeding , for all p.

Denote p_j the j-th prime number. Denote dp_n, j the number of k such that 1 ≤ k ≤ n, and all prime divisors of k are at least p_j (note that 1 is counted in all dp_n, j, since the set of its prime divisors is empty). dp_n, j satisfy a simple recurrence:

1. dp_n, 1 = n (since p₁ = 2)

2. dp_n, j = dp_{n, j + 1} + dp_{⌊ n / pj⌋, j}, hence dp_{n, j + 1} = dp_n, j - dp_{⌊ n / pj⌋, j}

Let p_k be the smallest prime greater than . Then π(n) = dp_n, k + k - 1 (by definition, the first summand accounts for all the primes not less than k).

If we evaluate the recurrence dp_n, k straightforwardly, all the reachable states will be of the form dp_{⌊ n / i⌋, j}. We can also note that if p_j and p_k are both greater than , then dp_n, j + j = dp_n, k + k. Thus, for each ⌊ n / i⌋ it makes sense to keep only values of dp_{⌊ n / i⌋, j}.

Instead of evaluating all DP states straightforwardly, we perform a two-step process:

1. Choose K.

2. Run recursive evaluation of dp_n, k. If we want to compute a state with n < K, memorize the query ``count the numbers not exceeding n with all prime divisors at least k''.

3. Answer all the queries off-line: compute the sieve for numbers up to K, then sort all numbers by the smallest prime divisor. Now all queries can be answered using RSQ structure. Store all the answers globally.

4. Run recurisive evaluation of dp_n, k yet again. If we want to compute a state with n < K, then we must have preprocessed a query for this state, so take it from the global set of answers.

The performance of this approach relies heavily on Q~--- the number of queries we have to preprocess.

Statement. .

Proof:

Each state we have to preprocess is obtained by following a dp_{⌊ n / pj⌋, j} transition from some greater state. It follows that Q doesn't exceed the total number of states for n > K.

[Q \leq \sum_{j = 1}^{n / K} \pi(\sqrt{n / j}) \sim \sum_{j = 1}^{n / K} \sqrt{n / j} / \log n \sim \frac{1}{\log n}\int_1^{n / K} \sqrt{n / x} dx \sim \frac{n}{\sqrt{K}\log n}]

The preprocessing of Q queries can be done in , and it is the heaviest part of the computation. Choosing optimal , we obtain the complexity .

li n;

bool read() {
	return !!(cin >> n);
}

const int K = 10 * 1000 * 1000;
const int N = K;
const int P = 700100;
const int Q = 25 * 1000 * 1000;

int szp, p[P];
int mind[N];

void prepare() {
	szp = 0;
	forn(i, N) mind[i] = -1;
	fore(i, 2, N) {
		if (mind[i] == -1) {
			assert(szp < P);
			mind[i] = szp;
			p[szp++] = i;
		}
		for (int j = 0; j < szp && j <= mind[i] && i * p[j] < N; j++)
			mind[i * p[j]] = j;
	}
}

inline int getk(li n) {
	int lf = 0, rg = szp - 1;
	while (lf != rg) {
		int md = (lf + rg) >> 1;
		if (p[md] * 1ll * p[md] > n) rg = md;
		else lf = md + 1;
	}
	assert(p[lf] * 1ll * p[lf] > n);
	return lf;
}

int t[K];
void inc(int i, int val) {
	for ( ; i < K; i |= i + 1)
		t[i] += val;
}
int sum(int i) {
	int ans = 0;
	for ( ; i >= 0; i = (i & (i + 1)) - 1)
		ans += t[i];
	return ans;
}

int szq;
pti q[Q];
int ans[Q];
vector<int> vs[P];

void process() {
	sort(q, q + szq, greater<pti> ());
	memset(t, 0, sizeof(t));

	forn(i, szp) vs[i].clear();
	fore(i, 2, K)
		vs[mind[i]].pb(i);

	inc(1, +1);

	int p = szp - 1;
	for (int i = 0, j = 0; i < szq; i = j) {
		while (p >= q[i].x) {
			for (auto v : vs[p])
				inc(v, +1);
			p--;
		}

		while (j < szq && q[j].x == q[i].x) {
			ans[j] = sum(q[j].y);
			j++;
		}
	}
}

map<pair<li, int>, li> z;

li solve(li n, int jj, bool fs) {
	if (!n) return 0;
	int j = min(jj, getk(n));
	if (!j) return n + j - jj;

	li ans = 0;
	if (n < K) {
		pti p(j, (int) n);
		if (fs) {
			assert(szq < Q);
			q[szq++] = p;
			ans = 0;
			return 0;
		} else {
			int idx = int(lower_bound(q, q + szq, p, greater<pti> ()) - q);
			assert(idx < szq && q[idx] == p);
			ans = ::ans[idx];
		}
	} else {
		if (!z.count(mp(n, j))) {
			ans = solve(n, j - 1, fs);
			ans -= solve(n / p[j - 1], j - 1, fs);
			z[mp(n, j)] = ans;
		} else {
			ans = z[mp(n, j)];
		}
	}

	ans += j - jj;

	return ans;
}

inline li pi(li n, bool fs) {
	int k = szp - 1;
	return solve(n, k, fs) + k - 1;
}

void solve() {
	szq = 0;
	z.clear();
	for (int j = 0; p[j] * 1ll * p[j] <= n; j++) {
		li nn = n / p[j];
		if (nn > p[j]) {
			pi(n / p[j], true);
		}
	}

	process();

	z.clear();
	li ans = 0;
	for (int j = 0; p[j] * 1ll * p[j] <= n; j++) {
		li nn = n / p[j];
		if (nn > p[j]) {
			ans += pi(n / p[j], false);
			ans -= j + 1;
		}
	}

	for (int i = 0; i < szp && p[i] * 1ll * p[i] * 1ll * p[i] <= n; i++)
		ans++;

	cout << ans << endl;
}

Complexity: .