Bell Number Function

Versions of combinat[bell](n) recently inspected leave something to be desired. Under Maple 9 performance was adequate for $n \sim 2500;$ under Maple 10 however, it plummeted inexplicably to the point where progress past $n = 168$ defeated this user's patience. The Maple procedures below evaluate a generalisation of the Bell numbers, the Tail Coefficients $R(n, k)$; they are accompanied by a table and an abbreviated application.

$R(n+1, k-1)$ counts sequences $\( a_0, a_1, \dots, a_n\)$ satisfying the restricted growth constraint $$0 \leq a_i \leq 1 + \max (a_0, \dots, a_{i-1})\quad \text{for}\quad i = 1, \dots, n;$$ where additionally $a_0 = k-1$; these numbers arise initially in connection with indexing partitions of a set in lexicographical order [see Williamson (2002), sect. 3.33 pp 98--99.]

For our investigation the original choice of origin proves inconvenient, so we instead translate the argument to $$G_k(n) = R(n+1, k-1).$$ With this notation the original Bell numbers appear both as $B(n) = G_0(n)$ and as $B(n+1) = G_1(n)$, for $n \geq 0.$

For $n \leq 40$ or so, a simple recursion costing time order $n^2$ and implemented as genbell(n, k) is recommended:

• $$G_k(n+1) = G_{k+1}(n) + k G_k(n)\quad \text{for}\quad n \geq 0 ; \quad G_k(0) = 1 .$$

For larger $n$, a more elaborate but faster alternative is dobinski(n, k) based on the Dobinski summation

• $$G_k(n) = \frac 1{\rm e} \sum_{i = 0}^\infty \frac {(i+k)^n}{i!} .$$

By the Stirling formula, the logarithm of term $x$ of this sum is approximately

• $$ f_1(x) = n\log(x+k) - ((x+1/2)\log x - x + c)$$

where the constant $c = \log(2\pi)/2 + 1$ (roughly 2) may quietly be dropped.

An interior maximum of $f_1(x)$ occurs at $x = i_0,$ the root of

• $$ f_0(x) = (d/dx)f_1(x) = n/(x+k) - (\log x + 1/2x).$$

Happily this function is monotonic in the intervals of interest, with no root in $0 <= x < -k$ and at most one root in $-k <= x < 2n + |k|.$ The only complication is that for $k < 0$, the maximum at $x = i_0$ may be smaller than at the left endpoint: the overall maximum value of $f_1$ therefore equals either $f_1(i_0)$ or $f_1(0).$

The point $i_1$ at which summation may be truncated occurs once the terms are smaller than unity: it is easily seen that they decrease exponentially from here on. Thus $i_1$ equals the (largest) root of $f_1$ in the interval $i_0 < x < 2n + |k|;$ it is advisable to avoid the less well-behaved region around $x = -k,$ where there is actually a zero term when $n > 0$ and $k \leq 0.$

Procedure bin_chop() employs a simple-minded binary subdivision algorithm to locate roots, terminating when the error is less than unity. During summation, the precision will be increased to $(f_1(i_0) + \log i_1)/\log 10$ decimal digits (or a little more), and the sum continued to $i_1$ terms, then further as necessary until the terms drop below unity.

[The overhead of solving logarithmic equations for $i_0, i_1$ probably accounts for the relatively poor performance of the algorithm when $n$ is small. The code is obliged to fudge the estimated precision by a linear adjustment, in an attempt to compensate for the embarrassing fact that as $n$ increases, it falls somewhat short of that actually demanded by the ensuing summation. For $k < 0$ the sign of $G_k(n)$ initially alternates with $n$, negative terms with $i < k$ dominating the sum; as $n$ increases, it finally becomes positive somewhere apparently very close to $n = 4|k|\log|k| + 5.$]

The theory behind the sample number-theoretic application largely follows the lines of the 1979 citation. A brief summary follows: umbral notation is employed, in which powers $B^n, G_k^n$ are to be substituted by $B(n), G_k(n)$ after an expression has been expanded.

Touchard's recurrence asserts that

• $$(B^p - B - 1) B^n \equiv 0 \pmod p $$

for prime $p$; it may be generalised in numerous directions. Our interest in restricted tail numbers arises in connection with the periodic behaviour of Bell numbers modulo a prime: the Touchard generalisation

• $$B^n = B^{\sum b_k p^k} \equiv \prod_k (B + k)^{b_k} \pmod p$$

