Add Lah numbers and clean up combinatorics section

[kwankyu]

According to the Wikipedia article https://en.wikipedia.org/wiki/Lah_number, a formula that computes the Lah number $L(n,k)$, also called the Stirling numbers of the third kind, is

$$\binom{n-1}{k-1}\frac{n!}{k!}$$

but the most efficient form of the formula seems to be

$$\binom{n}{k}\frac{k}{n}\binom{n}{k}(n-k)!$$

contributed by @user202729.

Long time ago, #17161 proposed yet another formula that gives $L(n,k)$ for all integers $n,k$. But that formula is extremely slow when implemented in Sage. Like existing functions for Stirling numbers, we only support the main case where $n$ and $k$ are nonnegative.

After this PR, we may close #17161, #17159, #17123 since we do not intend to support the negative $k$ case for combinatorial functions, which is I think only of minor interest and in controversy.

I made some cosmetic changes to the file src/sage/combinat/combinat.py, mostly for user friendliness.

While at it, I ended up making extensive changes to index pages of the combinatorics section, for tidiness. In particular,

• removed some of "forever" TODOs.
• removed excessive cross references.
• uniformized capitalization.

and also fixed #17421.

📝 Checklist

• The title is concise and informative.
• The description explains in detail what this PR is about.
• I have linked a relevant issue or discussion.
• I have created tests covering the changes.
• I have updated the documentation and checked the documentation preview.

⌛ Dependencies

[github-actions]

Documentation preview for this PR (built with commit 217b186; changes) is ready! 🎉
This preview will update shortly after each push to this PR.

[user202729]

While we're at it, maybe binomial(n, k)*(n-k).factorial() would be faster than n.factorial()//k.factorial()?

[kwankyu]

Not really, at least on my machine:

sage: def lah_numbers4(n, k):
....:     if k > n:
....:         return 0
....:     if k == 0:
....:         return 0 if n else 1
....:     n = Integer(n)
....:     k = Integer(k)
....:     a = n.factorial() // k.factorial()
....:     return a**2*k // (n*(n-k).factorial())
....: 
....: def lah_numbers5(n, k):
....:     if k > n:
....:         return 0
....:     if k == 0:
....:         return 0 if n else 1
....:     n = Integer(n)
....:     k = Integer(k)
....:     a = binomial(n, k)*(n-k).factorial()
....:     return a**2*k // (n*(n-k).factorial())
....: 
sage: timeit('T1 = matrix([[lah_numbers4(n,k) for n in range(100)] for k in rang
....: e(100)])')
25 loops, best of 3: 9.74 ms per loop
sage: timeit('T2 = matrix([[lah_numbers5(n,k) for n in range(100)] for k in rang
....: e(100)])')
5 loops, best of 3: 116 ms per loop
sage: T1 = matrix([[lah_numbers4(n,k) for k in range(10)] for n in range(10)])
sage: T2 = matrix([[lah_numbers5(n,k) for k in range(10)] for n in range(10)])
sage: T1 == T2
True

[user202729]

Apparently there's massive overhead in binomial(n, k) (symbolic version) over n.binomial(k) for small input. I only tested for large input.

Try this one

    a = n.binomial(k)
    return a * k // n * a * (n-k).factorial() 

(note that n.binomial(k) * k / n == (n-1).binomial(k-1) so the division is exact)

In my experiment, it has minimal overhead for small n and k (as in your matrix) and a few times faster for e.g. lah_numbers(10^6, 10^5) .

It looks impossible to have a version that calls binomial()/factorial() only twice and avoid division though. Best division-free version I got is (n-1).binomial(k-1) * n.binomial(k) * (n-k).factorial() but it takes 3 calls (slower).

[kwankyu]

Voilà!

sage: def lah_numbers4(n, k):
....:     if k > n:
....:         return 0
....:     if k == 0:
....:         return 0 if n else 1
....:     n = Integer(n)
....:     k = Integer(k)
....:     a = n.factorial() // k.factorial()
....:     return a**2*k // (n*(n-k).factorial())
....: 
sage: def lah_numbers6(n, k):
....:     if k > n:
....:         return 0
....:     if k == 0:
....:         return 0 if n else 1
....:     n = Integer(n)
....:     k = Integer(k)
....:     a = n.binomial(k)
....:     return a * k // n * a * (n-k).factorial()
....: 
sage: timeit('T1 = matrix([[lah_numbers4(n,k) for n in range(100)] for k in range(100)])')
25 loops, best of 3: 9.36 ms per loop
sage: timeit('T1 = matrix([[lah_numbers6(n,k) for n in range(100)] for k in range(100)])')
25 loops, best of 3: 11.1 ms per loop
sage: timeit('T1 = matrix([[lah_numbers4(10000+n,10000+k) for n in range(100)] for k in range(100)])')
5 loops, best of 3: 1.91 s per loop
sage: timeit('T2 = matrix([[lah_numbers6(10000+n,10000+k) for n in range(100)] for k in range(100)])')
25 loops, best of 3: 20.7 ms per loop
sage: timeit('T1 = matrix([[lah_numbers4(10000+n,k) for n in range(100)] for k in range(100)])')
5 loops, best of 3: 10.7 s per loop
sage: timeit('T2 = matrix([[lah_numbers6(10000+n,k) for n in range(100)] for k in range(100)])')
5 loops, best of 3: 2.19 s per loop

Nice!

[kwankyu]

@fchapoton As you currently seem to be in the "combinat" section, would you review this PR? This PR touches many files in the directory src/sage/combinat to improve the documentation, but most of the changes are quite small.

The new code for Lah numbers was, I think, already reviewed by @user202729.