The statement:
Given three integers $$$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 k$$$
Since the number can be big, output it under modulo $$$p$$$.
For convenient, you can assume $$$p$$$ is a large constant prime $$$10^9 + 7$$$
Notice that in this blog, we will solve in $$$O(n)$$$ for general $$$k$$$ from the original task
For harder variants you can see in this blog [Variants] An interesting counting problem related to square product
Solution for k = 1
The answer just simply be $$$n$$$
Solution for k = 2
Problem
Given $$$n (1 \leq n \leq 10^6)$$$, count the number of pair $$$(a, b)$$$ that satisfied
- $$$1 \leq a < b \leq n$$$
- $$$a \times b$$$ is a perfect square.
Since the number can be big, output it under modulo $$$10^9 + 7$$$.
You can submit here: https://lqdoj.edu.vn/problem/aicprtsp2.
Examples
Example 1Input 1:
1
Output 1:
0
Explanation 1:
There are no satisfied integer pair $$$(a, b)$$$ that $$$1 \leq a < b \leq 1$$$
Example 2Input 2:
10
Output 2:
4
Explanation 2:
There are $$$4$$$ satisfied pairs: {$$$1, 4$$$}, {$$$1, 9$$$}, {$$$2, 8$$$}, {$$$4, 9$$$}.
Example 3Input 3:
25
Output 3:
16
Explanation 3:
There are $$$16$$$ satisfied pairs: {$$$1, 4$$$}, {$$$1, 9$$$}, {$$$1, 16$$$}, {$$$1, 25$$$}, {$$$2, 8$$$}, {$$$2, 18$$$}, {$$$3, 12$$$}, {$$$4, 9$$$}, {$$$4, 16$$$}, {$$$4, 25$$$}, {$$$5, 20$$$}, {$$$6, 24$$$}, {$$$8, 18$$$}, {$$$9, 16$$$}, {$$$9, 25$$$}, {$$$16, 25$$$}.
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$$$.
Illustration
Squarefree Sieve 10x10 Table | Iteration 0 Squarefree Sieve 10x10 Table | Iteration 1 Squarefree Sieve 10x10 Table | Iteration 2 Squarefree Sieve 10x10 Table | Iteration 3 Squarefree Sieve 10x10 Table | Iteration 5 Squarefree Sieve 10x10 Table | Iteration 6 Squarefree Sieve 10x10 Table | Iteration 7 Squarefree Sieve 10x10 Table | Iteration 10 Squarefree Sieve 10x10 Table | Iteration 11 Squarefree Sieve 10x10 Table | Iteration 13 Squarefree Sieve 10x10 Table | Iteration 14 Squarefree Sieve 10x10 Table | Iteration 15 Squarefree Sieve 10x10 Table | Iteration 17 Squarefree Sieve 10x10 Table | Iteration 19 Squarefree Sieve 10x10 Table | Iteration 21 Squarefree Sieve 10x10 Table | Iteration 22 Squarefree Sieve 10x10 Table | Iteration 23 Squarefree Sieve 10x10 Table | Iteration 100 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;
}
Implementation related to Möbius function/// 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;
}
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)$$$.
For the last implementation, the complexity is Linear
Hint 1It depends on the sieve as the calculation is already linear.
Hint 2Each number can be represent as $$$x = lpf[x] \times k$$$ where $$$lpf[x]$$$ is lowest prime factor of $$$x$$$ and $$$k$$$ is an integer
ProofEach $$$x$$$ can be represented uniquely as the form $$$x = lpf[x] \times k$$$.
So you can make a sieve of linear, where each number is visisted once or atmost $$$O(1)$$$
Solution for general k
Problem
Given $$$k, n (1 \leq k \leq n \leq 10^6)$$$, 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 k$$$
Since the number can be big, output it under modulo $$$10^9 + 7$$$.
You can submit here: https://lqdoj.edu.vn/problem/aicprtspk.
For $$$k = 3$$$, you can subimit here: https://oj.vnoi.info/problem/dhbb2020_square.
Examples
Example 1Input 1:
2 1
Output 1:
2
Explanation 1:
There are $$$2$$$ satisfied array of size $$$1$$$: {$$$1$$$}, {$$$2$$$}.
Example 2Input 2:
10 2
Output 2:
4
Explanation 2:
There are $$$4$$$ satisfied array of size $$$2$$$: {$$$1, 4$$$}, {$$$1, 9$$$}, {$$$2, 8$$$}, {$$$4, 9$$$}.
Example 3Input 3:
27 3
Output 3:
12
Explanation 3:
There are $$$12$$$ satisfied array of size $$$3$$$: {$$$1, 4, 9$$$}, {$$$1, 4, 16$$$}, {$$$1, 4, 25$$$}, {$$$1, 9, 16$$$}, {$$$1, 9, 25$$$}, {$$$1, 16, 25$$$}, {$$$2, 8, 18$$$}, {$$$3, 12, 27$$$}, {$$$4, 9, 16$$$}, {$$$4, 9, 25$$$}, {$$$4, 16, 25$$$}, {$$$9, 16, 25$$$}.
Using the same logic above, we can easily solve the problem.
Now you face up with familliar binomial coefficient problem
This implementation here is using the assumption of $$$p$$$ prime and $$$p > max(n, k)$$$
You can still solve the problem for squarefree $$$p$$$ using lucas and CRT
Yet just let things simple as we only focus on the counting problem, we will assume $$$p$$$ is a large constant prime.
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_nck(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, int k)
{
memset(is_squarefree, true, sizeof(is_squarefree[0]) * (n + 1));
precal_nck(n);
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;
}
Extra Tasks
These are harder variants, and generalization from the original problem. You can see more detail here
A: Can we also use phi function or something similar to solve for $$$k = 2$$$ in $$$O(\sqrt{n})$$$ or faster ?
B: Can we also use phi function or something similar to solve for general $$$k$$$ in $$$O(\sqrt{n})$$$ or faster ?
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 \leq n$$$ ?
F: Can we solve for a specific range, whereas $$$L \leq a_1 < a_2 < \dots < a_k \leq 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$$$ ?
J: Can we also solve the problem where there are no order: Just simply $$$1 \leq a_i \leq n$$$ ?
K: Can we also solve the problem where there are no order: Just simply $$$0 \leq a_i \leq n$$$ ?
M: Can we solve for $$$q$$$-product $$$a_{i_1} \times a_{i_2} \times \dots \times a_{i_q} = x^q$$$ (for given constant $$$q$$$) ?
N: Given $$$0 \leq \delta \leq n$$$, can we also solve the problem when $$$1 \leq a_1 \leq a_1 + \delta + \leq a_2 \leq a_2 + \delta \leq \dots \leq a_k \leq n$$$ ?
O: What if the condition is just two nearby elements and not all pairs. Or you can say $$$a_i \times a_{i+1} \forall 1 \leq i < n$$$ is a perfect square ?
Auto comment: topic has been updated by SPyofcode (previous revision, new revision, compare).
Go get a life pls
agreed
SPyofgame Master when?
.
I've had only a quick look and haven't checked your references yet. Your definitions of $$$f$$$ and $$$g$$$ seem off. As defined, $$$f(49) = \mathrm{size} \{(1, 7)\} = 1$$$ and $$$g(49) = \mathrm{size} \{(1, 7), (49, 49)\} = 2$$$, and $$$g(n) = f(n) + 1$$$, not $$$f(n) - 1$$$. But from the rest of the post it seems likely intended that $$$f(49) = 6$$$ and $$$g(49) = 7$$$. Then in the solution derivation it seems that it should be that $$$\overset{n}{\underset{p=1}{\Large \Sigma}} f(p) = -n + \overset{n}{\underset{p=1}{\Large \Sigma}} g(p)$$$ rather than what is written.
It possible to solve this problem in $$$O(\sqrt{n} \log{\log{n}})$$$ for any $$$k$$$ using a similar idea as that for $$$k = 2$$$. Here's a quick sketch of how to generalize your approach to any $$$k$$$:
Let $$$F_k(n)$$$ be the answer, and let $$$f_k(n) := F_k(n) - F_k(n-1)$$$. Then, clearly, $$$f_k(n)$$$ is the number of valid sequences of length $$$k$$$ ending in $$$n$$$, which is equal to $$$\binom{g(n)-1}{k-1}$$$, where $$$g(n)$$$ is the largest integer whose square divides $$$n$$$, as in your argument. Letting $$$s_k(i) = \binom{i-1}{k-1}$$$, we have by Möbius inversion that $$$f_k(n) = s_k(g(n)) = \sum_{d^2 | n} (\mu * s_k)(d)$$$. Then, as you argue above, we have that $$$F_k(n) = F_k(0) + \sum_{d=1}^{\lfloor \sqrt{n} \rfloor} (\mu * s_k)(d) \cdot \lfloor \frac{n}{d^2} \rfloor$$$. Since Möbius inversion can be done in $$$O(n \log{\log{n}})$$$ for a general (not necessarily multiplicative) sequence of length $$$n$$$, this formula takes the claimed $$$O(\sqrt{n} \log{\log{n}})$$$ when implemented directly.
(EDIT: minor notation cleanup)
It should also be possible to perform the transpose of the Möbius inversion on the sequence $$$(d \mapsto \lfloor \frac{n}{d^2} \rfloor)$$$ and use prefix sums on the sequence $$$s_k$$$ to calculate the answer, which I would intuitively expect to work in $$$\tilde{O}(n^{1/3})$$$. (That transpose looks similar to the DP described in this blog.)
The solution is great. But I still find it hard to understand that the part using Dirichlet Convolution. Can you eleborated more ?
I'm not sure which part you're referring to. Here are my guesses:
If you mean the equality $$$s_k(g(n)) = \sum_{d^2 | n} (\mu * s_k)(d)$$$, this is just the Möbius inversion formula combined with the fact that $$$d|g(n)$$$ if and only if $$$d^2|n$$$.
If you mean performing Möbius inversion on a sequence of length $$$n$$$ in $$$O(n \log{\log{n}})$$$, this is done by factoring the Möbius function. In terms of Dirichlet series/generating functions, this is dividing by the Euler product formula for the zeta function. Since for each prime $$$p$$$ the relevant term can be convolved with in $$$\lfloor \frac{n}{p} \rfloor$$$ steps, this works out to a total of around $$$n \log{\log{n}}$$$ steps.
(EDIT: fixed broken link)
If you mean the part where I wanted to perform the transpose of a Dirichlet convolution on $$$(d \mapsto \lfloor \frac{n}{d^2} \rfloor)$$$, I didn't provide any detail because I hadn't thought through exactly how to do this efficiently yet! It turns out that my intuition was too optimistic, but this can be done in $$$O(n^{2/5})$$$. However, it was easier for me to come up with this algorithm un-transposed, where it turns out to just be a simple modification of the (standard?) $$$O(n^{2/3})$$$ algorithm for calculating prefix sums of a Dirichlet convolution. In this prefix sums of $$$(\mu * s_k)$$$ are needed up to all indices $$$i$$$ for which $$$\lfloor \frac{n}{i^2} \rfloor \gt \lfloor \frac{n}{(i+1)^2} \rfloor$$$. Start by calculating prefix sums of $$$\mu$$$ and $$$s_k$$$ up to these same indices. Then, $$$\mu(i) \cdot s_k(j)$$$ contributes to the value of $$$(\mu * s_k)(i \cdot j)$$$; break these contributions up into three cases:
In the first case iterate over $$$i \cdot j$$$; in the second case iterate over $$$i$$$ and then $$$j$$$; in the third case iterate over $$$j$$$ and then $$$i$$$. That it is possible to do all of these using only the prefix sums I mentioned earlier and that these three cases each take $$$O(n^{2/5})$$$ time to consider, I leave as an exercise. (Transposing this algorithm is still somewhat helpful for extra task H, since the result can be shared between queries. The resulting complexity is $$$O(n^{2/5})$$$ plus the cost of evaluating $$$q$$$ degree-$$$k$$$ polynomials at $$$O(n^{1/3})$$$ points each.)
Regarding the extra tasks:
A, B,: Solved.
C: Redefine $$$s_k(i) = \binom{i+k-1}{k} - \binom{i+k-2}{k} = \binom{i+k-2}{k-1}$$$ and the solution goes through with no other changes.
D: Pretty much solved. The cases where $$$k$$$ is larger than around $$$p$$$ or $$$n^{2/5}$$$ are more annoying, but not really more difficult, assuming $$$p$$$ is prime or at least square-free.
E: Since all perfect squares are positive, we must have that in any valid sequence either all terms are at least 0, or all terms are at most 0. So, this reduces to the positive-only problem.
G: For $$$k > 3$$$ a very similar solution using $$$d^3|n$$$ instead of $$$d^2|n$$$ should just work. For $$$k = 3$$$ I'd need to think about it more.
H: The result at the end of the third spoiler is better than solving the whole problem $$$q$$$ times, but still not that great.
(EDIT: corrected wrong formula from misreading task C)
Great solution. But can I ask that why $$$p \approx n^{2/5}$$$ is also complicated ? Is it for general $$$k$$$ but with $$$O(n^{2/5})$$$ solution ?
It's an issue with claiming $$$O(n^{2/5})$$$ indeed. That's not enough time to naively pre-calculate the factorials and inverse factorials to make calculating the necessary binomial coefficients fast, and how messy the fallback strategies are depends on how large $$$k$$$ is.
*Actually, it's a bit worse than this with FFT-unfriendly primes. But not by enough to really matter for this discussion.
No offensive :( I just want to say that I did not expect it is that hard when $$$p \approx n^{2/5}$$$.