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

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 pair $$$(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 pair $$$(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)$$$.

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

Yet this paper also takes you to something similar.

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

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

The complexity is depended on euler totient phi sieve and calculating for answer.

But as you can see, both can be done in $$$O(\sqrt n)$$$.

As what clyring decribed here.

We can generalize those definition above into $$$k$$$-term.

Let $$$f_k(n)$$$ is the number of set $$$(a_1, a_2, \dots, a_k, n)$$$ that $$$1 \leq a_1 < a_2 < \dots < a_k \leq n$$$ and $$$(a_1, a_2, \dots, a_k, n)$$$ is a $$$(k+1)$$$-term geometric progression.

Let $$$g_k(n)$$$ is the number of set $$$(a_1, a_2, \dots, a_k, n)$$$ that $$$1 \leq a_1 \leq a_2 \leq \dots \leq a_k \leq n$$$ and $$$(a_1, a_2, \dots, a_k, n)$$$ is a $$$(k+1)$$$-term geometric progression.

Let $$$F_k(n) = \overset{n}{\underset{p=1}{\Large \Sigma}} f_k(p)$$$.

Let $$$s_k(n)$$$ is the number of way to choose $$$p^2$$$ among those $$$k$$$ numbers when you fix squarefree $$$u$$$ (though we are doing in reverse).

$$$s_k(n) = C^{n-1}_{k-1}$$$

$$$f_k(n) = F_k(n) - F_k(n-1) = s_k(g_k(n)) = \underset{d^2 | n}{\Large \Sigma} (\mu \times s_k)(d) = \underset{d^2 | n}{\Large \Sigma} \underset{p | d}{\Large \Sigma} \mu(p) \times s_k\left(\frac{d}{p}\right)$$$

$$$F_k(n) = \overset{n}{\underset{p=1}{\Large \Sigma}} f_k(p)$$$

$$$= F_k(0) + \overset{\left \lfloor \sqrt{n} \right \rfloor}{\underset{d = 1}{\Large \Sigma}}(\mu \times s_k)(d) \times \left \lfloor \frac{n}{d^2} \right \rfloor$$$

$$$= \overset{\left \lfloor \sqrt{n} \right \rfloor}{\underset{d = 1}{\Large \Sigma}} \left \lfloor \frac{n}{d^2} \right \rfloor \times \left ( \underset{p | d}{\Large \Sigma} \mu(p) \times s_k\left(\frac{d}{p}\right) \right)$$$

const int LIM = 5e6 + 56;
const int SQRT_LIM = ceil(sqrt(LIM) + 1) + 1;
const int MOD = 1e9 + 7;

/// Precalculating factorials under prime modulo
int fact[SQRT_LIM + 10]; /// fact[n] = n!
int invs[SQRT_LIM + 10]; /// invs[n] = n^(-1)
int tcaf[SQRT_LIM + 10]; /// tcaf[n] = (n!)^(-1)
void precal_nck(int n = SQRT_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;
    }
}

/// Calculating binomial coefficient queries
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;
}

/// Linear Sieve
vector<int> prime;           /// prime list              = A000040
bool isPrime[SQRT_LIM + 10]; /// characteristic function = A010051
int lpf[SQRT_LIM + 10];      /// lowest prime factor     = A020639
int mu[SQRT_LIM + 10];       /// mobius                  = A008683
void linear_sieve(int n)
{
    if (n < 1) return ;
    /// Extension Sieve || You can add something more
    memset(lpf, 0, sizeof(lpf[0]) * (n + 1));
    fill_n(mu, n + 1, 1);
    /// Main Sieve || Without this, you barely able to achive linear complexity
    prime.clear();
    prime.reserve(n / log(n - 1));
    memset(isPrime, true, sizeof(isPrime[0]) * (n + 1));
    isPrime[0] = isPrime[1] = false;
    for (int x = 2; x <= n; ++x) /// For each number
    {
        if (isPrime[x]) /// Func[Prime]
        {
            mu[x] = -1;
            lpf[x] = x;
            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;
            mu[x * p] = (lpf[x] == p) ? 0 : -mu[x];
        }
    }
}

/// Divisor sieve
vector<int> divisors[SQRT_LIM];
void precal_div(int n) /// O(n log n)
{
    for (int u = n; u >= 1; --u)
    {
        divisors[u].clear();
        for (int v = u; v <= n; v += u)
            divisors[v].push_back(u);
    }
}

/// Solving for n, k
long long solve(int n, int k)
{
    /// We only care for d that 1 <= d <= sqrt(n)
    int t = ceil(sqrt(n) + 1) + 1;
    linear_sieve(t);
    precal_nck(t);
    precal_div(t);

    long long res = 0;
    for (int d = 1; d * d <= n; ++d) /// For each fixed p^2
    {
        long long sum = 0;
        for (int p : divisors[d]) /// For each (p | d)
            sum += mu[d / p] * nck(p - 1, k - 1);

        sum %= MOD;
        res += sum * (n / (d * d));
    }

    res %= MOD;
    return res;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(NULL);

    /// Assumming constant p = 10^9 + 7
    int n, k;
    cin >> n >> k;
    cout << solve(n, k);
    return 0;
}

