Number of array of size k that (1 <= a1 < a2 < ... < ak <= n) and each pair having a perfect square product

#	User	Rating
1	tourist	4009
2	jiangly	3823
3	Benq	3738
4	Radewoosh	3633
5	jqdai0815	3620
6	orzdevinwang	3529
7	ecnerwala	3446
8	Um_nik	3396
9	ksun48	3390
10	gamegame	3386

#	User	Contrib.
1	cry	167
2	Um_nik	163
3	maomao90	162
3	atcoder_official	162
5	adamant	159
6	-is-this-fft-	158
7	awoo	157
8	TheScrasse	154
9	Dominater069	153
9	nor	153

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

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

Implementation using factorization

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;
}

Implementation 1

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;
}

Implementation 2

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;
}

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

For the 2 implementations below: The complexity is $$$O\left(\underset{squarefree p \leq n}{\Large \Sigma} \left \lfloor \frac{n}{p^2} \right \rfloor \right)$$$

Solve for general k

Using the same logic above, we can easily solve the problem with combinatorics

O(n) solution

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;
}

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:

Swap for loop implementation

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;
}

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

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 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 also takes you to something similar.

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;
}

My question

A: Can we also use phi function or something similar to solve for general $$$k$$$ in $$$O(\sqrt{n})$$$ ?

B: Can we also solve the problem when $$$a_i \leq a_j (\forall\ i < j)$$$ and no longer $$$a_i < a_j (\forall\ i < j)$$$ ?

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

D: Can we solve for negative number, whereas $$$-n \leq a_1 < a_2 < \dots < a_k < n$$$

E: Can we solve for a specific range, wehreas $$$L \leq a_1 < a_2 < \dots < a_k < R$$$

F: Can we solve for cube product effectively ?

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