INVBINCF - EditorialEditorial (CONSIDERABLY IMPROVED!!!)
#PROBLEM LINK:
[Practice][111]
[Contest][222]
[Practice][practice]
[Contest][contest]
**Author:** [Anton Lunyov][6666] Lunyov][author]
**Tester:** [Hiroto Sekido][5555] Sekido][tester]
**Editorialist:** [Anton Lunyov][6666]
Lunyov][author]
#DIFFICULTY:
HARD
#PREREQUISITES:
Binomial coefficients modulo prime powers, 128-bit arithmetic
#PROBLEM:
For the given **n** and **R** such that **0 ≤ R < 2<sup>n</sup>** you need to find the smallest **K** such that **C(2<sup>n</sup> − 1, K) mod 2<sup>n</sup> = R**,
where **C(N, K) = N! / K! / (N − K)!**, or report that such **K** does not exist. Note, that **n** could be up to **120**.
#QUICK EXPLANATION:
#EXPLANATION:
At first note that the problem requires to code 128-bit arithmetic. This can be made by representing each number as a pair of 64-bit values (lowest 64 bits and highest 64 bits) and then some routine job should be made to code all necessary operations. The only non-trivial one is multiplication. It can be made using only six 64-bit multiplications (refer to the [tester's solution][tester_sol]). Actually only five is enough (refer to the [author's solution][author_sol2]).
The first thing to noticed is notice is:
**C(2<sup>n</sup> − 1, K) = F<sub>2</sub>(2<sup>n</sup> fact2(2<sup>n</sup> − 1) / F<sub>2</sub>(K) fact2(K) / F<sub>2</sub>(2<sup>n</sup> fact2(2<sup>n</sup> − 1 − k) K) * C(2<sup>n−1</sup> − 1, K div 2)**, (1)
where **F<sub>2</sub>(N)** **fact2(N)** is the product of all odd numbers up to **N**. (If some one want to see the proof ask me in comment I and will provide it when I will have a free time.)
It immediately follows from here that **C(2<sup>n</sup> − 1, K)** is always odd. So if **R** is even then the answer is **−1**.
We will prove that for each odd **R** in the range **{0, 1, 2, ..., 2<sup>n</sup> − 1}**,
there exists unique **K** in the range **{0, 1, 2, ..., 2<sup>n−1</sup> − 1}** such that **C(2<sup>n</sup> − 1, K) mod 2<sup>n</sup> = R**.
At first we infer the following congruence from the above equality: (1):
**C(2<sup>n</sup> − 1, 2 * X K + Y) X) ≡ (−1)<sup>X (−1)<sup>K + Y</sup> X</sup> * C(2<sup>n−1</sup> − 1, X) K) (mod 2<sup>n</sup>)**,
where **X = **X** is **0** or **1**.
(Again if some one want to see the proof ask me in comment I and will provide it when I will have a free time.)
Also we will need the following property:
**C(2<sup>n</sup> − 1, 2 * K) + C(2<sup>n</sup> − 1, 2 * K div + 1) ≡ 2<sup>n</sup> (mod 2<sup>n+1</sup>)**, for **n ≥ 2** and **Y = K mod 2**.
**0 ≤ S < 2<sup>n−1</sup>**.
Then after some non-trivial analysis we get the following simple iterative routine that will find the answer:
invbincf(n, INVBINCF(n, R)
K = (R mod 4) div 2
for j = 3 to n do
par = K mod 2
if C(2^(j-1) - 1, K) mod 2^j != R mod 2^j
K = K xor 1
K = 2 * K + par
return K
So The proof in Russian is contained [here (pages 57-58)][sbornik] and contains a lot of formulas. If you really want to see the proof in English, ask me in comment and maybe some day I will translate it :)
Thus, the problem boils down to effective calculation of **C(2<sup>n</sup> − 1, K)**. K)** modulo powers of two.
And here we ask for help of British mathematician [Andrew Granville][555], Granville][granville], the number theory expert :)
We refer to Theorem 2 of [his article][666] article][granville_paper] about binomial coefficients modulo prime powers (see end of the page 3).
For **p = 2** it sounds as follows
**F<sub>2</sub>(2 follows:
**fact2(2 * u) ≡ sgn * Product[ F<sub>2</sub>(2 fact2(2 * j)<sup>b(r, j, u)</sup> : 1 ≤ j ≤ r] (mod 2<sup>2 * r + 1</sup>)**, for **r ≥ 2**, (2) (2)
where
**b(r, j, u) = u / j * Product[ (u<sup>2</sup> − i<sup>2</sup>) / (j<sup>2</sup> − i<sup>2</sup>) : 1 ≤ i ≤ r, i ≠ j]**,
and **sgn = 1** **sgn** is **1** or **−1** and can be found by comparing residues of LHS and RHS modulo **4**.
Clearly **b(r, j, u)** could be calculated modulo **2<sup>2 * r − 1</sup>** in view of the property:
**a<sup>2<sup>n − 2</sup></sup> = 1 (mod 2<sup>n</sup>)** for **n ≥ 3** and odd **a**.
The proof of (2) is based on 2-adic exponent and logarithm and is like magic for me.
But straightforward implementation of **C(2<sup>n</sup> − 1, K)** using formulas (1) and (2), together with `invbincf` `INVBINCF` routine above leads to quite slow solution.
Note, that actually we can maintain value **C(2<sup>j</sup> − 1, K) mod 2<sup>n</sup>** in the loop of `invbincf` `INVBINCF` and recalculate it using formula (1). Namely, if **K = K xor 1** is performed we actually replace **K** by its neighbor adjacent integer and can update **C(2<sup>j</sup> − 1, K)** straightforwardly using one multiplication and one division.
While when we do **K = 2 * K + par** we should transform **C(2<sup>j</sup> − 1, K)** to **C(2<sup>j + 1</sup> − 1, 2 * K + par)** and this exactly can be made by formula (1) by calculating three values of **F<sub>2</sub>(X)** **fact2(X)** by formula (2).
The next step is to (2). Also we can maintain numerator and denominator of **C(2<sup>j</sup> − 1, K)** separately and get rid of divisions by rewriting check in `invbincf` `INVBINCF` in obvious way.
way. Refer to routine `invbin` in [author's solution][author_sol2] for clear implementation of this.
So `invbincf` `INVBINCF` is ready to be used but we need to figure it out how to calculate **F<sub>2</sub>(X) **fact2(X) mod 2<sup>n</sup>** efficiently. efficiently.
At first to deal with we should calculate **b(j, r, u)** we should maintain both numerator and denumerator u) mod 2<sup>n−2</sup>** efficiently. Here main complication is that denominator could have even factors and we can't directly apply technique of inverse residues. To handle this we will use extended representation of integers in the form **2<sup>a</sup> * b**...
b**, where **b** is odd. So we could create a data structure that deals with this representation. Involving negative **a** and inverse residue technique for odd part of the number, we can handle non-integer numbers as well. Thus, to calculate **b(r, j, u)** we use the following plan:
* Before processing test data we precalculate values **pseudo_fact[r][j] = 1 / j * Product[ 1 / (j<sup>2</sup> − i<sup>2</sup>) : 1 ≤ i ≤ r, i ≠ j]** in extended format. They equal to inverse of denominators of **b(r, j, u)**.
* Inside `fact2` routine we calculate values **pref[k] = u * Product[ (u<sup>2</sup> − i<sup>2</sup>) : 1 ≤ i ≤ k]** and **suf[k] = Product[ (u<sup>2</sup> − i<sup>2</sup>) : k ≤ i ≤ r]** in two simple **O(r)** loops (again in extended format).
* Then **b(r, j, u) = pref[j - 1] * suf[j + 1] * pseudo_fact[r][j]** and after this, we can transform it to usual integer representation.
Thus, the complexity of precalc is **O(N<sup>3</sup>)** (we calculate **O(N<sup>2</sup>)** values **pseudo_fact[r][j]** and each step requires finding inverse residue, which can be made in **O(N)** time) and calculation of all values **b(r, j, u)** for the given **u** is only **O(N)**. Actually, saving inverses for numbers up to **N** in **O(N<sup>2</sup>)** time we can speed up the precalc part to **O(N<sup>2</sup>)**, but even **O(N<sup>3</sup>)** precalc is many times faster then some tricky precalc below.
Now having values **b(r, j, u)** calculated we can find **fact2(u)** using **r** `PowerMod` operations, where `PowerMod[a, b, n]` returns **a<sup>b</sup> mod 2<sup>n</sup>**. But even using exponentiation by squaring this will lead to **O(N<sup>2</sup>)** complexity for the `fact2` routine and thus to **O(N<sup>3</sup>)** complexity for the `INVBINCF` routine, which is still too slow.
To speed up this part we can do the following. Note that we need to calculate powers of small fact2-values in `fact2` routine. So we can do some clever precalc for them. Namely, if we want to set up fast calculation of powers of some number **A** we can precalc the following its powers in advance: **powA[k][j] = A<sup>j * 2<sup>H * k</sup></sup>** for all **0 ≤ k < N/H** and **0 ≤ j < 2<sup>H</sup>**, where **H** is some constant (divisor of **N = 120** will be most convenient). Then, when we want to calculate **A<sup>B</sup>** for some **B < 2<sup>N</sup>**, we represent **B** in **2<sup>H</sup>**-ary system as
**B = B<sub>0</sub> + B<sub>1</sub> * 2<sup>H</sup> + ... + B<sub>K</sub> * 2<sup>H * K</sup>**
and calculate **A<sup>B</sup>** as the product of known powers **powA[0][B<sub>0</sub>] * powA[1][B<sub>1</sub>] * ... * powA[K][B<sub>K</sub>]** in **O(N / H)** time.
Thus precalculating such powers for each **fact2(2 * j)** (**1 ≤ j ≤ 60**) with constant **H = 15**, will speed up our solution in **H = 15** times to **O(N * N * N / H)**. But we can do even better. We can find prime factorization of **fact2(u)** and set up powers precalc only for odd primes up to **120** (there are only **29** of them). Then after calculation of **b(r, j, u)** we can do some simple **O(N)** loop to calculate degrees in factorization of **fact2(u)** and then calculate the answer as the product of prime powers in at most **29 * 8** multiplications. This is quite cumbersome stuff and I will provide explanation later. Currently refer to the [tester's solution][444].
But overall trick also speeds up powers precalc in two times.
Overall complexity of the solution assuming constant time for each 128bit 128-bit operation is
**O(π(N) **O(N<sup>3</sup> + π(N) * 2<sup>H</sup> * N/H + N<sup>3</sup> T * N * (N + T * N * π(N) * N/H)**, N/H))**,
where I choose **N = 120**, **H = 16** 15** and **π(N)** is the number of primes **≤ N**.
So as you may guess, I Here:
* term **N<sup>3</sup>** corresponds to precalc powers of the form **p<sup>2<sup>H * k</sup> * j</sup>**
of **pseudo_fact[r][j]**;
* term **π(N) * 2<sup>H</sup> * N/H** corresponds to precalc of prime powers;
* term **T * N * (N + π(N) * N/H)** means that for all primes **p ≤ N**, **0 ≤ k < N/H**, **0 ≤ j < 2<sup>H</sup>**.
This allows to each of **T** tests we calculate **F<sub>2</sub>(X) mod 2<sup>n</sup>** in **O(π(n) * n/H)** **O(N)** values **fatc2(u)** in **O(N + π(N) * N/H)** time after we calculate all **b(j, r, u)**. For this we find the factorization of **F<sub>2</sub>(X)** and then calculate each prime power in **O(n/H)** using precalculated prime powers.
Also the problem requires to code 128-bit arithmetic. This can be made by representing each number as pair of 64-bit values (lowest 64 bits and highest 64 bits) and then some routine job should be made to code all necessary operations. The only non-trivial one is multiplication. It can be made using only six 64-bit multiplications (refer to the [tester's solution][444]). Actually only five is enough (it will be demonstrated in the [author's solution]).
each.
You may also ask "Why did I annoy all by crappy 128-bit arithmetic?". This is because I almost invented more elementary solution with **O(2<sup>n/3</sup>)** **O(2<sup>N/3</sup>)** precalc and **poly(n)** **poly(N)** routine for each test. So I decided that staying under 64-bit arithmetic could be dangerous. Actually I spent the whole day to code the working 128-bit version having working 64-bit version. So this annoyed me too :)
#EXPLANATION:
Will be provided soon.
#ALTERNATIVE SOLUTION:
@winger solution is very interesting and differ differs from described one.
It would be nice if he will describe his solution.
Solution by @mugurelionut is based on more elementary math than Granville formula. You can read about it [here](http://discuss.ww2.codechef.com/answer_link/8420/).
#AUTHOR'S AND TESTER'S SOLUTIONS:
Author's solution will can be provided soon. found [here][author_sol2].
Tester's solution can be found [here][444]. [here][tester_sol].
#RELATED PROBLEMS:
[e-olimp - 322 - Binomial Coefficients 5][1111] 5][eolimp322]
The current problem is extreme version of this one.
I've set this problem (for **n ≤ 40**) at Kharkiv Winter Programming School 2009 and editorial in Russian is available [here][999] [here][sbornik] (see page 57). So those who knew about this problem could have some advantage against others, but it seems that all who solved the problem was unaware of this link :)
The good problems to practice your skills on direct problems on calculating reduced factorials modulo prime powers:
[Mathalon - 141 - Factorial trailing digits2][777] digits2][mathalon141]
[Mathalon - 144 - Binomial Coefficient 1][888] 1][mathalon144] (this one is mine :P)
[111]: [practice]: http://www.codechef.com/problems/INVBINCF
[222]: [contest]: http://www.codechef.com/APRIL13/problems/INVBINCF
[333]:
[author]: http://www.codechef.com/users/anton_lunyov
[author_sol]: http://www.codechef.com/download/Solutions/2013/April/Setter/INVBINCF.cpp
[444]: [author_sol2]: http://ww2.codechef.com/viewsolution/2071520
[tester]: http://www.codechef.com/users/laycurse
[tester_sol]: http://www.codechef.com/download/Solutions/2013/April/Tester/INVBINCF.cpp
[555]:
[granville]: http://en.wikipedia.org/wiki/Andrew_Granville
[666]: [granville_paper]: http://www.dms.umontreal.ca/~andrew/PDF/BinCoeff.pdf
[777]: [mathalon141]: http://mathalon.in/Mathalon/?page=show_problem.php&pid=141
[888]: [mathalon144]: http://mathalon.in/Mathalon/?page=show_problem.php&pid=144
[999]: [sbornik]: http://ejudge.kture.kharkov.ua/media/sbornik/Sbornik2009.pdf
[1111]: [eolimp322]: http://www.e-olimp.com.ua/en/problems/322
[5555]: http://www.codechef.com/users/laycurse
[6666]: http://www.codechef.com/users/anton_lunyov