Loading [MathJax]/jax/element/mml/optable/BasicLatin.js

Computing n! modulo pᵏ for small p

Revision en8, by adamant, 2023-05-23 03:07:15

Hi everyone!

There is an article on cp-algorithms about how to compute n! modulo prime number p. Today I was wondering, how to do it if we need to compute modulo pk rather than just p. Turns out one can do it efficiently if p is small.

Motivational example

Xiaoxu Guo Contest 3 — Binomial Coefficient. Given 1n,k1018, find (nk) modulo 232.

There are also some projecteuler problems that require it, including with several queries of distinct n and k.

Task formulation and outline of result

To clarify the task a bit, our ultimate goal here is to be able to compute e.g. binomial coefficients modulo pk. Thus, what we really need is to represent n!=pta, where gcd, and then report t and a modulo p^k. This is sufficient to also compute \binom{n}{r} modulo p^k.

We will show that, assuming that polynomial multiplication of size n requires O(n \log n) operations, and assuming that arithmetic operations modulo p^k take O(1), we can find t and a in O(d^2 + dk\log k), where d=\log_p n. It requires O(pk\log^2 k) pre-computation that takes O(pk \log k) memory.

Algorithm with polynomials

Great thanks to Endagorion for sharing this algorithm with me!

Let k=\lfloor \frac{n}{p} \rfloor, then we can represent n! as

n! = p^k k! \prod\limits_{\substack{1 \leq j \leq n \\ \gcd(j, n)=1}}j.

So, there is a part contributed by numbers divisible by p, which can be reduced to computing k!, and a part contributed by numbers not divisible by p. Let n = a_0 + a_1 p + \dots + a_d p^d, then the later can be represented as follows:

\prod\limits_{\substack{1 \leq j \leq n \\ \gcd(j, n)=1}}j = \prod\limits_{i=0}^d \prod\limits_{\substack{1 \leq j \leq a_i p^i \\\gcd(j, p)=1}} \left(\left\lfloor\frac{n}{p^{i+1}} \right\rfloor p^{i+1}+j\right).

To compute it quickly, we can define a family of polynomial P_{i,a}(x) for 0 \leq a \leq p such that

P_{i, a}(x) = \prod\limits_{\substack{1 \leq j\leq a p^{i-1} \\ \gcd(j, p)=1}}(xp^{i}+j),

so the value of the factorial would be represented as

n! = p^k k! \prod\limits_{i=0}^d P_{i+1,a_i}\left(\left\lfloor\frac{n}{p^{i+1}} \right\rfloor\right),

and expanding k! it then rewrites into

n! = \prod\limits_{i=0}^{d} p^{\left\lfloor \frac{n}{p^{i+1}} \right\rfloor} \prod\limits_{j=0}^{d-i} P_{j+1,a_{i+j}}\left(\left\lfloor\frac{n}{p^{i+j+1}} \right\rfloor\right),

which simplifies as

\boxed{n! = \prod\limits_{i=0}^{d} p^{\left\lfloor \frac{n}{p^{i+1}} \right\rfloor} \prod\limits_{j=0}^{i} P_{j+1,a_{i}}\left(\left\lfloor\frac{n}{p^{i+1}} \right\rfloor\right)}

Now, what would it take us to use this setup? First of all, notice that P_{i,a}(x) can be computed from one another:

P_{i,a}(x) = \prod\limits_{\substack{1 \leq j\leq a p^{i-1} \\ \gcd(j, p)=1}}(xp^{i}+j) = \prod\limits_{t=0}^{a-1}\prod\limits_{\substack{1 \leq j\leq p^{i-1} \\ \gcd(j, p)=1}}(xp^{i}+tp^{i-1}+j)=\prod\limits_{t=0}^{a-1} P_{i-1,p}(px+t).

Note that for shorter and more consistent implementation, this recurrent formula also mostly works for i=1 if we put P_{0, p}(x) = x+1, but for P_{1,p} we should go up to p-2 instead of p-1. We should also note that for P_{i,a}(x), we only care for coefficients up to x^{\lfloor x/i \rfloor}, as the larger ones are divisible by p^k.

This allows to compute P_{i,a}(x) in O(\frac{pk \log k}{i}) for all a for a given i. Over all i from 1 to k it sums up to O(pk \log^2 k) time and O(p k \log k) memory. Then, evaluating all the polynomials for any specific a requires O(d^2 + dk \log k) operations, where d = \log_p n.

As it requires some manipulations with polynomials, I implemented it in Python with sympy just as a proof of concept:

Code

Algorithm with p-adic logarithms and exponents

The algorithm above requires some heavy machinery on polynomials. I also found out an alternative approach that allows to compute t and a for any given n divisible by p in O(pk). It doesn't rely on polynomial operations, except for Lagrange interpolation.

