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: | sage-6.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 follow-ups: ↓ 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 ; follow-up: ↓ 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 ; follow-up: ↓ 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.6e-05 2.68e-05 200 cpu/wall(s): 8.2e-05 8.57e-05 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 remote-tracking 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_clifford-15300' 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_clifford-15300' 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