vector<int> prime;           /// prime list              = A000040
bool isPrime[SQRT_LIM + 10]; /// characteristic function = A010051
int lpf[SQRT_LIM + 10];      /// lowest prime factor     = A020639
void linear_sieve(int n)
{
    if (n < 1) return ;
    prime.clear();
    prime.reserve(n / log(n - 1));
    memset(lpf, 0, sizeof(lpf[0]) * (n + 1));
    memset(isPrime, true, sizeof(isPrime[0]) * (n + 1));
    isPrime[0] = isPrime[1] = false;
    for (int x = 2; x <= n; ++x)
    {
        if (isPrime[x]) /// Func[Prime]
        {
            lpf[x] = x;
            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;
        }
    }
}

long long res[SQRT_LIM + 10];
long long solve(int n, int k)
{
    int t = ceil(sqrt(n) + 1) + 1;
    linear_sieve(t);
    precal_nck(t);

    memset(res, 0, sizeof(res[0]) * (t + 1));
    for (int d = 1; d * d <= n; ++d) 
        res[d] = nck(d - 1, k - 1);

    for (int p : prime)
        for (int d = t / p; d > 0; --d)
            res[d * p] -= res[d];

    long long ans = 0;
    for (int d = 1; d * d <= n; ++d)
        ans += res[d] * (n / (d * d));

    ans %= MOD;
    return ans;
}

The complexity of the first implementation is $$$O(\sqrt{n} \log \sqrt{n})$$$.

The complexity of the second implementation is $$$O(\sqrt{n} \log \log \sqrt{n})$$$.

It is no hard to prove that we can use the same algorithm as described in task A, B or in original task.

The core of calculating is to find out the number of non-decreasing integer sequence of size $$$k$$$ where numbers are in $$$[1, n]$$$.

Can you prove it ?

Now it is done, just that it !

The idea is the same as what clyring described here but represented in the other way.

int fact[SQRT_LIM + 10];
int invs[SQRT_LIM + 10];
int tcaf[SQRT_LIM + 10];
void precal_nck(int n = SQRT_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;
}

bool is_squarefree[LIM];
int solve(int n, int k)
{
    memset(is_squarefree, true, sizeof(is_squarefree[0]) * (n + 1));
    precal_nck(2 * 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(k + j - 2, k);
    }

    res %= MOD;
    return res;
}

const int LIM = 5e6 + 56;
const int SQRT_LIM = ceil(sqrt(LIM) + 1) + 1;
const int MOD = 1e9 + 7;

/// Precalculating factorials under prime modulo
int fact[SQRT_LIM + 10]; /// fact[n] = n!
int invs[SQRT_LIM + 10]; /// invs[n] = n^(-1)
int tcaf[SQRT_LIM + 10]; /// tcaf[n] = (n!)^(-1)
void precal_nck(int n = SQRT_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;
    }
}

/// Calculating binomial coefficient queries
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;
}

/// Linear Sieve
vector<int> prime;           /// prime list              = A000040
bool isPrime[SQRT_LIM + 10]; /// characteristic function = A010051
int lpf[SQRT_LIM + 10];      /// lowest prime factor     = A020639
int mu[SQRT_LIM + 10];       /// mobius                  = A008683
void linear_sieve(int n)
{
    if (n < 1) return ;
    /// Extension Sieve || You can add something more
    memset(lpf, 0, sizeof(lpf[0]) * (n + 1));
    fill_n(mu, n + 1, 1);
    /// Main Sieve || Without this, you barely able to achive linear complexity
    prime.clear();
    prime.reserve(n / log(n - 1));
    memset(isPrime, true, sizeof(isPrime[0]) * (n + 1));
    isPrime[0] = isPrime[1] = false;
    for (int x = 2; x <= n; ++x) /// For each number
    {
        if (isPrime[x]) /// Func[Prime]
        {
            mu[x] = -1;
            lpf[x] = x;
            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;
            mu[x * p] = (lpf[x] == p) ? 0 : -mu[x];
        }
    }
}

/// Divisor sieve
vector<int> divisors[SQRT_LIM];
void precal_div(int n) /// O(n log n)
{
    for (int u = n; u >= 1; --u)
    {
        divisors[u].clear();
        for (int v = u; v <= n; v += u)
            divisors[v].push_back(u);
    }
}

/// Solving for n, k
long long solve(int n, int k)
{
    /// We only care for d that 1 <= d <= sqrt(n)
    int t = ceil(sqrt(n) + 1) + 1;
    linear_sieve(t);
    precal_nck(t);
    precal_div(t);

    long long res = 0;
    for (int d = 1; d * d <= n; ++d) /// For each fixed p^2
    {
        long long sum = 0;
        for (int p : divisors[d]) /// For each (p | d)
            sum += mu[d / p] * nck(d + k - 2, k - 1);

        sum %= MOD;
        res += sum * (n / (d * d));
    }

    res %= MOD;
    return res;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(NULL);

    /// Assumming constant p = 10^9 + 7
    int n, k;
    cin >> n >> k;
    cout << solve(n, k);
    return 0;
}