Let n=pt+b, where 0 \leq b < p, then we can represent its factorial as

n! = p^t t! \prod\limits_{i=1}^{b} \prod\limits_{j=0}^{t} \left(pj+i\right)\prod\limits_{i=b+1}^{p-1} \prod\limits_{j=0}^{t-1} \left(pj+i\right).

We can further rewrite it as

n! = p^t t! (p-1)!^t b! \prod\limits_{i=1}^{b} \prod\limits_{j=0}^{t} \left(1+\frac{j}{i}p\right)\prod\limits_{i=b+1}^{p-1} \prod\limits_{j=0}^{t-1} \left(1+\frac{j}{i}p\right)

Let's learn how to compute the product

A(b, t) = \prod\limits_{i=1}^b \prod\limits_{j=0}^t \left(1+\frac{j}{i}p\right),

as with it we can represent the factorial as

n! = p^t t! (p-1)!^t b! A(b, t)\frac{A(p-1,t-1)}{A(b,t-1)}.

We can take a p-adic logarithm (please refer to this article for definitions and properties) from the product to get

\boxed{\log A(b, t) = \sum\limits_{i=1}^{b} \sum\limits_{j=0}^{t} \log\left(1+\frac{j}{i}p\right)=\sum\limits_{z=1}^\infty \frac{(-1)^{z+1}p^z}{z} \sum\limits_{i=1}^{b} i^{-z} \sum\limits_{j=0}^{t} j^z}

We can precompute sums of i^{-z} for each z up to k and b up to p-1 in O(pk), and we can find the sum of j^z for j from 0 to t in O(z) via Lagrange interpolation. Therefore, with O(pk) precomputation, we can compute \log A(b, t) in O(k^2), which then allows to find n! in O(dk^2), where d = \log_p n. I implemented it in Python, as a proof of concept, but I didn't bother with faster sums over i and j:

Code

Note: Due to inherent properties of p-adic logarithms and exponents, the algorithm only works with p>2, and it might return -n! instead with p=2. This is because \exp \log z = -z with p=2 when z has remainder 3 modulo 4. It is accounted for in the code by tracking the number of integers below n that have remainder 3 modulo 4.

Some other approaches

There is also a blog by prabowo, based apparently on another blog by Min_25, and a scholarly article by Andrew Granville. I'm not sure whether they talk about the same algorithms as described here.

Tags tutorial, p-adic numbers, factorial

History

 
 
 
 
Revisions
 
 
  Rev. Lang. By When Δ Comment
en27 English adamant 2023-05-31 14:22:05 2 Tiny change: '\ \gcd(j, n)=1}}j = \' -> '\ \gcd(j, p)=1}}j = \'
en26 English adamant 2023-05-31 14:21:17 2 Tiny change: '\ \gcd(j, n)=1}}j.\n$' -> '\ \gcd(j, p)=1}}j.\n$'
en25 English adamant 2023-05-31 14:10:41 60
en24 English adamant 2023-05-31 14:10:03 531
en23 English adamant 2023-05-31 14:04:45 256
en22 English adamant 2023-05-31 13:58:15 2269
en21 English adamant 2023-05-31 02:45:28 41
en20 English adamant 2023-05-31 02:29:49 308
en19 English adamant 2023-05-31 02:26:05 308
en18 English adamant 2023-05-31 02:20:22 1037 Tiny change: 'lgorithm w' -> 'lgorithm when $p^k$ is also small\n\n\n\n#### Algorithm w'
en17 English adamant 2023-05-23 18:48:21 57
en16 English adamant 2023-05-23 17:46:45 41
en15 English adamant 2023-05-23 17:45:16 42
en14 English adamant 2023-05-23 17:44:45 1268
en13 English adamant 2023-05-23 03:22:28 22 Tiny change: 'omputation. It doesn' -> 'omputation, where $d = \log_p n$. It doesn'
en12 English adamant 2023-05-23 03:21:59 30 Tiny change: '$p$ in $O(pk)$. It doesn' -> '$p$ in $O(dk^2)$ with $O(pk)$ precomputation. It doesn'
en11 English adamant 2023-05-23 03:20:54 216
en10 English adamant 2023-05-23 03:19:19 175
en9 English adamant 2023-05-23 03:14:42 10
en8 English adamant 2023-05-23 03:07:15 184
en7 English adamant 2023-05-23 02:57:31 5 Tiny change: 'k, a = fct_test(t)\n ' -> 'k, a = fct(t)\n '
en6 English adamant 2023-05-23 00:40:56 153
en5 English adamant 2023-05-23 00:04:31 106
en4 English adamant 2023-05-22 23:23:37 625
en3 English adamant 2023-05-22 23:23:03 62
en2 English adamant 2023-05-22 23:22:24 339
en1 English adamant 2023-05-22 23:16:26 7960 Initial revision (published)