Efficient way to reduce a binomial coefficient as a fraction

Question

Here is the full problem.

You need to calculate Euler's totient function of a binomial coefficient $C_n^k$.

Input

The first line contains two integers: $n$ and $k$ $(0 \le k \le n \le 500000)$.

Output

Print one number $\varphi (C_n^k)$ modulo $10^9+7$.

My thoughts:

It is known that $$\varphi(a)=a \prod_{p|a}(1-\frac{1}{p}) $$ where $p$ are prime numbers divide $a$.

Hence, if we can obtain somehow vector<int> multipliers that contains divisors of $C_n^k$ then we can easily do the following steps in order to calculate $\varphi(C_n^k)$:

1. Multiply all elements of that vector modulo $10^9+7$. Let's call the result by result
2. Then we can iterate through all prime numbers that divide any element of multipliers(these prime numbers can be obtained by a minor modification of sieve of Eratosthenes). Since $1-\frac{1}{p}=\frac{p-1}{p}$ we can update the result by:

result = divideMod(multiplyMod(result, p-1), p)

where divideMod and multiplyMod are functions doing corresponding operations modulo $10^9+7$.

And yes, we can do modulus division since $10^9+7$ is prime.

By doing all that stuff we get what we needed: $\varphi(C_n^k)$ modulo $10^9+7$. This all idea now requires just a vector multipliers. Here is my attempt to get it:

I need to write a function calculates the combinations number $C_n^k$. The function shouldn't return the total result of the operation(because it can be too large since $(0 \le k \le n \le 500000)$). It should return the vector<int> which contains divisors of that number. Let's do some math:

$$ C_n^k = \frac{n!}{(n-k)! k!} \\ =\frac{n(n-1)(n-2)...(n-k+1)}{k(k-1)(k-2)...1} $$

So now I need to reduce this fraction. And the question is: what is the most efficient way to do this(in terms of time)?

I've tried the following. Consider the numerator and denominator are represented by vector<int> numerator={n, n-1, ..., n-k+1} and vector<int> denominator={k, k-1, ..., 1} respectively.

    vector<long> numerator(k);
    vector<long> denominator(k);
    for (int i = 0; i<k; i++) {
        numerator[i] = n-i;
        denominator[i] = k-i;
    }

    vector<long> multipliers;
    for (int i = 0; i < k; i++) {
        for (int j = 0; j < k; j++) {
            if (numerator[i] == 1) 
                break;
            
            long greatest_common_divisor = gcd(numerator[i], denominator[j]);
            numerator[i] /= greatest_common_divisor;
            denominator[j] /= greatest_common_divisor;
        }
        if (numerator[i] != 1) 
            multipliers.push_back(numerator[i]);
    }

As you can see I just go through all numbers in numerator and denominator and divide them by their greatest common divisor.

Time complexity of this algorithm is $O( k^2 log(nk) )$

It's too big and for this solution contest system returns time limit exceeded.($0 \le k \le n \le 500000$)

Does there exist more efficient way?

Answer 1

The formula $\varphi(a)=a \prod_{p|a}(1-\frac{1}{p})$ tells us to approach the problem by the prime factors.

Here is another useful formula.

Legendre's formula. For any prime number $p$ and any positive integer $n$, let $\nu _{p}(n)$ be the exponent of the largest power of $p$ that divides $n!$, i.e, $p^{\nu_{p}(n)}$ divides $n!$ but $p^{\nu _{p}(n)+1}$ does not divide $n!$. We have, $$\nu _{p}(n)=\lfloor\frac np\rfloor + \lfloor\frac n{p^2}\rfloor + \lfloor\frac n{p^3}\rfloor + \cdots,$$ where the ellipsis means the addition go on until the term becomes 0. (For a proof, check that Wikipedia link.)

Applying Legendre's formula, we see that $n!=\prod_{\text{prime } p\le n} p^{\nu_p(n)}$ for all $n$. Since $C_n^k = \frac{n!}{(n-k)! k!}$, we have $C_n^k=\prod_{\text{prime } p\le n}p^{\nu_p(n)-\nu_p(k)-\nu_p(n-k)}$.

---

Here is the outline of the algorithm.

1. Find all primes numbers not larger than $n$.
2. Initialize $answer$ to 1. Iterate over all prime p not larger than $n$.
   1. Compute $e = \nu_p(n)-\nu_p(k)-\nu_p(n-k)$.
   2. If $e\ge1$, replace $answer$ by $answer * p^{e-1} * (p-1) \pmod{10^9+7}$
3. return $answer$.

The complexity of the algorithm is $O(n\log n)$ time and $O(n)$ space.

---

Here is code in Java (which is almost valid c/c++ code as well). It takes less than a hundredth of a second to compute $C_{500000}^{250000}$ on my computer.

final static int MOD = 1000000007;

static long totientOfBinomialCoefficient(final int n, final int k) {
    boolean[] isComposite = new boolean[n + 1];
    for (int i = 2; i <= n; i++) {
        if (!isComposite[i]) {
            for (int j = 2 * i; j <= n; j += i) {
                isComposite[j] = true;
            }
        }
    }

    long answer = 1;
    for (int i = 2; i <= n; i++) {
        if (!isComposite[i]) {
            int exp = exponentInFactorial(i, n) - exponentInFactorial(i, k) - exponentInFactorial(i, n - k);
            if (exp > 0) {
                answer = answer * powerMod(i, exp - 1) % MOD * (i - 1) % MOD;
            }
        }
    }

    return answer;
}

static int exponentInFactorial(int prime, int n) {
    int answer = 0;
    while (n >= prime) {
        n /= prime;
        answer += n;
    }
    return answer;
}

static long powerMod(int base, int exp) {
    long answer = 1;
    while (exp >= 1) {
        answer = answer * base % MOD;
        exp--;
    }

    return answer;
}

Some further micro-optimization can be done. However, they are not needed for this problem.

Answer 2

You store the numbers n-k+1 to n in an array.

Then for each prime number p ≤ k: Find which power of p is a factor of k! (Thats k/p + k/p^2 + k/p^3 ... ) Then remove that power from the array: Find the first number divisible by p (that would be the number at index 0 if n-k+1 is divisible by p, otherwise at index p - ((n-k+1) modulo p)). That number is divisible by p, possible by p^2 etc., and the next number divisible by p is at the index p higher).

You are now left with an array of k numbers that you need to multiply.