int fact[SQRT_LIM + 10];
int invs[SQRT_LIM + 10];
int tcaf[SQRT_LIM + 10];
void precal_nck(int n = SQRT_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;
}

vector<int> prime;           /// prime list              = A000040
bool isPrime[SQRT_LIM + 10]; /// characteristic function = A010051
int lpf[SQRT_LIM + 10];      /// lowest prime factor     = A020639
void linear_sieve(int n)
{
    if (n < 1) return ;
    prime.clear();
    prime.reserve(n / log(n - 1));
    memset(lpf, 0, sizeof(lpf[0]) * (n + 1));
    memset(isPrime, true, sizeof(isPrime[0]) * (n + 1));
    isPrime[0] = isPrime[1] = false;
    for (int x = 2; x <= n; ++x)
    {
        if (isPrime[x]) /// Func[Prime]
        {
            lpf[x] = x;
            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;
        }
    }
}

long long res[SQRT_LIM + 10];
long long solve(int n, int k)
{
    int t = ceil(sqrt(n) + 1) + 1;
    linear_sieve(t + k);
    precal_nck(t + k);

    memset(res, 0, sizeof(res[0]) * (t + 1));
    for (int d = 1; d * d <= n; ++d) 
        res[d] = nck(d + k - 2, k - 1);

    for (int p : prime)
        for (int d = t / p; d > 0; --d)
            res[d * p] -= res[d];

    long long ans = 0;
    for (int d = 1; d * d <= n; ++d)
        ans += res[d] * (n / (d * d));

    ans %= MOD;
    return ans;
}

In the first implementation it is obviously linear.

The second and third implementation is also easy to show its complexity

Sadly, since $$$k \leq n$$$. We also conclude that the complexity is $$$O(n)$$$, and even worse it also contains large constant factor compared to that in the first implementation.

But it is still effecient enough to solve problem where $$$k$$$ is small.

So first of all, the result do depend on how you calculate binomial coefficient but they are calculated independently even if you can somehow manage to use the for loop of binomial coeffient go first.

Therefore even if there is no restriction between $$$k, n, p$$$, the counting part and the algorithm doesnt change.

You just need to change how you calculate binomial coefficient, and that is all for this task.

Let just ignore the fact that though this need more detail, but as the blog is not about nck problem I will just make it quick.

For large prime $$$p > max(n, k)$$$

• Just using normal combinatorics related to factorial (since $$$p > max(n, k)$$$ nothing will affect the result)
• For taking divides under modulo you can just take modular inversion (as a prime always exist such number)
• Yet this is standard problem, just becareful of the overflow part
• You can also optimize by precalculating factorial, inversion number and inversion factorial in linear too

For general prime $$$p$$$

• We can just ignore factors $$$p$$$ in calculating $$$n!$$$.
• You also need to know how many times factor $$$p$$$ appears in $$$1 \dots n$$$
• Then combining it back when calculating for the answer.
• If we dont do this $$$n!$$$ become might divides some factors of $$$p$$$.
• By precalculation you can answer queries in $$$O(1)$$$

For squarefree $$$p$$$

• Factorize $$$p = p_1 \times p_2 \times p_q$$$ that all $$$p_i$$$ is prime.
• Ignore all factors $$$p_i$$$ when calculate $$$n!$$$.
• Remember to calculate how many times factors $$$p_i$$$ appear in $$$1 \dots n$$$.
• When query for the answer we just combine all those part back.
• Remember you can just take modulo upto $$$\phi(p)$$$ which you can also calculate while factorizing $$$p$$$.
• Remember that $$$n!$$$ must not divides any factor $$$p_i$$$ otherwise you will get wrong answer.
• By precalculation you can answer queries in $$$O(\log p)$$$

For general positive modulo $$$p$$$

• Factorize $$$p = p_1^{f_1} \times p_2^{f_2} \times p_q^{f_q}$$$ that all $$$p_i$$$ is unique prime.
• We calculate $$$C(n, k)$$$ modulo $$$p_i^{f_i}$$$ for each $$$i = 1 \dots q$$$.
• To do that, we need to calculate $$$n!$$$ modulo $$$p_i^{f_i}$$$ which is described here.
• To get the final answer we can use CRT.
• Yet this is kinda hard to code and debug also easy to make mistake so you must becareful
• I will let the implementation for you lovely readers.
• Yet depends on how you calculate stuffs that might increase your query complexity
• There are few (effective or atleast fully correct) papers about this but you can read the one written here

