I am trying to understand how to find binomial coefficients modulo a power of a prime. I am reading the paper by Andrew Granville for this. But I am unable to understand it completely. More specifically, I am unable to work out how each of $(N_j!)_p$ is computed efficiently. It would be really awesome if someone could show a small hand-worked example too – say $\binom{16}{5}$ mod $3^3$ or any other small example. Thanks in advance!
Edit Note: Earlier I had mentioned that I am unable to work out an example of the theorem by hand, but found out I was making a mistake. So have edited this question to understand how these coefficients are computed efficiently.
Best Answer
In section 3 of the paper, he gives an overview of his proposed algorithm. I have read (that section of) the paper and implemented most of it, but I still have one remaining problem. I'm no mathematician, but here goes. I am going to leave out equations as they appear in the paper, and refer to variables and notations as they are used there, so you will need to refer to the paper itself.
Rewrite the binomial coefficient using Theorem 1 (note that there is a trivial typo in the statement: $m = n + r$ should be $n = m + r$). For $\binom{16}{5}$ mod $3^3$, you get N = (16, 5, 1), M = (5, 1, 0), R = (11, 3, 1); n = (1, 2, 1), m = (2, 1, 0), r = (2, 0, 1). Here's the first problem: it seems to me the method described for calculating the "carry" array $e_j$ only works if you use the greater of m and r. Consider the case $\binom{100}{2}$ mod $5^3$: n = (0, 0, 4), m = (2, 0, 0), r = (3, 4, 3). If you use the suggested method, you get e = (1, 0, 0). If you swap m for r (they are interchangeable), you get (2, 1, 0), which is correct (try adding m to r a "digit" at a time and see where the carries happen; e is a cumulative reverse sum). In your example, it comes out the same both ways as e = (1, 0, 0) (comparing each digit of n and m, 1 < 2, so set $e_0 = 1$. 2 >= 1 and 1 >= 0, so set the rest of e to 0. Now go back through the array in reverse and set each element to the cumulative sum of itself and the 1s that follow. In this case, we don't change anything.).
Now you have a bunch of numbers expressed as $(N_j!)_p$. Read the section on page 9 under equation 21. We convert this to the form $(up + v)!_p$: $u = N_j / p$, $v = N_j$ mod $p$. Rewrite this to $v!(up!)_p\binom{up+v}{v}_p$. Theorem 2 handles $(up!)_p$, Theorem 3 handles $\binom{up+v}{v}_p$. Let's concentrate on $N_0 = 16$: $u = 16/3 = 5, v = 16\%3 = 1$
Let's go to Theorem 3 first. First we calculate $\alpha_j$:
$\frac{u}{j} \prod_{1 \le i \le q - 1, i \ne j} (\frac{u - i}{j - i})$
$\alpha_1 = 5 \frac{5 - 2}{1 - 2} = -15$, $\alpha_2 = \frac{5}{2} \frac{5 - 1}{2 - 1} = 10$
Then the main portion of equation 21 becomes
j = 1: $\binom{3 + 1}{1}^{a_1 = -15}_3$, j = 2: $\binom{6 + 1}{1}^{a_2 = 10}_3$
The notation $\binom{n}{m}_p$ means $\frac{(n!)_p}{(m!)_p((n-m)!)}_p$, where $(n!)_p$ says to calculate $n!$ leaving out factors divisible by $p$. You can calculate those directly, as the numbers you're dealing with are on the order of $p$ and $q$. You can also take the opportunity to apply some straightforward optimizations to avoid having to multiply the same sequence more than once, but that is an implementation detail. Remember to use the multiplicative inverse mod $p^q$ for division. Working out the above:
j = 1: $(\frac{4\cdot2\cdot1}{2\cdot1\cdot1})^{-15} = 4^{-15} \equiv 7^{15}$ (mod 27) $\equiv 10$ (mod 27)
j = 2: $(\frac{7\cdot5\cdot4\cdot2\cdot1}{5\cdot4\cdot2\cdot1\cdot1})^{10} = 7^{10} \equiv 7$ (mod 27)
Multiply those together and you get 16 (mod 27). If we just calculate the other terms by hand for now, we get $(up!)_p \equiv 14$ (mod 27) and $v! = 1$, giving us
$(16!)_p \equiv 1\cdot 14 \cdot 16$ (mod 27) $\equiv 8$ (mod 27). Checking this by calculating $(16!)_p$ directly verifies the result.
Theorem 2 works pretty much the same way, except where the author says it "allows us to express any $(up!)_p($ mod $p^q)$ in terms of $(jp!)_p$ with $j \le [q/2]$", I'm not quite sure that it does. It looks to me like you can only calculate for numbers mod $p^{2r + 1}$, which are all odd powers (and, only in limited cases, mod $p^{2r}$). I do not understand his assertion, nor how to calculate for, say, mod $5^8$. If anyone can clear this up, that would be great.
Finishing the example, here are the other $(n!)_p$ numbers:
Theorem 1 calculation:
$(1/p^{e_0=1})\binom{16}{5} \equiv \frac{8 \cdot 13}{13 \cdot 25 \cdot 2}$ (mod 27) $\equiv 25$ (mod 27). Multiply by $p^1$ on both sides, and you get $\binom{16}{5} \equiv 21$ (mod $3^3$), which Google calculator informs me is correct.
There are some details which I've glossed over, such as special cases for signs that change under certain conditions, and only needing to calculate $\alpha$ and $\beta$ mod $q - 1$ and $2r$ respectively. The main obstacle is my lack of understanding of how Theorem 2 can be made to apply in all cases.