FOMBINATORIAL - Editorial

PROBLEM LINK

contest

practice

Author: Devendra Agarwal

Tester: Sergey Kulik

Editorialist: Florin Chirica

DIFFICULTY

easy-medium

PREREQUISITES

number theory

QUICK EXPLANATION

We need to calculate f(n) / f(r) * f(n - r) mod M. If f(r) and f(n - r) would be coprime with M, then the problem could be solved by modular inverses (extended Euclid’s algorithm). Another approach is to take each prime factor P from 2 to n and for each prime factor to calculate x = maximum x such as P^x is a factor of f(n), y = maximum y such as P^y is a factor of f(r) and z = maximum z such as P^z is a factor of f(n - r). The result will be product of P ^ (x - y - z), for each P from 2 to n. The solution is to calculate this only on prime divisors of M. If we ignore them and calculate the expression without them, it will be coprime with M, so we can use again modular inverse.

EXPLANATION

Particular case: f(x) is coprime with M

Let’s assume for a moment that, for every 1 <= i <= N, f(i) is coprime with M. This means every f(i) allows a modular inverse. A modular inverse can be calculated by Extended Euclid’s algorithm. A modular inverse is like division for modulo M. This means, a / b mod M is equivalent by a * ModularInverse(b) mod M. We can calculate f(N) / (f(r) * f(N - r)) by calculating (f(N) * ModularInverse(f(r)) * ModularInverse(f(N - r))) mod M. Modular inverse of a number x can be calculated in O(log(x)). We can precalculate f(i) in O(N * log(N)) time: f(i) = f(i - 1) * i ^ i mod M. i ^ i can be calculated by exponentiation by squaring, since the multiplication operation is associative.

General case

The approach works only if gcd(f(i), M) = 1. Gcd between two number is defined as common factors of both numbers at the smallest power. Let M = p1 ^ q1 * p2 ^ q2 * … for gcd(f(i), M) = 1 it means that all factors p1, p2, … need to appear at power 0 in f(i). Otherwise, it a factor pi has a power at least 1, it will also have a power at least 1 in M (otherwise it won’t be in its prime factorization) so gcd(f(i), M) is at least pi.

From here we get the following approach: we define g(i) = i^i from which we erase all factors of M. For example, let M = 6 and N = 10

g(1) = 1

g(2) = 1 ( Reason removed all factors of 2)

g(3) = 1 (Reason removed all factors of 3)

g(4) = 1 (removed all factors of 2)

g(5) = 5^5

g(6) = 1

g(7) = 7^7

g(8) = 1

g(9) = 1

g(10) = 5^10 (NOTE THIS)

Let’s define g’(x) = g’(x - 1) * g(x). Now g’(x) is coprime with M. This means you can apply the approach from above.

However, the result is not complete yet. So far we have g’(N) * ModularInverse(g’(N - r)) * ModularInverse(g’(r)), which is the answer, EXCEPT all prime factors for M. We need to multiply this number by p1 ^ (x1 - y1 - z1) * p2 ^ (x2 - y2 - z2) * …

xi is the largest power such as pi ^ xi is a divisor of f(N).
yi is the largest power such as pi ^ yi is a divisor of f(r).
zi is the largest power such as pi ^ zi is a divisor of f(N - r).

By doing this, all factors of M are multiplied to the result and we’re done. The number of factors of M is small (worst case 2 * 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23 * 27 > 10^9). So we can afford to process each factor of M one by one. For a factor P, for each f(i), we need to calculate maximum x such as P ^ x is a divisor of f(i). The power P ^ x modulo M can be calculated then by exponentiation by squaring, too.

Number x for 1^1 * 2^2 * … * i^i is calculated in the following way: we can write each y^y number as some P^y * another number. The whole product will be P^y1 * foo * P^y2 * bar * … = P^y1 * P^y2 * … * P^yi * some_number = P^(y1+y2+…+yi) * some_number. By yi I denote maximum number such as P^yi divides only i ^ i term.

How to get P^yi such as yi is the power of P in prime decomposition of i ^ i? i can be written as a product of prime factors. We’re interested only in power of factor P. So i = (P ^ power * other_numbers) ^ i. We get that i = P ^ (power * i) * other_numbers. Hence, yi = power * i. To get power, we simply divide i by P as much as possible (of course, we need to have remainder equal to zero) and we increment power each time.

This problem allows multiple approaches. If you have a different solution, please share it with us!

Complexity

Let cnt the number of prime factors of M. Precalculation step takes O(cnt * n). For a query, modular inverse stuff takes O(log(n)) and using precalculated info by exponentiation by squaring also takes O(log(n)) time.

AUTHOR’S AND TESTER’S SOLUTIONS:

To be updated soon
To be updated soon

11 Likes

Latex support should be added to this platform.

16 Likes

I had a lot of “fun” with that problem.

I thought that for 40 points I just need to implement modular inversion. I read question on forum, that condition about M beeing prime is not enough to use it, but my tests shown different…

