#14982 fix q_binomial

Author: videlec

src/sage/combinat/q_analogues.py

@@ -18,6 +18,7 @@
 
 from sage.misc.cachefunc import cached_function
 from sage.misc.all import prod
+from sage.structure.element import parent
 from sage.rings.all import ZZ
 from sage.combinat.dyck_word import DyckWords
 from sage.combinat.partition import Partition
@@ -215,7 +216,7 @@ def q_binomial(n, k, q=None, algorithm='auto'):
         sage: q_binomial(1/2,1)
         Traceback (most recent call last):
         ...
-        ValueError: arguments (1/2, 1) must be integers
+        TypeError: no conversion of this rational to integer
 
     One checks that `n` is nonnegative::
 
@@ -232,8 +233,13 @@ def q_binomial(n, k, q=None, algorithm='auto'):
 
     This also works for complex roots of unity::
 
-        sage: q_binomial(6,1,I)
-        1 + I
+        sage: q_binomial(6, 1, QQbar(I))
+        I + 1
+
+    Note that the symbolic computation works (see :trac:`14982`)::
+
+        sage: q_binomial(6, 1, I)
+        I + 1
 
     Check that the algorithm does not matter::
 
@@ -256,14 +262,10 @@ def q_binomial(n, k, q=None, algorithm='auto'):
     - Frederic Chapoton, David Joyner and William Stein
     """
     # sanity checks
-    if not( n in ZZ and k in ZZ ):
-        raise ValueError("arguments (%s, %s) must be integers" % (n, k))
+    n = ZZ(n)
+    k = ZZ(k)
     if n < 0:
         raise ValueError('n must be nonnegative')
-    if not(0 <= k and k <= n):
-        return 0
-
-    k = min(n-k,k) # Pick the smallest k
 
     # polynomiality test
     if q is None:
@@ -273,10 +275,19 @@ def q_binomial(n, k, q=None, algorithm='auto'):
     else:
         from sage.rings.polynomial.polynomial_element import Polynomial
         is_polynomial = isinstance(q, Polynomial)
-    from sage.symbolic.ring import SR
+
+    R = parent(q)
+    zero = R(0)
+    one = R(1)
+
+    if not(0 <= k and k <= n):
+        return zero
+
+    k = min(n-k,k) # Pick the smallest k
 
     # heuristic choice of the fastest algorithm
     if algorithm == 'auto':
+        from sage.symbolic.ring import SR
         if is_polynomial:
             if n <= 70 or k <= n/4:
                 algorithm = 'naive'
@@ -295,31 +306,29 @@ def q_binomial(n, k, q=None, algorithm='auto'):
         raise ValueError("invalid algorithm choice")
 
     # the algorithms
-    try:
-        if algorithm == 'naive':
-            denomin = prod([1 - q**i for i in range(1, k+1)])
-            if denomin == 0: # q is a root of unity, use the cyclotomic algorithm
-                algorithm = 'cyclo_generic'
-            else:
-                numerat = prod([1 - q**i for i in range(n-k+1, n+1)])
+    if algorithm == 'naive':
+        denom = prod(one - q**i for i in range(1, k+1))
+        if not denom: # q is a root of unity, use the cyclotomic algorithm
+            return cyclotomic_value(n, k, q, algorithm='cyclotomic')
+        else:
+            num = prod(one - q**i for i in range(n-k+1, n+1))
+            try:
+                return num//denom
+            except TypeError:
                 try:
-                    return numerat//denomin
-                except TypeError:
-                    return numerat/denomin
-        from sage.functions.all import floor
-        if algorithm == 'cyclo_generic':
-            from sage.rings.polynomial.cyclotomic import cyclotomic_value
-            return prod(cyclotomic_value(d,q)
-                        for d in range(2,n+1)
-                        if floor(n/d) != floor(k/d) + floor((n-k)/d))
-        if algorithm == 'cyclo_polynomial':
-            R = q.parent()
-            return prod(R.cyclotomic_polynomial(d)
-                        for d in range(2,n+1)
-                        if floor(n/d) != floor(k/d) + floor((n-k)/d))
-    except (ZeroDivisionError, TypeError):
-        # As a last attempt, do the computation formally and then substitute
-        return q_binomial(n, k)(q)
+                    return num/denom
+                except (TypeError,ZeroDivisionError):
+                    #try a substitution
+                    return q_binomial(n,k)(q)
+    elif algorithm == 'cyclo_generic':
+        from sage.rings.polynomial.cyclotomic import cyclotomic_value
+        return prod(cyclotomic_value(d,q)
+                    for d in range(2,n+1)
+                    if (n/d).floor() != (k/d).floor() + ((n-k)/d).floor())
+    elif algorithm == 'cyclo_polynomial':
+        return prod(R.cyclotomic_polynomial(d)
+                    for d in range(2,n+1)
+                    if (n/d).floor() != (k/d).floor() + ((n-k)/d).floor())
 
 def gaussian_binomial(n, k, q=None, algorithm='auto'):
     r"""

src/sage/rings/polynomial/cyclotomic.pyx

@@ -30,8 +30,11 @@ include "sage/ext/stdsage.pxi"
 include "sage/ext/interrupt.pxi"
 include "sage/ext/cdefs.pxi"
 
+from sage.structure.element cimport parent_c
+
 from sage.rings.arith import factor
 from sage.rings.infinity import infinity
+from sage.rings.integer_ring import ZZ
 from sage.misc.all import prod, subsets
 from sage.rings.integer cimport Integer
 from sage.rings.rational cimport Rational
@@ -276,19 +279,29 @@ def cyclotomic_value(n, x):
         -t^7 - t^6 - t^5 + t^2 + t + 1
         sage: cyclotomic_value(10,mod(3,4))
         1
+
+    Check that the issue with symbolic element in :trac:`14982` is fixed::
+
+        sage: a = cyclotomic_value(3, I)
+        sage: a.pyobject()
+        I
+        sage: parent(_)
+        Number Field in I with defining polynomial x^2 + 1
     """
-    n = int(n)
+    n = ZZ(n)
     if n < 3:
         if n == 1:
-            return x - 1
+            return x - ZZ.one()
         if n == 2:
-            return x + 1
+            return x + ZZ.one()
         raise ValueError("n must be positive")
 
+    P = parent_c(x)
     try:
-        return x.parent()(pari.polcyclo_eval(n, x))
+        return P(pari.polcyclo_eval(n, x).sage())
     except Exception:
         pass
+    one = P(1)
 
     # The following is modeled on the implementation in PARI and is
     # used for cases for which PARI doesn't work. These are in
@@ -297,11 +310,11 @@ def cyclotomic_value(n, x):
     # - x is some Sage type which cannot be converted to PARI;
     # - PARI's algorithm encounters a zero-divisor which is not zero.
 
-    factors = factor(n)
+    factors = n.factor()
     cdef Py_ssize_t i, j, ti, L, root_of_unity = -1
     primes = [p for p, e in factors]
     L = len(primes)
-    if any([e != 1 for p, e in factors]):
+    if any(e != 1 for p, e in factors):
         # If there are primes that occur in the factorization with multiplicity
         # greater than one we use the fact that Phi_ar(x) = Phi_r(x^a) when all
         # primes dividing a divide r.
@@ -318,7 +331,6 @@ def cyclotomic_value(n, x):
     xd = [x] # the x^d for d | n
     cdef char mu
     cdef char* md = <char*>sage_malloc(sizeof(char) * (1 << L)) # the mu(d) for d | n
-    one = x.parent()(1)
     try:
         md[0] = 1
         if L & 1: