unify the three implementations of gaussian q-binomial coefficients

[fchapoton]

In sage 5.8, one can find the gaussian q-binomial coefficients in (at least) three places :

• sage.combinat.sf.kfpoly.q_bin
• sage.combinat.q_analogues.q_binomial
• sage.rings.arith.gaussian_binomial

The syntax is not quite the same for q_bin as for the two others.

Some timings:

sage: timeit("sage.combinat.sf.kfpoly.q_bin(7,5)")
625 loops, best of 3: 819 µs per loop
sage: timeit("sage.combinat.q_analogues.q_binomial(12,5)")
625 loops, best of 3: 1.12 ms per loop
sage: timeit("sage.rings.arith.gaussian_binomial(12,5)")  
625 loops, best of 3: 294 µs per loop

The parents :

sage: sage.combinat.sf.kfpoly.q_bin(7,5).parent()
Fraction Field of Univariate Polynomial Ring in t over Integer Ring
sage: sage.combinat.q_analogues.q_binomial(12,5).parent()
Univariate Polynomial Ring in q over Integer Ring
sage: sage.rings.arith.gaussian_binomial(12,5).parent()
Fraction Field of Univariate Polynomial Ring in q over Integer Ring

The behaviour with an added integer parameter

sage: sage.combinat.sf.kfpoly.q_bin(7,5,2)               
114429029715
sage: sage.combinat.q_analogues.q_binomial(12,5,2)
DOES NOT WORK
sage: sage.rings.arith.gaussian_binomial(12,5,2)    
114429029715

The behaviour with a polynomial parameter

sage: w=polygen(ZZ,'w')
sage: sage.combinat.sf.kfpoly.q_bin(4,2,w)
w^8 + w^7 + 2*w^6 + 2*w^5 + 3*w^4 + 2*w^3 + 2*w^2 + w + 1
sage: sage.combinat.q_analogues.q_binomial(6,2,w)   
w^8 + w^7 + 2*w^6 + 2*w^5 + 3*w^4 + 2*w^3 + 2*w^2 + w + 1
sage: sage.rings.arith.gaussian_binomial(6,2,w) 
w^8 + w^7 + 2*w^6 + 2*w^5 + 3*w^4 + 2*w^3 + 2*w^2 + w + 1
sage: sage.combinat.sf.kfpoly.q_bin(4,2,w).parent()             
Univariate Polynomial Ring in w over Integer Ring
sage: sage.combinat.q_analogues.q_binomial(6,2,w).parent()      
Univariate Polynomial Ring in w over Integer Ring
sage: sage.rings.arith.gaussian_binomial(6,2,w).parent()        
Fraction Field of Univariate Polynomial Ring in w over Integer Ring

Maybe it would be better to have a single function ?

CC: @sagetrac-fwclarke @AndrewAtLarge

Component: combinatorics

Keywords: gaussian binomial

Author: Frédéric Chapoton

Reviewer: Francis Clarke, Travis Scrimshaw

Merged: sage-5.10.beta2

Issue created by migration from https://trac.sagemath.org/ticket/14496

[fchapoton]

comment:1

Here is already a seemingly faster function. There remains to deprecate the two others. The most annoying one is q_bin, with a different syntax and with a global variable inside.

[sagetrac-fwclarke]

comment:3

This looks good. It's certainly nice to have integer parameters giving integer results.

But n//d can causes problems if q takes a real value, and similarly in some other cases. Perhaps this should be a "try:", with the exception giving n/d or even gaussian_binomial(n, k)(q). [It is surely not necessary to introduce var when this occurs before.]

I also don't like the way (inherited from the existing code) that n plays two very different rôles; the line

        n = prod([1 - q**i for i in range(n-k+1,n+1)])

is particularly bad style. It would make the code more readable to call them "numerator" and "denominator" (possibly abbreviated)

[fchapoton]

comment:4

Here is a new version of the patch. I have modified the names d and n and introduced a try statement. I am not sure what kind of exceptions I need to catch.

Still it remains to choose if one keeps the gaussian_binomial or rather the q_binomial function.

[fchapoton]

Author: Frédéric Chapoton

[sagetrac-fwclarke]

comment:6

1. There's a doctest failure at line 341 of sage/combinat/sf/kfpoly.py. I don't know what's being computed here, so I don't know what's going wrong.

2. I thought the standard spacing was as

def gaussian_binomial(n, k, q=None): 

3. I think TypeError is enough; I never saw anything else when I was testing your code.

4. Is it intentional not to alter sage/combinat/q_analogues.py ?

So it's nearly ready for a positive review.

[fchapoton]

comment:7

Here is a new patch.

Point 1 is indeed a problem. I do not understand either what is computed in the weight function. The procedure q_bin was returning 1 for any negative values of n and any k. This is not a very natural behaviour.

Point 4 : I have chosen that the name to be kept would be q_binomial.

There is another problem remaining. The existing algorithm in sage.combinat.q_analogues uses cyclotomic polynomials and become faster than the naive one when the arguments are big (and of similar size ?), but is much slower in some cases (when n is small or k is small ?).

sage: timeit("sage.combinat.q_analogues.q_binomial(85,3)")
625 loops, best of 3: 246 µs per loop
sage: timeit("sage.combinat.q_analogues.q_binomial_via_cyclotomic(85,3)")
125 loops, best of 3: 2.78 ms per loop
sage: timeit("sage.combinat.q_analogues.q_binomial(85,43)")               
25 loops, best of 3: 9.21 ms per loop
sage: timeit("sage.combinat.q_analogues.q_binomial_via_cyclotomic(85,43)")
25 loops, best of 3: 8.18 ms per loop

So one needs to keep both and find a strategy of choice..

[tscrim]

comment:8

I'll take a look at why the weight function in SF fails; but the best way to fix it is in the weight, have a check for negative n or k and just circumvent the call to q_binomial() and return 1.

As for the speed of the naive vs. cyclotomic polys, I would pick a cutoff value which you think is about the breakeven point(s) for the two, and call which ever function is faster based on that.

Also, have you tried implementing it by computing the q-binomial, then substituting in a user defined q value to the result (in particular, have you compared the speed)? I'm thinking this might be slower than just computing it outright if q is in GF(3) for example... If you haven't, I can write it up and run the tests.

Thanks for noticing this and working on it,

Travis

[fchapoton]

comment:9

Good idea, I have done just that in the new patch, so there is no doctest failure now.

There remains to choose a strategy..

[fchapoton]

comment:10

Some timings were given in the ticket #13166

According to my own timings, it seems that for n <= ~ 78 the naive algorithm is faster than the cyclotomic algorithm for computing the gaussian polynomial.

[sagetrac-fwclarke]

comment:11

Replying to @fchapoton:

According to my own timings, it seems that for n <= ~ 78 the naive algorithm is faster than the cyclotomic algorithm for computing the gaussian polynomial.

