Number of array of size k that (1 <= a1 < a2 < ... < ak <= n) and each pair having a perfect square product
Difference between en2 and en3, changed 748 character(s)
-------------------------↵
## 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$ [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
: T, the complexity is $O\left(\underset{squarefree p \leq n}:↵

$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} \left \lfloor \frac{n1}{p^2} \right \rfloor \right) = O\left(n \times \frac{\pi^2}{6} \right) = O(n)$↵

-------------------------↵

## 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})$ ?↵

BC: Can we also solve the problem whenre there can be duplicate: $a_i \leq a_j (\forall\ i < j)$ and no longer $a_i < a_j (\forall\ i < j)$ ?↵

CD: Can we solve the problem where there is no restriction between $k, n, p$ ?↵

DE: Can we solve for negative numbintegers, whereas $-n \leq a_1 < a_2 < \dots < a_k < n$↵

EF: Can we solve for a specific range, wehhereas $L \leq a_1 < a_2 < \dots < a_k < R$↵

FG: Can we solve for cube product effectively ?↵

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$ ?

History

 
 
 
 
Revisions
 
 
  Rev. Lang. By When Δ Comment
en60 English SPyofcode 2021-11-09 12:21:19 4
en59 English SPyofcode 2021-11-09 12:20:38 6 Tiny change: 'n subimit [https://' -> 'n subimit here: [https://'
en58 English SPyofcode 2021-11-09 12:19:51 265
en57 English SPyofcode 2021-11-09 12:12:22 3178
en56 English SPyofcode 2021-11-05 19:03:46 3619
en55 English SPyofcode 2021-11-01 08:53:10 64
en54 English SPyofcode 2021-11-01 08:43:50 158
en53 English SPyofcode 2021-11-01 06:55:47 2 Tiny change: ' for $k = 3$ in $O(\s' -> ' for $k = 2$ in $O(\s'
en52 English SPyofcode 2021-11-01 06:44:45 71
en51 English SPyofcode 2021-11-01 06:40:09 16954 (published)
en50 English SPyofcode 2021-11-01 06:22:49 49881 (saved to drafts)
en49 English SPyofcode 2021-11-01 05:27:15 49
en48 English SPyofcode 2021-11-01 05:18:00 184
en47 English SPyofcode 2021-10-31 19:17:30 40
en46 English SPyofcode 2021-10-31 19:10:27 5305
en45 English SPyofcode 2021-10-31 12:29:39 196
en44 English SPyofcode 2021-10-31 12:15:21 229
en43 English SPyofcode 2021-10-31 12:12:51 5175
en42 English SPyofcode 2021-10-31 12:04:20 4012 (published)
en41 English SPyofcode 2021-10-31 07:41:25 0 (saved to drafts)
en40 English SPyofcode 2021-10-31 07:10:39 916 Tiny change: ') solution">\n\n```c' -> ') solution for k = 3">\n\n```c'
en39 English SPyofcode 2021-10-31 06:59:31 6996
en38 English SPyofcode 2021-10-30 20:35:28 436
en37 English SPyofcode 2021-10-30 19:51:26 8
en36 English SPyofcode 2021-10-30 19:50:57 168
en35 English SPyofcode 2021-10-30 19:49:17 38
en34 English SPyofcode 2021-10-30 19:48:18 710
en33 English SPyofcode 2021-10-30 19:43:42 6
en32 English SPyofcode 2021-10-30 19:42:58 714
en31 English SPyofcode 2021-10-30 19:36:20 3825
en30 English SPyofcode 2021-10-30 16:20:57 8
en29 English SPyofcode 2021-10-30 16:19:47 2877
en28 English SPyofcode 2021-10-30 16:18:30 56
en27 English SPyofcode 2021-10-30 16:16:24 6543
en26 English SPyofcode 2021-10-30 11:38:14 769
en25 English SPyofcode 2021-10-30 06:48:18 76
en24 English SPyofcode 2021-10-30 06:45:41 23 Reverted to en22
en23 English SPyofcode 2021-10-30 06:41:21 23
en22 English SPyofcode 2021-10-30 05:34:13 2034 Tiny change: ' summary="Hint 2">\n\nThe ' -> ' summary="Proof">\n\nThe '
en21 English SPyofcode 2021-10-30 04:47:01 229
en20 English SPyofcode 2021-10-29 20:05:05 68 Tiny change: '---\n\n## Tasks\n\n' -> '---\n\n## Extra Tasks\n\n'
en19 English SPyofcode 2021-10-29 19:57:01 603
en18 English SPyofcode 2021-10-29 19:50:27 5436 Tiny change: '------\n\n' -> '------\n\n-------------------------\n\n-------------------------\n\n'
en17 English SPyofcode 2021-10-29 17:43:10 26
en16 English SPyofcode 2021-10-29 12:50:12 3892
en15 English SPyofcode 2021-10-29 12:16:06 1739
en14 English SPyofcode 2021-10-29 03:52:02 2
en13 English SPyofcode 2021-10-28 18:40:23 452 Tiny change: '## Statement:\' -> '## The statement:\'
en12 English SPyofcode 2021-10-28 18:31:38 1487
en11 English SPyofcode 2021-10-28 18:09:33 748
en10 English SPyofcode 2021-10-28 13:42:36 278
en9 English SPyofcode 2021-10-28 13:31:37 99
en8 English SPyofcode 2021-10-28 13:23:42 78
en7 English SPyofcode 2021-10-28 13:21:55 1383
en6 English SPyofcode 2021-10-28 13:08:47 885 Tiny change: 'e product effective' -> 'e product $a_i \times a_j \times a_k$ effective'
en5 English SPyofcode 2021-10-28 13:01:32 1 Tiny change: 'n $O\left(sqrt{n} \l' -> 'n $O\left(\sqrt{n} \l'
en4 English SPyofcode 2021-10-28 13:00:39 478
en3 English SPyofcode 2021-10-28 12:53:55 748 (published)
en2 English SPyofcode 2021-10-28 12:42:55 7656
en1 English SPyofcode 2021-10-28 12:05:23 3517 Initial revision (saved to drafts)