/// SPyofgame linear template for precalculating factorials under large prime modulo
int fact[SQRT_LIM + 10]; /// fact[n] = n!
int invs[SQRT_LIM + 10]; /// invs[n] = n^(-1)
int tcaf[SQRT_LIM + 10]; /// tcaf[n] = (n!)^(-1)
void precal_nck(int n = SQRT_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;
    }
}

/// Calculating binomial coefficient queries
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;
}

vector<int> factor;
int factorize(int n) /// Calculating phi(n) while factorizing (n) in O(sqrt n)
{
    factor.clear();
    int phi = n;

    if (!(n & 1))
    {
        n >>= __builtin_ctz(n);
        factor.push_back(2);
        phi -= phi / 2;
    }

    for (int x = 3; x * x <= n; x += 2)
    {
        if (n % x == 0)
        {
            do n /= x; while (n % x == 0);
            factor.push_back(x);
            phi -= phi / x;
        }
    }

    if (n > 1)
    {
        factor.push_back(n);
        phi -= phi / n;
    }

    return phi;
}

int f[LIM];    /// f[x] = nck(n, x)
int fact[LIM]; /// n! 
int tcaf[LIM]; /// n!^(-1)
int divp[LIM]; /// x but ignore all factors p[i]
int cntp[LIM][LOG_LIM]; /// cntp[x][i] = Number of time factor p[i] appear in 1..x
void precal(int MOD) /// Calculate f[x] for all x = 1 -> n in O(n log mod + sqrt mod)
{
    int PHIMOD = factorize(MOD);
    for (int x = 1; x <= n; ++x) /// For each part x in n!
    {
        int &t = divp[x] = x;
        for (int i = 0; i < factor.size(); ++i) /// Ignore all factor p[i] of p
        {
            cntp[x][i] = cntp[x - 1][i];
            for (; t % factor[i] == 0; t /= factor[i]) /// Count how many times p[i] appears in 1..n
                ++cntp[x][i];
        }
    }

    fact[0] = fact[1] = 1;
    tcaf[0] = tcaf[1] = 1;
    for (int x = 2; x <= n; ++x) /// Finding n! and n!^(-1)
    {
        fact[x] = (1LL * fact[x - 1] * divp[x]) % MOD;
        tcaf[x] = powMOD(fact[x], PHIMOD - 1, MOD);
    }

    memset(f, 0, sizeof(f[0]) * k);
    for (int x = k; x <= n; ++x)
    {
        /// Calculate nck % p normally
        f[x] = fact[x];
        mulMOD(f[x], tcaf[k], MOD);
        mulMOD(f[x], tcaf[x - k], MOD);
        for (int i = 0; i < factor.size(); ++i) /// Bringing those factors back
        {
            int p = cntp[x][i] - cntp[k][i] - cntp[x - k][i];
            f[x] = 1LL * f[x] * powMOD(factor[i], p, MOD) % MOD;
        }
    }
}

In the first implementation it is obviously linear.

And for the second implementation.

So you got $$$O(n \times \log p + \sqrt{p})$$$ in final.

Though you can still optimize this but by doing that why dont you just go straight up to solve for non squarefree $$$p$$$ too ?

Yet this is the same as extra task C where only the counting part should be changed.

As we only care about integer therefore let not use complex math into this problem.

If there exist a negative number and a positive number, the product will be negative thus the sequence will not satisfied.

Becareful, there are the zeros too.

When the numbers are all unique, or $$$-n \leq a_1 < a_2 < \dots < a_k \leq n$$$

Thus give us the formula of $$$task_E(n, k) = 2 \times task_B(n, k) + 2 \times task_B(n, k - 1)$$$.

Remember that when $$$k = 0$$$ the answer is $$$0$$$ otherwise you might somewhat having wrong result for negative number in binomial coefficients formula

So what if I mix the problem with task C too ?

When the numbers can have duplicates, or $$$-n \leq a_1 \leq a_2 \leq \dots \leq a_k \leq n$$$

Yet once again you can simplified it with less cases for easier calculation.

Thus give us the formula of $$$task_E(n, k) = 1 + 2 \times \overset{k}{\underset{t = 1}{\Large \Sigma}} task_B(n, t)$$$.

But this give you a $$$O(k)$$$ solution.

You can do better with math

long long res[SQRT_LIM + 10];
long long solve(int n, int k)
{
    int t = ceil(sqrt(n) + 1) + 1;
    linear_sieve(t);
    precal_nck(t);

    memset(res, 0, sizeof(res[0]) * (t + 1));
    for (int d = 1; d * d <= n; ++d) 
    {
        res[d] += (k >= 1) * nck(d - 1, k - 1) * 2;
        res[d] += (k >= 2) * nck(d - 1, k - 2) * 2;  
    }

    for (int p : prime)
        for (int d = t / p; d > 0; --d)
            res[d * p] -= res[d];

    long long ans = 0;
    for (int d = 1; d * d <= n; ++d)
        ans += res[d] * (n / (d * d));

    ans %= MOD;
    return ans;
}

