Opened 6 years ago
Closed 6 years ago
#17157 closed enhancement (fixed)
Improve formula for Bell numbers
Reported by:  jdemeyer  Owned by:  

Priority:  minor  Milestone:  sage6.4 
Component:  combinatorics  Keywords:  
Cc:  Merged in:  
Authors:  Jeroen Demeyer  Reviewers:  Travis Scrimshaw 
Report Upstream:  N/A  Work issues:  
Branch:  330617c (Commits)  Commit:  330617caa6e52c1cc5d7b67ae41fd5cf30ef8dc9 
Dependencies:  #16428, #15300  Stopgaps: 
Description (last modified by )
Improve the error analysis for the implementation of Dobinski's formula to compute the Bell numbers leading to
 code which is actually consistent with the mathematical formulas
 code which is cleaner and simpler than before
 a complete error estimate including the final real approximation
 faster code:
Before:
sage: timeit('bell_number(200)', number=2000) 2000 loops, best of 3: 1.45 ms per loop
After:
sage: timeit('bell_number(200)', number=2000) 2000 loops, best of 3: 441 µs per loop
Change History (20)
comment:1 Changed 6 years ago by
 Description modified (diff)
comment:2 Changed 6 years ago by
 Description modified (diff)
 Summary changed from Cleanup formula for Bell numbers to Improve formula for Bell numbers
comment:3 Changed 6 years ago by
 Description modified (diff)
 Status changed from new to needs_review
comment:4 Changed 6 years ago by
 Branch set to u/jdemeyer/ticket/17157
 Created changed from 10/15/14 08:35:38 to 10/15/14 08:35:38
 Modified changed from 10/15/14 21:56:27 to 10/15/14 21:56:27
comment:5 Changed 6 years ago by
 Commit set to 0266cef42895b095b8c7eca5718a1466300e1304
 Description modified (diff)
comment:6 followups: ↓ 7 ↓ 8 Changed 6 years ago by
You might also consider wrapping flint's arith_bell_number. On this computer it seems to be about twice as fast as Sage's current bell_number(n, algorithm='dobinski') for n = 1000, 2000, 4000, 10000, 20000, 40000.
comment:7 in reply to: ↑ 6 ; followup: ↓ 9 Changed 6 years ago by
Replying to fredrik.johansson:
Sage's current bell_number(n, algorithm='dobinski')
With or without this patch?
comment:8 in reply to: ↑ 6 ; followup: ↓ 13 Changed 6 years ago by
Replying to fredrik.johansson:
You might also consider wrapping flint's arith_bell_number.
If I do so, will you review this ticket?
comment:9 in reply to: ↑ 7 Changed 6 years ago by
Replying to jdemeyer:
Replying to fredrik.johansson:
Sage's current bell_number(n, algorithm='dobinski')
With or without this patch?
Without. The difference is bigger for smaller n of course. Here is a more accurate comparison:
flint:
100 cpu/wall(s): 2.6e05 2.68e05 200 cpu/wall(s): 8.2e05 8.57e05 400 cpu/wall(s): 0.00039 0.000403 1000 cpu/wall(s): 0.004 0.00441 2000 cpu/wall(s): 0.023 0.026 4000 cpu/wall(s): 0.16 0.172 10000 cpu/wall(s): 0.91 0.949 20000 cpu/wall(s): 4.02 4.173 40000 cpu/wall(s): 18.12 21.109
sage (without patch):
100 5 loops, best of 3: 1.05 ms per loop 200 625 loops, best of 3: 1.52 ms per loop 400 125 loops, best of 3: 1.96 ms per loop 1000 125 loops, best of 3: 7.39 ms per loop 2000 5 loops, best of 3: 38 ms per loop 4000 5 loops, best of 3: 216 ms per loop 10000 5 loops, best of 3: 2.17 s per loop 20000 1 loops, best of 3: 13 s per loop 40000 1 loops, best of 3: 68 s per loop
comment:10 Changed 6 years ago by
You can just compute the sum in Dobinski's formula as an exact fraction (keeping the numerator and denominator separate, not actually dividing rational numbers by successive factorials). Indeed this is what flint does for n < 5000 (for larger n it uses a multimodular algorithm). This simplifies the error analysis.
comment:11 Changed 6 years ago by
Which reminds me that flint is missing two tricks: firstly, batch computing all the powers 1^n, 3^n, 5^n, ..., N^n
using sieving (one multiplication for every odd index, then one shift for every even index), and secondly, only working with approximate values of the powers near the end of that list.
The first trick would use much more memory, of course, but that's fine as one can fall back to the multimodular algorithm for very large n anyway.
comment:12 Changed 6 years ago by
 Commit changed from 0266cef42895b095b8c7eca5718a1466300e1304 to 58984af74dcf78a93b3751d7d25862168f8528a9