allows $B(n) \pmod p$ to be computed efficiently for large $n$ and small $p$ in terms of the base-$p$ representation of $n$; this leads in turn to the investigation of the numbers

• $$G_k^n = (B+k)^n = B(B-1)(B-2)...(B-k+1) B^n$$

for $n,k \geq 0$. [The second equation breaks down for $k < 0$; somewhat arcane considerations then lead to the conclusion that the first equation provides a more natural extension in this region.]

Setting $n = 0$ in Touchard's recurrence, we have $B(p) \equiv 2 \pmod p .$ David Wilson asked whether conversely satisfying the criterion $B(n) \equiv 2 \pmod n$ necessarily implies that $n$ is prime. The details of the ensuing search are a little technical; instead, the script below first merely verifies directly (at the cost of three hours of Macintosh CPU time) that the smallest known counter-example, the Bell pseudo-prime $n = 21361 = 41\times 521,$ also passes the criterion.

More generally, for prime $p$ and any $k$ it may be shown that

• $$G_k(p) \equiv k+2 \pmod p .$$

One might enquire whether $n$ above also passes this criterion when $k \neq 0$. The script now utilises a more sophisticated but less expensive approach to establish that this $n$ fails to satisfy case $k = 1$, that is

• $$G_1(n) = B(n+1) \equiv 3 \pmod n $$

is untrue.

Easily (or by the Chinese remainder theorem) if the congruence were satisfied, then $G_1(n) \equiv 3 \pmod p$ would hold for every factor $p$ of $n$. Now another generalisation of the Touchard recurrence takes the form

• $$G_k(p q) \equiv (G_k + 1)^q \pmod p ;$$

and an equation generalising the 'traditional' Bell number recursion states

• $$G_k(q + 1) = (G_k + 1)^q + k G_k^q .$$

Combining these observations, to satisfy the criterion necessarily

• $$G_k(q + 1) - k G_k^q \equiv 3 \pmod p $$

