Ticket #9541: trac_9541-add_fmpq.patch

module_list.py

@@ -468 +468 @@
               extra_compile_args=["-std=c99", "-D_XPG6"],
               depends = [SAGE_ROOT + "/local/include/FLINT/flint.h"]),
 
+    Extension('sage.libs.flint.fmpq_poly_example',
+              sources = ["sage/libs/flint/fmpq_poly_example.pyx", "sage/libs/flint/fmpq_poly.c"],
+              libraries = ["csage", "flint", "gmp", "gmpxx", "m", "stdc++"],
+              include_dirs = [SAGE_ROOT+'/local/include/FLINT/'],
+              extra_compile_args=["-std=c99", "-D_XPG6"],
+              depends = [SAGE_ROOT + "/local/include/FLINT/flint.h"]),
+
     Extension('sage.libs.flint.fmpz_poly',
               sources = ["sage/libs/flint/fmpz_poly.pyx"],
               libraries = ["csage", "flint", "gmp", "gmpxx", "m", "stdc++"],

new file sage/libs/flint/fmpq_poly.c

@@ +1 @@
+///////////////////////////////////////////////////////////////////////////////
+// Copyright (C) 2010 Sebastian Pancratz                                     //
+//                                                                           //
+// Distributed under the terms of the GNU General Public License as          //
+// published by the Free Software Foundation; either version 2 of the        //
+// License, or (at your option) any later version.                           //
+//                                                                           //
+// http://www.gnu.org/licenses/                                              //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "fmpq_poly.h"
+
+/**
+ * \file     fmpq_poly.c
+ * \brief    Fast implementation of the rational function field
+ * \author   Sebastian Pancratz
+ * \date     Mar 2010 -- July 2010
+ * \version  0.1.3
+ * 
+ * \mainpage 
+ * 
+ * <center>
+ * A FLINT-based implementation of the rational polynomial ring.
+ * </center>
+ * 
+ * \section Overview
+ * 
+ * The module <tt>fmpq_poly</tt> provides functions for performing 
+ * arithmetic on rational polynomials in \f$\mathbf{Q}[t]\f$, represented as 
+ * quotients of integer polynomials of type <tt>fmpz_poly_t</tt> and 
+ * denominators of type <tt>fmpz_t</tt>.  These functions start with the 
+ * prefix <tt>fmpq_poly_</tt>.
+ * 
+ * Rational polynomials are stored in objects of type #fmpq_poly_t, 
+ * which is an array of #fmpq_poly_struct's of length one.  This 
+ * permits passing parameters of type #fmpq_poly_t by reference.  
+ * We also define the type #fmpq_poly_ptr to be a pointer to 
+ * #fmpq_poly_struct's.
+ * 
+ * The representation of a rational polynomial as the quotient of an integer 
+ * polynomial and an integer denominator can be made canonical by demanding 
+ * the numerator and denominator to be coprime and the denominator to be 
+ * positive.  As the only special case, we represent the zero function as 
+ * \f$0/1\f$.  All arithmetic functions assume that the operands are in this 
+ * canonical form, and canonicalize their result.  If the numerator or 
+ * denominator is modified individually, for example using the methods in the 
+ * group \ref NumDen, it is the user's responsibility to canonicalize the 
+ * rational function using the function #fmpq_poly_canonicalize() if necessary.
+ * 
+ * All methods support aliasing of their inputs and outputs \e unless 
+ * explicitly stated otherwise, subject to the following caveat.  If 
+ * different rational polynomials (as objects in memory, not necessarily in 
+ * the mathematical sense) share some of the underlying integer polynomials  
+ * or integers, the behaviour is undefined.
+ * 
+ * \section Changelog  Version history
+ * - 0.1.3
+ *   - Changed #fmpq_poly_inv()
+ *   - Another sign check to <tt>fmpz_abs</tt> in #fmpq_poly_content()
+ * - 0.1.2
+ *   - Introduce a function #fmpq_poly_monic() and use this to simplify the 
+ *     code for the gcd and xgcd functions
+ *   - Make further use of #fmpq_poly_is_zero()
+ * - 0.1.1
+ *   - Replaced a few sign checks and negations by <tt>fmpz_abs</tt>
+ *   - Small changes to comments and the documentation
+ *   - Moved some function bodies from #fmpq_poly.h to #fmpq_poly.c
+ * - 0.1.0
+ *   - First draft, based on the author's Sage code
+ * 
+ * \bug  The method #fmpq_poly_to_string_pretty() leaks memory because the 
+ *       underlying FLINT 1.5.0 method does so.
+ */
+
+/**
+ * \defgroup Definitions        Type definitions
+ * \defgroup MemoryManagement   Memory management
+ * \defgroup NumDen             Accessing numerator and denominator
+ * \defgroup Assignment         Assignment and basic manipulation
+ * \defgroup Coefficients       Setting and retrieving individual coefficients
+ * \defgroup PolyParameters     Polynomial parameters
+ * \defgroup Comparison         Comparison
+ * \defgroup Addition           Addition and subtraction
+ * \defgroup ScalarMul          Scalar multiplication and division
+ * \defgroup Multiplication     Multiplication
+ * \defgroup Division           Euclidean division
+ * \defgroup Powering           Powering
+ * \defgroup GCD                Greatest common divisor
+ * \defgroup Derivative         Derivative
+ * \defgroup Evaluation         Evaluation
+ * \defgroup Content            Gaussian content
+ * \defgroup Resultant          Resultant
+ * \defgroup Composition        Composition
+ * \defgroup Squarefree         Square-free test
+ * \defgroup Subpolynomials     Subpolynomials
+ * \defgroup StringConversions  String conversions and I/O
+ */
+
+///////////////////////////////////////////////////////////////////////////////
+// Auxiliary functions                                                       //
+///////////////////////////////////////////////////////////////////////////////
+
+/**
+ * Returns the number of digits in the decimal representation of \c n.
+ */
+unsigned int fmpq_poly_places(unsigned long n)
+{
+    unsigned int count;
+    if (n == 0)
+        return 1u;
+    count = 0;
+    while (n > 0)
+    {
+        n /= 10ul;
+        count++;
+    }
+    return count;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Implementation                                                            //
+///////////////////////////////////////////////////////////////////////////////
+
+/**
+ * \brief    Puts <tt>f</tt> into canonical form.
+ * \ingroup  Definitions
+ * 
+ * This method ensures that the denominator is coprime to the content 
+ * of the numerator polynomial, and that the denominator is positive.
+ * If <tt>f-\>den == NULL</tt>, this method does not change this.
+ * 
+ * This methods assumes that the denominator, if initialised, is non-zero.
+ * However, it is allowed to be negative.
+ * 
+ * The optional parameter <tt>temp</tt> is the temporary variable that this 
+ * method would otherwise need to allocate.  If <tt>temp</tt> is provided as 
+ * an initialized <tt>fmpz_t</tt>, it is assumed \e without further checks 
+ * that it is large enough.  For example, 
+ * <tt>max(f-\>num-\>limbs, fmpz_size(f-\>den))</tt> limbs will certainly 
+ * suffice.
+ */
+void fmpq_poly_canonicalize(fmpq_poly_ptr f, fmpz_t temp)
+{
+    int tcheck;                         /* True if temp == NULL              */
+    fmpz_t tempgcd;                     /* Another temporary variable        */
+    
+    if (_fmpq_poly_den_is_one(f))       /* If the denominator is one..       */
+        return;
+    
+    if (fmpq_poly_is_zero(f))           /* If f is zero..                    */
+        fmpz_set_si(f->den, 1);
+    else if (fmpz_is_m1(f->den))        /* If f != 0 and den(f) == -1..      */
+    {
+        fmpz_poly_neg(f->num, f->num);
+        fmpz_set_si(f->den, 1);
+    }
+    else                                /* If f != 0 and den(f) != +- 1..    */
+    {
+        if (temp == NULL)
+        {
+            tcheck = 1;
+            temp = fmpz_init(FLINT_MAX(f->num->limbs, fmpz_size(f->den)));
+        }
+        else
+            tcheck = 0;
+        
+        fmpz_poly_content(temp, f->num);
+        fmpz_abs(temp, temp);
+        
+        if (fmpz_is_one(temp))          /* If the content of f is 1..        */
+        {
+            if (fmpz_sgn(f->den) < 0)
+            {
+                fmpz_poly_neg(f->num, f->num);
+                fmpz_neg(f->den, f->den);
+            }
+        }
+        else                            /* If the content exceeds 1..        */
+        {
+            tempgcd = fmpz_init(FLINT_MAX(f->num->limbs, fmpz_size(f->den)));
+            
+            fmpz_gcd(tempgcd, temp, f->den);
+            fmpz_abs(tempgcd, tempgcd);
+            
+            if (fmpz_is_one(tempgcd))
+            {
+                if (fmpz_sgn(f->den) < 0)
+                {
+                    fmpz_poly_neg(f->num, f->num);
+                    fmpz_neg(f->den, f->den);
+                }
+            }
+            else
+            {
+                if (fmpz_sgn(f->den) < 0)
+                    fmpz_neg(tempgcd, tempgcd);
+                fmpz_poly_scalar_div_fmpz(f->num, f->num, tempgcd);
+                fmpz_tdiv(temp, f->den, tempgcd);
+                fmpz_set(f->den, temp);
+            }
+            
+            fmpz_clear(tempgcd);
+        }
+        
+        if (tcheck)
+            fmpz_clear(temp);
+    }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Assignment and basic manipulation
+
+/**
+ * \ingroup  Assignment
+ * 
+ * Sets the element \c rop to the additive inverse of \c op.
+ */
+void fmpq_poly_neg(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
+{
+    if (rop == op)
+    {
+        fmpz_poly_neg(rop->num, op->num);
+    }
+    else
+    {
+        fmpz_poly_neg(rop->num, op->num);
+        fmpq_poly_set_den(rop, op->den);
+    }
+}
+
+/**
+ * \ingroup  Assignment
+ *
+ * Sets the element \c rop to the multiplicative inverse of \c op.
+ * 
+ * Assumes that the element \c op is a unit.  Otherwise, an exception 
+ * is raised in the form of an <tt>abort</tt> statement.
+ */
+void fmpq_poly_inv(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
+{
+    fmpz_t t;
+    
+    /* Assertion! */
+    if (fmpz_poly_degree(op->num) != 0)
+    {
+        printf("ERROR (fmpq_poly_inv).  Element is not a unit.\n");
+        abort();
+    }
+    
+    /* Case 1.  The denominator of op is one. */
+    if (_fmpq_poly_den_is_one(op))
+    {
+        fmpq_poly_set_den(rop, op->num->coeffs);
+        fmpz_poly_zero(rop->num);
+        fmpz_poly_set_coeff_si(rop->num, 0, 1);
+        if (rop->den != NULL && fmpz_sgn(rop->den) < 0)
+        {
+            fmpz_poly_neg(rop->num, rop->num);
+            fmpz_neg(rop->den, rop->den);
+        }
+    }
+    /* Case 2.  Now the denominator of op is not one. */
+    else
+    {
+        if (rop == op)
+        {
+            t = fmpz_init(rop->num->limbs);
+            fmpz_set(t, rop->num->coeffs);
+            fmpz_poly_zero(rop->num);
+            if (fmpz_sgn(t) > 0)
+            {
+                fmpz_poly_set_coeff_fmpz(rop->num, 0, rop->den);
+                rop->den = fmpz_realloc(rop->den, fmpz_size(t));
+                fmpz_set(rop->den, t);
+            }
+            else
+            {
+                fmpz_neg(rop->den, rop->den);
+                fmpz_poly_set_coeff_fmpz(rop->num, 0, rop->den);
+                rop->den = fmpz_realloc(rop->den, fmpz_size(t));
+                fmpz_neg(rop->den, t);
+            }
+            fmpz_clear(t);
+        }
+        else
+        {
+            fmpz_poly_zero(rop->num);
+            fmpz_poly_set_coeff_fmpz(rop->num, 0, op->den);
+            fmpq_poly_set_den(rop, op->num->coeffs);
+            if (rop->den != NULL && fmpz_sgn(rop->den) < 0)
+            {
+                fmpz_poly_neg(rop->num, rop->num);
+                fmpz_neg(rop->den, rop->den);
+            }
+        }
+    }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Setting/ retrieving individual coefficients
+
+/**
+ * \ingroup  Coefficients
+ * 
+ * Obtains the <tt>i</tt>th coefficient from the polynomial <tt>op</tt>.
+ */
+void fmpq_poly_get_coeff_mpq(mpq_t rop, const fmpq_poly_ptr op, unsigned long i)
+{
+    fmpz_poly_get_coeff_mpz(mpq_numref(rop), op->num, i);
+    if (op->den == NULL)
+        mpz_set_si(mpq_denref(rop), 1);
+    else
+        fmpz_to_mpz(mpq_denref(rop), op->den);
+    mpq_canonicalize(rop);
+}
+
+/**
+ * \ingroup  Coefficients
+ * 
+ * Sets the <tt>i</tt>th coefficient of the polynomial <tt>op</tt> to \c x.
+ */
+void fmpq_poly_set_coeff_fmpz(fmpq_poly_ptr rop, const unsigned long i, const fmpz_t x)
+{
+    fmpz_t t;
+    int canonicalize;
+    
+    /* Is the denominator 1?  This includes the case when rop == 0. */
+    if (_fmpq_poly_den_is_one(rop))
+    {
+        fmpz_poly_set_coeff_fmpz(rop->num, i, x);
+        return;
+    }
+    
+    t = fmpz_poly_get_coeff_ptr(rop->num, i);
+    canonicalize = !(t == NULL || fmpz_is_zero(t));
+    
+    t = fmpz_init(fmpz_size(x) + fmpz_size(rop->den));
+    fmpz_mul(t, x, rop->den);
+    fmpz_poly_set_coeff_fmpz(rop->num, i, t);
+    fmpz_clear(t);
+    if (canonicalize)
+        fmpq_poly_canonicalize(rop, NULL);
+}
+
+/**
+ * \ingroup  Coefficients
+ * 
+ * Sets the <tt>i</tt>th coefficient of the polynomial <tt>op</tt> to \c x.
+ */
+void fmpq_poly_set_coeff_mpz(fmpq_poly_ptr rop, const unsigned long i, const mpz_t x)
+{
+    mpz_t z;
+    fmpz_t t;
+    int canonicalize;
+    
+    /* Is the denominator 1?  This includes the case when rop == 0. */
+    if (_fmpq_poly_den_is_one(rop))
+    {
+        fmpz_poly_set_coeff_mpz(rop->num, i, x);
+        return;
+    }
+    
+    t = fmpz_poly_get_coeff_ptr(rop->num, i);
+    canonicalize = !(t == NULL || fmpz_is_zero(t));
+    
+    mpz_init(z);
+    fmpz_to_mpz(z, rop->den);
+    mpz_mul(z, z, x);
+    fmpz_poly_set_coeff_mpz(rop->num, i, z);
+    if (canonicalize)
+        fmpq_poly_canonicalize(rop, NULL);
+    mpz_clear(z);
+}
+
+/**
+ * \ingroup  Coefficients
+ * 
+ * Sets the <tt>i</tt>th coefficient of the polynomial <tt>op</tt> to \c x.
+ */
+void fmpq_poly_set_coeff_mpq(fmpq_poly_ptr rop, const unsigned long i, const mpq_t x)
+{
+    fmpz_t oldc;
+    mpz_t den, gcd, t;
+    int canonicalize;
+    
+    if (rop->den == NULL)
+        mpz_init_set_si(den, 1);
+    else
+    {
+        mpz_init(den);
+        fmpz_to_mpz(den, rop->den);
+    }
+    mpz_init(gcd);
+    mpz_gcd(gcd, den, mpq_denref(x));
+    
+    oldc = fmpz_poly_get_coeff_ptr(rop->num, i);
+    canonicalize = !(oldc == NULL || fmpz_is_zero(oldc));
+    
+    if (mpz_cmp(mpq_denref(x), gcd) == 0)
+    {
+        mpz_divexact(den, den, gcd);
+        mpz_mul(gcd, den, mpq_numref(x));
+        fmpz_poly_set_coeff_mpz(rop->num, i, gcd);
+        if (canonicalize)
+            fmpq_poly_canonicalize(rop, NULL);
+    }
+    else
+    {
+        mpz_init(t);
+        mpz_divexact(t, mpq_denref(x), gcd);
+        fmpz_poly_scalar_mul_mpz(rop->num, rop->num, t);
+        
+        mpz_divexact(gcd, den, gcd);
+        mpz_mul(gcd, gcd, mpq_numref(x));
+        
+        fmpz_poly_set_coeff_mpz(rop->num, i, gcd);
+        
+        mpz_mul(den, den, t);
+        _fmpq_poly_den_fit_limbs(rop, mpz_size(den));
+        mpz_to_fmpz(rop->den, den);
+        
+        if (canonicalize)
+            fmpq_poly_canonicalize(rop, NULL);
+        mpz_clear(t);
+    }
+    
+    /* Clean-up */
+    mpz_clear(den);
+    mpz_clear(gcd);
+}
+
+/**
+ * \ingroup  Coefficients
+ * 
+ * Sets the <tt>i</tt>th coefficient of the polynomial <tt>op</tt> to \c x.
+ */
+void fmpq_poly_set_coeff_si(fmpq_poly_ptr rop, const unsigned long i, long x)
+{
+    mpz_t z;
+    fmpz_t t;
+    int canonicalize;
+    
+    /* Is the denominator 1?  This includes the case when rop == 0. */
+    if (_fmpq_poly_den_is_one(rop))
+    {
+        fmpz_poly_set_coeff_si(rop->num, i, x);
+        return;
+    }
+    
+    t = fmpz_poly_get_coeff_ptr(rop->num, i);
+    canonicalize = !(t == NULL || fmpz_is_zero(t));
+    
+    mpz_init(z);
+    fmpz_to_mpz(z, rop->den);
+    mpz_mul_si(z, z, x);
+    fmpz_poly_set_coeff_mpz(rop->num, i, z);
+    if (canonicalize)
+        fmpq_poly_canonicalize(rop, NULL);
+    mpz_clear(z);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Comparison
+
+/**
+ * \brief    Returns whether \c op1 and \c op2 are equal.
+ * \ingroup  Comparison
+ * 
+ * Returns whether the two elements \c op1 and \c op2 are equal.
+ */
+int fmpq_poly_equal(const fmpq_poly_ptr op1, const fmpq_poly_ptr op2)
+{
+    if (_fmpq_poly_den_is_one(op1))
+    {
+        if (_fmpq_poly_den_is_one(op2))
+            return fmpz_poly_equal(op1->num, op2->num);
+        else
+            return 0;
+    }
+    else
+    {
+        if (_fmpq_poly_den_is_one(op2))
+            return 0;
+        else
+            return fmpz_poly_equal(op1->num, op2->num);
+    }
+}
+
+/**
+ * \brief   Compares the two polynomials <tt>left</tt> and <tt>right</tt>.
+ * \ingroup Comparison
+ * 
+ * Compares the two polnomials <tt>left</tt> and <tt>right</tt>, returning 
+ * <tt>-1</tt>, <tt>0</tt>, or <tt>1</tt> as <tt>left</tt> is less than, 
+ * equal to, or greater than <tt>right</tt>.
+ * 
+ * The comparison is first done by degree and then, in case of a tie, by 
+ * the individual coefficients, beginning with the highest.
+ */
+int fmpq_poly_cmp(const fmpq_poly_ptr left, const fmpq_poly_ptr right)
+{
+    int ans;
+    long degdiff, i;
+    unsigned long limbs;
+    fmpz_t lcoeff, rcoeff;
+    
+    /* Quick check whether left and right are the same object. */
+    if (left == right)
+        return 0;
+    
+    i = fmpz_poly_degree(left->num);
+    degdiff = i - fmpz_poly_degree(right->num);
+    
+    if (degdiff < 0)
+        return -1;
+    else if (degdiff > 0)
+        return 1;
+    else
+    {
+        if (_fmpq_poly_den_equal(left, right))
+        {
+            if (i == -1)  /* left and right are both zero */
+                return 0;
+            limbs = FLINT_MAX(left->num->limbs, right->num->limbs) + 1;
+            lcoeff = fmpz_init(limbs);
+            while (fmpz_equal(fmpz_poly_get_coeff_ptr(left->num, i), 
+                fmpz_poly_get_coeff_ptr(right->num, i)) && i > 0)
+                i--;
+            fmpz_sub(lcoeff, fmpz_poly_get_coeff_ptr(left->num, i), 
+                fmpz_poly_get_coeff_ptr(right->num, i));
+            ans = fmpz_sgn(lcoeff);
+            fmpz_clear(lcoeff);
+            return ans;
+        }
+        else if (_fmpq_poly_den_is_one(left))  /* Also right->den > 1 */
+        {
+            limbs = FLINT_MAX(left->num->limbs + fmpz_size(right->den), 
+                right->num->limbs) + 1;
+            lcoeff = fmpz_init(limbs);
+            fmpz_mul(lcoeff, fmpz_poly_get_coeff_ptr(left->num, i), right->den);
+            while (fmpz_equal(lcoeff, fmpz_poly_get_coeff_ptr(right->num, i)) 
+                && i > 0)
+            {
+                i--;
+                fmpz_mul(lcoeff, fmpz_poly_get_coeff_ptr(left->num, i), 
+                    right->den);
+            }
+            fmpz_sub(lcoeff, lcoeff, fmpz_poly_get_coeff_ptr(right->num, i));
+            ans = fmpz_sgn(lcoeff);
+            fmpz_clear(lcoeff);
+            return ans;
+        }
+        else if(_fmpq_poly_den_is_one(right))  /* Also left->den > 1 */
+        {
+            limbs = FLINT_MAX(left->num->limbs, 
+                right->num->limbs + fmpz_size(left->den)) + 1;
+            rcoeff = fmpz_init(limbs);
+            fmpz_mul(rcoeff, fmpz_poly_get_coeff_ptr(right->num, i), left->den);
+            while (fmpz_equal(fmpz_poly_get_coeff_ptr(left->num, i), rcoeff) 
+                && i > 0)
+            {
+                i--;
+                fmpz_mul(rcoeff, fmpz_poly_get_coeff_ptr(right->num, i), 
+                    left->den);
+            }
+            fmpz_sub(rcoeff, fmpz_poly_get_coeff_ptr(left->num, i), rcoeff);
+            ans = fmpz_sgn(rcoeff);
+            fmpz_clear(rcoeff);
+            return ans;
+        }
+        else
+        {
+            limbs = FLINT_MAX(left->num->limbs + fmpz_size(right->den), 
+                right->num->limbs + fmpz_size(left->den)) + 1;
+            lcoeff = fmpz_init(limbs);
+            rcoeff = fmpz_init(right->num->limbs + fmpz_size(left->den));
+            fmpz_mul(lcoeff, fmpz_poly_get_coeff_ptr(left->num, i), right->den);
+            fmpz_mul(rcoeff, fmpz_poly_get_coeff_ptr(right->num, i), left->den);
+            while (fmpz_equal(lcoeff, rcoeff) && i > 0)
+            {
+                i--;
+                fmpz_mul(lcoeff, fmpz_poly_get_coeff_ptr(left->num, i), 
+                    right->den);
+                fmpz_mul(rcoeff, fmpz_poly_get_coeff_ptr(right->num, i), 
+                    left->den);
+            }
+            fmpz_sub(lcoeff, lcoeff, rcoeff);
+            ans = fmpz_sgn(lcoeff);
+            fmpz_clear(lcoeff);
+            fmpz_clear(rcoeff);
+            return ans;
+        }
+    }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Addition/ subtraction
+
+/**
+ * \ingroup  Addition
+ * 
+ * Sets \c rop to the sum of \c rop and \c op.
+ * 
+ * \todo  This is currently implemented by creating a copy!
+ */
+void _fmpq_poly_add_in_place(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
+{
+    if (rop == op)               /* If rop == op, return 2 * op              */
+    {
+        fmpq_poly_scalar_mul_si(rop, rop, 2l);
+        return;
+    }
+    
+    if (fmpq_poly_is_zero(rop))  /* Is one of the polynomials zero?          */
+    {
+        fmpq_poly_set(rop, op);
+        return;
+    }
+    if (fmpq_poly_is_zero(op))
+    {
+        return;
+    }
+    
+    /* Now we may assume that rop and op refer to distinct objects in        */
+    /* memory and that both polynomials are non-zero.                        */
+    fmpq_poly_t t;
+    fmpq_poly_init(t);
+    fmpq_poly_add(t, rop, op);
+    fmpq_poly_swap(rop, t);
+    fmpq_poly_clear(t);
+}
+
+/**
+ * \ingroup  Addition
+ * 
+ * Sets \c rop to the sum of \c op1 and \c op2.
+ */
+void fmpq_poly_add(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2)
+{
+    fmpz_poly_t tpoly;
+    fmpz_t tfmpz;
+    unsigned long limbs;
+    
+    if (op1 == op2)                         /* If op1 == op2, return 2 * op1 */
+    {
+        fmpq_poly_scalar_mul_si(rop, op1, 2l);
+        return;
+    }
+    
+    if (rop == op1)                         /* Is rop == op1 or rop == op2?  */
+    {
+        _fmpq_poly_add_in_place(rop, op2);
+        return;
+    }
+    if (rop == op2)
+    {
+        _fmpq_poly_add_in_place(rop, op1);
+        return;
+    }
+    
+    /* From here on, we may assume that rop, op1 and op2 all refer to        */
+    /* distinct objects in memory, although they may still be equal.         */
+    
+    if (_fmpq_poly_den_is_one(op1))
+    {
+        if (_fmpq_poly_den_is_one(op2))
+        {
+            fmpz_poly_add(rop->num, op1->num, op2->num);
+            if (rop->den != NULL)
+                fmpz_set_si(rop->den, 1);
+        }
+        else
+        {
+            fmpz_poly_scalar_mul_fmpz(rop->num, op1->num, op2->den);
+            fmpz_poly_add(rop->num, rop->num, op2->num);
+            fmpq_poly_set_den(rop, op2->den);
+            fmpq_poly_canonicalize(rop, NULL);
+        }
+    }
+    else
+    {
+        if (_fmpq_poly_den_is_one(op2))
+        {
+            fmpz_poly_scalar_mul_fmpz(rop->num, op2->num, op1->den);
+            fmpz_poly_add(rop->num, op1->num, rop->num);
+            fmpq_poly_set_den(rop, op1->den);
+            fmpq_poly_canonicalize(rop, NULL);
+        }
+        else
+        {
+            fmpz_poly_init(tpoly);
+            
+            limbs = fmpz_size(op1->den) + fmpz_size(op2->den);
+            _fmpq_poly_den_fit_limbs(rop, limbs);
+            
+            limbs = FLINT_MAX(limbs, fmpz_size(op2->den) + op1->num->limbs);
+            limbs = FLINT_MAX(limbs, fmpz_size(op1->den) + op2->num->limbs);
+            tfmpz = fmpz_init(limbs);
+            
+            fmpz_gcd(rop->den, op1->den, op2->den);
+            fmpz_tdiv(tfmpz, op2->den, rop->den);
+            fmpz_poly_scalar_mul_fmpz(rop->num, op1->num, tfmpz);
+            fmpz_tdiv(tfmpz, op1->den, rop->den);
+            fmpz_poly_scalar_mul_fmpz(tpoly, op2->num, tfmpz);
+            fmpz_poly_add(rop->num, rop->num, tpoly);
+            fmpz_mul(rop->den, tfmpz, op2->den);
+            
+            fmpq_poly_canonicalize(rop, tfmpz);
+            
+            fmpz_poly_clear(tpoly);
+            fmpz_clear(tfmpz);
+        }
+    }
+}
+
+/**
+ * \ingroup  Addition
+ * 
+ * Sets \c rop to the difference of \c rop and \c op.
+ * 
+ * \note  This is implemented using the methods #fmpq_poly_neg() and 
+ *        #_fmpq_poly_add_in_place().
+ */
+void _fmpq_poly_sub_in_place(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
+{
+    if (rop == op)
+    {
+        fmpq_poly_zero(rop);
+        return;
+    }
+    
+    fmpq_poly_neg(rop, rop);
+    _fmpq_poly_add_in_place(rop, op);
+    fmpq_poly_neg(rop, rop);
+}
+
+/**
+ * \ingroup  Addition
+ * 
+ * Sets \c rop to the difference of \c op1 and \c op2.
+ */
+void fmpq_poly_sub(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2)
+{
+    fmpz_poly_t tpoly;
+    fmpz_t tfmpz;
+    unsigned long limbs;
+    
+    if (op1 == op2)
+    {
+        fmpq_poly_zero(rop);
+        return;
+    }
+    if (rop == op1)
+    {
+        _fmpq_poly_sub_in_place(rop, op2);
+        return;
+    }
+    if (rop == op2)
+    {
+        _fmpq_poly_sub_in_place(rop, op1);
+        fmpq_poly_neg(rop, rop);
+        return;
+    }
+    
+    /* From here on, we know that rop, op1 and op2 refer to distinct objects */
+    /* in memory, although as rational functions they may still be equal     */
+    
+    if (_fmpq_poly_den_is_one(op1))
+    {
+        if (_fmpq_poly_den_is_one(op2))
+        {
+            fmpz_poly_sub(rop->num, op1->num, op2->num);
+            if (rop->den != NULL)
+                fmpz_set_si(rop->den, 1);
+        }
+        else
+        {
+            fmpz_poly_scalar_mul_fmpz(rop->num, op1->num, op2->den);
+            fmpz_poly_sub(rop->num, rop->num, op2->num);
+            fmpq_poly_set_den(rop, op2->den);
+            fmpq_poly_canonicalize(rop, NULL);
+        }
+    }
+    else
+    {
+        if (_fmpq_poly_den_is_one(op2))
+        {
+            fmpz_poly_scalar_mul_fmpz(rop->num, op2->num, op1->den);
+            fmpz_poly_sub(rop->num, op1->num, rop->num);
+            fmpq_poly_set_den(rop, op1->den);
+            fmpq_poly_canonicalize(rop, NULL);
+        }
+        else
+        {
+            fmpz_poly_init(tpoly);
+            
+            limbs = fmpz_size(op1->den) + fmpz_size(op2->den);
+            _fmpq_poly_den_fit_limbs(rop, limbs);
+            
+            limbs = FLINT_MAX(limbs, fmpz_size(op2->den) + op1->num->limbs);
+            limbs = FLINT_MAX(limbs, fmpz_size(op1->den) + op2->num->limbs);
+            tfmpz = fmpz_init(limbs);
+            
+            fmpz_gcd(rop->den, op1->den, op2->den);
+            fmpz_tdiv(tfmpz, op2->den, rop->den);
+            fmpz_poly_scalar_mul_fmpz(rop->num, op1->num, tfmpz);
+            fmpz_tdiv(tfmpz, op1->den, rop->den);
+            fmpz_poly_scalar_mul_fmpz(tpoly, op2->num, tfmpz);
+            fmpz_poly_sub(rop->num, rop->num, tpoly);
+            fmpz_mul(rop->den, tfmpz, op2->den);
+            
+            fmpq_poly_canonicalize(rop, tfmpz);
+            
+            fmpz_poly_clear(tpoly);
+            fmpz_clear(tfmpz);
+        }
+    }
+}
+
+/**
+ * \ingroup  Addition
+ * 
+ * Sets \c rop to <tt>rop + op1 * op2</tt>.
+ * 
+ * Currently, this method refers to the methods #fmpq_poly_mul() and 
+ * #fmpq_poly_add() to form the result in the naive way.
+ * 
+ * \todo  Implement this method more efficiently.
+ */
+void fmpq_poly_addmul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2)
+{
+    fmpq_poly_t t;
+    fmpq_poly_init(t);
+    fmpq_poly_mul(t, op1, op2);
+    fmpq_poly_add(rop, rop, t);
+    fmpq_poly_clear(t);
+}
+
+/**
+ * \ingroup  Addition
+ * 
+ * Sets \c rop to <tt>rop - op1 * op2</tt>.
+ * 
+ * Currently, this method refers to the methods #fmpq_poly_mul() and 
+ * #fmpq_poly_sub() to form the result in the naive way.
+ * 
+ * \todo  Implement this method more efficiently.
+ */
+void fmpq_poly_submul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2)
+{
+    fmpq_poly_t t;
+    fmpq_poly_init(t);
+    fmpq_poly_mul(t, op1, op2);
+    fmpq_poly_sub(rop, rop, t);
+    fmpq_poly_clear(t);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Scalar multiplication and division
+
+/**
+ * \ingroup  ScalarMul
+ * 
+ * Sets \c rop to the scalar product of \c op and the integer \c x.
+ */
+void fmpq_poly_scalar_mul_si(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long x)
+{
+    fmpz_t fx, g, t;
+    
+    if (_fmpq_poly_den_is_one(op))
+    {
+        fmpz_poly_scalar_mul_si(rop->num, op->num, x);
+        if (rop->den != NULL)
+            fmpz_set_si(rop->den, 1);
+    }
+    else
+    {
+        g  = fmpz_init(fmpz_size(op->den));
+        fx = fmpz_init(1);
+        fmpz_set_si(fx, x);
+        fmpz_gcd(g, op->den, fx);
+        fmpz_abs(g, g);
+        if (fmpz_is_one(g))
+        {
+            fmpz_poly_scalar_mul_si(rop->num, op->num, x);
+            fmpq_poly_set_den(rop, op->den);
+        }
+        else
+        {
+            if (rop == op)
+            {
+                t = fmpz_init(fmpz_size(op->den));
+                fmpz_tdiv(t, fx, g);
+                fmpz_poly_scalar_mul_fmpz(rop->num, op->num, t);
+                fmpz_tdiv(t, op->den, g);
+                fmpz_clear(rop->den);
+                rop->den = t;
+            }
+            else
+            {
+                _fmpq_poly_den_fit_limbs(rop, fmpz_size(op->den));
+                fmpz_tdiv(rop->den, fx, g);
+                fmpz_poly_scalar_mul_fmpz(rop->num, op->num, rop->den);
+                fmpz_tdiv(rop->den, op->den, g);
+            }
+        }
+        fmpz_clear(g);
+        fmpz_clear(fx);
+    }
+}
+
+/**
+ * \ingroup  ScalarMul
+ * 
+ * Sets \c rop to the scalar multiple of \c op with the \c mpz_t \c x.
+ */
+void fmpq_poly_scalar_mul_mpz(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpz_t x)
+{
+    fmpz_t fx, g, t;
+    
+    if (_fmpq_poly_den_is_one(op))
+    {
+        fmpz_poly_scalar_mul_mpz(rop->num, op->num, x);
+        if (rop->den != NULL)
+            fmpz_set_si(rop->den, 1);
+    }
+    else
+    {
+        g  = fmpz_init(fmpz_size(op->den));
+        fx = fmpz_init(mpz_size(x));
+        mpz_to_fmpz(fx, x);
+        fmpz_gcd(g, op->den, fx);
+        fmpz_abs(g, g);
+        if (fmpz_is_one(g))
+        {
+            fmpz_poly_scalar_mul_mpz(rop->num, op->num, x);
+            fmpq_poly_set_den(rop, op->den);
+        }
+        else
+        {
+            if (rop == op)
+            {
+                t = fmpz_init(fmpz_size(op->den));
+                fmpz_tdiv(t, fx, g);
+                fmpz_poly_scalar_mul_fmpz(rop->num, op->num, t);
+                fmpz_tdiv(t, op->den, g);
+                fmpz_clear(rop->den);
+                rop->den = t;
+            }
+            else
+            {
+                _fmpq_poly_den_fit_limbs(rop, fmpz_size(op->den));
+                fmpz_tdiv(rop->den, fx, g);
+                fmpz_poly_scalar_mul_fmpz(rop->num, op->num, rop->den);
+                fmpz_tdiv(rop->den, op->den, g);
+            }
+        }
+        fmpz_clear(g);
+        fmpz_clear(fx);
+    }
+}
+
+/**
+ * \ingroup  ScalarMul
+ * 
+ * Sets \c rop to the scalar multiple of \c op with the \c mpq_t \c x.
+ */
+void fmpq_poly_scalar_mul_mpq(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x)
+{
+    fmpz_t s, t;
+    
+    fmpz_poly_scalar_mul_mpz(rop->num, op->num, mpq_numref(x));
+    if (_fmpq_poly_den_is_one(op))
+    {
+        _fmpq_poly_den_fit_limbs(rop, mpz_size(mpq_denref(x)));
+        mpz_to_fmpz(rop->den, mpq_denref(x));
+    }
+    else
+    {
+        s = fmpz_init(mpz_size(mpq_denref(x)));
+        t = fmpz_init(fmpz_size(op->den));
+        mpz_to_fmpz(s, mpq_denref(x));
+        fmpz_set(t, op->den);
+        _fmpq_poly_den_fit_limbs(rop, fmpz_size(s) + fmpz_size(t));
+        fmpz_mul(rop->den, s, t);
+        fmpz_clear(s);
+        fmpz_clear(t);
+    }
+    fmpq_poly_canonicalize(rop, NULL);
+}
+
+/**
+ * /ingroup  ScalarMul
+ * 
+ * Divides \c rop by the integer \c x.  
+ * 
+ * Assumes that \c x is non-zero.  Otherwise, an exception is raised in the 
+ * form of an <tt>abort</tt> statement.
+ */
+void _fmpq_poly_scalar_div_si_in_place(fmpq_poly_ptr rop, long x)
+{
+    fmpz_t cont, fx, gcd, t;
+    
+    /* Assertion! */
+    if (x == 0l)
+    {
+        printf("ERROR (_fmpq_poly_scalar_div_si_in_place).  Division by zero.\n");
+        abort();
+    }
+    
+    if (x == 1l)
+    {
+        return;
+    }
+    
+    cont = fmpz_init(rop->num->limbs);
+    fmpz_poly_content(cont, rop->num);
+    fmpz_abs(cont, cont);
+    
+    if (fmpz_is_one(cont))
+    {
+        if (x > 0l)
+        {
+            if (rop->den == NULL)
+            {
+                rop->den = fmpz_init(1);
+                fmpz_set_si(rop->den, x);
+            }
+            else
+            {
+                t = fmpz_init(fmpz_size(rop->den) + 1);
+                fmpz_mul_ui(t, rop->den, (unsigned long) x);
+                fmpz_clear(rop->den);
+                rop->den = t;
+            }
+        }
+        else
+        {
+            fmpz_poly_neg(rop->num, rop->num);
+            if (rop->den == NULL)
+            {
+                rop->den = fmpz_init(1);
+                fmpz_set_si(rop->den, -x);
+            }
+            else
+            {
+                t = fmpz_init(fmpz_size(rop->den) + 1);
+                fmpz_mul_ui(t, rop->den, (unsigned long) -x);
+                fmpz_clear(rop->den);
+                rop->den = t;
+            }
+        }
+        fmpz_clear(cont);
+        return;
+    }
+    
+    fx = fmpz_init(1);
+    fmpz_set_si(fx, x);
+    
+    gcd = fmpz_init(FLINT_MAX(rop->num->limbs, fmpz_size(fx)));
+    fmpz_gcd(gcd, cont, fx);
+    fmpz_abs(gcd, gcd);
+    
+    if (fmpz_is_one(gcd))
+    {
+        if (x > 0l)
+        {
+            if (rop->den == NULL)
+            {
+                rop->den = fmpz_init(1);
+                fmpz_set_si(rop->den, x);
+            }
+            else
+            {
+                t = fmpz_init(fmpz_size(rop->den) + 1);
+                fmpz_mul_ui(t, rop->den, (unsigned long) x);
+                fmpz_clear(rop->den);
+                rop->den = t;
+            }
+        }
+        else
+        {
+            fmpz_poly_neg(rop->num, rop->num);
+            if (rop->den == NULL)
+            {
+                rop->den = fmpz_init(1);
+                fmpz_set_si(rop->den, -x);
+            }
+            else
+            {
+                t = fmpz_init(fmpz_size(rop->den) + 1);
+                fmpz_mul_ui(t, rop->den, (unsigned long) -x);
+                fmpz_clear(rop->den);
+                rop->den = t;
+            }
+        }
+    }
+    else
+    {
+        fmpz_poly_scalar_div_fmpz(rop->num, rop->num, gcd);
+        if (rop->den == NULL)
+        {
+            rop->den = fmpz_init(1);
+            fmpz_tdiv(rop->den, fx, gcd);
+        }
+        else
+        {
+            fmpz_tdiv(cont, fx, gcd);
+            fx = fmpz_realloc(fx, fmpz_size(rop->den));
+            fmpz_set(fx, rop->den);
+            rop->den = fmpz_realloc(rop->den, fmpz_size(rop->den) + 1);
+            fmpz_mul(rop->den, fx, cont);
+        }
+        if (x < 0l)
+        {
+            fmpz_poly_neg(rop->num, rop->num);
+            fmpz_neg(rop->den, rop->den);
+        }
+    }
+    
+    fmpz_clear(cont);
+    fmpz_clear(fx);
+    fmpz_clear(gcd);
+}
+
+
+/**
+ * \ingroup  ScalarMul
+ * 
+ * Sets \c rop to the scalar multiple of \c op with the multiplicative inverse 
+ * of the integer \c x.
+ *
+ * Assumes that \c x is non-zero.  Otherwise, an exception is raised in the 
+ * form of an <tt>abort</tt> statement.
+ */
+void fmpq_poly_scalar_div_si(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long x)
+{
+    fmpz_t cont, fx, gcd;
+    
+    /* Assertion! */
+    if (x == 0l)
+    {
+        printf("ERROR (fmpq_poly_scalar_div_si).  Division by zero.\n");
+        abort();
+    }
+    
+    if (rop == op)
+    {
+        _fmpq_poly_scalar_div_si_in_place(rop, x);
+        return;
+    }
+    
+    /* From here on, we may assume that rop and op denote two different      */
+    /* rational polynomials (as objects in memory).                          */
+    
+    if (x == 1l)
+    {
+        fmpq_poly_set(rop, op);
+        return;
+    }
+    
+    cont = fmpz_init(op->num->limbs);
+    fmpz_poly_content(cont, op->num);
+    fmpz_abs(cont, cont);
+    
+    if (fmpz_is_one(cont))
+    {
+        if (x > 0l)
+        {
+            fmpz_poly_set(rop->num, op->num);
+            if (op->den == NULL)
+            {
+                if (rop->den == NULL)
+                    rop->den = fmpz_init(1);
+                fmpz_set_si(rop->den, x);
+            }
+            else
+            {
+                _fmpq_poly_den_fit_limbs(rop, fmpz_size(op->den) + 1);
+                fmpz_mul_ui(rop->den, op->den, (unsigned long) x);
+            }
+        }
+        else
+        {
+            fmpz_poly_neg(rop->num, op->num);
+            if (op->den == NULL)
+            {
+                if (rop->den == NULL)
+                    rop->den = fmpz_init(1);
+                fmpz_set_si(rop->den, -x);
+            }
+            else
+            {
+                _fmpq_poly_den_fit_limbs(rop, fmpz_size(op->den) + 1);
+                fmpz_mul_ui(rop->den, op->den, (unsigned long) -x);
+            }
+        }
+        fmpz_clear(cont);
+        return;
+    }
+    
+    fx = fmpz_init(1);
+    fmpz_set_si(fx, x);
+    
+    gcd = fmpz_init(FLINT_MAX(op->num->limbs, fmpz_size(fx)));
+    fmpz_gcd(gcd, cont, fx);
+    fmpz_abs(gcd, gcd);
+    
+    if (fmpz_is_one(gcd))
+    {
+        if (x > 0l)
+        {
+            fmpz_poly_set(rop->num, op->num);
+            if (op->den == NULL)
+            {
+                if (rop->den == NULL)
+                    rop->den = fmpz_init(1);
+                fmpz_set_si(rop->den, x);
+            }
+            else
+            {
+                _fmpq_poly_den_fit_limbs(rop, fmpz_size(op->den) + 1);
+                fmpz_mul_ui(rop->den, op->den, (unsigned long) x);
+            }
+        }
+        else
+        {
+            fmpz_poly_neg(rop->num, op->num);
+            if (op->den == NULL)
+            {
+                if (rop->den == NULL)
+                    rop->den = fmpz_init(1);
+                fmpz_set_si(rop->den, -x);
+            }
+            else
+            {
+                _fmpq_poly_den_fit_limbs(rop, fmpz_size(op->den) + 1);
+                fmpz_mul_ui(rop->den, op->den, (unsigned long) -x);
+            }
+        }
+    }
+    else
+    {
+        fmpz_poly_scalar_div_fmpz(rop->num, op->num, gcd);
+        if (_fmpq_poly_den_is_one(op))
+        {
+            _fmpq_poly_den_fit_limbs(rop, fmpz_size(fx));
+            fmpz_tdiv(rop->den, fx, gcd);
+        }
+        else
+        {
+            _fmpq_poly_den_fit_limbs(rop, fmpz_size(op->den) + fmpz_size(fx));
+            fmpz_tdiv(cont, fx, gcd);  /* fx and gcd are word-sized */
+            fmpz_mul(rop->den, op->den, cont);
+        }
+        if (x < 0l)
+        {
+            fmpz_poly_neg(rop->num, rop->num);
+            fmpz_neg(rop->den, rop->den);
+        }
+    }
+    
+    fmpz_clear(cont);
+    fmpz_clear(fx);
+    fmpz_clear(gcd);
+}
+
+/**
+ * \ingroup  ScalarMul
+ * 
+ * Sets \c rop to the scalar multiple of \c op with the multiplicative inverse 
+ * of the integer \c x.
+ *
+ * Assumes that \c x is non-zero.  Otherwise, an exception is raised in the 
+ * form of an <tt>abort</tt> statement.
+ */
+void _fmpq_poly_scalar_div_mpz_in_place(fmpq_poly_ptr rop, const mpz_t x)
+{
+    fmpz_t cont, fx, gcd, t;
+    
+    /* Assertion! */
+    if (mpz_sgn(x) == 0)
+    {
+        printf("ERROR (_fmpq_poly_scalar_div_mpz_in_place).  Division by zero.\n");
+        abort();
+    }
+    
+    if (mpz_cmp_si(x, 1) == 0)
+        return;
+    
+    cont = fmpz_init(rop->num->limbs);
+    fmpz_poly_content(cont, rop->num);
+    fmpz_abs(cont, cont);
+    
+    if (fmpz_is_one(cont))
+    {
+        if (rop->den == NULL)
+        {
+            rop->den = fmpz_init(mpz_size(x));
+            mpz_to_fmpz(rop->den, x);
+        }
+        else
+        {
+            fx = fmpz_init(mpz_size(x));
+            mpz_to_fmpz(fx, x);
+            t = fmpz_init(fmpz_size(rop->den) + mpz_size(x));
+            fmpz_mul(t, rop->den, fx);
+            fmpz_clear(rop->den);
+            rop->den = t;
+            fmpz_clear(fx);
+        }
+        if (mpz_sgn(x) < 0)
+        {
+            fmpz_poly_neg(rop->num, rop->num);
+            fmpz_neg(rop->den, rop->den);
+        }
+        fmpz_clear(cont);
+        return;
+    }
+    
+    fx = fmpz_init(mpz_size(x));
+    mpz_to_fmpz(fx, x);
+    
+    gcd = fmpz_init(FLINT_MAX(rop->num->limbs, fmpz_size(fx)));
+    fmpz_gcd(gcd, cont, fx);
+    fmpz_abs(gcd, gcd);
+    
+    if (fmpz_is_one(gcd))
+    {
+        if (rop->den == NULL)
+        {
+            rop->den = fmpz_init(mpz_size(x));
+            mpz_to_fmpz(rop->den, x);
+        }
+        else
+        {
+            t = fmpz_init(fmpz_size(rop->den) + mpz_size(x));
+            fmpz_mul(t, rop->den, fx);
+            fmpz_clear(rop->den);
+            rop->den = t;
+        }
+        if (mpz_sgn(x) < 0)
+        {
+            fmpz_poly_neg(rop->num, rop->num);
+            fmpz_neg(rop->den, rop->den);
+        }
+    }
+    else
+    {
+        fmpz_poly_scalar_div_fmpz(rop->num, rop->num, gcd);
+        if (rop->den == NULL)
+        {
+            rop->den = fmpz_init(fmpz_size(fx));
+            fmpz_tdiv(rop->den, fx, gcd);
+        }
+        else
+        {
+            cont = fmpz_realloc(cont, fmpz_size(fx));
+            fmpz_tdiv(cont, fx, gcd);
+            fx = fmpz_realloc(fx, fmpz_size(rop->den));
+            fmpz_set(fx, rop->den);
+            _fmpq_poly_den_fit_limbs(rop, fmpz_size(rop->den) + fmpz_size(cont));
+            fmpz_mul(rop->den, fx, cont);
+        }
+        if (mpz_sgn(x) < 0)
+        {
+            fmpz_poly_neg(rop->num, rop->num);
+            fmpz_neg(rop->den, rop->den);
+        }
+    }
+    
+    fmpz_clear(cont);
+    fmpz_clear(fx);
+    fmpz_clear(gcd);
+}
+
+/**
+ * \ingroup  ScalarMul
+ * 
+ * Sets \c rop to the scalar multiple of \c op with the multiplicative inverse 
+ * of the integer \c x.
+ *
+ * Assumes that \c x is non-zero.  Otherwise, an exception is raised in the 
+ * form of an <tt>abort</tt> statement.
+ */
+void fmpq_poly_scalar_div_mpz(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpz_t x)
+{
+    fmpz_t cont, fx, gcd, t;
+    
+    /* Assertion! */
+    if (mpz_sgn(x) == 0)
+    {
+        printf("ERROR (fmpq_poly_scalar_div_mpz).  Division by zero.\n");
+        abort();
+    }
+    
+    if (rop == op)
+    {
+        _fmpq_poly_scalar_div_mpz_in_place(rop, x);
+        return;
+    }
+    
+    /* From here on, we may assume that rop and op denote two different      */
+    /* rational polynomials (as objects in memory).                          */
+    
+    if (mpz_cmp_si(x, 1) == 0)
+    {
+        fmpq_poly_set(rop, op);
+        return;
+    }
+    
+    cont = fmpz_init(op->num->limbs);
+    fmpz_poly_content(cont, op->num);
+    fmpz_abs(cont, cont);
+    
+    if (fmpz_is_one(cont))
+    {
+        if (op->den == NULL)
+        {
+            if (rop->den == NULL)
+                rop->den = fmpz_init(mpz_size(x));
+            mpz_to_fmpz(rop->den, x);
+        }
+        else
+        {
+            t = fmpz_init(mpz_size(x));
+            mpz_to_fmpz(t, x);
+            _fmpq_poly_den_fit_limbs(rop, fmpz_size(op->den) + fmpz_size(t));
+            fmpz_mul(rop->den, op->den, t);
+            fmpz_clear(t);
+        }
+        if (mpz_sgn(x) > 0)
+        {
+            fmpz_poly_set(rop->num, op->num);
+        }
+        else
+        {
+            fmpz_poly_neg(rop->num, op->num);
+            fmpz_neg(rop->den, rop->den);
+        }
+        fmpz_clear(cont);
+        return;
+    }
+    
+    fx = fmpz_init(mpz_size(x));
+    mpz_to_fmpz(fx, x);
+    
+    gcd = fmpz_init(FLINT_MAX(op->num->limbs, fmpz_size(fx)));
+    fmpz_gcd(gcd, cont, fx);
+    fmpz_abs(gcd, gcd);
+    
+    if (fmpz_is_one(gcd))
+    {
+        if (op->den == NULL)
+        {
+            if (rop->den == NULL)
+                rop->den = fmpz_init(mpz_size(x));
+            mpz_to_fmpz(rop->den, x);
+        }
+        else
+        {
+            t = fmpz_init(mpz_size(x));
+            mpz_to_fmpz(t, x);
+            _fmpq_poly_den_fit_limbs(rop, fmpz_size(op->den) + fmpz_size(t));
+            fmpz_mul(rop->den, op->den, t);
+            fmpz_clear(t);
+        }
+        if (mpz_sgn(x) > 0)
+        {
+            fmpz_poly_set(rop->num, op->num);
+        }
+        else
+        {
+            fmpz_poly_neg(rop->num, op->num);
+            fmpz_neg(rop->den, rop->den);
+        }
+    }
+    else
+    {
+        fmpz_poly_scalar_div_fmpz(rop->num, op->num, gcd);
+        if (op->den == NULL)
+        {
+            if (rop->den == NULL)
+            {
+                rop->den = fmpz_init(fmpz_size(fx));
+                fmpz_tdiv(rop->den, fx, gcd);
+            }
+            else
+            {
+                rop->den = fmpz_realloc(rop->den, fmpz_size(fx));
+                fmpz_tdiv(rop->den, fx, gcd);
+            }
+        }
+        else
+        {
+            cont = fmpz_realloc(cont, fmpz_size(fx));
+            fmpz_tdiv(cont, fx, gcd);
+            _fmpq_poly_den_fit_limbs(rop, fmpz_size(op->den) + fmpz_size(cont));
+            fmpz_mul(rop->den, op->den, cont);
+        }
+        if (mpz_sgn(x) < 0)
+        {
+            fmpz_poly_neg(rop->num, rop->num);
+            fmpz_neg(rop->den, rop->den);
+        }
+    }
+    
+    fmpz_clear(cont);
+    fmpz_clear(fx);
+    fmpz_clear(gcd);
+}
+
+/**
+ * \ingroup  ScalarMul
+ * 
+ * Sets \c rop to the scalar multiple of \c op with the multiplicative inverse 
+ * of the rational \c x.
+ *
+ * Assumes that the rational \c x is in lowest terms and non-zero.  If the 
+ * rational is not in lowest terms, the resulting value of <tt>rop</tt> is 
+ * undefined.  If <tt>x</tt> is zero, an exception is raised in the form 
+ * of an <tt>abort</tt> statement.
+ */
+void fmpq_poly_scalar_div_mpq(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x)
+{
+    fmpz_t s, t;
+    
+    /* Assertion! */
+    if (mpz_sgn(mpq_numref(x)) == 0)
+    {
+        printf("ERROR (fmpq_poly_scalar_div_mpq).  Division by zero.\n");
+        abort();
+    }
+   
+    fmpz_poly_scalar_mul_mpz(rop->num, op->num, mpq_denref(x));
+    if (_fmpq_poly_den_is_one(op))
+    {
+        _fmpq_poly_den_fit_limbs(rop, mpz_size(mpq_numref(x)));
+        mpz_to_fmpz(rop->den, mpq_numref(x));
+    }
+    else
+    {
+        s = fmpz_init(mpz_size(mpq_numref(x)));
+        t = fmpz_init(fmpz_size(op->den));
+        mpz_to_fmpz(s, mpq_numref(x));
+        fmpz_set(t, op->den);
+        _fmpq_poly_den_fit_limbs(rop, fmpz_size(s) + fmpz_size(t));
+        fmpz_mul(rop->den, s, t);
+        fmpz_clear(s);
+        fmpz_clear(t);
+    }
+    fmpq_poly_canonicalize(rop, NULL);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Multiplication
+
+/**
+ * \ingroup  Multiplication
+ * 
+ * Multiplies <tt>rop</tt> by <tt>op</tt>.
+ */
+void _fmpq_poly_mul_in_place(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
+{
+    fmpq_poly_t t;
+    
+    if (rop == op)
+    {
+        fmpz_poly_power(rop->num, op->num, 2ul);
+        if (rop->den != NULL)
+        {
+            rop->den = fmpz_realloc(rop->den, 2 * fmpz_size(rop->den));
+            fmpz_pow_ui(rop->den, op->den, 2ul);
+        }
+        return;
+    }
+    
+    if (fmpq_poly_is_zero(rop) || fmpq_poly_is_zero(op))
+    {
+        fmpq_poly_zero(rop);
+        return;
+    }
+    
+    /* From here on, rop and op point to two different objects in memory,    */
+    /* and these are both non-zero rational polynomials.                     */
+    
+    fmpq_poly_init(t);
+    fmpq_poly_mul(t, rop, op);
+    fmpq_poly_swap(rop, t);
+    fmpq_poly_clear(t);
+}
+
+/**
+ * \ingroup  Multiplication
+ * 
+ * Sets \c rop to the product of \c op1 and \c op2.
+ */
+void fmpq_poly_mul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2)
+{
+    unsigned long limbs;
+    fmpz_t t, t1, t2;
+    
+    if (op1 == op2)
+    {
+        fmpz_poly_power(rop->num, op1->num, 2ul);
+        if (_fmpq_poly_den_is_one(op1))
+        {
+            if (rop->den != NULL)
+                fmpz_set_si(rop->den, 1);
+        }
+        else
+        {
+            if (rop == op1)
+            {
+                t = fmpz_init(2 * fmpz_size(op1->den));
+                fmpz_pow_ui(t, op1->den, 2ul);
+                fmpz_clear(rop->den);
+                rop->den = t;
+            }
+            else
+            {
+                _fmpq_poly_den_fit_limbs(rop, 2 * fmpz_size(op1->den));
+                fmpz_pow_ui(rop->den, op1->den, 2ul);
+            }
+        }
+        return;
+    }
+    
+    if (rop == op1)
+    {
+        _fmpq_poly_mul_in_place(rop, op2);
+        return;
+    }
+    if (rop == op2)
+    {
+        _fmpq_poly_mul_in_place(rop, op1);
+        return;
+    }
+    
+    if (_fmpq_poly_den_is_one(op1))
+    {
+        /* Case 1.  a.den == b.den == 1 */
+        if (_fmpq_poly_den_is_one(op2))
+        {
+            fmpz_poly_mul(rop->num, op1->num, op2->num);
+            if (rop->den != NULL)
+                fmpz_set_si(rop->den, 1);
+        }
+        /* Case 2.  a.den == 1, b.den > 1 */
+        else
+        {
+            _fmpq_poly_den_fit_limbs(rop, FLINT_MAX(op1->num->limbs, 
+                fmpz_size(op2->den)));
+            
+            fmpz_poly_content(rop->den, op1->num);
+            fmpz_abs(rop->den, rop->den);
+            
+            if (fmpz_is_one(rop->den))  /* The content of a.num is 1         */
+            {
+                fmpz_poly_mul(rop->num, op1->num, op2->num);
+                fmpz_set(rop->den, op2->den);
+            }
+            else                        /* The content of a.num exceeds 1    */
+            {
+                t1 = fmpz_init(FLINT_MAX(op1->num->limbs, fmpz_size(op2->den)));
+                fmpz_gcd(t1, rop->den, op2->den);
+                fmpz_abs(t1, t1);
+                if (fmpz_is_one(t1))
+                {
+                    fmpz_poly_mul(rop->num, op1->num, op2->num);
+                    fmpz_set(rop->den, op2->den);
+                }
+                else
+                {
+                    fmpz_poly_scalar_div_fmpz(rop->num, op1->num, t1);
+                    fmpz_poly_mul(rop->num, rop->num, op2->num);
+                    fmpz_tdiv(rop->den, op2->den, t1);
+                }
+                fmpz_clear(t1);
+            }
+        }
+    }
+    else
+    {
+        /* Case 3.  a.den > 1, b.den == 1 */
+        if (_fmpq_poly_den_is_one(op2))
+            fmpq_poly_mul(rop, op2, op1);
+        
+        /* Case 4.  a.den > 1, b.den > 1 */
+        else
+        {
+            limbs = FLINT_MAX(op1->num->limbs, op2->num->limbs);
+            limbs = FLINT_MAX(limbs, fmpz_size(op1->den) + fmpz_size(op2->den));
+            _fmpq_poly_den_fit_limbs(rop, limbs);
+            
+            fmpz_poly_content(rop->den, op1->num);
+            fmpz_abs(rop->den, rop->den);
+            
+            if (fmpz_is_one(rop->den))
+            {
+                fmpz_poly_content(rop->den, op2->num);
+                fmpz_abs(rop->den, rop->den);
+                
+                /* Case 4.A.  c(a.num) == c(b.num) == 1 */
+                if (fmpz_is_one(rop->den))
+                {
+                    fmpz_poly_mul(rop->num, op1->num, op2->num);
+                    fmpz_mul(rop->den, op1->den, op2->den);
+                }
+                
+                /* Case 4.B.  c(a.num) == 1, c(b.num) > 1 */
+                else
+                {
+                    t1 = fmpz_init(FLINT_MAX(op2->num->limbs, fmpz_size(op1->den)));
+                    fmpz_gcd(t1, rop->den, op1->den);
+                    fmpz_abs(t1, t1);
+                    if (fmpz_is_one(t1))
+                    {
+                        fmpz_poly_mul(rop->num, op1->num, op2->num);
+                        fmpz_mul(rop->den, op1->den, op2->den);
+                    }
+                    else
+                    {
+                        fmpz_poly_scalar_div_fmpz(rop->num, op2->num, t1);
+                        fmpz_poly_mul(rop->num, op1->num, rop->num);
+                        fmpz_tdiv(rop->den, op1->den, t1);
+                        fmpz_mul(t1, rop->den, op2->den);
+                        fmpz_clear(rop->den);
+                        rop->den = t1;
+                        t1 = NULL;
+                    }
+                    if (t1 != NULL)
+                        fmpz_clear(t1);
+                }
+            }
+            else  /* The content of a.num exceeds 1 */
+            {
+                limbs = FLINT_MAX(op1->num->limbs, op2->num->limbs);
+                limbs = FLINT_MAX(limbs, fmpz_size(op1->den) + fmpz_size(op2->den));
+                t1 = fmpz_init(limbs);
+                
+                fmpz_poly_content(t1, op2->num);
+                fmpz_abs(t1, t1);
+                
+                /* Case 4.C.  c(a.num) > 1, c(b.num) == 1 */
+                if (fmpz_is_one(t1))
+                {
+                    fmpz_gcd(t1, rop->den, op2->den);
+                    fmpz_abs(t1, t1);
+                    if (fmpz_is_one(t1))
+                    {
+                        fmpz_poly_mul(rop->num, op1->num, op2->num);
+                        fmpz_mul(rop->den, op1->den, op2->den);
+                    }
+                    else
+                    {
+                        fmpz_poly_scalar_div_fmpz(rop->num, op1->num, t1);
+                        fmpz_poly_mul(rop->num, rop->num, op2->num);
+                        fmpz_tdiv(rop->den, op2->den, t1);
+                        fmpz_mul(t1, op1->den, rop->den);
+                        fmpz_clear(rop->den);
+                        rop->den = t1;
+                        t1 = NULL;
+                    }
+                }
+                
+                /* Case 4.D.  This is the general case:  positive            */
+                /* denominators and positive contents.  At this point,       */
+                /* res.den == c(a.num) and t1 == c(b.num).                   */
+                else
+                {
+                    limbs = FLINT_MAX(op1->num->limbs, op2->num->limbs);
+                    limbs = FLINT_MAX(limbs, fmpz_size(op1->den) + fmpz_size(op2->den));
+                    t2 = fmpz_init(limbs);
+                    
+                    fmpz_gcd(t2, rop->den, op2->den);
+                    fmpz_abs(t2, t2);
+                    
+                    if (fmpz_is_one(t2))
+                    {
+                        fmpz_gcd(t2, t1, op1->den);
+                        fmpz_abs(t2, t2);
+                        
+                        /* Case 4.D.1.  c(a.num) and b.den are coprime,      */
+                        /*              c(b.num) and a.den are coprime       */
+                        if (fmpz_is_one(t2))
+                        {
+                            fmpz_poly_mul(rop->num, op1->num, op2->num);
+                            fmpz_mul(rop->den, op1->den, op2->den);
+                        }
+                        /* Case 4.D.2.  c(a.num) and b.den are coprime,      */
+                        /*              c(b.num) and a.den are *not* coprime */
+                        else
+                        {
+                            fmpz_poly_scalar_div_fmpz(rop->num, op2->num, t2);
+                            fmpz_poly_mul(rop->num, op1->num, rop->num);
+                            fmpz_tdiv(t1, op1->den, t2);
+                            fmpz_mul(rop->den, t1, op2->den);
+                        }
+                    }
+                    else
+                    {
+                        fmpz_gcd(rop->den, t1, op1->den);
+                        fmpz_abs(rop->den, rop->den);
+                        
+                        /* Case 4.D.3.  c(a.num) and b.den are *not* coprime,*/
+                        //              c(b.num) and a.den are coprime       */
+                        if (fmpz_is_one(rop->den))
+                        {
+                            fmpz_poly_scalar_div_fmpz(rop->num, op1->num, t2);
+                            fmpz_poly_mul(rop->num, rop->num, op2->num);
+                            fmpz_tdiv(t1, op2->den, t2);
+                            fmpz_mul(rop->den, op1->den, t1);
+                        }
+                        
+                        /* Case 4.D.4.  c(a.num) and b.den are *not* coprime,*/
+                        /*              c(b.num) and a.den are *not* coprime */
+                        /*                                                   */
+                        /* NOTES:  How do we best handle this case?  The     */
+                        /* problem is that we can't store both               */
+                        /* ``a.num / res.den`` and ``b.num / t2`` in the one */
+                        /* polynomial ``res.num`` simultaneously.            */
+                        /* Suggestion:  In order to avoid allocating another */
+                        /* polynomial, only cancel the common factor in one  */
+                        /* polynomial.  Choose the larger one.               */
+                        else
+                        {
+                            if (fmpz_cmpabs(t2, rop->den) < 0)
+                            {
+                                fmpz_poly_scalar_div_fmpz(rop->num, op2->num, rop->den);
+                                fmpz_poly_mul(rop->num, op1->num, rop->num);
+                                fmpz_poly_scalar_div_fmpz(rop->num, rop->num, t2);
+                            }
+                            else
+                            {
+                                fmpz_poly_scalar_div_fmpz(rop->num, op1->num, t2);
+                                fmpz_poly_mul(rop->num, rop->num, op2->num);
+                                fmpz_poly_scalar_div_fmpz(rop->num, rop->num, rop->den);
+                            }
+                            fmpz_tdiv(t1, op1->den, rop->den);
+                            fmpz_tdiv(rop->den, op2->den, t2);
+                            fmpz_mul(t2, t1, rop->den);
+                            fmpz_clear(rop->den);
+                            rop->den = t2;
+                            t2 = NULL;
+                        }
+                    }
+                    if (t2 != NULL)
+                        fmpz_clear(t2);
+                }
+                if (t1 != NULL)
+                    fmpz_clear(t1);
+            }
+        }
+    }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Division
+
+/**
+ * \ingroup  Division
+ * 
+ * Returns the quotient of the Euclidean division of \c a by \c b.
+ * 
+ * Assumes that \c b is non-zero.  Otherwise, an exception is raised in the 
+ * form of an <tt>abort</tt> statement.
+ */
+void fmpq_poly_floordiv(fmpq_poly_ptr q, const fmpq_poly_ptr a, const fmpq_poly_ptr b)
+{
+    unsigned long limbs, m;
+    fmpz_t lead, t;
+    fmpq_poly_t tpoly;
+    
+    /* Assertion! */
+    if (fmpq_poly_is_zero(b))
+    {
+        printf("ERROR (fmpq_poly_floordiv).  Division by zero.\n");
+        abort();
+    }
+    
+    /* Catch the case when a and b are the same objects in memory.           */
+    /* As of FLINT version 1.5.1, there is a bug in this case.               */
+    if (a == b)
+    {
+        fmpq_poly_set_si(q, 1);
+        return;
+    }
+    
+    /* Deal with the various other cases of aliasing.                        */
+    if (q == a | q == b)
+    {
+        fmpq_poly_init(tpoly);
+        fmpq_poly_floordiv(tpoly, a, b);
+        fmpq_poly_swap(q, tpoly);
+        fmpq_poly_clear(tpoly);
+        return;
+    }
+    
+    /* Deal separately with the case deg(b) = 0.                             */
+    if (fmpq_poly_degree(b) == 0)
+    {
+        lead = fmpz_poly_get_coeff_ptr(b->num, 0);
+        if (_fmpq_poly_den_is_one(a))
+        {
+            if (_fmpq_poly_den_is_one(b))  /* a->den == b->den == 1          */
+            {
+                fmpz_poly_set(q->num, a->num);
+                fmpq_poly_set_den(q, lead);
+                fmpq_poly_canonicalize(q, NULL);
+            }
+            else                           /* a->den == 1, b->den > 1        */
+            {
+                fmpz_poly_scalar_mul_fmpz(q->num, a->num, b->den);
+                fmpq_poly_set_den(q, lead);
+                fmpq_poly_canonicalize(q, NULL);
+            }
+        }
+        else
+        {
+            if (_fmpq_poly_den_is_one(b))  /* a->den > 1, b->den == 1        */
+            {
+                fmpz_poly_set(q->num, a->num);
+                _fmpq_poly_den_fit_limbs(q, fmpz_size(a->den)+fmpz_size(lead));
+                fmpz_mul(q->den, a->den, lead);
+                fmpq_poly_canonicalize(q, NULL);
+            }
+            else                           /* a->den, b->den > 1             */
+            {
+                fmpz_poly_scalar_mul_fmpz(q->num, a->num, b->den);
+                _fmpq_poly_den_fit_limbs(q, fmpz_size(a->den)+fmpz_size(lead));
+                fmpz_mul(q->den, a->den, lead);
+                fmpq_poly_canonicalize(q, NULL);
+            }
+        }
+        return;
+    }
+    
+    /* General case..                                                        */
+    /* Set q to b->den q->num / (a->den lead^m).                             */
+    
+    fmpz_poly_pseudo_div(q->num, &m, a->num, b->num);
+    
+    lead = fmpz_poly_lead(b->num);
+    
+    /* Case 1:  lead^m is 1 */
+    if (fmpz_is_one(lead) || m == 0 || (fmpz_is_m1(lead) & m % 2 == 0))
+    {
+        if (!_fmpq_poly_den_is_one(b))
+            fmpz_poly_scalar_mul_fmpz(q->num, q->num, b->den);
+        fmpq_poly_set_den(q, a->den);
+        fmpq_poly_canonicalize(q, NULL);
+    }
+    /* Case 2:  lead^m is -1 */
+    else if (fmpz_is_m1(lead) & m % 2)
+    {
+        if (!_fmpq_poly_den_is_one(b))
+            fmpz_poly_scalar_mul_fmpz(q->num, q->num, b->den);
+        fmpz_poly_neg(q->num, q->num);
+        fmpq_poly_set_den(q, a->den);
+        fmpq_poly_canonicalize(q, NULL);
+    }
+    /* Case 3:  lead^m is not +-1 */
+    else
+    {
+        if (!_fmpq_poly_den_is_one(b))
+            fmpz_poly_scalar_mul_fmpz(q->num, q->num, b->den);
+        
+        limbs = m * fmpz_size(lead);
+        
+        if (_fmpq_poly_den_is_one(a))
+        {
+            _fmpq_poly_den_fit_limbs(q, limbs);
+            fmpz_pow_ui(q->den, lead, m);
+        }
+        else
+        {
+            t = fmpz_init(limbs);
+            _fmpq_poly_den_fit_limbs(q, limbs + fmpz_size(a->den));
+            fmpz_pow_ui(t, lead, m);
+            fmpz_mul(q->den, t, a->den);
+            fmpz_clear(t);
+        }
+        fmpq_poly_canonicalize(q, NULL);
+    }
+}
+
+/**
+ * \ingroup  Division
+ * 
+ * Sets \c r to the remainder of the Euclidean division of \c a by \c b.
+ * 
+ * Assumes that \c b is non-zero.  Otherwise, an exception is raised in the 
+ * form of an <tt>abort</tt> statement.
+ */
+void fmpq_poly_mod(fmpq_poly_ptr r, const fmpq_poly_ptr a, const fmpq_poly_ptr b)
+{
+    unsigned long limbs, m;
+    fmpz_t lead, t;
+    fmpq_poly_t tpoly;
+    
+    /* Assertion! */
+    if (fmpq_poly_is_zero(b))
+    {
+        printf("ERROR (fmpq_poly_mod).  Division by zero.\n");
+        abort();
+    }
+    
+    /* Catch the case when a and b are the same objects in memory.           */
+    /* As of FLINT version 1.5.1, there is a bug in this case.               */
+    if (a == b)
+    {
+        fmpq_poly_set_si(r, 0);
+        return;
+    }
+    
+    /* Deal with the various other cases of aliasing.                        */
+    if (r == a | r == b)
+    {
+        fmpq_poly_init(tpoly);
+        fmpq_poly_mod(tpoly, a, b);
+        fmpq_poly_swap(r, tpoly);
+        fmpq_poly_clear(tpoly);
+        return;
+    }
+    
+    /* Deal separately with the case deg(b) = 0.                             */
+    if (fmpq_poly_degree(b) == 0)
+    {
+        fmpz_poly_zero(r->num);
+        if (r->den != NULL)
+            fmpz_set_si(r->den, 1);
+        return;
+    }
+    
+    /* General case..                                                        */
+    /* Set r to r->num / (a->den lead^m).                                    */
+    
+    fmpz_poly_pseudo_rem(r->num, &m, a->num, b->num);
+    
+    lead = fmpz_poly_lead(b->num);
+    
+    /* Case 1:  lead^m is 1 */
+    if (fmpz_is_one(lead) || m == 0 || (fmpz_is_m1(lead) & m % 2 == 0))
+    {
+        fmpq_poly_set_den(r, a->den);
+        fmpq_poly_canonicalize(r, NULL);
+    }
+    /* Case 2:  lead^m is -1 */
+    else if (fmpz_is_m1(lead) & m % 2)
+    {
+        fmpq_poly_set_den(r, a->den);
+        fmpz_neg(r->den, r->den);
+        fmpq_poly_canonicalize(r, NULL);
+    }
+    /* Case 3:  lead^m is not +-1 */
+    else
+    {
+        limbs = m * fmpz_size(lead);
+        
+        if (_fmpq_poly_den_is_one(a))
+        {
+            _fmpq_poly_den_fit_limbs(r, limbs);
+            fmpz_pow_ui(r->den, lead, m);
+        }
+        else
+        {
+            t = fmpz_init(limbs);
+            _fmpq_poly_den_fit_limbs(r, limbs + fmpz_size(a->den));
+            fmpz_pow_ui(t, lead, m);
+            fmpz_mul(r->den, t, a->den);
+            fmpz_clear(t);
+        }
+        fmpq_poly_canonicalize(r, NULL);
+    }
+}
+
+/**
+ * \ingroup  Division
+ * 
+ * Sets \c q and \c r to the quotient and remainder of the Euclidean 
+ * division of \c a by \c b.
+ * 
+ * Assumes that \c b is non-zero, and that \c q and \c r refer to distinct 
+ * objects in memory.
+ */
+void fmpq_poly_divrem(fmpq_poly_ptr q, fmpq_poly_ptr r, const fmpq_poly_ptr a, const fmpq_poly_ptr b)
+{
+    unsigned long limbs, m;
+    fmpz_t lead, t;
+    fmpq_poly_t tempq, tempr;
+    
+    /* Assertion! */
+    if (fmpq_poly_degree(b))
+    {
+        printf("ERROR (fmpq_poly_divrem).  Division by zero.\n");
+        abort();
+    }
+    if (q == r)
+    {
+        printf("ERROR (fmpq_poly_divrem).  Output arguments aliased.\n");
+        abort();
+    }
+    
+    /* Catch the case when a and b are the same objects in memory.           */
+    /* As of FLINT version 1.5.1, there is a bug in this case.               */
+    if (a == b)
+    {
+        fmpq_poly_set_si(q, 1);
+        fmpq_poly_set_si(r, 0);
+        return;
+    }
+    
+    /* Deal with the various other cases of aliasing.                        */
+    if (r == a | r == b)
+    {
+        if (q == a | q == b)
+        {
+            fmpq_poly_init(tempq);
+            fmpq_poly_init(tempr);
+            fmpq_poly_divrem(tempq, tempr, a, b);
+            fmpq_poly_swap(q, tempq);
+            fmpq_poly_swap(r, tempr);
+            fmpq_poly_clear(tempq);
+            fmpq_poly_clear(tempr);
+            return;
+        }
+        else
+        {
+            fmpq_poly_init(tempr);
+            fmpq_poly_divrem(q, tempr, a, b);
+            fmpq_poly_swap(r, tempr);
+            fmpq_poly_clear(tempr);
+            return;
+        }
+    }
+    else
+    {
+        if (q == a | q == b)
+        {
+            fmpq_poly_init(tempq);
+            fmpq_poly_divrem(tempq, r, a, b);
+            fmpq_poly_swap(q, tempq);
+            fmpq_poly_clear(tempq);
+            return;
+        }
+    }
+    
+    // TODO: Implement the case `\deg(b) = 0` more efficiently!
+    
+    /* General case..                                                        */
+    /* Set q to b->den q->num / (a->den lead^m)                              */
+    /* and r to r->num / (a->den lead^m).                                    */
+    
+    fmpz_poly_pseudo_divrem(q->num, r->num, &m, a->num, b->num);
+    
+    lead = fmpz_poly_lead(b->num);
+    
+    /* Case 1.  lead^m is 1 */
+    if (fmpz_is_one(lead) || m == 0 || (fmpz_is_m1(lead) & m % 2 == 0))
+    {
+        if (!_fmpq_poly_den_is_one(b))
+            fmpz_poly_scalar_mul_fmpz(q->num, q->num, b->den);
+        fmpq_poly_set_den(q, a->den);
+        fmpq_poly_canonicalize(q, NULL);
+        
+        fmpq_poly_set_den(r, a->den);
+        fmpq_poly_canonicalize(r, NULL);
+    }
+    /* Case 2.  lead^m is -1 */
+    else if (fmpz_is_m1(lead) & m % 2)
+    {
+        if (!_fmpq_poly_den_is_one(b))
+            fmpz_poly_scalar_mul_fmpz(q->num, q->num, b->den);
+        fmpz_poly_neg(q->num, q->num);
+        fmpq_poly_set_den(q, a->den);
+        fmpq_poly_canonicalize(q, NULL);
+        
+        fmpz_poly_neg(r->num, r->num);
+        fmpq_poly_set_den(r, a->den);
+        fmpq_poly_canonicalize(r, NULL);
+    }
+    /* Case 3.  lead^m is not +-1 */
+    else
+    {
+        if (!_fmpq_poly_den_is_one(b))
+            fmpz_poly_scalar_mul_fmpz(q->num, q->num, b->den);
+        
+        limbs = m * fmpz_size(lead);
+        
+        if (_fmpq_poly_den_is_one(a))
+        {
+            _fmpq_poly_den_fit_limbs(q, limbs);
+            fmpz_pow_ui(q->den, lead, m);
+            
+            fmpq_poly_set_den(r, q->den);
+        }
+        else
+        {
+            t = fmpz_init(limbs);
+            _fmpq_poly_den_fit_limbs(q, limbs + fmpz_size(a->den));
+            fmpz_pow_ui(t, lead, m);
+            fmpz_mul(q->den, t, a->den);
+            fmpz_clear(t);
+            
+            fmpq_poly_set_den(r, q->den);
+        }
+        fmpq_poly_canonicalize(q, NULL);
+        fmpq_poly_canonicalize(r, NULL);
+    }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Powering
+
+/**
+ * \ingroup  Powering
+ * 
+ * Sets \c rop to the <tt>exp</tt>th power of \c op.
+ *
+ * The corner case of <tt>exp == 0</tt> is handled by setting \c rop to the 
+ * constant function \f$1\f$.  Note that this includes the case \f$0^0 = 1\f$.
+ */
+void fmpq_poly_power(fmpq_poly_ptr rop, const fmpq_poly_ptr op, unsigned long exp)
+{
+    fmpz_t t;
+    
+    if (exp == 0ul)
+    {
+        fmpq_poly_one(rop);
+    }
+    else
+    {
+        fmpz_poly_power(rop->num, op->num, exp);
+        if (_fmpq_poly_den_is_one(op))
+        {
+            if (rop->den != NULL)
+                fmpz_set_si(rop->den, 1);
+        }
+        else
+        {
+            if (rop == op)        /* op->den != NULL, hence rop->den != NULL */
+            {
+                t = fmpz_init(exp * fmpz_size(op->den));
+                fmpz_pow_ui(t, op->den, exp);
+                fmpz_clear(rop->den);
+                rop->den = t;
+            }
+            else
+            {
+                _fmpq_poly_den_fit_limbs(rop, exp * fmpz_size(op->den));
+                fmpz_pow_ui(rop->den, op->den, exp);
+            }
+        }
+    }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Greatest common divisor
+
+/**
+ * \ingroup  GCD
+ * 
+ * Returns the (monic) greatest common divisor \c res of \c a and \c b.
+ * 
+ * Corner cases:  If \c a and \c b are both zero, returns zero.  If only 
+ * one of them is zero, returns the other polynomial, up to normalisation.
+ */
+void fmpq_poly_gcd(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b)
+{
+    fmpz_t lead, t;
+    fmpz_poly_t num;
+    fmpq_poly_t tpoly;
+    
+    /* Deal with aliasing */
+    if (rop == a | rop == b)
+    {
+        fmpq_poly_init(tpoly);
+        fmpq_poly_gcd(tpoly, a, b);
+        fmpq_poly_swap(rop, tpoly);
+        fmpq_poly_clear(tpoly);
+        return;
+    }
+    
+    /* Deal with corner cases */
+    if (fmpq_poly_is_zero(a))
+    {
+        if (fmpq_poly_is_zero(b))  /* a == b == 0                            */
+            fmpq_poly_zero(rop);
+        else                       /* a == 0, b != 0                         */
+            fmpq_poly_monic(rop, b);
+        return;
+    }
+    else
+    {
+        if (fmpq_poly_is_zero(b))  /* a != 0, b == 0                         */
+        {
+            fmpq_poly_monic(rop, a);
+            return;
+        }
+    }
+    
+    /* General case..                                                        */
+    fmpz_poly_init(num);
+    fmpz_poly_primitive_part(rop->num, a->num);
+    fmpz_poly_primitive_part(num, b->num);
+    
+    /* Since rop.num is the greatest common divisor of the primitive parts   */
+    /* of a.num and b.num, it is also primitive.  But as of FLINT 1.4.0, the */
+    /* leading term *might* be negative.                                     */
+    fmpz_poly_gcd(rop->num, rop->num, num);
+    
+    lead = fmpz_poly_lead(rop->num);
+    if (fmpz_sgn(lead) < 0)
+        fmpz_poly_neg(rop->num, rop->num);
+    fmpq_poly_set_den(rop, lead);
+    
+    fmpz_poly_clear(num);
+}
+
+/**
+ * \ingroup  GCD
+ * 
+ * Returns polynomials \c s, \c t and \c rop such that \c rop is 
+ * (monic) greatest common divisor of \c a and \c b, and such that 
+ * <tt>rop = s a + t b</tt>.
+ * 
+ * Corner cases:  If \c a and \c b are zero, returns zero polynomials. 
+ * Otherwise, if only \c a is zero, returns <tt>(res, s, t) = (b, 0, 1)</tt> 
+ * up to normalisation, and similarly if only \c b is zero.
+ * 
+ * Assumes that the output parameters \c rop, \c s, and \c t refer to 
+ * distinct objects in memory.  Otherwise, an exception is raised in the 
+ * form of an <tt>abort</tt> statement.
+ */
+void fmpq_poly_xgcd(fmpq_poly_ptr rop, fmpq_poly_ptr s, fmpq_poly_ptr t, const fmpq_poly_ptr a, const fmpq_poly_ptr b)
+{
+    fmpz_t lead, temp;
+    unsigned long bound, limbs;
+    fmpq_poly_t tempg, temps, tempt;
+    
+    /* Assertion! */
+    if (rop == s | rop == t | s == t)
+    {
+        printf("ERROR (fmpq_poly_xgcd).  Output arguments aliased.\n");
+        abort();
+    }
+    
+    /* Deal with aliasing                                                    */
+    if (rop == a | rop == b)
+    {
+        if (s == a | s == b)
+        {
+            /* We know that t does not coincide with a or b, since otherwise */
+            /* two of rop, s, and t coincide, too.                           */
+            fmpq_poly_init(tempg);
+            fmpq_poly_init(temps);
+            fmpq_poly_xgcd(tempg, temps, t, a, b);
+            fmpq_poly_swap(rop, tempg);
+            fmpq_poly_swap(s, temps);
+            fmpq_poly_clear(tempg);
+            fmpq_poly_clear(temps);
+            return;
+        }
+        else
+        {
+            if (t == a | t == b)
+            {
+                fmpq_poly_init(tempg);
+                fmpq_poly_init(tempt);
+                fmpq_poly_xgcd(tempg, s, tempt, a, b);
+                fmpq_poly_swap(rop, tempg);
+                fmpq_poly_swap(t, tempt);
+                fmpq_poly_clear(tempg);
+                fmpq_poly_clear(tempt);
+                return;
+            }
+            else
+            {
+                fmpq_poly_init(tempg);
+                fmpq_poly_xgcd(tempg, s, t, a, b);
+                fmpq_poly_swap(rop, tempg);
+                fmpq_poly_clear(tempg);
+                return;
+            }
+        }
+    }
+    else
+    {
+        if (s == a | s == b)
+        {
+            if (t == a | t == b)
+            {
+                fmpq_poly_init(temps);
+                fmpq_poly_init(tempt);
+                fmpq_poly_xgcd(rop, temps, tempt, a, b);
+                fmpq_poly_swap(s, temps);
+                fmpq_poly_swap(t, tempt);
+                fmpq_poly_clear(temps);
+                fmpq_poly_clear(tempt);
+                return;
+            }
+            else
+            {
+                fmpq_poly_init(temps);
+                fmpq_poly_xgcd(rop, temps, t, a, b);
+                fmpq_poly_swap(s, temps);
+                fmpq_poly_clear(temps);
+                return;
+            }
+        }
+        else
+        {
+            if (t == a | t == b)
+            {
+                fmpq_poly_init(tempt);
+                fmpq_poly_xgcd(rop, s, tempt, a, b);
+                fmpq_poly_swap(t, tempt);
+                fmpq_poly_clear(tempt);
+                return;
+            }
+        }
+    }
+    
+    /* From here on, we may assume that none of the output variables are     */
+    /* aliases for the input variables.                                      */
+    
+    /* Deal with the following three corner cases:                           */
+    /*   a == 0, b == 0                                                      */
+    /*   a == 0, b =! 0                                                      */
+    /*   a != 0, b == 0                                                      */
+    if (fmpq_poly_is_zero(a))
+    {
+        if (fmpq_poly_is_zero(b))          /* Case 1.  a == b == 0           */
+        {
+            fmpq_poly_zero(rop);
+            fmpq_poly_zero(s);
+            fmpq_poly_zero(t);
+            return;
+        }
+        else                               /* Case 2.  a == 0, b != 0        */
+        {
+            fmpq_poly_monic(rop, b);
+            fmpq_poly_zero(s);
+            
+            lead = fmpz_poly_lead(b->num);
+            fmpz_poly_zero(t->num);
+            if (_fmpq_poly_den_is_one(b))
+            {
+                fmpz_poly_set_coeff_si(t->num, 0, 1);
+                fmpq_poly_set_den(t, lead);
+            }
+            else
+            {
+                temp = fmpz_init(FLINT_MAX(fmpz_size(b->den), fmpz_size(lead)));
+                fmpz_gcd(temp, b->den, lead);
+                fmpz_poly_set_coeff_fmpz(t->num, 0, b->den);
+                fmpz_poly_scalar_div_fmpz(t->num, t->num, temp);
+                _fmpq_poly_den_fit_limbs(t, fmpz_size(lead));
+                fmpz_tdiv(t->den, lead, temp);
+                fmpz_clear(temp);
+            }
+            if (t->den != NULL && fmpz_sgn(t->den) < 0)
+            {
+                fmpz_poly_neg(t->num, t->num);
+                fmpz_neg(t->den, t->den);
+            }
+            return;
+        }
+    }
+    else
+    {
+        if (fmpq_poly_is_zero(b))          /* Case 3.  a != 0, b == 0        */
+        {
+            fmpq_poly_xgcd(rop, t, s, b, a);
+            return;
+        }
+    }
+    
+    /* We are now in the general case where a and b are non-zero.            */
+    
+    _fmpq_poly_den_fit_limbs(s, a->num->limbs);
+    _fmpq_poly_den_fit_limbs(t, b->num->limbs);
+    fmpz_poly_content(s->den, a->num);
+    fmpz_poly_content(t->den, b->num);
+    fmpz_poly_scalar_div_fmpz(s->num, a->num, s->den);
+    fmpz_poly_scalar_div_fmpz(t->num, b->num, t->den);
+    
+    /* Note that, since s->num and t->num are primitive, rop->num is         */
+    /* primitive, too.  In fact, it is the rational greatest common divisor  */
+    /* of a and b.  As of FLINT 1.4.0, the leading coefficient of res.num    */
+    /* *might* be negative.                                                  */
+    
+    fmpz_poly_gcd(rop->num, s->num, t->num);
+    if (fmpz_sgn(fmpz_poly_lead(rop->num)) < 0)
+        fmpz_poly_neg(rop->num, rop->num);
+    lead = fmpz_poly_lead(rop->num);
+    
+    /* Now rop->num is a (primitive) rational greatest common divisor of     */
+    /* a and b.                                                              */
+    
+    if (fmpz_poly_degree(rop->num) > 0)
+    {
+        fmpz_poly_div(s->num, s->num, rop->num);
+        fmpz_poly_div(t->num, t->num, rop->num);
+    }
+    
+    bound = fmpz_poly_resultant_bound(s->num, t->num);
+    if (fmpz_is_one(lead))
+        _fmpq_poly_den_fit_limbs(rop, bound/FLINT_BITS + 2);
+    else
+        _fmpq_poly_den_fit_limbs(rop, FLINT_MAX(bound/FLINT_BITS + 2, 
+            fmpz_size(lead)));
+    fmpz_poly_xgcd(rop->den, s->num, t->num, s->num, t->num);
+    
+    /* Now the following equation holds:                                     */
+    /*   rop->den rop->num ==                                                */
+    /*             (s->num a->den / s->den) a +  (t->num b->den / t->den) b. */
+    
+    limbs = FLINT_MAX(s->num->limbs, t->num->limbs);
+    limbs = FLINT_MAX(limbs, fmpz_size(s->den));
+    limbs = FLINT_MAX(limbs, fmpz_size(t->den) + fmpz_size(rop->den) + fmpz_size(lead));
+    temp = fmpz_init(limbs);
+    
+    _fmpq_poly_den_fit_limbs(s, fmpz_size(s->den) + fmpz_size(rop->den) 
+        + fmpz_size(lead));
+    if (!_fmpq_poly_den_is_one(a))
+        fmpz_poly_scalar_mul_fmpz(s->num, s->num, a->den);
+    fmpz_mul(temp, s->den, rop->den);
+    fmpz_mul(s->den, temp, lead);
+    
+    _fmpq_poly_den_fit_limbs(t, fmpz_size(t->den) + fmpz_size(rop->den) 
+        + fmpz_size(lead));
+    if (!_fmpq_poly_den_is_one(b))
+        fmpz_poly_scalar_mul_fmpz(t->num, t->num, b->den);
+    fmpz_mul(temp, t->den, rop->den);
+    fmpz_mul(t->den, temp, lead);
+    
+    fmpq_poly_canonicalize(s, temp);
+    fmpq_poly_canonicalize(t, temp);
+    
+    fmpz_set(rop->den, lead);
+    
+    fmpz_clear(temp);
+}
+
+/**
+ * \ingroup  GCD
+ * 
+ * Computes the monic (or zero) least common multiple of \c a and \c b.
+ * 
+ * If either of \c a and \c b is zero, returns zero.  This behaviour ensures 
+ * that the relation 
+ * \f[
+ * \text{lcm}(a,b) \gcd(a,b) \sim a b
+ * \f]
+ * holds, where \f$\sim\f$ denotes equality up to units.
+ */
+void fmpq_poly_lcm(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b)
+{
+    fmpz_poly_t prod;
+    fmpz_t lead;
+    fmpq_poly_t tpoly;
+    
+    /* Handle aliasing */
+    if (rop == a | rop == b)
+    {
+        fmpq_poly_init(tpoly);
+        fmpq_poly_lcm(tpoly, a, b);
+        fmpq_poly_swap(rop, tpoly);
+        fmpq_poly_clear(tpoly);
+        return;
+    }
+    
+    if (fmpq_poly_is_zero(a))
+        fmpq_poly_zero(rop);
+    else if (fmpq_poly_is_zero(b))
+        fmpq_poly_zero(rop);
+    else
+    {
+        fmpz_poly_init(prod);
+        fmpq_poly_gcd(rop, a, b);
+        fmpz_poly_mul(prod, a->num, b->num);
+        fmpz_poly_primitive_part(prod, prod);
+        fmpz_poly_div(rop->num, prod, rop->num);
+        
+        /* Note that GCD returns a monic rop and so a primitive rop.num.     */
+        /* Dividing the primitive prod by this yields a primitive quotient   */
+        /* rop->num.                                                         */
+        
+        lead = fmpz_poly_lead(rop->num);
+        if (fmpz_is_one(lead))
+        {
+            if (rop->den != NULL)
+                fmpz_set_si(rop->den, 1);
+        }
+        else if (fmpz_is_m1(lead))
+        {
+            fmpz_poly_neg(rop->num, rop->num);
+            if (rop->den != NULL)
+                fmpz_set_si(rop->den, 1);
+        }
+        else if (fmpz_sgn(lead) < 0)
+        {
+            fmpq_poly_set_den(rop, lead);
+            fmpz_poly_neg(rop->num, rop->num);
+            fmpz_neg(rop->den, rop->den);
+        }
+        else
+            fmpq_poly_set_den(rop, lead);
+        fmpz_poly_clear(prod);
+    }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Derivative
+
+/**
+ * \ingroup  Derivative
+ * 
+ * Sets \c rop to the derivative of \c op.
+ * 
+ * \todo  The second argument should be declared \c const, but as of 
+ *        FLINT 1.5.0 this generates a compile-time warning.
+ */
+void fmpq_poly_derivative(fmpq_poly_ptr rop, fmpq_poly_ptr op)
+{
+    if (fmpq_poly_degree(op) < 1)
+        fmpq_poly_zero(rop);
+    else
+    {
+        fmpz_poly_derivative(rop->num, op->num);
+        fmpq_poly_set_den(rop, op->den);
+        fmpq_poly_canonicalize(rop, NULL);
+    }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Evaluation
+
+/**
+ * \ingroup  Evaluation
+ * 
+ * Evaluates the integer polynomial \c f at the rational \c a using Horner's 
+ * method.
+ */
+void _fmpz_poly_evaluate_mpq_horner(mpq_t rop, const fmpz_poly_t f, const mpq_t a)
+{
+    mpq_t temp;
+    mpq_t tempr;      /* Temporary variable to handle aliasing               */
+    unsigned long n;  /* Indexing variable                                   */
+    
+    /* Handle aliasing                                                       */
+    if (rop == a)
+    {
+        mpq_init(tempr);
+        _fmpz_poly_evaluate_mpq_horner(tempr, f, a);
+        mpq_swap(rop, tempr);
+        mpq_clear(tempr);
+        return;
+    }
+    
+    n = fmpz_poly_length(f);
+    
+    if (n == 0)
+    {
+        mpq_set_si(rop, 0, 1);
+    }
+    else if (n == 1)
+    {
+        fmpz_poly_get_coeff_mpz(mpq_numref(rop), f, 0);
+        mpz_set_si(mpq_denref(rop), 1);
+    }
+    else
+    {
+        n--;
+        fmpz_poly_get_coeff_mpz(mpq_numref(rop), f, n);
+        mpz_set_si(mpq_denref(rop), 1);
+        mpq_init(temp);
+        do {
+            n--;
+            mpq_mul(temp, rop, a);
+            fmpz_poly_get_coeff_mpz(mpq_numref(rop), f, n);
+            mpz_set_si(mpq_denref(rop), 1);
+            mpq_add(rop, rop, temp);
+        } while (n);
+        mpq_clear(temp);
+    }
+}
+
+/**
+ * \ingroup  Evaluation
+ * 
+ * Evaluates the rational polynomial \c f at the integer \c a.
+ * 
+ * Assumes that the numerator and denominator of the <tt>mpq_t</tt>
+ * \c rop are distinct (as objects in memory) from the <tt>mpz_t</tt> 
+ * \c a.
+ */
+void fmpq_poly_evaluate_mpz(mpq_t rop, fmpq_poly_ptr f, const mpz_t a)
+{
+    fmpz_t num, t;
+    unsigned long limbs, max;
+    
+    if (fmpq_poly_is_zero(f))
+    {
+        mpq_set_si(rop, 0, 1);
+        return;
+    }
+    
+    /* Establish a bound on the size of f->num evaluated at a                */
+    max = (f->num->length) * (f->num->limbs + f->num->length * mpz_size(a));
+    
+    /* Compute the result                                                    */
+    num = fmpz_init(max);
+    t   = fmpz_init(mpz_size(a));
+    mpz_to_fmpz(t, a);
+    fmpz_poly_evaluate(num, f->num, t);
+    fmpz_to_mpz(mpq_numref(rop), num);
+    if (f->den == NULL)
+    {
+        mpz_set_si(mpq_denref(rop), 1);
+    }
+    else
+    {
+        fmpz_to_mpz(mpq_denref(rop), f->den);
+        mpq_canonicalize(rop);
+    }
+    
+    /* Clean-up                                                              */
+    fmpz_clear(num);
+    fmpz_clear(t);
+}
+
+/**
+ * \ingroup  Evaluation
+ * 
+ * Evaluates the rational polynomial \c f at the rational \c a.
+ */
+void fmpq_poly_evaluate_mpq(mpq_t rop, fmpq_poly_ptr f, const mpq_t a)
+{
+    mpq_t tempr;
+    mpz_t den;
+    
+    if (rop == a)
+    {
+        mpq_init(tempr);
+        fmpq_poly_evaluate_mpq(tempr, f, a);
+        mpq_swap(rop, tempr);
+        mpq_clear(tempr);
+        return;
+    }
+    
+    _fmpz_poly_evaluate_mpq_horner(rop, f->num, a);
+    if (!_fmpq_poly_den_is_one(f))
+    {
+        mpz_init(den);
+        fmpz_to_mpz(den, f->den);
+        mpz_mul(mpq_denref(rop), mpq_denref(rop), den);
+        mpq_canonicalize(rop);
+        mpz_clear(den);
+    }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Gaussian content
+
+/**
+ * \ingroup  Content
+ * 
+ * Returns the non-negative content of \c op.
+ * 
+ * The content of \f$0\f$ is defined to be \f$0\f$.
+ */
+void fmpq_poly_content(mpq_t rop, const fmpq_poly_ptr op)
+{
+    fmpz_t numc;
+    
+    numc = fmpz_init(op->num->limbs);
+    fmpz_poly_content(numc, op->num);
+    fmpz_abs(numc, numc);
+    fmpz_to_mpz(mpq_numref(rop), numc);
+    if (op->den == NULL)
+        mpz_set_si(mpq_denref(rop), 1);
+    else
+        fmpz_to_mpz(mpq_denref(rop), op->den);
+    fmpz_clear(numc);
+}
+
+/**
+ * \ingroup  Content
+ * 
+ * Returns the primitive part (with non-negative leading coefficient) of 
+ * \c op as an element of type #fmpq_poly_t.
+ */
+void fmpq_poly_primitive_part(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
+{
+    if (fmpq_poly_is_zero(op))
+        fmpq_poly_zero(rop);
+    else
+    {
+        fmpz_poly_primitive_part(rop->num, op->num);
+        if (fmpz_sgn(fmpz_poly_lead(rop->num)) < 0)
+            fmpz_poly_neg(rop->num, rop->num);
+        if (rop->den != NULL)
+            fmpz_set_si(rop->den, 1);
+    }
+}
+
+/**
+ * \brief    Returns whether \c op is monic.
+ * \ingroup  Content
+ * 
+ * Returns whether \c op is monic.
+ * 
+ * By definition, the zero polynomial is \e not monic.
+ */
+int fmpq_poly_is_monic(const fmpq_poly_ptr op)
+{
+    fmpz_t lead;
+    
+    if (fmpq_poly_is_zero(op))
+        return 0;
+    else
+    {
+        lead = fmpz_poly_lead(op->num);
+        if (_fmpq_poly_den_is_one(op))
+            return fmpz_is_one(lead);
+        else
+            return (fmpz_sgn(lead) > 0) && (fmpz_cmpabs(lead, op->den) == 0);
+    }
+}
+
+/**
+ * Sets \c rop to the unique monic scalar multiple of \c op.
+ * 
+ * As the only special case, if \c op is the zero polynomial, \c rop is set 
+ * to zero, too.
+ */
+void fmpq_poly_monic(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
+{
+    fmpz_t lead;
+    
+    if (fmpq_poly_is_zero(op))
+    {
+        fmpq_poly_zero(rop);
+        return;
+    }
+    
+    fmpz_poly_primitive_part(rop->num, op->num);
+    lead = fmpz_poly_lead(rop->num);
+    if (fmpz_sgn(lead) < 0)
+        fmpz_poly_neg(rop->num, rop->num);
+    fmpq_poly_set_den(rop, lead);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Resultant
+
+/**
+ * \brief    Returns the resultant of \c a and \c b.
+ * \ingroup  Resultant
+ * 
+ * Returns the resultant of \c a and \c b.
+ * 
+ * Enumerating the roots of \c a and \c b over \f$\bar{\mathbf{Q}}\f$ as 
+ * \f$r_1, \dotsc, r_m\f$ and \f$s_1, \dotsc, s_n\f$, respectively, and 
+ * letting \f$x\f$ and \f$y\f$ denote the leading coefficients, the resultant 
+ * is defined as 
+ * \f[
+ * x^{\deg(b)} y^{\deg(a)} \prod_{1 \leq i, j \leq n} (r_i - s_j).
+ * \f]
+ * 
+ * We handle special cases as follows:  if one of the polynomials is zero, 
+ * the resultant is zero.  Note that otherwise if one of the polynomials is 
+ * constant, the last term in the above expression is the empty product.
+ */
+void fmpq_poly_resultant(mpq_t rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b)
+{
+    fmpz_t rest, t1, t2;
+    fmpz_poly_t anum, bnum;
+    fmpz_t acont, bcont;
+    long d1, d2;
+    unsigned long bound, denbound, numbound;
+    fmpz_poly_t g;
+    
+    d1 = fmpz_poly_degree(a->num);
+    d2 = fmpz_poly_degree(b->num);
+    
+    /* We first handle special cases.                                        */
+    if (d1 < 0 | d2 < 0)
+    {
+        mpq_set_si(rop, 0, 1);
+        return;
+    }
+    if (d1 == 0)
+    {
+        if (d2 == 0)
+            mpq_set_si(rop, 1, 1);
+        else if (d2 == 1)
+        {
+            fmpz_to_mpz(mpq_numref(rop), fmpz_poly_lead(a->num));
+            if (a->den == NULL)
+                mpz_set_si(mpq_denref(rop), 1);
+            else
+                fmpz_to_mpz(mpq_denref(rop), a->den);
+        }
+        else
+        {
+            if (a->den == NULL)
+                bound = a->num->limbs;
+            else
+                bound = FLINT_MAX(a->num->limbs, fmpz_size(a->den));
+            t1 = fmpz_init(d2 * bound);
+            fmpz_pow_ui(t1, fmpz_poly_lead(a->num), d2);
+            fmpz_to_mpz(mpq_numref(rop), t1);
+            if (a->den == NULL)
+                mpz_set_si(mpq_denref(rop), 1);
+            else
+            {
+                fmpz_pow_ui(t1, a->den, d2);
+                fmpz_to_mpz(mpq_denref(rop), t1);
+            }
+            fmpz_clear(t1);
+        }
+        return;
+    }
+    if (d2 == 0)
+    {
+        fmpq_poly_resultant(rop, b, a);
+        return;
+    }
+    
+    /* We are now in the general case, with both polynomials of degree at    */
+    /* least 1.                                                              */
+    
+    /* We set a->num =: acont anum with acont > 0 and anum primitive.        */
+    acont = fmpz_init(a->num->limbs);
+    fmpz_poly_content(acont, a->num);
+    fmpz_abs(acont, acont);
+    fmpz_poly_init(anum);
+    fmpz_poly_scalar_div_fmpz(anum, a->num, acont);
+    
+    /* We set b->num =: bcont bnum with bcont > 0 and bnum primitive.        */
+    bcont = fmpz_init(b->num->limbs);
+    fmpz_poly_content(bcont, b->num);
+    fmpz_abs(bcont, bcont);
+    fmpz_poly_init(bnum);
+    fmpz_poly_scalar_div_fmpz(bnum, b->num, bcont);
+    
+    fmpz_poly_init(g);
+    fmpz_poly_gcd(g, anum, bnum);
+    
+    /* If the gcd has positive degree, the resultant is zero.                */
+    if (fmpz_poly_degree(g) > 0)
+    {
+        mpq_set_si(rop, 0, 1);
+        
+        /* Clean up                                                          */
+        fmpz_clear(acont);
+        fmpz_clear(bcont);
+        fmpz_poly_clear(anum);
+        fmpz_poly_clear(bnum);
+        fmpz_poly_clear(g);
+        return;
+    }
+    
+    /* Set some bounds                                                       */
+    if (a->den == NULL)
+    {
+        if (b->den == NULL)
+        {
+            numbound = FLINT_MAX(d2 * fmpz_size(acont), d1 * fmpz_size(bcont));
+            denbound = 1;
+        }
+        else
+        {
+            numbound = FLINT_MAX(d2 * fmpz_size(acont), d1 * fmpz_size(bcont));
+            denbound = d1 * fmpz_size(b->den);
+        }
+    }
+    else
+    {
+        if (b->den == NULL)
+        {
+            numbound = FLINT_MAX(d2 * fmpz_size(acont), d1 * fmpz_size(bcont));
+            denbound = d2 * fmpz_size(a->den);
+        }
+        else
+        {
+            numbound = FLINT_MAX(d2 * fmpz_size(acont), d1 * fmpz_size(bcont));
+            denbound = FLINT_MAX(d2 * fmpz_size(a->den), d1 * fmpz_size(b->den));
+        }
+    }
+    
+    /* Now anum and bnum are coprime and we compute their resultant using    */
+    /* the method from the fmpz_poly module.                                 */
+    bound = fmpz_poly_resultant_bound(anum, bnum);
+    bound = bound/FLINT_BITS + 2 + d2 * fmpz_size(acont) 
+        + d1 * fmpz_size(bcont);
+    rest = fmpz_init(FLINT_MAX(bound, denbound));
+    fmpz_poly_resultant(rest, anum, bnum);
+    
+    /* Finally, w take into account the factors acont/a.den and bcont/b.den. */
+    t1 = fmpz_init(FLINT_MAX(bound, denbound));
+    t2 = fmpz_init(FLINT_MAX(numbound, denbound));
+    
+    if (!fmpz_is_one(acont))
+    {
+        fmpz_pow_ui(t2, acont, d2);
+        fmpz_set(t1, rest);
+        fmpz_mul(rest, t1, t2);
+    }
+    if (!fmpz_is_one(bcont))
+    {
+        fmpz_pow_ui(t2, bcont, d1);
+        fmpz_set(t1, rest);
+        fmpz_mul(rest, t1, t2);
+    }
+    
+    fmpz_to_mpz(mpq_numref(rop), rest);
+    
+    if (_fmpq_poly_den_is_one(a))
+    {
+        if (_fmpq_poly_den_is_one(b))
+            fmpz_set_si(rest, 1);
+        else
+            fmpz_pow_ui(rest, b->den, d1);
+    }
+    else
+    {
+        if (_fmpq_poly_den_is_one(b))
+            fmpz_pow_ui(rest, a->den, d2);
+        else
+        {
+            fmpz_pow_ui(t1, a->den, d2);
+            fmpz_pow_ui(t2, b->den, d1);
+            fmpz_mul(rest, t1, t2);
+        }
+    }
+    
+    fmpz_to_mpz(mpq_denref(rop), rest);
+    mpq_canonicalize(rop);
+    
+    /* Clean up                                                              */
+    fmpz_clear(rest);
+    fmpz_clear(t1);
+    fmpz_clear(t2);
+    fmpz_poly_clear(anum);
+    fmpz_poly_clear(bnum);
+    fmpz_clear(acont);
+    fmpz_clear(bcont);
+    fmpz_poly_clear(g);
+}
+
+/**
+ * \brief    Computes the discriminant of \c a.
+ * \ingroup  Discriminant
+ * 
+ * Computes the discriminant of the polynomial \f$a\f$.
+ * 
+ * The discriminant \f$R_n\f$ is defined as 
+ * \f[
+ * R_n = a_n^{2 n-2} \prod_{1 \le i < j \le n} (r_i - r_j)^2,
+ * \f]
+ * where \f$n\f$ is the degree of this polynomial, \f$a_n\f$ is the leading 
+ * coefficient and \f$r_1, \ldots, r_n\f$ are the roots over 
+ * \f$\bar{\mathbf{Q}}\f$ are.
+ * 
+ * The discriminant of constant polynomials is defined to be \f$0\f$.
+ * 
+ * This implementation uses the identity 
+ * \f[
+ * R_n(f) := (-1)^(n (n-1)/2) R(f,f') / a_n,
+ * \f]
+ * where \f$n\f$ is the degree of this polynomial, \f$a_n\f$ is the leading 
+ * coefficient and \f$f'\f$ is the derivative of \f$f\f$.
+ * 
+ * \see  #fmpq_poly_resultant()
+ */
+void fmpq_poly_discriminant(mpq_t disc, fmpq_poly_t a)
+{
+    fmpq_poly_t der;
+    mpq_t t;
+    long n;
+    
+    n = fmpq_poly_degree(a);
+    if (n < 1)
+    {
+        mpq_set_si(disc, 0, 1);
+        return;
+    }
+    
+    fmpq_poly_init(der);
+    fmpq_poly_derivative(der, a);
+    fmpq_poly_resultant(disc, a, der);
+    mpq_init(t);
+    mpq_set_si(t, 1, 1);
+    
+    if (a->den != NULL)
+        fmpz_to_mpz(mpq_numref(t), a->den);
+    if (!fmpz_is_one(fmpz_poly_lead(a->num)))
+        fmpz_to_mpz(mpq_denref(t), fmpz_poly_lead(a->num));
+    mpq_canonicalize(t);
+    mpq_mul(disc, disc, t);
+    
+    if (n % 4 > 1)
+        mpz_neg(mpq_numref(disc), mpq_numref(disc));
+    
+    fmpq_poly_clear(der);
+    mpq_clear(t);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Composition
+
+/**
+ * \brief    Returns the composition of \c a and \c b.
+ * \ingroup  Composition
+ * 
+ * Returns the composition of \c a and \c b.
+ * 
+ * To be clear about the order of the composition, denoting the polynomials 
+ * \f$a(T)\f$ and \f$b(T)\f$, returns the polynomial \f$a(b(T))\f$.
+ */
+void fmpq_poly_compose(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b)
+{
+    mpq_t x;            /* Rational to hold the inverse of b->den            */
+    fmpq_poly_t tempr;  /* Temporary variable to handle aliasing             */
+    
+    if (_fmpq_poly_den_is_one(b))
+    {
+        fmpz_poly_compose(rop->num, a->num, b->num);
+        fmpq_poly_set_den(rop, a->den);
+        fmpq_poly_canonicalize(rop, NULL);
+        return;
+    }
+    
+    /* Aliasing.                                                             */
+    /*                                                                       */
+    /* Note that rop and a, as well as a and b may be aliased, but rop and   */
+    /* b may not be aliased.                                                 */
+    if (rop == b)
+    {
+        fmpq_poly_init(tempr);
+        fmpq_poly_compose(tempr, a, b);
+        fmpq_poly_swap(rop, tempr);
+        fmpq_poly_clear(tempr);
+        return;
+    }
+    
+    /* Set x = 1/b.den, and note this is already in canonical form.          */
+    mpq_init(x);
+    mpz_set_si(mpq_numref(x), 1);
+    fmpz_to_mpz(mpq_denref(x), b->den);
+    
+    /* First set rop = a(T / b.den) and then use FLINT's composition to      */
+    /* set rop->num = rop->num(b->num).                                      */
+    fmpq_poly_rescale(rop, a, x);
+    fmpz_poly_compose(rop->num, rop->num, b->num);
+    fmpq_poly_canonicalize(rop, NULL);
+    
+    mpq_clear(x);
+}
+
+/**
+ * \brief    Rescales the co-ordinate.
+ * \ingroup  Composition
+ * 
+ * Denoting this polynomial \f$a(T)\f$, computes the polynomial \f$a(x T)\f$.
+ */
+void fmpq_poly_rescale(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x)
+{
+    unsigned long numsize, densize;
+    unsigned long limbs;
+    fmpz_t num, den;
+    fmpz_t coeff, power, t;
+    long i, n;
+    
+    fmpq_poly_set(rop, op);
+    
+    if (fmpq_poly_degree(rop) < 1)
+        return;
+    
+    num = fmpz_init(mpz_size(mpq_numref(x)));
+    den = fmpz_init(mpz_size(mpq_denref(x)));
+    mpz_to_fmpz(num, mpq_numref(x));
+    mpz_to_fmpz(den, mpq_denref(x));
+    numsize = fmpz_size(num);
+    densize = fmpz_size(den);
+    
+    n = fmpz_poly_degree(rop->num);
+    
+    if (fmpz_is_one(den))
+    {
+        coeff = fmpz_init(rop->num->limbs + n * numsize);
+        power = fmpz_init(n * numsize);
+        t = fmpz_init(rop->num->limbs + n * numsize);
+        
+        fmpz_set(power, num);
+        
+        fmpz_poly_get_coeff_fmpz(t, rop->num, 1);
+        fmpz_mul(coeff, t, power);
+        fmpz_poly_set_coeff_fmpz(rop->num, 1, coeff);
+        
+        for (i = 2; i <= n; i++)
+        {
+            fmpz_set(t, power);
+            fmpz_mul(power, t, num);
+            fmpz_poly_get_coeff_fmpz(t, rop->num, i);
+            fmpz_mul(coeff, t, power);
+            fmpz_poly_set_coeff_fmpz(rop->num, i, coeff);
+        }
+    }
+    else
+    {
+        coeff = fmpz_init(rop->num->limbs + n * (numsize + densize));
+        power = fmpz_init(n * (numsize + densize));
+        limbs = rop->num->limbs + n * (numsize + densize);
+        if (rop->den != NULL)
+            limbs = FLINT_MAX(limbs, fmpz_size(rop->den));
+        t = fmpz_init(limbs);
+        
+        fmpz_pow_ui(power, den, n);
+        
+        if (rop->den == NULL)
+        {
+            rop->den = fmpz_init(n * densize);
+            fmpz_set(rop->den, power);
+        }
+        else
+        {
+            fmpz_set(t, rop->den);
+            limbs = n * densize + fmpz_size(rop->den);
+            if (fmpz_size(rop->den) < limbs)
+                rop->den = fmpz_realloc(rop->den, limbs);
+            fmpz_mul(rop->den, power, t);
+        }
+        
+        fmpz_set_si(power, 1);
+        for (i = n-1; i >= 0; i--)
+        {
+            fmpz_set(t, power);
+            fmpz_mul(power, t, den);
+            fmpz_poly_get_coeff_fmpz(t, rop->num, i);
+            fmpz_mul(coeff, t, power);
+            fmpz_poly_set_coeff_fmpz(rop->num, i, coeff);
+        }
+        
+        fmpz_set_si(power, 1);
+        for (i = 1; i <= n; i++)
+        {
+            fmpz_set(t, power);
+            fmpz_mul(power, t, num);
+            fmpz_poly_get_coeff_fmpz(t, rop->num, i);
+            fmpz_mul(coeff, t, power);
+            fmpz_poly_set_coeff_fmpz(rop->num, i, coeff);
+        }
+    }
+    fmpq_poly_canonicalize(rop, NULL);
+    fmpz_clear(num);
+    fmpz_clear(den);
+    fmpz_clear(coeff);
+    fmpz_clear(power);
+    fmpz_clear(t);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Square-free
+
+/**
+ * \brief    Returns whether \c op is squarefree.
+ * \ingroup  Squarefree
+ * 
+ * Returns whether \c op is squarefree.
+ * 
+ * By definition, a polynomial is square-free if it is not a multiple of a 
+ * non-unit factor.  In particular, polynomials up to degree 1 (including) 
+ * are square-free.
+ */
+int fmpq_poly_is_squarefree(const fmpq_poly_ptr op)
+{
+    fmpz_poly_t prim;
+    int ans;
+    if (fmpq_poly_degree(op) < 2)
+        return 1;
+    else
+    {
+        fmpz_poly_init(prim);
+        fmpz_poly_primitive_part(prim, op->num);
+        ans = fmpz_poly_is_squarefree(prim);
+        fmpz_poly_clear(prim);
+        return ans;
+    }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Subpolynomials
+
+/**
+ * \brief    Returns a contiguous subpolynomial.
+ * \ingroup  Subpolynomials
+ * 
+ * Returns the slice with coefficients from \f$x^i\f$ (including) to 
+ * \f$x^j\f$ (excluding).
+ */
+void fmpq_poly_getslice(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long i, long j)
+{
+    long k;
+    
+    i = FLINT_MAX(0, i);
+    j = FLINT_MIN(fmpq_poly_degree(op) + 1, j);
+    
+    /* Aliasing */
+    if (rop == op)
+    {
+        for (k = 0; k < i; k++)
+            fmpz_poly_set_coeff_si(rop->num, k, 0);
+        for (k = j; k <= fmpz_poly_degree(rop->num); k++)
+            fmpz_poly_set_coeff_si(rop->num, k, 0);
+        fmpq_poly_canonicalize(rop, NULL);
+        return;
+    }
+    
+    fmpz_poly_zero(rop->num);
+    if (i < j)
+    {
+        for (k = i; k < j; k++)
+            fmpz_poly_set_coeff_fmpz(rop->num, k, 
+                fmpz_poly_get_coeff_ptr(op->num, k));
+        fmpq_poly_set_den(rop, op->den);
+        fmpq_poly_canonicalize(rop, NULL);
+    }
+    else
+    {
+        if (rop->den != NULL)
+            fmpz_set_si(rop->den, 1);
+    }
+}
+
+/**
+ * \brief    Shifts this polynomial.
+ * \ingroup  Subpolynomials
+ * 
+ * Multiplies \c op by \f$x^n\f$.  If \f$n\f$ is negative, terms below 
+ * \f$x^n\f$ are simply discarded.
+ */
+void fmpq_poly_shift(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long n)
+{
+    /* XXX:  As a workaround for a bug in FLINT 1.5.1, we need to handle  */
+    /* the zero polynomial separately.                                    */
+    if (fmpq_poly_is_zero(op))
+    {
+        fmpq_poly_zero(rop);
+        return;
+    }
+    
+    if (n > 0)
+    {
+        fmpz_poly_left_shift(rop->num, op->num, n);
+        fmpq_poly_set_den(rop, op->den);
+    }
+    else if (n < 0)
+    {
+        fmpz_poly_right_shift(rop->num, op->num, -n);
+        fmpq_poly_set_den(rop, op->den);
+        fmpq_poly_canonicalize(rop, NULL);
+    }
+    else
+    {
+        fmpq_poly_set(rop, op);
+    }
+}
+
+/**
+ * \brief    Truncates this polynomials.
+ * \ingroup  Subpolynomials
+ * 
+ * Truncates <tt>op</tt>  modulo \f$x^n\f$ whenever \f$n\f$ is positive.  
+ * Returns zero otherwise.
+ */
+void fmpq_poly_truncate(fmpq_poly_ptr rop, const fmpq_poly_ptr op, unsigned long n)
+{
+    fmpq_poly_set(rop, op);
+    if (fmpq_poly_degree(rop) >= n)
+    {
+        fmpz_poly_truncate(rop->num, n);
+        fmpq_poly_canonicalize(rop, NULL);
+    }
+}
+
+/**
+ * \brief    Reverses this polynomial.
+ * \ingroup  Subpolynomials
+ * 
+ * Reverses the coefficients of \c op and places the result in <tt>rop</tt>.
+ */
+void fmpq_poly_reverse(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
+{
+    if (fmpz_poly_degree(op->num) < 1)
+        fmpq_poly_set(rop, op);
+    else
+    {
+        fmpz_poly_reverse(rop->num, op->num, fmpz_poly_length(op->num));
+        fmpq_poly_set_den(rop, op->den);
+    }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// String conversion
+
+/**
+ * \addtogroup StringConversions
+ * 
+ * The following three methods enable users to construct elements of type 
+ * \c fmpq_poly_ptr from strings or to obtain string representations of such 
+ * elements.
+ *
+ * The format used is based on the FLINT format for integer polynomials of 
+ * type <tt>fmpz_poly_t</tt>, which we recall first: 
+ * 
+ * A non-zero polynomial \f$a_0 + a_1 X + \dotsb + a_n X^n\f$ of length 
+ * \f$n + 1\f$ is represented by the string <tt>n+1  a_0 a_1 ... a_n</tt>, 
+ * where there are two space characters following the length and single space 
+ * characters separating the individual coefficients.  There is no leading or 
+ * trailing white-space.  In contrast, the zero polynomial is represented by 
+ * <tt>0</tt>.
+ * 
+ * We adapt this notation for rational polynomials by using the <tt>mpq_t</tt> 
+ * notation for the coefficients without any additional white-space.
+ * 
+ * There is also a <tt>_pretty</tt> variant available.
+ * 
+ * Note that currently these functions are not optimized for performance and 
+ * are intended to be used only for debugging purposes or one-off input and 
+ * output, rather than as a low-level parser.
+ */
+
+/**
+ * \ingroup  StringConversions
+ * 
+ * Sets the rational polynomial \c rop to the value specified by the 
+ * null-terminated string \c str.
+ * 
+ * The behaviour is undefined if the format of the string \c str does not 
+ * conform to the specification.
+ * 
+ * \todo  In a future version, it would be nice to have this function return 
+ *        <tt>0</tt> or <tt>1</tt> depending on whether the input format 
+ *        was correct.  Currently, the method always returns <tt>1</tt>.
+ */
+int fmpq_poly_from_string(fmpq_poly_ptr rop, const char * str)
+{
+    mpq_t * a;
+    mpz_t den, t;
+    char * strcopy;
+    unsigned long strcopy_len;
+    unsigned long i, j, k, n;
+    
+    n = atoi(str);
+    
+    /* Handle the case of the zero polynomial                                */
+    if (n == 0)
+    {
+        fmpq_poly_zero(rop);
+        return 1;
+    }
+    
+    /* Compute the maximum length that the copy buffer has to be             */
+    strcopy_len = 0;
+    for (j = 0; str[j] != ' '; j++);
+    j += 2;
+    for (i = 0; i < n; i++)
+    {
+        for (k = j; !(str[k] == ' ' || str[k] == '\0'); k++);
+        strcopy_len = FLINT_MAX(strcopy_len, k - j + 1);
+        j = k + 1;
+    }
+    strcopy = (char *) calloc(strcopy_len, sizeof(char));
+    
+    /* Read the data into the array a of mpq_t's                             */
+    mpz_init_set_si(den, 1);
+    a = (mpq_t *) calloc(n, sizeof(mpq_t));
+    for (j = 0; str[j] != ' '; j++);
+    j += 2;
+    for (i = 0; i < n; i++)
+    {
+        for (k = j; !(str[k] == ' ' || str[k] == '\0'); k++);
+        memcpy(strcopy, str + j, k - j + 1);
+        strcopy[k - j] = '\0';
+        
+        mpq_init(a[i]);
+        mpq_set_str(a[i], strcopy, 10);
+        mpq_canonicalize(a[i]);
+        j = k + 1;
+        
+        mpz_lcm(den, den, mpq_denref(a[i]));
+    }
+    
+    /* Re-adjust the coefficients to the common denominator                  */
+    mpz_init(t);
+    for (i = 0; i < n; i++)
+    {
+        mpz_divexact(t, den, mpq_denref(a[i]));
+        mpz_mul(mpq_numref(a[i]), mpq_numref(a[i]), t);
+    }
+    
+    /* Transfer the data to the fmpq_poly_t rop                              */
+    fmpz_poly_realloc(rop->num, n);
+    i = n;
+    do
+    {
+        i--;
+        fmpz_poly_set_coeff_mpz(rop->num, i, mpq_numref(a[i]));
+    } while (i != 0);
+    if (mpz_cmp_si(den, 1) == 0)
+    {
+        if (rop->den != NULL)
+        {
+            fmpz_clear(rop->den);
+            rop->den = NULL;
+        }
+    }
+    else
+    {
+        _fmpq_poly_den_fit_limbs(rop, mpz_size(den));
+        mpz_to_fmpz(rop->den, den);
+    }
+    
+    /* Clean-up                                                              */
+    free(strcopy);
+    for (i = 0; i < n; i++)
+        mpq_clear(a[i]);
+    free(a);
+    mpz_clear(den);
+    mpz_clear(t);
+    
+    return 1;
+}
+
+/**
+ * \ingroup  StringConversions
+ * 
+ * Returns the string representation of the rational polynomial \c op.
+ */
+char * fmpq_poly_to_string(const fmpq_poly_ptr op)
+{
+    unsigned long i, j;
+    unsigned long len;     /* Upper bound on the length                      */
+    unsigned long denlen;  /* Size of the denominator in base 10             */
+    mpz_t z;
+    mpq_t q;
+    
+    char * str;
+    
+    if (fmpq_poly_is_zero(op))
+    {
+        str = calloc(2, sizeof(char));
+        str[0] = '0';
+        str[1] = '\0';
+        return str;
+    }
+    
+    mpz_init(z);
+    if (_fmpq_poly_den_is_one(op))
+    {
+        denlen = 0;
+    }
+    else
+    {
+        fmpz_to_mpz(z, op->den);
+        denlen = mpz_sizeinbase(z, 10);
+    }
+    len = fmpq_poly_places(fmpq_poly_length(op)) + 2;
+    for (i = 0; i < fmpq_poly_length(op); i++)
+    {
+        fmpz_poly_get_coeff_mpz(z, op->num, i);
+        len += mpz_sizeinbase(z, 10) + 1;
+        if (mpz_sgn(z) != 0)
+            len += 2 + denlen;
+    }
+    
+    mpq_init(q);
+    str = calloc(len, sizeof(char));
+    sprintf(str, "%lu", fmpq_poly_length(op));
+    for (j = 0; str[j] != '\0'; j++);
+    str[j++] = ' ';
+    for (i = 0; i < fmpq_poly_length(op); i++)
+    {
+        str[j++] = ' ';
+        fmpz_poly_get_coeff_mpz(mpq_numref(q), op->num, i);
+        if (op->den == NULL)
+            mpz_set_si(mpq_denref(q), 1);
+        else
+            fmpz_to_mpz(mpq_denref(q), op->den);
+        mpq_canonicalize(q);
+        mpq_get_str(str + j, 10, q);
+        for ( ; str[j] != '\0'; j++);
+    }
+    
+    mpq_clear(q);
+    mpz_clear(z);
+    
+    return str;
+}
+
+/**
+ * \ingroup  StringConversions
+ * 
+ * Returns the pretty string representation of \c op.
+ * 
+ * Returns the pretty string representation of the rational polynomial \c op, 
+ * using the string \c var as the variable name.
+ */
+char * fmpq_poly_to_string_pretty(const fmpq_poly_ptr op, const char * var)
+{
+    long i;
+    unsigned long j;
+    unsigned long len;     /* Upper bound on the length                      */
+    unsigned long denlen;  /* Size of the denominator in base 10             */
+    unsigned long varlen;  /* Length of the variable name                    */
+    mpz_t z;               /* op->den (if this is not 1)                     */
+    mpq_t q;
+    char * str;
+    
+    if (fmpq_poly_is_zero(op))     /* Zero polynomial                        */
+    {
+        str = calloc(2, sizeof(char));
+        str[0] = '0';
+        str[1] = '\0';
+        return str;
+    }
+    if (fmpq_poly_degree(op) == 0) /* Constant polynomials                   */
+    {
+        mpq_init(q);
+        fmpz_to_mpz(mpq_numref(q), fmpz_poly_lead(op->num));
+        if (op->den == NULL)
+            mpz_set_si(mpq_denref(q), 1);
+        else
+        {
+            fmpz_to_mpz(mpq_denref(q), op->den);
+            mpq_canonicalize(q);
+        }
+        str = mpq_get_str(NULL, 10, q);
+        mpq_clear(q);
+        return str;
+    }
+    
+    varlen = strlen(var);
+    
+    /* Copy the denominator into an mpz_t                                    */
+    mpz_init(z);
+    if (_fmpq_poly_den_is_one(op))
+    {
+        denlen = 0;
+    }
+    else
+    {
+        fmpz_to_mpz(z, op->den);
+        denlen = mpz_sizeinbase(z, 10);
+    }
+    
+    /* Estimate the length                                                   */
+    len = 0;
+    for (i = 0; i < fmpq_poly_length(op); i++)
+    {
+        fmpz_poly_get_coeff_mpz(z, op->num, i);
+        len += mpz_sizeinbase(z, 10);            /* Numerator                */
+        if (mpz_sgn(z) != 0)
+            len += 1 + denlen;                   /* Denominator and /        */
+        len += 3;                                /* Operator and whitespace  */
+        len += 1 + varlen + 1;                   /* *, x and ^               */
+        len += fmpq_poly_places(i);              /* Exponent                 */
+    }
+    
+    mpq_init(q);
+    str = calloc(len, sizeof(char));
+    j = 0;
+    
+    /* Print the leading term                                                */
+    fmpz_to_mpz(mpq_numref(q), fmpz_poly_lead(op->num));
+    if (op->den == NULL)
+        mpz_set_si(mpq_denref(q), 1);
+    else
+        fmpz_to_mpz(mpq_denref(q), op->den);
+    mpq_canonicalize(q);
+    
+    if (mpq_cmp_si(q, 1, 1) != 0)
+    {
+        if (mpq_cmp_si(q, -1, 1) == 0)
+            str[j++] = '-';
+        else
+        {
+            mpq_get_str(str, 10, q);
+            for ( ; str[j] != '\0'; j++);
+            str[j++] = '*';
+        }
+    }
+    sprintf(str + j, "%s", var);
+    j += varlen;
+    if (fmpz_poly_degree(op->num) > 1)
+    {
+        str[j++] = '^';
+        sprintf(str + j, "%li", fmpz_poly_degree(op->num));
+        for ( ; str[j] != '\0'; j++);
+    }
+    
+    for (i = fmpq_poly_degree(op) - 1; i >= 0; i--)
+    {
+        if (fmpz_is_zero(fmpz_poly_get_coeff_ptr(op->num, i)))
+            continue;
+        
+        fmpz_poly_get_coeff_mpz(mpq_numref(q), op->num, i);
+        if (op->den == NULL)
+            mpz_set_si(mpq_denref(q), 1);
+        else
+            fmpz_to_mpz(mpq_denref(q), op->den);
+        mpq_canonicalize(q);
+        
+        str[j++] = ' ';
+        if (mpq_sgn(q) < 0)
+        {
+            mpq_abs(q, q);
+            str[j++] = '-';
+        }
+        else
+            str[j++] = '+';
+        str[j++] = ' ';
+        
+        mpq_get_str(str + j, 10, q);
+        for ( ; str[j] != '\0'; j++);
+        
+        if (i > 0)
+        {
+            str[j++] = '*';
+            sprintf(str + j, "%s", var);
+            j += varlen;
+            if (i > 1)
+            {
+                str[j++] = '^';
+                sprintf(str + j, "%li", i);
+                for ( ; str[j] != '\0'; j++);
+            }
+        }
+    }
+    
+    mpq_clear(q);
+    mpz_clear(z);
+    return str;
+}
+

new file sage/libs/flint/fmpq_poly.h

@@ +1 @@
+///////////////////////////////////////////////////////////////////////////////
+// Copyright (C) 2010 Sebastian Pancratz                                     //
+//                                                                           //
+// Distributed under the terms of the GNU General Public License as          //
+// published by the Free Software Foundation; either version 2 of the        //
+// License, or (at your option) any later version.                           //
+//                                                                           //
+// http://www.gnu.org/licenses/                                              //
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef _FMPQ_POLY_H_
+#define _FMPQ_POLY_H_
+
+#include <stdio.h>
+#include <gmp.h>
+
+#include "fmpz.h"
+#include "fmpz_poly.h"
+
+/**
+ * \file     fmpq_poly.h
+ * \brief    Fast implementation of the rational polynomial ring
+ * \author   Sebastian Pancratz
+ * \date     Mar 2010 -- July 2010
+ * \version  0.1.2
+ */
+
+/**
+ * \ingroup  Definitions
+ * \brief    Data type for a rational polynomial.
+ * \details  We represent a rational polynomial as the quotient of an integer 
+ *           polynomial of type <tt>fmpz_poly_t</tt> and a denominator of type 
+ *           <tt>fmpz_t</tt>, enforcing coprimality and that the denominator 
+ *           is positive.  Moreover, the denominator may be <tt>NULL</tt>, 
+ *           which is interpreted as \f$1\f$.  The zero polynomial is 
+ *           represented as \f$0/1\f$.
+ */
+typedef struct
+{
+    fmpz_poly_p num; //!< Numerator
+    fmpz_t den;      //!< Denominator
+} fmpq_poly_struct;
+
+/**
+ * \ingroup  Definitions
+ * \brief    Array of #fmpq_poly_struct of length one.
+ * \details  This allows passing <em>by reference</em> without having to 
+ *           explicitly allocate memory for the structure, as one would have 
+ *           to when using pointers.
+ */
+typedef fmpq_poly_struct fmpq_poly_t[1];
+
+/**
+ * \ingroup  Definitions
+ * \brief    Pointer to #fmpq_poly_struct.
+ */
+typedef fmpq_poly_struct * fmpq_poly_ptr;
+
+void fmpq_poly_canonicalize(fmpq_poly_ptr f, fmpz_t temp);
+
+///////////////////////////////////////////////////////////////////////////////
+// _fmpq_poly_*                                                              //
+///////////////////////////////////////////////////////////////////////////////
+
+/**
+ * \ingroup  Helper
+ * 
+ * Ensures that <tt>op-\>den</tt> can hold <tt>limbs</tt> limbs (plus a 
+ * sign/ size limb).
+ * 
+ * This method never shrinks <tt>op-\>den</tt>.  Moreover, after a call to 
+ * this method, <tt>op-\>den</tt> is guaranteed to be initialised.
+ */
+static inline
+void _fmpq_poly_den_fit_limbs(fmpq_poly_ptr op, unsigned long limbs)
+{
+    if (op->den == NULL)
+    {
+        op->den = fmpz_init(limbs > 1 ? limbs : 1);
+        fmpz_set_si(op->den, 1);  /* XXX: FLINT, Valgrind problem */
+    }
+    else if (fmpz_size(op->den) < limbs)
+        op->den = fmpz_realloc(op->den, limbs);
+}
+
+/**
+ * \ingroup  Helper
+ * 
+ * Compares the denominators of <tt>op1</tt> and <tt>op2</tt>.
+ * 
+ * This method handles <tt>NULL</tt> values correctly.  This method only deals 
+ * with the <tt>den</tt> attributes of the arguments.  In particular, it does 
+ * not ensure that the contents of the numerators are coprime to the 
+ * denominators.
+ * 
+ * All arguments may be aliased.
+ */
+static inline
+int _fmpq_poly_den_equal(fmpq_poly_ptr op1, fmpq_poly_ptr op2)
+{
+    if (op1->den == NULL)
+        return op2->den == NULL || fmpz_is_one(op2->den);
+    else
+        return op2->den == NULL ? fmpz_is_one(op1->den) 
+                                : fmpz_equal(op1->den, op2->den);
+}
+
+/**
+ * \ingroup  Helper
+ * 
+ * Returns whether the denominator of <tt>op</tt> is one.
+ * 
+ * This method does not check that the denominator is coprime to the 
+ * numerator.
+ */
+static inline
+int _fmpq_poly_den_is_one(fmpq_poly_ptr op)
+{
+    return op->den == NULL || fmpz_is_one(op->den);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// fmpq_poly_*                                                               //
+///////////////////////////////////////////////////////////////////////////////
+
+///////////////////////////////////////////////////////////////////////////////
+// Accessing numerator and denominator
+
+/**
+ * \def      fmpq_poly_numref(op)
+ * \brief    Returns a reference to the numerator of \c op.
+ * \ingroup  NumDen
+ * \details  The <tt>fmpz_poly</tt> functions can be used on the result of 
+ *           this macro.
+ */
+#define fmpq_poly_numref(op)  ((op)->num)
+
+/**
+ * \def      fmpq_poly_denref(op)
+ * \brief    Returns a reference to the denominator of \c op.
+ * \ingroup  NumDen
+ * \details  The <tt>fmpz</tt> functions can be used on the result of 
+ *           this macro \e if the result is not <tt>NULL</tt>.
+ */
+#define fmpq_poly_denref(op)  ((op)->den)
+
+/**
+ * \brief    Sets \c rop to the numerator of \c op.
+ * \ingroup  NumDen
+ * 
+ * \note  The arguments <tt>rop</tt> and <tt>op-\>num</tt> may be aliased.
+ * 
+ * \param[out] rop  Integer polynomial
+ * \param[in] op    Rational polynomial
+ */
+static inline 
+void fmpq_poly_get_num(fmpz_poly_p rop, const fmpq_poly_ptr op)
+{
+    fmpz_poly_set(rop, op->num);
+}
+
+/**
+ * \brief    Sets \c rop to the denominator of \c op.
+ * \ingroup  NumDen
+ * 
+ * This method handles <tt>NULL</tt> values correctly.
+ * 
+ * \note  The arguments <tt>rop</tt> and <tt>op-\>den</tt> may be aliased.
+ * 
+ * \param[out] rop  Integer
+ * \param[in] op    Rational polynomial
+ */
+static inline 
+void fmpq_poly_get_den(fmpz_t rop, const fmpq_poly_ptr op)
+{
+    if (rop == op->den)  /* Aliasing */
+        return;
+    
+    if (op->den == NULL)
+        fmpz_set_si(rop, 1);
+    else
+    {
+        fmpz_set_si(rop, 1);
+        if (fmpz_size(rop) < fmpz_size(op->den))
+            rop = fmpz_realloc(rop, fmpz_size(op->den));
+        fmpz_set(rop, op->den);
+    }
+}
+
+/**
+ * \brief    Sets the numerator of \c rop to the polynomial \c op. 
+ * \ingroup  NumDen
+ * 
+ * \note  The arguments <tt>rop-\>num</tt> and <tt>op</tt> may be aliased.
+ * 
+ * \param[out] rop  Rational polynomial
+ * \param[in] op    Integer polynomial
+ */
+static inline 
+void fmpq_poly_set_num(fmpq_poly_ptr rop, const fmpz_poly_p op)
+{
+    fmpz_poly_set(rop->num, op);
+}
+
+/**
+ * \brief    Sets the denominator of \c rop to the integer \c op. 
+ * \ingroup  NumDen
+ * 
+ * This method handles <tt>NULL</tt> values correctly.  It does \e not ensure 
+ * that \c rop ends up in canonical form.
+ * 
+ * The arguments <tt>rop-\>den</tt> and <tt>op</tt> may be aliased.
+ * 
+ * \param[out] rop  Rational polynomial
+ * \param[in] op    Integer
+ */
+static inline 
+void fmpq_poly_set_den(fmpq_poly_ptr rop, const fmpz_t op)
+{
+    if (op == NULL)
+    {
+        if (rop->den != NULL)
+            fmpz_set_si(rop->den, 1);
+    }
+    else
+    {
+        if (rop->den == NULL)
+        {
+            if (!fmpz_is_one(op))
+            {
+                rop->den = fmpz_init(fmpz_size(op));
+                fmpz_set(rop->den, op);
+            }
+        }
+        else
+        {
+            if (fmpz_size(rop->den) < fmpz_size(op))
+                rop->den = fmpz_realloc(rop->den, fmpz_size(op));
+            fmpz_set(rop->den, op);
+        }
+    }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Memory management
+
+/**
+ * \ingroup  MemoryManagement
+ * 
+ * Initializes the element \c rop to zero.
+ * 
+ * This function should be called once for the #fmpq_poly_ptr \c rop, before 
+ * using it with other <tt>fmpq_poly</tt> functions, or following a 
+ * preceeding call to #fmpq_poly_clear().
+ */
+static inline 
+void fmpq_poly_init(fmpq_poly_ptr rop)
+{
+    rop->num = (fmpz_poly_p) malloc(sizeof(fmpz_poly_struct));
+    fmpz_poly_init(rop->num);
+    rop->den = NULL;
+}
+
+/**
+ * \ingroup  MemoryManagement
+ *
+ * Clears the element \c rop.
+ * 
+ * This function should only be called on an element which has been 
+ * initialised.
+ */
+static inline 
+void fmpq_poly_clear(fmpq_poly_ptr rop)
+{
+    if (rop->num != NULL)
+    {
+        fmpz_poly_clear(rop->num);
+        free(rop->num);
+        rop->num = NULL;
+    }
+    if (rop->den != NULL);
+    {
+        fmpz_clear(rop->den);
+        rop->den = NULL;
+    }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Assignment and basic manipulation
+
+/**
+ * \ingroup  Assignment
+ * 
+ * Sets the element \c rop to the same value as the element \c op.
+ */
+static inline 
+void fmpq_poly_set(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
+{
+    if (rop != op)
+    {
+        fmpz_poly_set(rop->num, op->num);
+        if (op->den == NULL)
+        {
+            if (rop->den != NULL)
+            {
+                fmpz_clear(rop->den);
+                rop->den = NULL;
+            }
+        }
+        else
+        {
+            _fmpq_poly_den_fit_limbs(rop, fmpz_size(op->den));
+            fmpz_set(rop->den, op->den);
+        }
+    }
+}
+
+/**
+ * \ingroup  Assignment
+ * 
+ * Sets the element \c rop to the value given by the \c int \c op.
+ */
+static inline  
+void fmpq_poly_set_si(fmpq_poly_ptr rop, const long op)
+{
+    fmpz_poly_zero(rop->num);
+    fmpz_poly_set_coeff_si(rop->num, 0, op);
+    if (rop->den != NULL)
+        fmpz_set_si(rop->den, 1);
+}
+
+/**
+ * \ingroup  Assignment
+ * 
+ * Sets the element \c rop to the integer \c x.
+ */
+static inline
+void fmpq_poly_set_fmpz(fmpq_poly_ptr rop, const fmpz_t x)
+{
+    fmpz_poly_zero(rop->num);
+    fmpz_poly_set_coeff_fmpz(rop->num, 0, x);
+    if (rop->den != NULL)
+        fmpz_set_si(rop->den, 1);
+}
+
+/**
+ * \ingroup  Assignment
+ * 
+ * Sets the element \c rop to the integer \c x.
+ */
+static inline
+void fmpq_poly_set_mpz(fmpq_poly_ptr rop, const mpz_t x)
+{
+    fmpz_poly_zero(rop->num);
+    fmpz_poly_set_coeff_mpz(rop->num, 0, x);
+    if (rop->den != NULL)
+        fmpz_set_si(rop->den, 1);
+}
+
+/**
+ * \ingroup  Assignment
+ * 
+ * Sets the element \c rop to the rational \c x.
+ */
+static inline
+void fmpq_poly_set_mpq(fmpq_poly_ptr rop, const mpq_t x)
+{
+    fmpz_poly_zero(rop->num);
+    fmpz_poly_set_coeff_mpz(rop->num, 0, mpq_numref(x));
+    if (mpz_cmp_si(mpq_denref(x), 1))  /* x->den > 1 */
+    {
+        _fmpq_poly_den_fit_limbs(rop, mpz_size(mpq_denref(x)));
+        mpz_to_fmpz(rop->den, mpq_denref(x));
+    }
+    else                               /* x.den == 1 */
+        if (rop->den != NULL)
+            fmpz_set_si(rop->den, 1);
+}
+
+/**
+ * \ingroup  Assignment
+ * 
+ * Swaps the elements \c op1 and \c op2.
+ * 
+ * This is done efficiently by swapping pointers.
+ */
+static inline 
+void fmpq_poly_swap(fmpq_poly_ptr op1, fmpq_poly_ptr op2)
+{
+    fmpz_t t;
+    if (op1 != op2)
+    {
+        fmpz_poly_swap(op1->num, op2->num);
+        t        = op1->den;
+        op1->den = op2->den;
+        op2->den = t;
+    }
+}
+
+/**
+ * \ingroup  Assignment
+ * 
+ * Sets the element \c rop to zero.
+ */
+static inline
+void fmpq_poly_zero(fmpq_poly_ptr rop)
+{
+    fmpz_poly_zero(rop->num);
+    if (rop->den != NULL)
+        fmpz_set_si(rop->den, 1);
+}
+
+/**
+ * \ingroup  Assignment
+ * 
+ * Sets the element \c rop to one.
+ */
+static inline
+void fmpq_poly_one(fmpq_poly_ptr rop)
+{
+    fmpz_poly_zero(rop->num);
+    fmpz_poly_set_coeff_si(rop->num, 0, 1);
+    if (rop->den != NULL)
+        fmpz_set_si(rop->den, 1);
+}
+
+void fmpq_poly_neg(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
+void fmpq_poly_inv(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
+
+///////////////////////////////////////////////////////////////////////////////
+// Setting/ retrieving individual coefficients
+
+void fmpq_poly_get_coeff_mpq(mpq_t rop, const fmpq_poly_ptr op, unsigned long i);
+
+void fmpq_poly_set_coeff_fmpz(fmpq_poly_ptr rop, const unsigned long i, const fmpz_t x);
+void fmpq_poly_set_coeff_mpz(fmpq_poly_ptr rop, const unsigned long i, const mpz_t x);
+void fmpq_poly_set_coeff_mpq(fmpq_poly_ptr rop, const unsigned long i, const mpq_t x);
+void fmpq_poly_set_coeff_si(fmpq_poly_ptr rop, const unsigned long i, long x);
+
+///////////////////////////////////////////////////////////////////////////////
+// Comparison
+
+/**
+ * \brief    Returns whether \c op is zero.
+ * \ingroup  Comparison
+ * 
+ * Returns whether the element \c op is zero.
+ */
+static inline
+int fmpq_poly_is_zero(const fmpq_poly_ptr op)
+{
+    return fmpz_poly_degree(op->num) < 0;
+}
+
+/**
+ * \brief    Returns whether \c op is equal to \f$1\f$.
+ * \ingroup  Comparison
+ * 
+ * Returns whether the element \c op is equal to the constant polynomial 
+ * \f$1\f$.
+ */
+static inline
+int fmpq_poly_is_one(const fmpq_poly_ptr op)
+{
+    return (op->num->length == 1) && fmpz_is_one(op->num->coeffs) 
+        && _fmpq_poly_den_is_one(op);
+}
+
+int fmpq_poly_equal(const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
+int fmpq_poly_cmp(const fmpq_poly_ptr left, const fmpq_poly_ptr right);
+
+///////////////////////////////////////////////////////////////////////////////
+// Polynomial parameters
+
+/**
+ * \brief    Returns the length of \c op.
+ * \ingroup  PolyParameters
+ * \details  Returns the length of the polynomial \c op, which is one greater 
+ *           than its degree.
+ */
+static inline 
+unsigned long fmpq_poly_length(const fmpq_poly_ptr op)
+{
+    return fmpz_poly_length(op->num);
+}
+
+/**
+ * \brief    Returns the degree of \c op.
+ * \ingroup  Parameters
+ * \details  Returns the degree of the polynomial \c op.
+ */
+static inline
+long fmpq_poly_degree(const fmpq_poly_ptr op)
+{
+    return fmpz_poly_degree(op->num);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Addition/ subtraction
+
+void fmpq_poly_add(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
+void fmpq_poly_sub(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
+
+void fmpq_poly_addmul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
+void fmpq_poly_submul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
+
+///////////////////////////////////////////////////////////////////////////////
+// Scalar multiplication and division
+
+void fmpq_poly_scalar_mul_si(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long x);
+void fmpq_poly_scalar_mul_mpz(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpz_t x);
+void fmpq_poly_scalar_mul_mpq(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x);
+
+void fmpq_poly_scalar_div_si(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long x);
+void fmpq_poly_scalar_div_mpz(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpz_t x);
+void fmpq_poly_scalar_div_mpq(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x);
+
+///////////////////////////////////////////////////////////////////////////////
+// Multiplication
+
+void fmpq_poly_mul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
+
+///////////////////////////////////////////////////////////////////////////////
+// Division
+
+void fmpq_poly_floordiv(fmpq_poly_ptr q, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
+void fmpq_poly_mod(fmpq_poly_ptr r, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
+void fmpq_poly_divrem(fmpq_poly_ptr q, fmpq_poly_ptr r, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
+
+///////////////////////////////////////////////////////////////////////////////
+// Powering
+
+void fmpq_poly_power(fmpq_poly_ptr rop, const fmpq_poly_ptr op, unsigned long exp);
+
+///////////////////////////////////////////////////////////////////////////////
+// Greatest common divisor
+
+void fmpq_poly_gcd(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
+void fmpq_poly_xgcd(fmpq_poly_ptr rop, fmpq_poly_ptr s, fmpq_poly_ptr t, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
+void fmpq_poly_lcm(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
+
+///////////////////////////////////////////////////////////////////////////////
+// Derivative
+
+void fmpq_poly_derivative(fmpq_poly_ptr rop, fmpq_poly_ptr op);
+
+///////////////////////////////////////////////////////////////////////////////
+// Evaluation
+
+void fmpq_poly_evaluate_mpz(mpq_t rop, fmpq_poly_ptr f, const mpz_t a);
+void fmpq_poly_evaluate_mpq(mpq_t rop, fmpq_poly_ptr f, const mpq_t a);
+
+///////////////////////////////////////////////////////////////////////////////
+// Gaussian content
+
+void fmpq_poly_content(mpq_t rop, const fmpq_poly_ptr op);
+void fmpq_poly_primitive_part(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
+int fmpq_poly_is_monic(const fmpq_poly_ptr op);
+void fmpq_poly_monic(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
+
+///////////////////////////////////////////////////////////////////////////////
+// Resultant
+
+void fmpq_poly_resultant(mpq_t rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
+void fmpq_poly_discriminant(mpq_t disc, fmpq_poly_t a);
+
+///////////////////////////////////////////////////////////////////////////////
+// Composition
+
+void fmpq_poly_compose(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
+void fmpq_poly_rescale(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x);
+
+///////////////////////////////////////////////////////////////////////////////
+// Square-free
+
+int fmpq_poly_is_squarefree(const fmpq_poly_ptr op);
+
+///////////////////////////////////////////////////////////////////////////////
+// Subpolynomials
+
+void fmpq_poly_getslice(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long i, long j);
+void fmpq_poly_shift(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long n);
+void fmpq_poly_truncate(fmpq_poly_ptr rop, const fmpq_poly_ptr op, unsigned long n);
+void fmpq_poly_reverse(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
+
+///////////////////////////////////////////////////////////////////////////////
+// String conversion
+
+int fmpq_poly_from_string(fmpq_poly_ptr rop, const char * s);
+char * fmpq_poly_to_string(const fmpq_poly_ptr op);
+char * fmpq_poly_to_string_pretty(const fmpq_poly_ptr op, const char * x);
+
+#endif
+

new file sage/libs/flint/fmpq_poly.pxd

@@ +1 @@
+#include "fmpz_poly.pxi"
+#include "fmpz.pxi"
+
+cdef extern from "fmpq_poly.h":
+    ctypedef void* fmpz_t
+    ctypedef void* fmpz_poly_p
+    struct fmpq_poly:
+        fmpz_poly_p num
+        fmpz_t den
+        
+    ctypedef fmpq_poly fmpq_poly_struct
+    ctypedef fmpq_poly_struct fmpq_poly_t[1]
+
+    void fmpq_poly_init(fmpq_poly_t)
+    void fmpq_poly_from_string(fmpq_poly_t, char*)
+    char* fmpq_poly_to_string_pretty(fmpq_poly_t, char*)
+    void fmpq_poly_mul(fmpq_poly_t, fmpq_poly_t, fmpq_poly_t)
+    void fmpq_poly_clear(fmpq_poly_t)

new file sage/libs/flint/fmpq_poly_example.pyx

@@ +1 @@
+###############################################################################
+# Copyright (C) 2010 Sebastian Pancratz                                       #
+#                                                                             #
+# Distributed under the terms of the GNU General Public License as            #
+# published by the Free Software Foundation; either version 2 of the          #
+# License, or (at your option) any later version.                             #
+#                                                                             #
+# http:#www.gnu.org/licenses/                                                 #
+###############################################################################
+
+# This file contains a simple example for how to use fmpq_poly package.
+
+include "fmpq_poly.pxd"
+
+cdef extern from "stdlib.h":
+    void free(void*)
+
+def example():
+    """
+    This is an example of how to use fmpq from Cython.
+    
+    EXAMPLES::
+    
+        sage: import sage.libs.flint.fmpq_poly_example as e
+        sage: e.example()
+        f     = 3/5*t + 1/2
+        g     = 2/7*t^2 - 3/2*t - 1
+        f * g = 6/35*t^3 - 53/70*t^2 - 27/20*t - 1/2
+    """
+    cdef char *str, *strf, *strg
+    cdef fmpq_poly_t f, g
+    fmpq_poly_init(f)
+    fmpq_poly_init(g)
+    fmpq_poly_from_string(f, "2  1/2 3/5")
+    fmpq_poly_from_string(g, "3  -1 -3/2 2/7")
+    strf = fmpq_poly_to_string_pretty(f, "t")
+    strg = fmpq_poly_to_string_pretty(g, "t")
+    print("f     = %s"%strf)
+    print("g     = %s"%strg)
+    fmpq_poly_mul(f, f, g)
+    str  = fmpq_poly_to_string_pretty(f, "t")
+    print("f * g = %s"%str);
+    free(str)
+    free(strf)
+    free(strg)
+    fmpq_poly_clear(f)
+    fmpq_poly_clear(g)
+
+