I get much the same. For gaussian_binomial(2*k, k) I found the naive method faster for k < 30, and the cyclotomic method faster for `k > 30'.

One point about the naive method:

sage: %timeit gaussian_binomial(50, 4)
1000 loops, best of 3: 287 µs per loop
sage: %timeit gaussian_binomial(50, 46)
100 loops, best of 3: 3.41 ms per loop

Since the code contains two loops of length k, it's obviously better to evaluate gaussian_binomial(n, k) as gaussian_binomial(n, n-k) if k > n/2.

[fchapoton]

comment:12

Here is a new patch

• using the minimum of k and n-k in the naive algo

• implementing the following strategy :

if n <= 70 or k <= N/4 or q is not a polynomial, use the naive algo

otherwise use the cyclotomic algo

I have found this strategy using some timings (I only timed the polynomial case).

[fchapoton]

Reviewer: Francis Clarke, Travis Scrimshaw

[sagetrac-fwclarke]

comment:15

Two points remain, I think:

1. On speed, I see that you are not using the cyclotomic method when q is not a polynomial. But by using cyclotomic_value the faster method can be efficient for "scalar" values too.

2. I don't see any reason why gaussian_binomial should be deprecated. This terminology is in common use, and a mathematician who knows only this name should not be told that it is "wrong". She should be able to use it without necessarily ever knowing that some people use a different name. Doesn't deprecation mean that eventually the function will disappear from Sage? This would only cause confusion. Surely the two should be synonyms, as are, for example, the number-field methods ring_of_integers and maximal_order.

[tscrim]

comment:18

Hey Frederic and Francis,

I've uploaded a review patch which does two things I wanted, an input for the algorithm and handling objects which may not have division (by computing it using a formal variable and then substituting in q). If you're happy with my changes, feel free to set this to a positive review.

Best,

Travis

[fchapoton]

comment:19

hello,

it seems to me that at the end you can write simply

return q_binomial(n, k)(q) 

because then the polygen import will be done inside the function

[tscrim]

comment:20

Attachment: trac_14496-qbinomial-review-ts.patch.gz

Yep, that works. I've uploaded the updated patch.

[sagetrac-fwclarke]

comment:21

I was about to say positive review, when I noticed two things:

1. gaussian_polynomial(5, 2, 7r) fails. But this caused by a feature of cyclotomic_value, so should probably be ignored.

2. gaussian_binomial(4, 2, Zmod(6)(2), algorithm='naive') fails. This can be very simply rectified by excepting a ZeroDivisionError as well as a TypeError just before the "last attempt".

[fchapoton]

comment:22

Point 2 is done. I vote to ignore point 1 (anyway, who cares about that ?)

[sagetrac-fwclarke]

comment:23

Replying to @fchapoton:

Point 2 is done. I vote to ignore point 1 (anyway, who cares about that ?)

Agreed. Positive review.

[jdemeyer]

comment:24

dochtml.log:[combinat ] /mazur/release/merger/sage-5.10.beta2/local/lib/python2.7/site-packages/sage/combinat/q_analogues.py:docstring of sage.combinat.q_analogues.q_binomial:141: WARNING: Duplicate explicit target name: "ch2006".
dochtml.log:[combinat ] /mazur/release/merger/sage-5.10.beta2/local/lib/python2.7/site-packages/sage/combinat/q_analogues.py:docstring of sage.combinat.q_analogues.q_binomial:141: WARNING: duplicate citation CH2006, other instance in /mazur/release/merger/sage-5.10.beta2/devel/sage/doc/en/reference/combinat/sage/combinat/q_analogues.rst

[fchapoton]

comment:25

ah, well, sorry.

This is because I defined an alias in an incorrect way, using
gaussian_binomial = q_binomial

But what is the correct way ?

[fchapoton]

comment:26

Hum, it seems that I have used the correct way to create an alias. But as the doc is duplicated (instead of saying that one function is an alias for the other), it will necessarily create a conflict any time it contains a reference. I am surprised that this problem has never been met before !

One can compare the situation to the definition of charpoly and characteristic_polynomial in sage.misc.functional, which does not contain a reference, hence dos not raise this problem.

[tscrim]

comment:27

There are two workarounds I can see:

1. Put the reference in a place that isn't duplicated such as at the module level.
2. Write in gaussian_binomial() which calls q_binomial() and it's doc just says "see q_binomial()"

Option 2 is more work and less clean IMO. Option 1 isn't perfect to me either but is the lesser of 2 evils. Up to you; either way would be okay with me.

[fchapoton]

comment:29

ok, this should be good now. Back to positive review

[tscrim]

comment:30

Hey Frederic,

Two minor things on the doc for gaussian_binomial():

• There should be periods '.' at the end of the sentences.
• The 'this function' is slightly ambiguous to me, I would explicitly put :func:`q_binomial()` to avoid this.

Otherwise I'm also happy with the patch. Thanks.

[fchapoton]

Attachment: trac_14496-qbinomial-rev-rev-fc.patch.gz

[fchapoton]

comment:31

Done, and thanks for the review

[jdemeyer]

Merged: sage-5.10.beta2