for every prime divisor $p$ of $n = p q$. This method takes about 4 seconds to establish that $n = 21361$ satisfies at $k = 0$ but fails at $k = 1$, for both prime factors.

   1 ################################################################################
   2 
   3 # Generalised Bell number G_k(n) via order n^2 recursion 
   4 #   G_k(0) = 1; G_k(n+1) = G_{k+1}(n) + k G_k(n) for n >= 0. 
   5 #   Note G_0(n) = B(n) when k = 0. 
   6 genbell := proc (n, k) local i,j,bel,beln; 
   7   if n < 0 then beln := floor(1.0/0.0) # oo for n < 0 
   8   else bel := array(0..n); 
   9     for j from 0 to n do bel[j] := 1 od; 
  10     for i from 0 to n-1 do for j from 0 to n-1-i do 
  11         bel[j] := (j+k)*bel[j] + bel[j+1] od od; 
  12     beln := bel[0] fi; 
  13   beln end; 
  14 
  15 # Short table --- Bell numbers B(n) along middle row k = 0
  16 kmin := -6; kmax := 6; nmin := 0; nmax := 10; 
  17 matrix([seq([seq(genbell(n, k), n = nmin..nmax)], k = kmin..kmax)]); 
  18 
  19 # Solve f(x) in x0_ <= x <= x1_ to accuracy eps, via iterated subdivision; 
  20 #   if no solution, returns x0; max_iter must exceed log_2(relative error)  
  21 max_iter := 40; # loop iteration bound 
  22 max_err := 1; # absolute error bound 
  23 bin_chop := proc (f, x, x0_, x1_) 
  24   local f0, f1, f2, x0, x1, x2, i; 
  25   x0 := evalf(x0_); x1 := evalf(x1_); 
  26   f0 := evalf(subs(x = x0, f)); 
  27   f1 := evalf(subs(x = x1, f)); 
  28   i := 0; 
  29   while abs(x0 - x1) > max_err and i < max_iter do 
  30     i := i+1; 
  31     x2 := evalf((x0 + x1)/2); 
  32     f2 := evalf(subs(x = x2, f)); 
  33     if f0*f2 <= 0 then x1 := x2; f1 := f2; 
  34     elif f1*f2 <= 0 then x0 := x2; f0 := f2; 
  35     else x1 := x0 fi od; 
  36   x0 end; 
  37 
  38 # Generalised Bell number G_k(n) via Dobinski (1/e) \sum_i (i+k)^n/i! 
  39 #   Note G_0(n) = B(n) when k = 0. 
  40 dobinski := proc (n, k) 
  41   local i,i_0,i_1,beln,f_0,f_1,x,mtl,dig0,dig1,dig2,faci,teri; 
  42   if n < 0 then beln := floor(1.0/0.0) # oo for n < 0 
  43   else dig0 := Digits; Digits := 10; # low precision term indices i 
  44     x := 'x': f_0 := n/(x+k) - (log(x) + 1/2/x + 1); # = (d/dx) f_1 
  45     f_1 := n*log(x+k) - ((x+1/2)*log(x) - x); # omit -(log(2*Pi)/2 + 1) 
  46     i_0 := round(bin_chop(f_0, x, max(2,2-k), 2*n+abs(k))); # max term location 
  47     i_1 := round(bin_chop(f_1, x, i_0, 2*n+abs(k))); # term count 
  48     mtl := max(evalf(subs(x = i_0, f_1)), evalf(n*log(max(1,abs(k))))); 
  49       # compare logs of located maximum and first term 
  50     dig1 := round((mtl + log(i_1))/log(10.0) * 1.005 + 4); 
  51       # log(max term * length) 
  52     Digits := dig1; # high precision terms in sum 
  53     beln := 0; faci := 1.0; i := 0; teri := 1.0; 
  54     while i <= i_1 or abs(teri) > 0.1 do # loop till long and small enough 
  55       teri := evalf(i+k)^n/faci; beln := beln + teri; 
  56       i := i+1; faci := faci*i od; 
  57     beln := round(beln/exp(1.0)); 
  58     Digits := 10; 
  59     dig2 := round(log(1+abs(evalf(beln)))/log(10.0) + 1); # verify precision 
  60     if dig2 >= dig1 then 
  61       print("dobinski: precision loss", n, k, i_0, i_1, mtl, dig1, dig2) fi; 
  62     Digits := dig0 fi; # restore precision 
  63   beln end; 
  64 
  65 # Compare with Q&D version --- all zeros in 42 sec 
  66 kmin := -15; kmax := 15; nmin := 0; nmax := 40; 
  67 matrix([seq([seq(dobinski(n, k) - genbell(n, k), 
  68   n = nmin..nmax)], k = kmin..kmax)]); 
  69 
  70 # A Bell pseudo-prime --- direct approach 
  71 n := 21361; ifactor(n); # n = 41*521 composite
  72 tim := time(): 
  73 beln := dobinski(n, 0): beln mod n; # B(n) = 2 (mod n) 
  74 tim := time() - tim; # elapsed time in seconds 
  75 ceil(log(evalf(beln))/log(10.0)); # number of digits 
  76 # B(n) has 65205 digits in 10800 sec 
  77 
  78 # Faster method: B(n) = 2 (mod n), but B(n+1) <> 3 (mod n) 
  79 tim := time(): 
  80 p := 41; q := 521; n := p*q; 
  81 k := 0; 
  82   (dobinski(p, k) - (k+2)) mod p; 
  83   (dobinski(q, k) - (k+2)) mod q; # 0, 0 
  84   (dobinski(p+1, k) - k*dobinski(p, k) - (k+2)) mod q; 
  85   (dobinski(q+1, k) - k*dobinski(q, k) - (k+2)) mod p; # 0, 0 
  86 k := 1; 
  87   (dobinski(p, k) - (k+2)) mod p; 
  88   (dobinski(q, k) - (k+2)) mod q; # 0, 0 
  89   (dobinski(p+1, k) - k*dobinski(p, k) - (k+2)) mod q; 
  90   (dobinski(q+1, k) - k*dobinski(q, k) - (k+2)) mod p; # 278, 21 
  91 tim := time() - tim; # in 4 sec 
  92 
  93 ################################################################################

References:

• S. G. Williamson Ranking Algorithms for Lists of Partitions

  • SIAM J. Comput. 5(4): 602-617 (1976);

  S. G. Williamson Combinatorics for Computer Science

  • Computer Science Press, Inc (1985), Dover (2002);

  F. Ruskey Combinatorial Generation draft, sect 4.7;

  • _
  W. F. Lunnon, P. A. B Pleasants, N. M. Stephens

  • Arithmetic Properties of Bell Numbers to Composite Modulus I Acta Arith 35: 1-16 (1979).

  N. Sloane OEIS A108087 and associated rows and columns;

Posted by FredLunnon, Nov 2008.

See also MapleGrumbles by FredLunnon.

CategoryOEIS CategoryMaple