Opened 7 years ago

Closed 7 years ago

#14496 closed enhancement (fixed)

unify the three implementations of gaussian q-binomial coefficients

Reported by: chapoton Owned by: tbd
Priority: major Milestone: sage-5.10
Component: combinatorics Keywords: gaussian binomial
Cc: fwclarke, andrew.mathas Merged in: sage-5.10.beta2
Authors: Frédéric Chapoton Reviewers: Francis Clarke, Travis Scrimshaw
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description

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 ?

Attachments (3)

trac_14496_qbinomial-fc.patch (15.8 KB) - added by chapoton 7 years ago.
trac_14496-qbinomial-review-ts.patch (7.8 KB) - added by tscrim 7 years ago.
trac_14496-qbinomial-rev-rev-fc.patch (2.4 KB) - added by chapoton 7 years ago.

Download all attachments as: .zip

Change History (35)

comment:1 Changed 7 years ago by chapoton

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.

comment:2 Changed 7 years ago by chapoton

  • Cc fwclarke added
  • Status changed from new to needs_review

comment:3 Changed 7 years ago by fwclarke

  • Status changed from needs_review to needs_work

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)

comment:4 Changed 7 years ago by chapoton

  • Status changed from needs_work to needs_review

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.

comment:5 Changed 7 years ago by chapoton

  • Authors set to Frédéric Chapoton
  • Component changed from PLEASE CHANGE to combinatorics

comment:6 Changed 7 years ago by fwclarke

  • Status changed from needs_review to needs_work
  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.
  1. I thought the standard spacing was as
    def gaussian_binomial(n, k, q=None): 
    
  1. I think TypeError is enough; I never saw anything else when I was testing your code.
  1. Is it intentional not to alter sage/combinat/q_analogues.py ?

So it's nearly ready for a positive review.

comment:7 Changed 7 years ago by chapoton

  • Status changed from needs_work to needs_info

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..

comment:8 Changed 7 years ago by tscrim

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

comment:9 Changed 7 years ago by chapoton

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

There remains to choose a strategy..

comment:10 follow-up: Changed 7 years ago by chapoton

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.

Last edited 7 years ago by chapoton (previous) (diff)

comment:11 in reply to: ↑ 10 Changed 7 years ago by fwclarke

Replying to chapoton:

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.

comment:12 Changed 7 years ago by chapoton

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).

comment:13 Changed 7 years ago by chapoton

  • Reviewers set to Francis Clarke, Travis Scrimshaw
  • Status changed from needs_info to needs_review

comment:14 Changed 7 years ago by chapoton

  • Type changed from task to enhancement

comment:15 Changed 7 years ago by fwclarke

  • Status changed from needs_review to needs_work

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.
  1. 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.

Changed 7 years ago by chapoton

comment:16 Changed 7 years ago by chapoton

Here is a new patch.

For point 2, I have now defined gaussian_binomial as being an alias to q_binomial. The global namespace only contains gaussian_binomial, as before the patch.

For point 1, I have made some changes in the choice of algo, mainly about symbolic inputs. But I do not think that it is really worth spending more time on trying harder to enhance the speed for arbitrary inputs. This function is supposed to be used mainly without 3rd argument.

comment:17 Changed 7 years ago by chapoton

  • Status changed from needs_work to needs_review

comment:18 Changed 7 years ago by tscrim

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

comment:19 Changed 7 years ago by chapoton

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

Changed 7 years ago by tscrim

comment:20 Changed 7 years ago by tscrim

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

comment:21 Changed 7 years ago by fwclarke

  • Status changed from needs_review to needs_work

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.
  1. 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".

comment:22 follow-up: Changed 7 years ago by chapoton

  • Status changed from needs_work to needs_review

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

comment:23 in reply to: ↑ 22 Changed 7 years ago by fwclarke

  • Status changed from needs_review to positive_review

Replying to chapoton:

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

Agreed. Positive review.

comment:24 Changed 7 years ago by jdemeyer

  • Status changed from positive_review to needs_work
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

comment:25 Changed 7 years ago by chapoton

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 ?

comment:26 Changed 7 years ago by chapoton

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.

comment:27 Changed 7 years ago by tscrim

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.

comment:28 Changed 7 years ago by andrew.mathas

  • Cc andrew.mathas added

comment:29 Changed 7 years ago by chapoton

  • Status changed from needs_work to positive_review

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

comment:30 Changed 7 years ago by tscrim

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.

Changed 7 years ago by chapoton

comment:31 Changed 7 years ago by chapoton

Done, and thanks for the review

comment:32 Changed 7 years ago by jdemeyer

  • Merged in set to sage-5.10.beta2
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.