bool is_squarefree[LIM];
int brute(int n, int k)
{
    memset(is_squarefree, true, sizeof(is_squarefree[0]) * (n + 1));
    precal_nck(2 * 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(k + j - 2, k);
    }

    res %= MOD;
    return res;
}


long long solve(int n, int k)
{
    long long res = 1;
    for (int t = 1; t <= k; ++t)
        res += brute(n, t) * 2;

    res %= MOD;
    return res;
}

long long res[SQRT_LIM + 10];
long long brute(int n, int k)
{
    /// We only care for d that 1 <= d <= sqrt(n)
    int t = ceil(sqrt(n) + 1) + 1;
    linear_sieve(t);
    precal_nck(t + k);
    precal_div(t);

    long long res = 0;
    for (int d = 1; d * d <= n; ++d) /// For each fixed p^2
    {
        long long sum = 0;
        for (int p : divisors[d]) /// For each (p | d)
            sum += mu[d / p] * nck(d + k - 2, k - 1);

        sum %= MOD;
        res += sum * (n / (d * d));
    }

    res %= MOD;
    return res;
}

long long solve(int n, int k)
{
    long long res = 1;
    for (int t = 1; t <= k; ++t)
        res += brute(n, t) * 2;

    res %= MOD;
    return res;
}

long long res[SQRT_LIM + 10];
long long brute(int n, int k)
{
    int t = ceil(sqrt(n) + 1) + 1;
    linear_sieve(t);
    precal_nck(t + k);

    memset(res, 0, sizeof(res[0]) * (t + 1));
    for (int d = 1; d * d <= n; ++d) 
        res[d] = nck(d + k - 2, k - 1);

    for (int p : prime)
        for (int d = t / p; d > 0; --d)
            res[d * p] -= res[d];

    long long ans = 0;
    for (int d = 1; d * d <= n; ++d)
        ans += res[d] * (n / (d * d));

    ans %= MOD;
    return ans;
}

long long solve(int n, int k)
{
    long long res = 1;
    for (int t = 1; t <= k; ++t)
        res += brute(n, t) * 2;

    res %= MOD;
    return res;
}

long long res[SQRT_LIM + 10];
long long solve(int n, int k)
{
    int t = ceil(sqrt(n) + 1) + 1;
    linear_sieve(t);
    precal_nck(t + k);

    memset(res, 0, sizeof(res[0]) * (t + 1));
    for (int d = 1; d * d <= n; ++d) 
        for (int t = 1; t <= k; ++t)
            res[d] += nck(d + t - 2, t - 1) * 2;

    for (int p : prime)
        for (int d = t / p; d > 0; --d)
            res[d * p] -= res[d];

    long long ans = 1;
    for (int d = 1; d * d <= n; ++d)
        ans += res[d] * (n / (d * d));

    ans %= MOD;
    return ans;
}

long long res[SQRT_LIM + 10];
long long solve(int n, int k)
{
    int t = ceil(sqrt(n) + 1) + 1;
    linear_sieve(t);
    precal_nck(t + k);

    memset(res, 0, sizeof(res[0]) * (t + 1));
    for (int d = 1; d * d <= n; ++d) 
        res[d] += nck(d + k - 1, k - 1) * 2;

    for (int p : prime)
        for (int d = t / p; d > 0; --d)
            res[d * p] -= res[d];

    long long ans = 1;
    for (int d = 1; d * d <= n; ++d)
        ans += res[d] * (n / (d * d));

    ans %= MOD;
    return ans;
}

So.. you might be tired of calculating complexity again and again for those are too familar to you.

So I gonna skip this as the proof is as the same as what you can read above.

But as a bonus, you optimize the problem using this

Then you can solve the problem in $$$O(\sqrt{n} \log \log \sqrt{n} + \sqrt{p} \log^q \ p)$$$ ($$$q$$$ is depended on how you implement it).

For unique array, the complexity can be $$$O(\sqrt{n} \log \log \sqrt{n})$$$

For duplicate array, the complexity can be $$$O(\sqrt{n} \log \log \sqrt{n} + k)$$$

For duplicate array and small prime $$$p > max(n, k)$$$, the complexity can be $$$O(\sqrt{n} \log \log \sqrt{n} + \sqrt{p} \log^q \ p)$$$

We split into 4 cases

You can easilier solve for each cases in linear

Now assumming $$$1 \leq L \leq R$$$, here is how we solve it.

As a simple approach, we can just do like original problem for general $$$k$$$.

We just need to iterative for each fixed squarefree $$$u$$$ and count the number of way to select $$$p^2$$$ as usual but a bit of change.

By doing so trivially we can have the complexity of $$$O((R - L + 1) \log R)$$$ time and $$$O(R)$$$ space.

Is it over ? Cant we do better ?

So what is that the better solution and how do we implement it steps by steps ?

