Speed up the computation of Bell numbers

[gerbicz]

For large n values use the formula: bell(n) = exp(-1)*sum(k=0, infinity, k**n/k!)

Currently my code beats all known implementations, also sage's code which wraps GAP's Bell is slow and using lots of memory, and unable to compute bell_number(n) if n is large.

Some timings on Amd 2.1 GHz:

n time (Wall time)

300 0.01 sec.

1000 0.03 sec.

3000 0.18 sec.

10000 2.88 sec.

30000 46.79 sec.

100000 819.37 sec.

Compare these with another programs: ​http://fredrik-j.blogspot.com/2009/03/computing-generalized-bell-numbers.html

Also exposes mpmath's bell() as an alternative algorithm.

---

Apply: trac_10170-bell_number_improvements-ts.patch​

[jbandlow]

Hello and thanks for the patch! I'm curious, have you compared your solution with wrapping mpmath directly? For example,

sage: import mpmath
sage: mpmath.bell?

Also, this ticket is basically a dup of #7812, so one of these (maybe that one, since this has a patch) should be closed.

Finally, don't forget to set the status as 'needs_review' once you've added a patch!

Thanks, Jason

[gerbicz]

Replying to jbandlow:

Hello and thanks for the patch! I'm curious, have you compared your solution with wrapping mpmath directly? For example,

sage: import mpmath
sage: mpmath.bell?

Yes, but times are not saved I will rerun them: here are the timings for mpmath:

n time (Wall time)

300 0.02 sec.

1000 0.06 sec.

3000 0.41 sec.

10000 10.26 sec.

30000 168.27 sec.

100000 (still computing)

So here mpmath is about 2-4 times slower, and seems that using more memory than my code.

[gerbicz]

Replying to gerbicz:

100000 (still computing)

it took 2383 sec. for mpmath.

[allmar]

Patch will not apply to sage 4.5.3

[tscrim]

I've uploaded a new version of the patch which also adds mpmath (so #7812 can be closed as a duplicate after this patch is done). I've also come across an error in calling with the 'mpmath' agrument and started #14247.

For patchbot:

Apply: trac_10170-bell_number_improvements-ts.patch

[tscrim]

I've folded in the changes from #14247 and marked that as a duplicate since there's no need for a separate ticket.

For patchbot:

Apply: trac_10170-bell_number_improvements-ts.patch

[ncohen]

Hmmm... Looks nice, but shouldn't this be in that module ?

​http://www.sagemath.org/doc/reference/combinat/sage/combinat/expnums.html

Nathann

[tscrim]

From my 10 minutes of googling, no, it should not be. The exponential numbers seem to be a generalization of the Bell numbers and I don't know that this implementation can be generalized to that (and hence go into that module). My guess is that it's probably yes, but I'll leave that to someone with more expertise to answer.

[ncohen]

At the line where you define E_1 and E_2 they appear to be perfectly equal. Is that intended ? O_o

Nathann

[gerbicz]

That is not true, in fact: E_2=exp(-1)*E_1 as you can read in my original bell_proof.txt

About exponential numbers: that is a different problem. Here we compute only a single Bell number not the first n Bell numbers and there the algorithm is optimal (O(n*n) time), the same question for Bell numbers is an open problem.

[ncohen]

That is not true, in fact:

I have no idea, but this is what is written in the patch.

Nathann

[gerbicz]

Yes, there is a missing parentheses in line 255 in the patch.

[tscrim]

Folded in mpmath patch

[tscrim]

Fixed.