# Calculating large Bell number modulo a composite number

I have been trying to solve http://www.javaist.com/rosecode/problem-511-Bell-Numbers-Modulo-Factorial-askyear-2018

It is not an ongoing contest problem.

We can calculate $n$th Bell number modulo a prime number greater than $10^7$ in $O(n)$ arithmetic operations, using this formula $$B_n=\sum_{k=1}^{n}\frac{k^n}{k!}\sum_{j=0}^{n-k}\frac{(-1)^j}{j!}$$ We can store the partial inner sum for future use, without calculating again and again. For example $B_{10^7}$ $\text{mod}$ $1000000007$ $=29987405$.

But in the problem the modulus is $30!$, so we can't do modular inverse operations. I tried to reduce the formula to $$B_n=\sum_{k=1}^{n}\frac{k^n!(n-k)}{k!(n-k)!}=\frac{1}{n!}\sum_{k=1}^{n}\binom{n}{k}k^n!(n-k)$$ $!n$ denotes subfactorial function.

I am stuck here. Can anyone help?

For this problem, it's more convenient to use the formula for Bell numbers expressing them as the sum of Stirling numbers of the second kind: $$B_n=\sum_{k=0}^n \left\{ {n \atop k} \right\}.$$ This formula does not involve any divisions and can be evaluated modulo any given $m$, provided that we can compute $n$-th row of Stirling numbers modulo $m$.

Stirling numbers modulo $m$ can be computed row-by-row (hence keeping only two rows in the memory at once) using the recurrence: $$\left\{{n+1\atop k}\right\} = k \left\{{ n \atop k }\right\} + \left\{{n\atop k-1}\right\}\quad (k>0)$$ with the initial conditions $\left\{{ 0 \atop 0 }\right\} = 1$ and $\left\{{ n \atop 0 }\right\} = \left\{{ 0 \atop n }\right\} = 0$ for any $n>0$.

• We can't afford to do such huge calculations. To calculate $10^7$th Bell number, we need almost $O(10^{14})$ operations. I think $p$-adic arithmetic is the right approach. Commented Sep 17, 2018 at 15:11
• @Asif; Yes. Alternative to $p$-adic numbers is to make computations modulo larger power: e.g., if you need answer modulo $2^s$ and this involves division by $2^t$ (knowing that the result is integer), then make computation modulo $2^{s+t}$ so that division by $2^t$ give you the answer modulo $2^s$. Commented Sep 17, 2018 at 15:25
• I thought about the alternate approach. But that doesn't work either. We have to divide by $n!$ which is divisible by $2^{9999992}$. We can't store such huge numbers. Commented Sep 17, 2018 at 15:32
• $p$-adic number approach doesn't seem to work due to precision issue or maybe I am doing something wrong. Still stuck at this problem :( Commented Sep 18, 2018 at 13:04
• @Asif: Well, working modulo numbers like $2^{10^7}$ is quite manageable - e.g., GIMPS is now getting closer to $2^{10^8}$. Commented Sep 18, 2018 at 18:19

As mentioned in the comment, large Bell number modulo factorials can be calculated efficiently using following congruence.

$$\sum_{k=0}^{n}D_{n,k}B_{m+k}\equiv 0\ (\textrm{mod}\ n!)$$ $$D_{n+1,k}=D_{n,k-1}-(n+1)D_{n,k}-nD_{n-1,k}$$ $n \geq 0,k \geq 0$

$D_{0,0}=1$

$D_{n,k}=0,\hspace{5 mm}k>n,k<0$

It takes $O(MN)$ time to calculate $M^{th}$ Bell number modulo $N!$.

• Thought the question is "Calculating large Bell number modulo a composite number", $n!$ is not generic composite number.
– Sil
Commented Sep 20, 2018 at 6:28
• There is a linear recurrence of order $r$ for $B_n \pmod{p^e}$, where $r$ is minimum integer such that $p^e$ divides $r!$. So, there is a linear recurrence of order $m$ for $B_n \pmod{m!}$. So given a composite number, we have to factorize and compute linear recurrence relation for each prime power and use CRT for final answer. Linear recurrence can be found using Berlekamp-Massey algorithm as demonstrated at codeforces.com/blog/entry/61306. We can also use FindLinearRecurrence function available in Mathematica. Commented Sep 20, 2018 at 8:01
• This works only if $r$ is small. Also, computing a linear recurrence from scratch has a questionable benefit, given the explicit order $r$ recurrence by Gessel. Commented Sep 20, 2018 at 14:25
• Yeah, you are right. To get a linear recurrence of order $r$, we need $O(r^2)$ Bell numbers. So we can calculate up to $r=100$ in reasonable time. Anything beyond that needs other method, I don't think we have one. Otherwise Mathematica would have implemented it in its BellB function. Commented Sep 20, 2018 at 14:44
• Btw, via modular matrix exponentiation the computation of $B_M \bmod N!$ can be done in $O(N^3\log M)$, which may be preferable when $M$ is much larger than $N$. Commented Sep 20, 2018 at 16:51