bool is_squarefree[LIM + 10];
int solve(int l, int r, int k)
{
    if (l > r) return 0;
    if (r < 0) return solve(-r, -l, k);
    if (l <= 0 && 0 <= r)
    {
        long long res = 0LL + taskB(abs(l), k) + taskB(abs(l), k - 1) + taskB(abs(r), k) + taskB(abs(r), k - 1);
        while (res >= MOD) res -= MOD;
        return res;
    }

    memset(is_squarefree, true, sizeof(is_squarefree[0]) * (r + l - 1));
    precal_nck(r - l + 1);

    int tot = r - l + 1;
    long long res = 0;
    for (int i = 1; i <= r; ++i)
    {
        int j = sqrt(l / i);
        while (i * j * j > l) --j;
        while (i * j * j < l) ++j;
        if (is_squarefree[i * j * j - l])
        {
            int cnt = 0;
            for (; i * j * j <= r; ++j)
            {
                ++cnt;
                --tot;
                is_squarefree[i * j * j - l] = false;
            }
            res += nck(cnt, k);
            if (tot == 0) break;
        }
    }

    res %= MOD;
    return res;
}

long long res[SQRT_LIM + 10];
long long taskB(int n, int k)
{
    if (k == 0) return 0;
    if (k == 1) return n;

    int t = ceil(sqrt(n) + 1) + 1;
    linear_sieve(t);
    precal_nck(t);

    memset(res, 0, sizeof(res[0]) * (t + 1));
    for (int d = 1; d * d <= n; ++d) 
        res[d] = nck(d - 1, k - 1);

    for (int p : prime)
        for (int d = t / p; d > 0; --d)
            res[d * p] -= res[d];

    long long ans = 0;
    for (int d = 1; d * d <= n; ++d)
        ans += res[d] * (n / (d * d));

    ans %= MOD;
    return ans;
}

bool is_squarefree[LIM + 10];
int solve(int l, int r, int k)
{
    if (l > r) return 0;
    if (r < 0) return solve(-r, -l, k);
    if (l <= 0 && 0 <= r)
    {
        long long res = 0LL + taskB(abs(l), k) + taskB(abs(l), k - 1) + taskB(abs(r), k) + taskB(abs(r), k - 1);
        while (res >= MOD) res -= MOD;
        return res;
    }
    
    int t = ceil(sqrt(r + 1) + 1) + 1;
    memset(is_squarefree, true, sizeof(is_squarefree[0]) * (r + l - 1));
    precal_nck(r - l + 1);
    linear_sieve(t);

    long long res = 0;
    for (int x = l; x <= r; ++x) if (is_squarefree[x - l])
    {
        int u = x;
        int v = 1;
        for (int p : prime)
        {
            int q = p * p;
            if (q > u) break;
            while (u % q == 0)
            {
                u /= q;
                v *= p;
            }
        }

        int cnt = 0;
        for (; u * v * v <= r; ++v)
        {
            is_squarefree[u * v * v - l] = false;
            ++cnt;
        }

        res += nck(cnt, k);
    }

    res %= MOD;
    return res;
}

long long res[SQRT_LIM + 10];
long long taskB(int n, int k)
{
    if (k == 0) return 0;
    if (k == 1) return n;

    int t = ceil(sqrt(n) + 1) + 1;
    linear_sieve(t);
    precal_nck(t);

    memset(res, 0, sizeof(res[0]) * (t + 1));
    for (int d = 1; d * d <= n; ++d) 
        res[d] = nck(d - 1, k - 1);

    for (int p : prime)
        for (int d = t / p; d > 0; --d)
            res[d * p] -= res[d];

    long long ans = 0;
    for (int d = 1; d * d <= n; ++d)
        ans += res[d] * (n / (d * d));

    ans %= MOD;
    return ans;
}

#include <algorithm>
#define all(x) (x).begin(), (x).end()
bool is_squarefree[LIM + 10];
int squarefree[LIM + 10];
int sqf_factor[LIM + 10];
int solve(int l, int r, int k)
{
    if (l > r) return 0;
    if (r < 0) return solve(-r, -l, k);
    if (l <= 0 && 0 <= r)
    {
        long long res = 0LL + taskB(abs(l), k) + taskB(abs(l), k - 1) + taskB(abs(r), k) + taskB(abs(r), k - 1);
        while (res >= MOD) res -= MOD;
        return res;
    }
    
    int t = ceil(sqrt(r + 1) + 1) + 1;
    memset(is_squarefree, true, sizeof(is_squarefree[0]) * (r + l - 1));
    precal_nck(r - l + 1);
    linear_sieve(t);
    
    for (int x = l; x <= r; ++x)
    {
        squarefree[x - l] = x;
        sqf_factor[x - l] = 1;
    }

    long long res = 0;
    for (int p : prime)
    {
        for (int q = p * p, x = max(q, (l + q - 1) / q * q); x <= r; x += q)
        {
            while (squarefree[x - l] % q == 0)
            {
                squarefree[x - l] /= q;
                sqf_factor[x - l] *= p;
            }
        }
    }
    
    for (int x = l; x <= r; ++x)
    {
        int u = squarefree[x - l];
        int p = sqf_factor[x - l];
        if (is_squarefree[x - l])
        {
            int cnt = 0;
            for (; u * p * p <= r; ++p)
            {
                is_squarefree[u * p * p - l] = false;
                ++cnt;
            }
    
            res += nck(cnt, k);   
        }
    }
    
    res %= MOD;
    return res;
}

