Opened 5 years ago

Closed 5 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 jdemeyer)

Improve the error analysis for the implementation of Dobinski's formula to compute the Bell numbers leading to

  1. code which is actually consistent with the mathematical formulas
  2. code which is cleaner and simpler than before
  3. a complete error estimate including the final real approximation
  4. 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 5 years ago by jdemeyer

  • Description modified (diff)

comment:2 Changed 5 years ago by jdemeyer

  • Description modified (diff)
  • Summary changed from Cleanup formula for Bell numbers to Improve formula for Bell numbers

comment:3 Changed 5 years ago by jdemeyer

  • Description modified (diff)
  • Status changed from new to needs_review

comment:4 Changed 5 years ago by jdemeyer

  • 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 5 years ago by jdemeyer

  • Commit set to 0266cef42895b095b8c7eca5718a1466300e1304
  • Description modified (diff)

New commits:

0266cefImprove formula for Bell numbers

comment:6 follow-ups: Changed 5 years ago by fredrik.johansson

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: Changed 5 years ago by jdemeyer

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: Changed 5 years ago by jdemeyer

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 5 years ago by fredrik.johansson

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 5 years ago by fredrik.johansson

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 5 years ago by fredrik.johansson

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 5 years ago by git

  • Commit changed from 0266cef42895b095b8c7eca5718a1466300e1304 to 58984af74dcf78a93b3751d7d25862168f8528a9

Branch pushed to git repo; I updated commit sha1. New commits:

58984afRound down instead of up

comment:13 in reply to: ↑ 8 Changed 5 years ago by tscrim

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 5 years ago by git

  • Commit changed from 58984af74dcf78a93b3751d7d25862168f8528a9 to c93c4e92544c2cc2cca5119b95a31b0e62969797

Branch pushed to git repo; I updated commit sha1. New commits:

ce22602Cleanup/reorganization of FLINT imports
fc17740Merge remote-tracking branch 'trac/u/jdemeyer/ticket/16428' into mpir
a669ec7Fix test in cython.py for new FLINT layout.
b69be0dFix some ctypedef in FLINT pxd files.
5851a85Merge commit 'b69be0da26e43f9454e676867afb54de530d2a58' into HEAD
c93c4e9Add FLINT implementation of Bell numbers

comment:15 Changed 5 years ago by tscrim

  • Dependencies set to #16428

Waiting now for all my cython files to recompile.

comment:16 Changed 5 years ago by tscrim

  • 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 5 years ago by vbraun

  • Status changed from positive_review to needs_work

Conflicts with #15300 most likely

comment:18 Changed 5 years ago by git

  • Commit changed from c93c4e92544c2cc2cca5119b95a31b0e62969797 to 330617caa6e52c1cc5d7b67ae41fd5cf30ef8dc9

Branch pushed to git repo; I updated commit sha1. Last 10 new commits:

165f943Rewrite Weyl algebra and reviewer suggestions added.
34c9b4cChanges from the reviewers.
b47af6aA minor tweaking of the doc.
53a9d7eFamilies return output based on order of variable names. Fixed coercion issue.
339b71dFixed some typos.
3ce3f6aMerge branch 'public/algebras/weyl_clifford-15300' of git://trac.sagemath.org/sage into clifford
8698275lifting bilinear forms onto exterior algebras
a002f4bSome tweaks to lifted_bilinear_form and fixed doctest.
ff27bdcfurther fixes
330617cMerge commit 'ff27bdc5d87967d0e7fd5c1305cef4340db816c7' into ticket/17157

comment:19 Changed 5 years ago by jdemeyer

  • Dependencies changed from #16428 to #16428, #15300
  • Status changed from needs_work to positive_review

Last 10 new commits:

165f943Rewrite Weyl algebra and reviewer suggestions added.
34c9b4cChanges from the reviewers.
b47af6aA minor tweaking of the doc.
53a9d7eFamilies return output based on order of variable names. Fixed coercion issue.
339b71dFixed some typos.
3ce3f6aMerge branch 'public/algebras/weyl_clifford-15300' of git://trac.sagemath.org/sage into clifford
8698275lifting bilinear forms onto exterior algebras
a002f4bSome tweaks to lifted_bilinear_form and fixed doctest.
ff27bdcfurther fixes
330617cMerge commit 'ff27bdc5d87967d0e7fd5c1305cef4340db816c7' into ticket/17157

comment:20 Changed 5 years ago by vbraun

  • Branch changed from u/jdemeyer/ticket/17157 to 330617caa6e52c1cc5d7b67ae41fd5cf30ef8dc9
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.