Computing n! modulo pᵏ for small p

Revision en6, by adamant, 2023-05-23 00:40:56

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(a,p)=1, and then report t and a modulo pk. This is sufficient to also compute (nr) modulo pk.

We will show that, assuming that polynomial multiplication of size n requires O(nlogn) operations, and assuming that arithmetic operations modulo pk take O(1), we can find t and a in O(d2+dklogk), where d=logpn. It requires O(pklog2k) pre-computation that takes O(pklogk) memory.

Algorithm with polynomials

Great thanks to Endagorion for sharing this algorithm with me!

Let k=np, then we can represent n! as

n!=pkk!1jngcd(j,n)=1j.

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=a0+a1p++adpd, then the later can be represented as follows:

1jngcd(j,n)=1j=di=01jaipigcd(j,p)=1(npi+1pi+1+j).

To compute it quickly, we can define a family of polynomial Pi,a(x) for 0ap such that

Pi,a(x)=1japi1gcd(j,p)=1(xpi+j),

so the value of the factorial would be represented as

n!=pkk!di=0Pi+1,ai(npi+1),

and expanding k! it then rewrites into

n!=di=0pnpi+1dij=0Pj+1,ai+j(npi+j+1),

which simplifies as

n!=di=0pnpi+1ij=0Pj+1,ai(npi+1)

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

Pi,a(x)=1japi1gcd(j,p)=1(xpi+j)=a1t=01jpi1gcd(j,p)=1(xpi+tpi1+j)=a1t=0Pi1,p(px+t).

Note that for shorter and more consistent implementation, this recurrent formula also mostly works for i=1 if we put P0,p(x)=x+1, but for P1,p we should go up to p2 instead of p1. We should also note that for Pi,a(x), we only care for coefficients up to xx/i, as the larger ones are divisible by pk.

This allows to compute Pi,a(x) in O(pklogki) for all a for a given i. Over all i from 1 to k it sums up to O(pklog2k) time and O(pklogk) memory. Then, evaluating all the polynomials for any specific a requires O(d2+dklogk) operations, where d=logpn.

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 0b<p, then we can represent its factorial as

n!=ptt!bi=1tj=0(pj+i)p1i=b+1t1j=0(pj+i).

We can further rewrite it as

n!=ptt!(p1)!tb!bi=1tj=0(1+jip)p1i=b+1t1j=0(1+jip)

Let's learn how to compute the product

A(b,t)=bi=1tj=0(1+jip),

as with it we can represent the factorial as

n!=ptt!(p1)!tb!A(b,t)A(p1,t1)A(b,t1).

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

logA(b,t)=bi=1tj=0log(1+jip)=z=1(1)z+1pzzbi=1iztj=0jz

We can precompute sums of iz for each z up to k and b up to p1 in O(pk), and we can find the sum of jz for j from 0 to t in O(z) via Lagrange interpolation. Therefore, with O(pk) precomputation, we can compute logA(b,t) in O(k2), which then allows to find n! in O(dk2), where d=logpn. 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 explogz=z with p=2 when z has remainder 3 modulo 4, and it should be tracked separately.

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)