Well, you might feel bored of long chain of proving for just a simple complexity, now I just make it quick.

Which part you find hard to understand, you can comment below if you cant accept this as an exercise.

For the first implemenetation it is obviously linear time in $$$max(|L|, |R|)$$$ and linear space in $$$|L - R|$$$

For the second implementation it is harder to prove. Yet it is bounded by the complexity of factorizing each number in range $$$[L, R]$$$.

For the third implementation it is based on the fact that:

• The sum of inversion of prime squared is constant.
• Each number $$$x$$$ will not be factorized more than $$$O(\log(\sqrt{x}))$$$.
• You can precalculate binomial coefficient in linear and query for constant time.
• You can sieve in linear, and only sieve to the square root of $$$R$$$.
• For each number in the range we care, we just need to visit once.

Yet you can furthur optimize it to $$$O\left(\frac{\sqrt{R} \log \log \sqrt{R}}{\log \sqrt{R}} + \frac{R - L}{\pi^2}\right)$$$.

*I am working on a $$$O {\Large ( } Õ(\sqrt{R}) + Õ(\sqrt{R - L}) {\Large ) }$$$ solution right now as it is much harder to apply the same difficult math idea described in task B.

For $$$k < 3$$$, the formula is easy to shown as the second condition always satisfied.

For $$$k > 3$$$, you can prove that every number you selected must share same cubefree therefore just make a cubefree sieve in linear.

But still, you can apply the same idea used in extra task B to achive better complexity.

Instead of throwing a bunch of math using weird formulas with long chain of theorems and proving stuffs.

We can see the algorithm in the other way then come to the formula.

So now we come for the formula:

Yet in task $$$k = 3$$$ things got differently.

By the theorem used above that if $$$b_x \times b_y$$$ is a cube and $$$b_x, b_y$$$ are cubefree numbers. Then if you know $$$b_x$$$ you can always find one and only one $$$b_y$$$.

So instead of a brute-forces, you can fix $$$2$$$ numbers and then you factorize it to the form of cubefree number.

Yet you can know the number of ways to select $$$b_x \times b_y$$$ with precalculation.

Then you can just iterate all over $$$b_z$$$ and count the total of ways to select $$$b_x \times b_y\times b_z$$$.

Since they are ordered, you must becareful, you can update new $$$b_x \times b_y (x < y)$$$ exists whenever your $$$z > y$$$.

int cnt[LIM];
vector<int> appear[LIM];
int solveG(int n, int k)
{
    linear_sieve(n * n);
    for (int i = 1; i <= n; ++i)
        appear[i].clear();

    for (int i = 1; i <= n; ++i)
    {
        for (int j = i + 1; j <= n; ++j)
        {
            int u = 1;
            for (int x = i * j; x > 1; )
            {
                int a = lpf[x], b = a * a * a;
                for (; x % b == 0; x /= b);
                for (; x % a == 0; x /= a) u *= a;
            }
            appear[j + 1].push_back(u);
        }
    }

    long long res = 0;
    memset(cnt, 0, sizeof(cnt[0]) * (n * n + 1));
    for (int i = 1; i <= n; ++i)
    {
        for (int u : appear[i]) ++cnt[u];
        int v = 1;
        for (int x = i; x > 1; )
        {
            int a = lpf[x], b = a * a * a;
            for (; x % b == 0; x /= b);
            if (x % (a * a) == 0)
            {
                x /= a * a;
                v *= a;
            }
            else if (x % a == 0)
            {
                x /= a;
                v *= a * a;
            }
        }

        res += cnt[v];
    }    

    res %= MOD;
    return res;
}

bool is_cubefree[LIM + 10];
int solve(int n, int k)
{
    if (k == 0) return 0;
    if (k == 1) return n;
    if (k == 2) return (1LL * n * (n - 1) / 2) % MOD;
    if (k == 3) return solveG(n, k);

    memset(is_cubefree, true, sizeof(is_cubefree[0]) * (n + 1));
    precal_nck(n);

    long long res = 0;
    for (int i = 1, j; i <= n; ++i) if (is_cubefree[i]) 
    {
        for (j = 1; i * j * j * j<= n; ++j)
        {
            is_cubefree[i * j * j * j] = false;
        }

        res += nck(j - 1, k);
    }

    res %= MOD;
    return res;
}