…and guess what, my test were wrong… I implemented brute force using BitInteger in Java to calculate up to N = 300 testing for all rs, but still showing WA. The problem was, that for my random test - N = 10, 100 and 300 (all rs in 1…N-1 and primes 2, 3, 4 and 109+9) it was working, but not for N = 6 and r = 2 (I found that after 2 days), funny isn’t it…

For 40 points I implemented solution using polynomes.

Where Element class represents xn. There are all necessary functions like multiplication of two polinomes (merge in my code), polynomen pow and so on…

7 Likes

Well I tried a different segment tree approach But failed

Here’s my code CodeChef: Practical coding for everyone

Here is my logic, which does not require much maths -

Logic ANS(N,R) = [ ANS(N, R-1) * X ] / Y

where

X = (N-R+1)^(N-R+1)

Y = R^R

Now the ans at any stage is just the product of some prime powers.

With each r, at most log(VAL) primes have change in their powers

This change in power can be -ve [when considering denominator] or +ve

Hence we maintain a segment tree of powers of each prime

Change is just a range update

For each test case complexity then is just

->O(N) (For each r)

->logN for prime factorization

->logN for range update + log(Large number) for binary exponentiation

Hence net = N * logN * logX for each test case. Query is answered in O(1)

7 Likes

I couldn’t get the solution … can somebody explain it simply giving some example of a value of n,r (n larger than 500) .

1 Like

Don’t like the easy-medium tag of the problem. Should come under medium-hard.

6 Likes

Here is my solution:

First factorize M = p1 ^a1 * … * pk ^ ak using some short precoputed sieve. Each pi is prime number.

For each i = 1…k compute pairs (p^b[i][1], Fi(1)) … (p^b[i][N], Fi(N)) such that Fi(j) is taken modulo pi^ai, Fi(j) contains no pi factors and pi^b[i][j] are these missing factors.

Now suppose we have a query (n, r). We will compute G(n, r, i) = F(n) / (F® * F(n-r)) mod pi^ai for each i = 1…k. So:

G(n, r, i) = p^(b[i][N] - b[i][r] - b[i][N-r]) * Fi(N) * revMod(Fi®, pi^ai) * revMod(Fi(N-r), pi^ai).
Observe that now we can use modular reverse because each Fi(m) is coprime to pi^ai.

Having G(n, r, i) for all i = 1…k enables us compute F(n) / (F® * F(n-r)) mod M. Observe that each pi^ai are pairwise coprime so we can use chinese remainder theorem to compute F(n) / (F® * F(n-r)) mod M.

I used exponentiation by squaring for fast modular powering. I also used Fermat’s little theorem a^(p-1) = 1 mod p to reduce exponent when I could i.e a^k mod p = a^(k mod p-1) mod p provided p is prime.

3 Likes

A more intuitive approach:

If M had been prime-only, the problem was fairly easy. Just take the modulo inverse of the denominator as you do for nCr and multiply with the numerator. For this problem, you can keep two lists. the first list contains, for each i, those primes which appear in f[i] and are NOT co-prime with M. You also store the exponent of the prime along with it. For each i, this wont be more than 10. In the other list, you keep, for each i, the product of all the primes which appear in f[i] and are co-prime to M. when you take the product, you take the primes along with their exponents (ie, the number of times they appear in f[i]). Now, for the answer, you can use extended euclidean algorithm for finding modular inverse of the co-prime product. For those primes which ain’t co-prime with M, you just find their final powers in the expression and multiply with the previous term (the co-prime one). This works in NlogN.

Link: http://www.codechef.com/viewplaintext/5417817

4 Likes

I have implemented different segment tree approach for this problem .
accepted solution

@MarioYC thank you so much for your solution, I just read your code using segment trees, every time you update the segment tree using the Update_exp() method, the array exp is storing the exponents but why aren’t you initializing them back to zero after every function call to update_exp()??

if there are say 3 number 2, 8, 16 so for 16 exp[2]=1(for 2)+3(for 8)+4(for 16) implies= exp[2]=8, so you’re passing a call to mod_pow(2,8).

this is one part I am not able to understand!!

I have one question about the editorial. Does it have the situation of x - y - z < 0 for p^(x-y-z)?

Links to contest and practice are missing.

Please can you be clear on what do you want to convey by writing this ^ ?

Few things, like random tests (even when lot of them) can really miss some important cases. Do you know how many time I saw: “I tested for all inputs…” complains. And I shared different approach for getting 40 points, maybe someone will be interested. Not you? Ok, I have no problem with that…

I am sorry if I sounded rude and I am interested I just didn’t get what you were trying to say initially

1 Like

I agree. I read it on the first day, thought of using modular inverse using extended Euclidean but it had little problems when M is not prime, so used Python :wink:
It should be under "medium " tag.

2 Likes

No problem at all, I received several down votes, so it’s clear it was difficult to understand what I wanted to say…

If you describe your approach in bigger detail, someone will help you, I’m sure :wink:

I also made an implementation with segment trees (keeping exponents in an array and the products of the powers of primers in the segment tree) and got AC, you can check it:
http://www.codechef.com/viewsolution/5346679