↵
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$ [he](https://lqdoj.edu.vn/problem/square)[re](https://oj.vnoi.info/problem/dhbb2020_square)↵
↵
-------------------------↵
↵
## Solve for k = 1↵
↵
The answer just simply be $n$↵
↵
-------------------------↵
↵
## Solve for k = 2↵
↵
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$↵
↵
Therefore the answer just simply be↵
↵
<spoiler summary="Implementation using factorization">↵
↵
```cpp↵
vector<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;↵
}↵
```↵
</spoiler>↵
↵
<spoiler summary="Implementation 1">↵
↵
```cpp↵
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 += 1LL * (j - 1) * (j - 2) / 2;↵
}↵
↵
res %= MOD;↵
return res;↵
}↵
```↵
</spoiler>↵
↵
<spoiler summary="Implementation 2">↵
↵
```cpp↵
int 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;↵
}↵
```↵
</spoiler>↵
↵
So about the complexity....↵
↵
For the implementation using factorization. It is ofcourse $O(n \log \log n)$, you can easily proove it. As it just related to [the sum of inversion of primes](https://en.wikipedia.org/wiki/Divergence_of_the_sum_of_the_reciprocals_of_the_primes)↵
↵
For the 2 implementations below
↵
$O\left(\underset{\text{squarefree } p \leq n}{\Large \Sigma} \left \lfloor \frac{n}{p^2} \right \rfloor \right) = O\left(\underset{p \leq n}{\Large \Sigma} \frac{n}{p^2} \right) = O\left(n \times {\Large \Sigma}
↵
-------------------------↵
↵
## Solve for general k↵
↵
Using the same logic above, we can easily solve the problem. But we can count with combinatorics↵
↵
<spoiler summary="O(n) solution">↵
↵
```cpp↵
const 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;↵
}↵
```↵
</spoiler>↵
↵
-------------------------↵
↵
## A better solution for k = 2↵
↵
In this approach, instead of fixing $u$ as a squarefree and count $p^2$. We do the reverse, let count the number of way to select $u$ as we fix $p^2$.↵
↵
Normaly, it will still lead you to $O(n)$ solution:↵
↵
<spoiler summary="Swap for loop implementation">↵
↵
```cpp↵
int 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;↵
}↵
```↵
</spoiler>↵
↵
-------------------------↵
↵
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↵
↵
Why do we need this function, well. You see since $(a, b, n)$ form a three-term geometric progression. Therefore we have $b^2 = an$.↵
↵
Let $F(N) = \overset{n}{\underset{p=1}{\Large \Sigma}} f(p)$↵
↵
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$↵
↵
-------------------------↵
↵
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↵
↵
It is no hard to proove that $g(n) = f(n) - 1$↵
↵
But this interesting sequence $g(n)$ is [A000188](https://oeis.org/A000188)↵
↵
This sequence have many property, 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$.↵
- bla bla bla↵
↵
Well, actually to make the problem whole easier, I gonna skip all the proofs (still, you can use the link in the [sequence](https://oeis.org/A000188) to prove). And use this property↵
↵
$g(n) = \underset{d^2 | n}{\Large \Sigma} \phi(d)$↵
↵
Then we have↵
↵
$F(N)$↵
$= \overset{n}{\underset{p=1}{\Large \Sigma}} f(p)$↵
$= \overset{n}{\underset{p=2}{\Large \Sigma}} g(p)$↵
$= \overset{n}{\underset{p=2}{\Large \Sigma}} \underset{d^2 | p}{\Large \Sigma} \phi(d)$↵
$= \overset{\left \lfloor \sqrt{n} \right \rfloor}{\underset{p=2}{\Large \Sigma}} \underset{1 \leq u \times p^2 \leq n}{\Large \Sigma} \phi(p)$↵
$= \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](https://www.austms.org.au/wp-content/uploads/Gazette/2008/Jul08/TechPaperMyerson.pdf) also takes you to something similar.↵
↵
<spoiler summary="O(sqrt n log log sqrt n) solution">↵
↵
```cpp↵
#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;↵
}↵
```↵
</spoiler>↵
↵
<spoiler summary="O(sqrt) solution">↵
↵
```cpp↵
#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;↵
}↵
```↵
</spoiler>↵
↵
↵
-------------------------↵
↵
## My question↵
↵
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})$ ?↵
↵
↵
↵
↵
↵
↵
H: Can I solve if it is given $n$ and queries for $k$ ?↵
↵
I: Can I solve if it is given $k$ and queries for $n$ ?
