PROBLEM LINK:
Author: Anton Lunyov
Tester: Hiroto Sekido
Editorialist: Anton Lunyov
DIFFICULTY:
HARD
PREREQUISITES:
Binomial coefficients modulo prime powers, 128-bit arithmetic
PROBLEM:
For the given n and R such that 0 ≤ R < 2n you need to find the smallest K such that C(2n − 1, K) mod 2n = 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.
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). Actually only five is enough (refer to the author’s solution).
The first thing to notice is:
C(2n − 1, K) = fact2(2n − 1) / fact2(K) / fact2(2n − 1 − K) * C(2n−1 − 1, K div 2), (1)
where 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(2n − 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, …, 2n − 1},
there exists unique K in the range {0, 1, 2, …, 2n−1 − 1} such that C(2n − 1, K) mod 2n = R.
At first we infer the following congruence from (1):
C(2n − 1, 2 * K + X) ≡ (−1)K + X * C(2n−1 − 1, K) (mod 2n), where 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(2n − 1, 2 * K) + C(2n − 1, 2 * K + 1) ≡ 2n (mod 2n+1), for n ≥ 2 and 0 ≤ S < 2n−1.
Then after some non-trivial analysis we get the following simple iterative routine that will find the answer:
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
The proof in Russian is contained here (pages 57-58) 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(2n − 1, K) modulo powers of two.
And here we ask for help of British mathematician Andrew Granville, the number theory expert
We refer to Theorem 2 of his article about binomial coefficients modulo prime powers (see end of the page 3).
For p = 2 it sounds as follows:
fact2(2 * u) ≡ sgn * Product[ fact2(2 * j)b(r, j, u) : 1 ≤ j ≤ r] (mod 22 * r + 1), for r ≥ 2, (2)
where
b(r, j, u) = u / j * Product[ (u2 − i2) / (j2 − i2) : 1 ≤ i ≤ r, i ≠ j],
and 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 22 * r − 1 in view of the property:
a2n − 2 = 1 (mod 2n) 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(2n − 1, K) using formulas (1) and (2), together with INVBINCF routine above leads to quite slow solution.
Note, that actually we can maintain value C(2j − 1, K) mod 2n in the loop of INVBINCF and recalculate it using formula (1). Namely, if K = K xor 1 is performed we actually replace K by adjacent integer and can update C(2j − 1, K) straightforwardly using one multiplication and one division.
While when we do K = 2 * K + par we should transform C(2j − 1, K) to C(2j + 1 − 1, 2 * K + par) and this exactly can be made by formula (1) by calculating three values of fact2(X) by formula (2). Also we can maintain numerator and denominator of C(2j − 1, K) separately and get rid of divisions by rewriting check in INVBINCF in obvious way. Refer to routine invbin in author’s solution for clear implementation of this.
So INVBINCF is ready to be used but we need to figure it out how to calculate fact2(X) mod 2n efficiently.
At first we should calculate b(j, r, u) mod 2n−2 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 2a * 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 / (j2 − i2) : 1 ≤ i ≤ r, i ≠ j] in extended format. They equal to inverse of denominators of b(r, j, u).
- Inside
fact2routine we calculate values pref[k] = u * Product[ (u2 − i2) : 1 ≤ i ≤ k] and suf[k] = Product[ (u2 − i2) : k ≤ i ≤ r] in two simple O® 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(N3) (we calculate O(N2) 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(N2) time we can speed up the precalc part to O(N2), but even O(N3) 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 ab mod 2n. But even using exponentiation by squaring this will lead to O(N2) complexity for the fact2 routine and thus to O(N3) 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] = Aj * 2H * k for all 0 ≤ k < N/H and 0 ≤ j < 2H, where H is some constant (divisor of N = 120 will be most convenient). Then, when we want to calculate AB for some B < 2N, we represent B in 2H-ary system as
B = B0 + B1 * 2H + … + BK * 2H * K
and calculate AB 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 trick also speeds up powers precalc in two times.
Overall complexity of the solution assuming constant time for each 128-bit operation is
O(N3 + π(N) * 2H * N/H + T * N * (N + π(N) * N/H)),
where N = 120, H = 15 and π(N) is the number of primes ≤ N. Here:
- term N3 corresponds to precalc of pseudo_fact[r][j];
- term π(N) * 2H * N/H corresponds to precalc of prime powers;
- term T * N * (N + π(N) * N/H) means that for each of T tests we calculate O(N) values fatc2(u) in O(N + π(N) * N/H) time 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(2N/3) precalc and 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
ALTERNATIVE SOLUTION:
@winger solution is very interesting and 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.
AUTHOR’S AND TESTER’S SOLUTIONS:
Author’s solution can be found here.
Tester’s solution can be found here.
RELATED PROBLEMS:
e-olimp - 322 - Binomial Coefficients 5
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 (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
Mathalon - 144 - Binomial Coefficient 1 (this one is mine :P)
