Answers to: INVBINCF - Editorialhttps://discuss.codechef.com/questions/8154/invbincf-editorial<h1>PROBLEM LINK:</h1>
<p><a href="http://www.codechef.com/problems/INVBINCF">Practice</a><br>
<a href="http://www.codechef.com/APRIL13/problems/INVBINCF">Contest</a></p>
<p><strong>Author:</strong> <a href="http://www.codechef.com/users/anton_lunyov">Anton Lunyov</a> <br>
<strong>Tester:</strong> <a href="http://www.codechef.com/users/laycurse">Hiroto Sekido</a><br>
<strong>Editorialist:</strong> <a href="http://www.codechef.com/users/anton_lunyov">Anton Lunyov</a></p>
<h1>DIFFICULTY:</h1>
<p>HARD</p>
<h1>PREREQUISITES:</h1>
<p>Binomial coefficients modulo prime powers, 128-bit arithmetic</p>
<h1>PROBLEM:</h1>
<p>For the given <strong>n</strong> and <strong>R</strong> such that <strong>0 ≤ R < 2<sup>n</sup></strong> you need to find the smallest <strong>K</strong> such that <strong>C(2<sup>n</sup> − 1, K) mod 2<sup>n</sup> = R</strong>,<br>
where <strong>C(N, K) = N! / K! / (N − K)!</strong>, or report that such <strong>K</strong> does not exist. Note, that <strong>n</strong> could be up to <strong>120</strong>.</p>
<h1>EXPLANATION:</h1>
<p>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 <a href="http://www.codechef.com/download/Solutions/2013/April/Tester/INVBINCF.cpp">tester's solution</a>). Actually only five is enough (refer to the <a href="http://ww2.codechef.com/viewsolution/2071520">author's solution</a>).</p>
<p>The first thing to notice is:<br>
<strong>C(2<sup>n</sup> − 1, K) = fact2(2<sup>n</sup> − 1) / fact2(K) / fact2(2<sup>n</sup> − 1 − K) * C(2<sup>n−1</sup> − 1, K div 2)</strong>, (1)<br>
where <strong>fact2(N)</strong> is the product of all odd numbers up to <strong>N</strong>. (If some one want to see the proof ask me in comment I and will provide it when I will have a free time.)</p>
<p>It immediately follows from here that <strong>C(2<sup>n</sup> − 1, K)</strong> is always odd. So if <strong>R</strong> is even then the answer is <strong>−1</strong>.<br>
We will prove that for each odd <strong>R</strong> in the range <strong>{0, 1, 2, ..., 2<sup>n</sup> − 1}</strong>,<br>
there exists unique <strong>K</strong> in the range <strong>{0, 1, 2, ..., 2<sup>n−1</sup> − 1}</strong> such that <strong>C(2<sup>n</sup> − 1, K) mod 2<sup>n</sup> = R</strong>.</p>
<p>At first we infer the following congruence from (1):<br>
<strong>C(2<sup>n</sup> − 1, 2 * K + X) ≡ (−1)<sup>K + X</sup> * C(2<sup>n−1</sup> − 1, K) (mod 2<sup>n</sup>)</strong>, where <strong>X</strong> is <strong>0</strong> or <strong>1</strong>.<br>
(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.)</p>
<p>Also we will need the following property:<br>
<strong>C(2<sup>n</sup> − 1, 2 * K) + C(2<sup>n</sup> − 1, 2 * K + 1) ≡ 2<sup>n</sup> (mod 2<sup>n+1</sup>)</strong>, for <strong>n ≥ 2</strong> and <strong>0 ≤ S < 2<sup>n−1</sup></strong>.</p>
<p>Then after some non-trivial analysis we get the following simple iterative routine that will find the answer:</p>
<pre><code>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
</code></pre>
<p>The proof in Russian is contained <a href="http://ejudge.kture.kharkov.ua/media/sbornik/Sbornik2009.pdf">here (pages 57-58)</a> 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 :)</p>
<p>Thus, the problem boils down to effective calculation of <strong>C(2<sup>n</sup> − 1, K)</strong> modulo powers of two.<br>
And here we ask for help of British mathematician <a href="http://en.wikipedia.org/wiki/Andrew_Granville">Andrew Granville</a>, the number theory expert :)<br>
We refer to Theorem 2 of <a href="http://www.dms.umontreal.ca/~andrew/PDF/BinCoeff.pdf">his article</a> about binomial coefficients modulo prime powers (see end of the page 3).<br>
For <strong>p = 2</strong> it sounds as follows:<br>
<strong>fact2(2 * u) ≡ sgn * Product[ fact2(2 * j)<sup>b(r, j, u)</sup> : 1 ≤ j ≤ r] (mod 2<sup>2 * r + 1</sup>)</strong>, for <strong>r ≥ 2</strong>, (2)<br>
where<br>
<strong>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]</strong>,<br>
and <strong>sgn</strong> is <strong>1</strong> or <strong>−1</strong> and can be found by comparing residues of LHS and RHS modulo <strong>4</strong>.<br>
Clearly <strong>b(r, j, u)</strong> could be calculated modulo <strong>2<sup>2 * r − 1</sup></strong> in view of the property:<br>
<strong>a<sup>2<sup>n − 2</sup></sup> = 1 (mod 2<sup>n</sup>)</strong> for <strong>n ≥ 3</strong> and odd <strong>a</strong>.<br>
The proof of (2) is based on 2-adic exponent and logarithm and is like magic for me.</p>
<p>But straightforward implementation of <strong>C(2<sup>n</sup> − 1, K)</strong> using formulas (1) and (2), together with <code>INVBINCF</code> routine above leads to quite slow solution.</p>
<p>Note, that actually we can maintain value <strong>C(2<sup>j</sup> − 1, K) mod 2<sup>n</sup></strong> in the loop of <code>INVBINCF</code> and recalculate it using formula (1). Namely, if <strong>K = K xor 1</strong> is performed we actually replace <strong>K</strong> by adjacent integer and can update <strong>C(2<sup>j</sup> − 1, K)</strong> straightforwardly using one multiplication and one division.<br>
While when we do <strong>K = 2 * K + par</strong> we should transform <strong>C(2<sup>j</sup> − 1, K)</strong> to <strong>C(2<sup>j + 1</sup> − 1, 2 * K + par)</strong> and this exactly can be made by formula (1) by calculating three values of <strong>fact2(X)</strong> by formula (2). Also we can maintain numerator and denominator of <strong>C(2<sup>j</sup> − 1, K)</strong> separately and get rid of divisions by rewriting check in <code>INVBINCF</code> in obvious way. Refer to routine <code>invbin</code> in <a href="http://ww2.codechef.com/viewsolution/2071520">author's solution</a> for clear implementation of this.</p>
<p>So <code>INVBINCF</code> is ready to be used but we need to figure it out how to calculate <strong>fact2(X) mod 2<sup>n</sup></strong> efficiently.</p>
<p>At first we should calculate <strong>b(j, r, u) mod 2<sup>n−2</sup></strong> 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 <strong>2<sup>a</sup> * b</strong>, where <strong>b</strong> is odd. So we could create a data structure that deals with this representation. Involving negative <strong>a</strong> and inverse residue technique for odd part of the number, we can handle non-integer numbers as well. Thus, to calculate <strong>b(r, j, u)</strong> we use the following plan:</p>
<ul>
<li>Before processing test data we precalculate values <strong>pseudo_fact[r][j] = 1 / j * Product[ 1 / (j<sup>2</sup> − i<sup>2</sup>) : 1 ≤ i ≤ r, i ≠ j]</strong> in extended format. They equal to inverse of denominators of <strong>b(r, j, u)</strong>.</li>
<li>Inside <code>fact2</code> routine we calculate values <strong>pref[k] = u * Product[ (u<sup>2</sup> − i<sup>2</sup>) : 1 ≤ i ≤ k]</strong> and <strong>suf[k] = Product[ (u<sup>2</sup> − i<sup>2</sup>) : k ≤ i ≤ r]</strong> in two simple <strong>O(r)</strong> loops (again in extended format).</li>
<li>Then <strong>b(r, j, u) = pref[j - 1] * suf[j + 1] * pseudo_fact[r][j]</strong> and after this, we can transform it to usual integer representation.</li>
</ul>
<p>Thus, the complexity of precalc is <strong>O(N<sup>3</sup>)</strong> (we calculate <strong>O(N<sup>2</sup>)</strong> values <strong>pseudo_fact[r][j]</strong> and each step requires finding inverse residue, which can be made in <strong>O(N)</strong> time) and calculation of all values <strong>b(r, j, u)</strong> for the given <strong>u</strong> is only <strong>O(N)</strong>. Actually, saving inverses for numbers up to <strong>N</strong> in <strong>O(N<sup>2</sup>)</strong> time we can speed up the precalc part to <strong>O(N<sup>2</sup>)</strong>, but even <strong>O(N<sup>3</sup>)</strong> precalc is many times faster then some tricky precalc below.</p>
<p>Now having values <strong>b(r, j, u)</strong> calculated we can find <strong>fact2(u)</strong> using <strong>r</strong> <code>PowerMod</code> operations, where <code>PowerMod[a, b, n]</code> returns <strong>a<sup>b</sup> mod 2<sup>n</sup></strong>. But even using exponentiation by squaring this will lead to <strong>O(N<sup>2</sup>)</strong> complexity for the <code>fact2</code> routine and thus to <strong>O(N<sup>3</sup>)</strong> complexity for the <code>INVBINCF</code> routine, which is still too slow.</p>
<p>To speed up this part we can do the following. Note that we need to calculate powers of small fact2-values in <code>fact2</code> routine. So we can do some clever precalc for them. Namely, if we want to set up fast calculation of powers of some number <strong>A</strong> we can precalc the following its powers in advance: <strong>powA[k][j] = A<sup>j * 2<sup>H * k</sup></sup></strong> for all <strong>0 ≤ k < N/H</strong> and <strong>0 ≤ j < 2<sup>H</sup></strong>, where <strong>H</strong> is some constant (divisor of <strong>N = 120</strong> will be most convenient). Then, when we want to calculate <strong>A<sup>B</sup></strong> for some <strong>B < 2<sup>N</sup></strong>, we represent <strong>B</strong> in <strong>2<sup>H</sup></strong>-ary system as<br>
<strong>B = B<sub>0</sub> + B<sub>1</sub> * 2<sup>H</sup> + ... + B<sub>K</sub> * 2<sup>H * K</sup></strong><br>
and calculate <strong>A<sup>B</sup></strong> as the product of known powers <strong>powA[0][B<sub>0</sub>] * powA[1][B<sub>1</sub>] * ... * powA[K][B<sub>K</sub>]</strong> in <strong>O(N / H)</strong> time.</p>
<p>Thus precalculating such powers for each <strong>fact2(2 * j)</strong> (<strong>1 ≤ j ≤ 60</strong>) with constant <strong>H = 15</strong>, will speed up our solution in <strong>H = 15</strong> times to <strong>O(N * N * N / H)</strong>. But we can do even better. We can find prime factorization of <strong>fact2(u)</strong> and set up powers precalc only for odd primes up to <strong>120</strong> (there are only <strong>29</strong> of them). Then after calculation of <strong>b(r, j, u)</strong> we can do some simple <strong>O(N)</strong> loop to calculate degrees in factorization of <strong>fact2(u)</strong> and then calculate the answer as the product of prime powers in at most <strong>29 * 8</strong> multiplications. This trick also speeds up powers precalc in two times.</p>
<p>Overall complexity of the solution assuming constant time for each 128-bit operation is<br>
<strong>O(N<sup>3</sup> + π(N) * 2<sup>H</sup> * N/H + T * N * (N + π(N) * N/H))</strong>,<br>
where <strong>N = 120</strong>, <strong>H = 15</strong> and <strong>π(N)</strong> is the number of primes <strong>≤ N</strong>. Here:</p>
<ul>
<li>term <strong>N<sup>3</sup></strong> corresponds to precalc of <strong>pseudo_fact[r][j]</strong>;</li>
<li>term <strong>π(N) * 2<sup>H</sup> * N/H</strong> corresponds to precalc of prime powers;</li>
<li>term <strong>T * N * (N + π(N) * N/H)</strong> means that for each of <strong>T</strong> tests we calculate <strong>O(N)</strong> values <strong>fatc2(u)</strong> in <strong>O(N + π(N) * N/H)</strong> time each.</li>
</ul>
<p>You may also ask "Why did I annoy all by crappy 128-bit arithmetic?". This is because I almost invented more elementary solution with <strong>O(2<sup>N/3</sup>)</strong> precalc and <strong>poly(N)</strong> 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 :)</p>
<h1>ALTERNATIVE SOLUTION:</h1>
<p><a href="/users/10592/winger">@winger</a> solution is very interesting and differs from described one.<br>
It would be nice if he will describe his solution.</p>
<p>Solution by <a href="/users/7523/mugurelionut"><a href="/users/7523/mugurelionut">@mugurelionut</a></a> is based on more elementary math than Granville formula. You can read about it <a href="http://discuss.ww2.codechef.com/questions/8154/invbincf-editorial-considerably-improved/8420">here</a>.</p>
<h1>AUTHOR'S AND TESTER'S SOLUTIONS:</h1>
<p>Author's solution can be found <a href="http://ww2.codechef.com/viewsolution/2071520">here</a>.<br>
Tester's solution can be found <a href="http://www.codechef.com/download/Solutions/2013/April/Tester/INVBINCF.cpp">here</a>. </p>
<h1>RELATED PROBLEMS:</h1>
<p><a href="http://www.e-olimp.com.ua/en/problems/322">e-olimp - 322 - Binomial Coefficients 5</a><br>
The current problem is extreme version of this one.<br>
I've set this problem (for <strong>n ≤ 40</strong>) at Kharkiv Winter Programming School 2009 and editorial in Russian is available <a href="http://ejudge.kture.kharkov.ua/media/sbornik/Sbornik2009.pdf">here</a> (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 :)</p>
<p>The good problems to practice your skills on direct problems on calculating reduced factorials modulo prime powers:<br>
<a href="http://mathalon.in/Mathalon/?page=show_problem.php&pid=141">Mathalon - 141 - Factorial trailing digits2</a><br>
<a href="http://mathalon.in/Mathalon/?page=show_problem.php&pid=144">Mathalon - 144 - Binomial Coefficient 1</a> (this one is mine :P)</p>enThu, 18 Apr 2013 05:07:02 +0530Answer by mugurelionuthttps://discuss.codechef.com/questions/8154/invbincf-editorial/8420<p>I will try to explain the main ideas in my solution for this problem. Actually, my solution is quite similar to the general idea described in the editorial - the major difference consists of the way I compute <strong>F2(X)</strong> = the product of all the odd numbers <= X.</p>
<p>Let's consider the binary representation of <strong>X</strong> and let's traverse its bits from the most significant one to the least significant one. Thus, let's assume that <strong>X = 2^p(1) + 2^p(2) + ... + 2^p(q)</strong> (where <strong>p(1) > p(2) > ... > p(q)</strong>).</p>
<p>The first bit is handled in a special manner. We want to compute the product of all the odd numbers <strong>1 * 3 * ... * (2^p(1) - 1)</strong>. Actually, this value will be precomputed and I will denote it by <strong>SPODD(p(1),0)</strong> (I will explain later what <strong>SPODD(i,j)</strong> means).</p>
<p>Let's assume now that we reached the bit <strong>p(r >= 2)</strong> and let's denote by <strong>Y = 2^p(1) + 2^p(2) + ... + 2^p(r-1)</strong> (i.e. the sum of all the bits of X which are more significant than p(r)). We now want to compute the product of all the odd numbers in the range <strong>(Y + 1) ... (Y + 2^p(r))</strong>. If <strong>p(r)</strong> is <strong>1</strong> or <strong>0</strong> then the product consists of just one number: <strong>(Y+1)</strong>.</p>
<p>Otherwise (if <strong>p(r) >= 2</strong>) then the product will be <strong>(Y + 1) * (Y + 3) * ... * (Y + 2^p(r) - 1)</strong> (there are <strong>2^(p(r)-1)</strong> odd numbers in the product). This product can be expressed as:</p>
<ul>
<li><strong>Y^0</strong> * (the product of all the odd numbers <strong>1 * 3 * ... * (2^p(r) - 1))</strong> +</li>
<li><strong>Y^1</strong> * (the sum of all the products <strong>u(1) * ... * u(2^(p(r) - 1) - 1)</strong>, where <strong>{u(1), ..., u(2^(p(r) - 1) - 1)}</strong> goes over all the subsets of the odd numbers <strong>{1,3,...,2^p(r) - 1}</strong> having <strong>2^(p(r) - 1) - 1</strong> elements) +</li>
<li>... +</li>
<li><strong>Y^j</strong> * (the sum of all the products <strong>u(1) * ... * u(2^(p(r) - 1) - j)</strong>, where <strong>{u(1), ..., u(2^(p(r) - 1) - j}</strong> goes over all the subsets of the odd numbers <strong>{1, 3, ..., 2^p(r) - 1}</strong> having <strong>2^(p(r) - 1) - j</strong> elements) +</li>
<li>... (the sum contains <strong>2^(p(r) - 1) + 1</strong> terms overall).</li>
</ul>
<p>Now I can introduce the values <strong>SPODD(i,j)</strong> = the sum of all the products <strong>u(1) * ... * u(2^(i - 1) - j)</strong>, where <strong>{u(1), ..., u(2^(i - 1) - j}</strong> goes over all the subsets of the odd numbers <strong>{1, 3, ..., 2^i - 1}</strong> having <strong>2^(i - 1) - j</strong> elements. Thus, <strong>SPODD(i,j)</strong> means that we consider all the subsets of odd numbers < 2^i having <strong>2^(i-1) - j</strong> elements. Obviously, by the definition, <strong>SPODD(i,0) = 1 * 3 * ... * (2^i - 1)</strong>.</p>
<p>Our sum of terms mentioned earlier (Y^0 * ... + Y^1 * ... + ...) can be expressed now as <strong>Y^0 * SPODD(p(r), 0) + Y^1 * SPODD(p(r), 1) + ... + Y^j * SPODD(p(r), j) + ...</strong></p>
<p>Let's notice now that we only ever need at most <strong>120</strong> terms from this sum, when computing it modulo <strong>2^120</strong>. Note that <strong>Y</strong> is an even number. Thus, only the terms corresponding to <strong>Y^0, Y^1, ..., Y^119</strong> may have a non-zero remainder when divided by <strong>2^120</strong> (all the other terms are <strong>0</strong> modulo <strong>2^120</strong>). Actually, we can do even better. When we reached the bit <strong>p(r)</strong> (<strong>r>=2</strong>), <strong>Y</strong> is an even number which is a multiple of at least <strong>2^(p(r) + 1)</strong> (because <strong>Y=2^p(1) + ... + 2^p(r-1)</strong> and <strong>p(r-1) > p(r)</strong>). Thus, only the terms corresponding to <strong>Y^0, Y^1, ..., Y^k</strong> may be non-zero (mod <strong>2^120</strong>), where <strong>k = 119 / p(r-1)</strong>.</p>
<p>When we add things up, we notice that, in the worst case, we will need around <strong>119/119 + ... + 119/3 = O(n * log(n))</strong> 120-bit operations (<strong>n = 120</strong>) for computing <strong>F2(X)</strong> (I stopped at 119/3 because we only use this method for p(r)>=2 and, thus, p(r)+1>=3 ; for p(r)=1 or 0 I explained earlier that the result is easy to compute).</p>
<p>Since we evaluate <strong>O(n)</strong> different <strong>F2(X)</strong> values per test case, this adds up to an <strong>O(n^2 * log(n))</strong> time complexity per test case.</p>
<p>For now I left out a very important part - how to precompute all the numbers <strong>SPODD(i,j)</strong> (<strong>i,j <= 120</strong>) in <strong>O(120^3)</strong> time. I will explain this at another time, as the explanation is more complicated than anything I wrote here.</p>mugurelionutThu, 18 Apr 2013 05:07:02 +0530https://discuss.codechef.com/questions/8154/invbincf-editorial/8420Answer by mugurelionuthttps://discuss.codechef.com/questions/8154/invbincf-editorial/8183<p>1) Can you estimate the time complexity of the official solution ? (per test case, as a function of n)<br>
<strong>Editorialist reply: Done. Check the editorial.</strong></p>
<p>I am asking because at some point I had an <strong>O(n^3)</strong> per test case solution and was too slow (I mean O(n^3) 128-bit operations). I had to come up with extra ideas in order to cut the time complexity down to <strong>O(n^2 * log(n))</strong> per test case (note that I don't mean base 2 logarithm). Of course, the time complexity per test case is not all that matters - preprocessing should be counted, too (my solution had <em>O(n^3)</em> preprocessing).</p>
<p>2) During the contest (while searching the web in the hope of finding useful mathematics to use for solving the problem) I did come across the e-olimp problem you mentioned (with n <= 40). I didn't know there was any editorial available and I noticed that Anton was the only one with an accepted solution at that problem on e-olimp :) Still.. I though that Russian or Ukrainian competitors would have an advantage because of that problem.. it seems that I was wrong.</p>
<p>3) During the contest I did come across Granville's paper (among many other papers). But I was looking for simpler equations/formulas so I automatically skipped Theorem 2, which looked rather complicated to me (so I thought it wouldn't be useful) :) In the end, I used only elementary mathematics in my solution (and I had to prove everything myself from scratch, except for how to compute the modular multiplicative inverse of an odd number modulo 2^n). I think <a href="/users/14599/scli">@scli</a> used the property from Granville's paper (I browsed his solution earlier and I noticed a i * i - j * j somewhere : I didn't understand it then, but now it's starting to make sense)</p>mugurelionutMon, 15 Apr 2013 18:02:00 +0530https://discuss.codechef.com/questions/8154/invbincf-editorial/8183Answer by amitupadhyayhttps://discuss.codechef.com/questions/8154/invbincf-editorial/8178<p>though I was not able to solve the problem with some analysis I found out some facts</p>
<ol>
<li>For R=even no. The answer is always -1</li>
<li>R=1 then k=0</li>
<li>R=2<sup>n</sup>-1 then k=1</li>
<li>R=2<sup>n-1</sup>+1 then k=2</li>
<li>R=2<sup>n-1</sup>-1 then k=3</li>
<li>R=2<sup>n-2</sup>+1 then k=4</li>
<li>R=3.2<sup>n-2</sup>-1 then k=5</li>
<li>R=3.2<sup>n-2</sup>+1 then k=6</li>
<li>R=2<sup>n-2</sup>-1 then k=7</li>
<li>R=2<sup>n-3</sup>+1 then k=10</li>
<li>R=7.2<sup>n-3</sup>-1 then k=11</li>
<li>R=3.2<sup>n-3</sup>+1 then k=12</li>
</ol>
<p><a href="/users/2962/anton_lunyov">@anton_lunyov</a> correct me if I am wrong.</p>
<p>I was unable to form a series or any generalised formula. but my analysis gave me above relations. May be they help in making the program</p>amitupadhyayMon, 15 Apr 2013 17:47:42 +0530https://discuss.codechef.com/questions/8154/invbincf-editorial/8178Answer by argoshttps://discuss.codechef.com/questions/8154/invbincf-editorial/8173<p>I used the same mathematics based on Granville formula adapted for 2^n. But my solution need to be speed up to 100 times. So I decided, that you invented something cool for this this problem, but I was wrong =(</p>argosMon, 15 Apr 2013 17:21:34 +0530https://discuss.codechef.com/questions/8154/invbincf-editorial/8173