vector<int> divisor[LIM];
int solveG(int n, int k)
{
    linear_sieve(n);
    for (int i = 1; i <= n; ++i)
        divisor[i].clear();
 
    for (int j = 1; j <= n; ++j)
    {
        int x = 1;
        for (int t = j; t > 1; )
        {
            if (x > n) break;
            int a = lpf[t];
            ll b = 1LL * a * a * a;
            while (t % b == 0)
            {
                t /= b;
                if (1LL * x * a > n) goto skip;
                x *= a;
            }
 
            if (t % a == 0)
            {
                if (1LL * x * a > n) goto skip;
                x *= a;
                do t /= a; while (t % a == 0);
            }
        }
 
        for (int i = x; i <= n; i += x)
            divisor[i].push_back(j);
 
        skip:{};
    }
 
    int res = 0;
    for (int i = 1; i <= n; ++i)
    {
        ll t = 1LL * i * i * i;
        for (int x = 0; x < divisor[i].size(); ++x)
        {
            for (int y = x + 1; y < divisor[i].size(); ++y)
            {
                int a = divisor[i][x];
                int b = divisor[i][y];
                int c = t / (a * b);
                res += (n >= c && c > b && 1LL * a * b * c == 1LL * i * i * i);
            }
        }
    }
 
    return res;
}
 
int res[LIM + 10];
int solve(int n, int k)
{
    if (k == 0) return 0;
    if (k == 1) return n;
    if (k == 2) return (1LL * n * (n - 1) / 2) % MOD;
    if (k == 3) return solveG(n, k);
    int t = ceil(cbrt(n) + 1) + 1;
    linear_sieve(t);
    precal_nck(t);
    precal_div(t);

    long long res = 0;
    for (int d = 1; d * d * d <= n; ++d) /// For each fixed p^2
    {
        long long sum = 0;
        for (int p : divisors[d]) /// For each (p | d)
            sum += mu[d / p] * nck(p - 1, k - 1);

        sum %= MOD;
        res += sum * (n / (d * d * d));
    }

    res %= MOD;
}

int ceil_sqrt(ll x)
{
    int t = sqrt(x);
    while (1LL * t * t > x) --t;
    while (1LL * t * t < x) ++t;
    return t;
}

#define sz(x) int((x).size())
#define all(x) (x).begin(), (x).end()
#define rall(x) (x).rbegin(), (x).rend()
#define lb(x, v) lower_bound(all(x), v) - (x).begin()
#define ub(x, v) upper_bound(all(x), v) - (x).begin()

ll scm[LIM];
vector<int> divisor[LIM];
int solveG(int n, int k)
{
    linear_sieve(n);
    fill_n(scm, n + 1, 1);
    for (int i = 1; i <= n; ++i)
        divisor[i].clear();
        
    for (int p : prime)
    {
        ll q = 1LL * p * p * p;
        for (int t = p; ; t *= q)
        {
            for (int i = t; i <= n; i += t)
                scm[i] = (scm[i] > n / p) ? n + 1 : scm[i] * p;
            
            if (1LL * t > n / q) break;
        }
    }

    for (int j = 1; j <= n; ++j)
        for (int i = scm[j]; i <= n; i += scm[j])
            divisor[i].push_back(j);
    
    int res = 0;
    for (int i = 1; i <= n; ++i)
    {
        ll t = 1LL * i * i * i;
        for (int x = 0; x < divisor[i].size(); ++x)
        {
            int a = divisor[i][x];
            int v = ceil_sqrt(t / a);
            int y0 = (divisor[i].back() < v) ? sz(divisor[i]) : lb(divisor[i], v);
            for (int y = y0 - 1; y > x; --y)
            {
                int b = divisor[i][y];
                int c = t / (a * b);
                if (c > n) break;
                res += (1LL * a * b * c == t);
            }
        }
    }

    return res;
}

int res[LIM + 10];
int solve(int n, int k)
{
    if (k == 0) return 0;
    if (k == 1) return n;
    if (k == 2) return (1LL * n * (n - 1) / 2) % MOD;
    if (k == 3) return solveG(n, k);
    int t = ceil(sqrt(n) + 1) + 1;
    linear_sieve(t);
    precal_nck(t);

    memset(res, 0, sizeof(res[0]) * (t + 1));
    for (int d = 1; d * d * d <= n; ++d) 
        res[d] = nck(d - 1, k - 1);

    for (int p : prime)
        for (int d = t / p; d > 0; --d)
            res[d * p] -= res[d];

    long long ans = 0;
    for (int d = 1; d * d * d <= n; ++d)
        ans += res[d] * (n / (d * d * d));

    ans %= MOD;
    return ans;
}

The complexity is obviously $$$O(1)$$$ and we cant optimize further more lol.

In the first implementation, we use a cubefree sieve in $$$O(n)$$$ and precalculatings for bionomial coefficients is also $$$O(n)$$$.

In the second and third implementation. We just apply the idea we used in task B.

In the first implementation as you must using factorization, the cost is $$$O(\log n)$$$ for each number, hence you got that complexity.

Yet for the second implementations things gone crazy.

It has unprovable complexity, though I tried to search for papers, and blogs, even with the help of some GMs I cant find nothing good enough to claim its real complexity.

Yet it might be not the real complexity but also kind of an illusion assumptioning high constant and faking its real higher complexity.

And the third implementation is also hard to prove