The statement:
Given two integer $$$n, k, p$$$, $$$(1 \leq k \leq n < p)$$$.
Count the number of array $$$a[]$$$ of size $$$k$$$ that satisfied
- $$$1 \leq a_1 < a_2 < \dots < a_k \leq n$$$
- $$$a_i \times a_j$$$ is perfect square $$$\forall 1 \leq i < j \leq n$$$
Since the number can be big, output it under modulo $$$p$$$.
Yet you can submit the problem for $$$k = 3$$$ here.
Solution for k = 1
The answer just simply be $$$n$$$
Solution for k = 2
Algorithm
We need to count the number of pair $$$(a, b)$$$ that $$$1 \leq a < b \leq n$$$ and $$$a \times b$$$ is perfect square.
Every positive integer $$$x$$$ can be represent uniquely as $$$x = u \times p^2$$$ for some positive integer $$$u, p$$$ and $$$u$$$ as small as possible ($$$u$$$ is squarefree number).
Let represent $$$x = u \times p^2$$$ and $$$y = v \times q^2$$$ (still, minimum $$$u$$$, $$$v$$$ ofcourse).
We can easily proove that $$$x \times y$$$ is a perfect square if and if only $$$u = v$$$.
So for a fixed squarefree number $$$u$$$. You just need to count the number of ways to choose $$$p^2$$$.
The answer will be the sum of such ways for each fixed $$$u$$$.
Implementation
Implementation using factorizationvector<int> prime; /// prime list
vector<int> lpf; /// Lowest prime factor, lpf[x] is smallest prime divisor of x
void sieve(int lim = LIM) /// O(n)
{
prime.assign(1, 2);
lpf.assign(lim + 1, 2);
lpf[1] = 1; /// For easier calculation but can cause inf loops
for (int i = 3; i <= lim; i += 2) {
if (lpf[i] == 2) prime.push_back(lpf[i] = i);
for (int j = 0; j < sz(prime) && prime[j] <= lpf[i] && prime[j] * i <= lim; ++j)
lpf[prime[j] * i] = prime[j];
}
}
/// mask(x) is smallest positive number that mask(x) * x is a perfect square
int getMask(int x) /// O(log n)
{
int mask = 1;
while (x > 1) {
int p = lpf[x], f = 0;
do x /= p, f++; while (p == lpf[x]);
if (f & 1) mask *= p; /// if current power is odd, we mutiple mask with current prime
}
return mask;
}
int cnt[LIM];
int magic(int n) /// O(n log max(a))
{
memset(cnt, 0, sizeof(cnt[0]) * (n + 1));
ll res = 0;
for (int a = 1; a <= n; ++a) /// Check all cases of a
res += cnt[getMask(a)]++;
res %= MOD;
return res;
}
int main() /// O(n log max(a))
{
int n;
cin >> n;
sieve(n + 500);
cout << magic(n);
return 0;
}
Implementation 1int solve(int n)
{
memset(is_squarefree, true, sizeof(is_squarefree[0]) * (n + 1));
long long res = 0;
for (int i = 1, j; i <= n; ++i) if (is_squarefree[i])
{
for (j = 1; i * j * j <= n; ++j)
is_squarefree[i * j * j] = false;
res += 1LL * (j - 1) * (j - 2) / 2;
}
res %= MOD;
return res;
}
Implementation 2int solve(int n)
{
memset(is_squarefree, true, sizeof(is_squarefree[0]) * (n + 1));
long long res = 0;
for (int i = 1; i <= n; ++i) if (is_squarefree[i])
{
for (int j = 1; i * j * j <= n; ++j)
{
is_squarefree[i * j * j] = false;
res += j - 1;
}
}
res %= MOD;
return res;
}
Complexity
So about the complexity....
For the implementation using factorization, it is $$$O(n \log n)$$$.
Hint 2Precalculating smallest prime factor for each $$$n$$$.
ProofSince you known the smallest prime factor for each $$$n$$$, you can factorize any number in $$$O(\log n)$$$.
There are $$$n$$$ number total. Hence the calculating part take $$$O(n \log n)$$$.
The sieve only cost you $$$O(n)$$$ give the total complexity of $$$O(n \log n)$$$.
BonusThis just a clickbait. Have a nice day.
Real Bonus The true BonusYou can optimize it to solve in $$$O(n \log \log n)$$$ though it is not that necessary
For the 2 implementations below, the complexity is linear.
Hint 1Each number will be calculated once.
Hint 2Though some numbers might be visited many times. It is actually armotized to constant.
Hint 3Let remind you about something called the sum of inversion of squares.
Hint 4Just ignore those squarefree. A normal loop through it is still linear.
Proof$$$O\left(\underset{\text{squarefree } p \leq n}{\Large \Sigma} \left \lfloor {\Large \frac{n}{p^2}} \right \rfloor \right) = O\left(\underset{p \leq n}{\Large \Sigma} {\Large \frac{n}{p^2}} \right) = O\left(n \times {\Large \Sigma} {\Large \frac{1}{p^2}} \right) = O\left(n \times {\Large \frac{\pi^2}{6}} \right) = O(n)$$$.
Solution for general k
Using the same logic above, we can easily solve the problem.
But for a fixed $$$u$$$, we can improve the counting using combinatorics.
Yet nCk is a well known problem therefore I wont describe it here.
O(n) solutionconst int LIM = 1e7 + 17;
const int MOD = 1e9 + 7;
int fact[LIM + 1]; /// factorial: fact[n] = n!
int invs[LIM + 1]; /// inverse modular: invs[n] = n^(-1)
int tcaf[LIM + 1]; /// inverse factorial: tcaf[n] = (n!)^(-1)
void precal(int n = LIM)
{
fact[0] = fact[1] = 1;
invs[0] = invs[1] = 1;
tcaf[0] = tcaf[1] = 1;
for (int i = 2; i <= n; ++i)
{
fact[i] = (1LL * fact[i - 1] * i) % MOD;
invs[i] = MOD - 1LL * (MOD / i) * invs[MOD % i] % MOD;
tcaf[i] = (1LL * tcaf[i - 1] * invs[i]) % MOD;
}
}
int nCk(int n, int k)
{
k = min(k, n - k);
if (k < 0) return 0;
long long res = fact[n];
res *= tcaf[k]; res %= MOD;
res *= tcaf[n - k]; res %= MOD;
return res;
}
int solve(int n)
{
memset(is_squarefree, true, sizeof(is_squarefree[0]) * (n + 1));
long long res = 0;
for (int i = 1, j; i <= n; ++i) if (is_squarefree[i])
{
for (j = 1; i * j * j <= n; ++j)
is_squarefree[i * j * j] = false;
res += nCk(j - 1, k);
}
res %= MOD;
return res;
}
A better solution for k = 2
Idea
In the above approach, we fix $$$u$$$ as a squarefree and count $$$p^2$$$.
But what if I fix $$$p^2$$$ to count $$$u$$$ instead ?
Yet you can see that the first loop now is $$$O(\sqrt{n})$$$, but it will still $$$O(n)$$$ total because of the second loop
Swap for loop implementationint solve(int n)
{
memset(is_squarefree, true, sizeof(is_squarefree[0]) * (n + 1));
long long res = 0;
int t = sqrt(n);
while (t * t < n) ++t;
while (t * t > n) --t;
for (int j = t; j > 1; --j)
{
for (int i = 1; i * j * j <= n; ++i)
{
if (!used[i * j * j])
{
used[i * j * j] = true;
res += j - 1;
}
}
}
res %= MOD;
return res;
}
Approach
Let $$$f(n)$$$ is the number of $$$(a, b)$$$ that $$$1 \leq a < b \leq n$$$ and $$$(a, b, n)$$$ is a three-term geometric progression.
Let $$$g(n)$$$ is the number of $$$(a, b)$$$ that $$$1 \leq a \leq b \leq n$$$ and $$$(a, b, n)$$$ is a three-term geometric progression.
Let $$$F(N) = \overset{n}{\underset{p=1}{\Large \Sigma}} f(p)$$$.
But why do we need these functions anywayWell, you see, since $$$(a, b, n)$$$ form a three-term geometric progression then you have $$$b^2 = an$$$.
It is not hard to proove that $$$F(N)$$$ will be our answer as we count over all possible squarefree $$$u$$$ for every fixed $$$p^2$$$.
We use $$$g(n)$$$ because its beautiful property, which make the calculation a lot easier.
So it is no hard to prove that $$$g(n) = f(n) + 1$$$.
This interesting sequence $$$g(n)$$$ is A000188, having many properties, such as
- Number of solutions to $$$x^2 \equiv 0 \pmod n$$$.
- Square root of largest square dividing $$$n$$$.
- Max $$$gcd \left(d, \frac{n}{d}\right)$$$ for all divisor $$$d$$$.
Well, to make the problem whole easier, I gonna skip all the proofs to use this property (still, you can use the link in the sequence for references).
$$$g(n) = \underset{d^2 | n}{\Large \Sigma} \phi(d)$$$.
From this property, we can solve the problem in $$$O(\sqrt{n})$$$.
Hint 1$$$F(N) = \overset{n}{\underset{p=1}{\Large \Sigma}} f(p)$$$.
and
$$$f(n) = g(n) + 1$$$.
Hint 2We have
$$$F(n)$$$
$$$= \overset{n}{\underset{p=2}{\Large \Sigma}} g(p) - n$$$
$$$= \overset{n}{\underset{p=2}{\Large \Sigma}} \underset{d^2 | p}{\Large \Sigma} \phi(d) - n$$$
Hint 3Break the loop, make it simpler.
Hint 4Reduce a for loop to constant.
Solution$$$F(N)$$$
$$$= \overset{n}{\underset{p=1}{\Large \Sigma}} f(p)$$$
$$$= -n + \overset{n}{\underset{p=2}{\Large \Sigma}} g(p)$$$
$$$= -n + \overset{n}{\underset{p=2}{\Large \Sigma}} \underset{d^2 | p}{\Large \Sigma} \phi(d)$$$
$$$= -n + \overset{\left \lfloor \sqrt{n} \right \rfloor}{\underset{p=1}{\Large \Sigma}} \underset{1 \leq u \times p^2 \leq n}{\Large \Sigma} \phi(p)$$$
$$$= -n + \overset{\left \lfloor \sqrt{n} \right \rfloor}{\underset{p=1}{\Large \Sigma}} \phi(p) \times \left \lfloor \frac{n}{p^2} \right \rfloor$$$
$$$= \overset{\left \lfloor \sqrt{n} \right \rfloor}{\underset{p=2}{\Large \Sigma}} \phi(p) \times \left \lfloor \frac{n}{p^2} \right \rfloor$$$
Well, you can sieve $$$\phi(p)$$$ for all $$$2 \leq p \leq \sqrt{n}$$$ in $$$O\left(\sqrt{n} \log \log \sqrt{n} \right)$$$ or improve it with linear sieve to $$$O(\sqrt{n})$$$.
Therefore you can solve the whole problem in $$$O(\sqrt{n})$$$.
Yet this paper also takes you to something similar.
Implementation
O(sqrt n log log sqrt n) solution#include <iostream>
#include <cstring>
#include <numeric>
#include <cmath>
using namespace std;
const int MOD = 1e9 + 7;
const int LIM = 1e7 + 17;
const int SQRT_LIM = ceil(sqrt(LIM) + 1) + 1;
int euler[SQRT_LIM];
void sieve_phi(int n)
{
iota(euler, euler + n + 1, 0);
for (int x = 2; x <= n; x++) if (euler[x] == x)
for (int j = x; j <= n; j += x)
euler[j] -= euler[j] / x;
}
int solve(int n)
{
sieve_phi(ceil(sqrt(n) + 1) + 1);
long long res = 0;
for (int p = 2; p * p <= n; ++p)
res += 1LL * euler[p] * (n / (p * p));
res %= MOD;
return res;
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(NULL);
int n;
cin >> n;
cout << solve(n);
return 0;
}
O(sqrt) solution#include <iostream>
#include <cstring>
#include <numeric>
#include <vector>
#include <cmath>
using namespace std;
const int MOD = 1e9 + 7;
vector<int> lpf;
vector<int> prime;
vector<int> euler;
void linear_sieve_phi(int n)
{
lpf.assign(n + 1, 0);
euler.assign(n + 1, 1);
for (int x = 2; x <= n; ++x)
{
if (lpf[x] == 0)
{
prime.push_back(lpf[x] = x);
euler[x] = x - 1;
}
for (int i = 0; i < prime.size() && x * prime[i] <= n; ++i)
{
lpf[x * prime[i]] = prime[i];
if (x % prime[i] == 0) {
euler[x * prime[i]] = euler[x] * prime[i];
break;
}
euler[x * prime[i]] = euler[x] * euler[prime[i]];
}
}
}
int solve(int n)
{
linear_sieve_phi(ceil(sqrt(n) + 1) + 1);
long long res = 0;
for (int p = 2; p * p <= n; ++p)
res += 1LL * euler[p] * (n / (p * p));
res %= MOD;
return res;
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(NULL);
int n;
cin >> n;
cout << solve(n);
return 0;
}
Miscellaneous
In case you still wondering why is that function $$$\phi(n)$$$ related here, then there is a implement using Möbius function in linear too
Spoiler/// Linear Sieve
vector<bool> isPrime; /// Characteristic function = A010051
vector<int> prime; /// Prime list = A000040
vector<int> lpf; /// Lowest prime factor = A020639
vector<int> mu; /// Mobius = A008683
vector<int> phi; /// Euler totient phi = A000010
void sieve(int n)
{
if (n < 1) return ;
/// Extension part
mu.assign(n + 1, 1);
phi.assign(n + 1, 1);
/// Main part
prime.clear();
lpf.assign(n + 1, 0);
isPrime.assign(n + 1, true);
isPrime[0] = isPrime[1] = false;
for (int x = 2; x <= n; ++x)
{
if (isPrime[x]) /// Func[Prime]
{
lpf[x] = x;
mu[x] = -1;
phi[x] = x - 1;
prime.push_back(x);
}
for (int p : prime) /// Func[Prime * X] <- Func[Prime]
{
if (p > lpf[x] || x * p > n) break;
isPrime[x * p] = 0;
lpf[x * p] = p;
if (lpf[x] == p)
{
mu[x * p] = 0;
phi[x * p] = phi[x] * p;
}
else
{
mu[x * p] = -mu[x];
phi[x * p] = phi[x] * phi[p];
}
}
}
}
long long solve(int n)
{
sieve(n + 1);
long long res = 0;
for (int x = 1; x <= n; ++x) if (mu[x])
{
int t = sqrt(n / x);
res += 1LL * t * (t - 1) / 2;
}
res %= MOD;
return res;
}
Tasks
A: Can we also use phi function or something similar to solve for $$$k = 3$$$ in $$$O(\sqrt{n})$$$ ?
B: Can we also use phi function or something similar to solve for general $$$k$$$ in $$$O(\sqrt{n})$$$ ?
C: Can we also solve the problem where there can be duplicate: $$$a_i \leq a_j (\forall\ i < j)$$$ and no longer $$$a_i < a_j (\forall\ i < j)$$$ ?
D: Can we solve the problem where there is no restriction between $$$k, n, p$$$ ?
E: Can we solve for negative integers, whereas $$$-n \leq a_1 < a_2 < \dots < a_k < n$$$
F: Can we solve for a specific range, whereas $$$L \leq a_1 < a_2 < \dots < a_k < R$$$
G: Can we solve for cube product $$$a_i \times a_j \times a_k$$$ effectively ?
H: Can we solve if it is given $$$n$$$ and queries for $$$k$$$ ?
I: Can we solve if it is given $$$k$$$ and queries for $$$n$$$ ?
Contribution
Yurushia for pointing out the linear complexity of squarefree sieve.
clyring for fixing typos, and approach for tasks A, B