Branch pushed to git repo; I updated commit sha1. New commits:
58984af  Round down instead of up

comment:13 in reply to: ↑ 8 Changed 6 years ago by
Replying to jdemeyer:
Replying to fredrik.johansson:
You might also consider wrapping flint's arith_bell_number.
If I do so, will you review this ticket?
I will review it if you do so. I'm for including access to multiple implementations.
comment:14 Changed 6 years ago by
 Commit changed from 58984af74dcf78a93b3751d7d25862168f8528a9 to c93c4e92544c2cc2cca5119b95a31b0e62969797
Branch pushed to git repo; I updated commit sha1. New commits:
ce22602  Cleanup/reorganization of FLINT imports

fc17740  Merge remotetracking branch 'trac/u/jdemeyer/ticket/16428' into mpir

a669ec7  Fix test in cython.py for new FLINT layout.

b69be0d  Fix some ctypedef in FLINT pxd files.

5851a85  Merge commit 'b69be0da26e43f9454e676867afb54de530d2a58' into HEAD

c93c4e9  Add FLINT implementation of Bell numbers

comment:15 Changed 6 years ago by
 Dependencies set to #16428
Waiting now for all my cython files to recompile.
comment:16 Changed 6 years ago by
 Reviewers set to Travis Scrimshaw
 Status changed from needs_review to positive_review
Some timings:
sage: %timeit bell_number(200) 1000 loops, best of 3: 612 µs per loop sage: %timeit bell_number(200, algorithm='dobinski') 100 loops, best of 3: 2.2 ms per loop sage: %timeit bell_number(200, algorithm='gap') 1 loops, best of 3: 21 ms per loop sage: %timeit bell_number(200, algorithm='mpmath') 10 loops, best of 3: 21.2 ms per loop sage: %timeit bell_number(1000) 10 loops, best of 3: 40.1 ms per loop sage: %timeit bell_number(1000, algorithm='dobinski') 10 loops, best of 3: 49.5 ms per loop sage: %timeit bell_number(1000, algorithm='gap') 1 loops, best of 3: 1.58 s per loop sage: %timeit bell_number(1000, algorithm='mpmath') 1 loops, best of 3: 737 ms per loop sage: %timeit bell_number(10000, algorithm='flint') 1 loops, best of 3: 3.57 s per loop sage: %timeit bell_number(10000, algorithm='dobinski') 1 loops, best of 3: 20 s per loop
So LGTM.
comment:17 Changed 6 years ago by
 Status changed from positive_review to needs_work
Conflicts with #15300 most likely
comment:18 Changed 6 years ago by
 Commit changed from c93c4e92544c2cc2cca5119b95a31b0e62969797 to 330617caa6e52c1cc5d7b67ae41fd5cf30ef8dc9
Branch pushed to git repo; I updated commit sha1. Last 10 new commits:
165f943  Rewrite Weyl algebra and reviewer suggestions added.

34c9b4c  Changes from the reviewers.

b47af6a  A minor tweaking of the doc.

53a9d7e  Families return output based on order of variable names. Fixed coercion issue.

339b71d  Fixed some typos.

3ce3f6a  Merge branch 'public/algebras/weyl_clifford15300' of git://trac.sagemath.org/sage into clifford

8698275  lifting bilinear forms onto exterior algebras

a002f4b  Some tweaks to lifted_bilinear_form and fixed doctest.

ff27bdc  further fixes

330617c  Merge commit 'ff27bdc5d87967d0e7fd5c1305cef4340db816c7' into ticket/17157

comment:19 Changed 6 years ago by
 Dependencies changed from #16428 to #16428, #15300
 Status changed from needs_work to positive_review
Last 10 new commits:
165f943  Rewrite Weyl algebra and reviewer suggestions added.

34c9b4c  Changes from the reviewers.

b47af6a  A minor tweaking of the doc.

53a9d7e  Families return output based on order of variable names. Fixed coercion issue.

339b71d  Fixed some typos.

3ce3f6a  Merge branch 'public/algebras/weyl_clifford15300' of git://trac.sagemath.org/sage into clifford

8698275  lifting bilinear forms onto exterior algebras

a002f4b  Some tweaks to lifted_bilinear_form and fixed doctest.

ff27bdc  further fixes

330617c  Merge commit 'ff27bdc5d87967d0e7fd5c1305cef4340db816c7' into ticket/17157

comment:20 Changed 6 years ago by
 Branch changed from u/jdemeyer/ticket/17157 to 330617caa6e52c1cc5d7b67ae41fd5cf30ef8dc9
 Resolution set to fixed
 Status changed from positive_review to closed
New commits:
Improve formula for Bell numbers