Ticket #4000: trac_4000-combined.2.patch

File trac_4000-combined.2.patch, 219.2 KB (added by mpatel, 8 years ago)

Fix 32/64-bit number_field_ideal doctest. Apply only this patch.

  • module_list.py

    # HG changeset patch
    # User Mitesh Patel <qed777@gmail.com>
    # Date 1286533396 25200
    # Node ID c53a9e75930a099e4fc00fbd72a0d02c8ff14215
    # Parent  9776abf682fb70fad1c77fdb1b478e06af715e6d
    #4000: Implement QQ['x'] via FLINT's ZZ['x'] + denominator
    * * *
    Trac 4000.  Implements QQ[] via FLINT
    * * *
    Trac 4000.  C code and header files for FMPQ_POLY.
    * * *
    Trac 4000.  Adds fmpq_poly.pxd to provide an interface to the C code
    * * *
    trac 4000: Implement QQ['x'] via Flint ZZ['x'] + denominator -- 64-bit update
    * * *
    trac 4000: is doctest
    * * *
    #4000: Fix PARI error in hensel_lift()
    * * *
    4000: Adjust graph-vertex sorting doctest to new module names
    * * *
    #4000: QQ via FLINT ZZ, work around a Solaris gcc problem
    
    diff --git a/module_list.py b/module_list.py
    a b ext_modules = [ 
    13831383              language = 'c++',
    13841384              include_dirs = ['sage/libs/ntl/']),
    13851385
     1386    Extension('sage.rings.polynomial.polynomial_rational_flint',
     1387              sources = ['sage/rings/polynomial/polynomial_rational_flint.pyx', 'sage/libs/flint/fmpq_poly.c'],
     1388              language = 'c++',
     1389              extra_compile_args=["-std=c99"] + uname_specific('SunOS', [], ['-D_XPG6']),
     1390              libraries = ["csage", "flint", "ntl", "gmpxx", "gmp"],
     1391              include_dirs = [SAGE_ROOT + '/local/include/FLINT/', SAGE_ROOT + '/devel/sage/sage/libs/flint/'],
     1392              depends = [SAGE_ROOT + "/local/include/FLINT/flint.h"]),
     1393
    13861394    Extension('sage.rings.polynomial.polynomial_modn_dense_ntl',
    13871395              sources = ['sage/rings/polynomial/polynomial_modn_dense_ntl.pyx'],
    13881396              libraries = ['ntl', 'stdc++', 'gmp'],
  • sage/calculus/all.py

    diff --git a/sage/calculus/all.py b/sage/calculus/all.py
    a b def symbolic_expression(x): 
    4040        sage: type(a)
    4141        <type 'sage.symbolic.expression.Expression'>
    4242        sage: R.<x> = QQ[]; type(x)
    43         <class 'sage.rings.polynomial.polynomial_element_generic.Polynomial_rational_dense'>
     43        <type 'sage.rings.polynomial.polynomial_rational_flint.Polynomial_rational_flint'>
    4444        sage: a = symbolic_expression(2*x^2 + 3); a
    4545        2*x^2 + 3
    4646        sage: type(a)
  • sage/categories/primer.py

    diff --git a/sage/categories/primer.py b/sage/categories/primer.py
    a b Applying an operation is generally done  
    8080    sage: p.factor()
    8181    6*(x + 1)^2
    8282
    83     sage: pQ = QQ[x] ( p )
     83    sage: R.<x> = PolynomialRing(QQ, sparse=True)
     84    sage: pQ = R ( p )
    8485    sage: type(pQ)
    85     <class 'sage.rings.polynomial.polynomial_element_generic.Polynomial_rational_dense'>
     86    <class 'sage.rings.polynomial.polynomial_element_generic.Polynomial_generic_sparse_field'>
    8687    sage: pQ.factor()
    8788    (6) * (x + 1)^2
    8889
  • sage/combinat/sf/jack.py

    diff --git a/sage/combinat/sf/jack.py b/sage/combinat/sf/jack.py
    a b class JackPolynomials_p(JackPolynomials_ 
    495495            sage: l = lambda c: [ (i[0],[j for j in sorted(i[1].items())]) for i in sorted(c.items())]
    496496            sage: P._m_cache(2)
    497497            sage: l(P._self_to_m_cache[2])
    498             [([1, 1], [([1, 1], 1)]), ([2], [([1, 1], 2/(t + 1)), ([2], 1)])]
     498            [([1, 1], [([1, 1], 1)]), ([2], [([1, 1], -2/(-t - 1)), ([2], 1)])]
    499499            sage: l(P._m_to_self_cache[2])
    500             [([1, 1], [([1, 1], 1)]), ([2], [([1, 1], -2/(t + 1)), ([2], 1)])]
     500            [([1, 1], [([1, 1], 1)]), ([2], [([1, 1], 2/(-t - 1)), ([2], 1)])]
    501501            sage: P._m_cache(3)
    502502            sage: l(P._m_to_self_cache[3])
    503503            [([1, 1, 1], [([1, 1, 1], 1)]),
    504504             ([2, 1], [([1, 1, 1], -6/(t + 2)), ([2, 1], 1)]),
    505              ([3], [([1, 1, 1], 6/(t^2 + 3*t + 2)), ([2, 1], -3/(2*t + 1)), ([3], 1)])]
     505             ([3], [([1, 1, 1], -6/(-t^2 - 3*t - 2)), ([2, 1], -3/(2*t + 1)), ([3], 1)])]
    506506            sage: l(P._self_to_m_cache[3])
    507507            [([1, 1, 1], [([1, 1, 1], 1)]),
    508508             ([2, 1], [([1, 1, 1], 6/(t + 2)), ([2, 1], 1)]),
    509              ([3], [([1, 1, 1], 6/(2*t^2 + 3*t + 1)), ([2, 1], 3/(2*t + 1)), ([3], 1)])]
     509             ([3], [([1, 1, 1], -6/(-2*t^2 - 3*t - 1)), ([2, 1], 3/(2*t + 1)), ([3], 1)])]
    510510        """
    511511        if n in self._self_to_m_cache:
    512512            return
    class JackPolynomials_p(JackPolynomials_ 
    536536            sage: p11 = Partition([1,1])
    537537            sage: f = P._to_m(p2)
    538538            sage: f(p11)
    539             2/(t + 1)
     539            -2/(-t - 1)
    540540            sage: f(p2)
    541541            1
    542542            sage: f = P._to_m(p11)
    class JackPolynomials_p(JackPolynomials_ 
    554554       
    555555            sage: P = JackPolynomialsP(QQ)
    556556            sage: P([1])^2 # indirect doctest
    557             (2*t/(t+1))*JackP[1, 1] + JackP[2]
     557            (-2*t/(-t-1))*JackP[1, 1] + JackP[2]
    558558            sage: P._m(_)
    559559            2*m[1, 1] + m[2]
    560560            sage: P = JackPolynomialsP(QQ, 2)
    class JackPolynomials_j(JackPolynomials_ 
    685685       
    686686            sage: s = SFASchur(J.base_ring())
    687687            sage: J(s([3])) # indirect doctest
    688             ((t^2-3*t+2)/(6*t^2+18*t+12))*JackJ[1, 1, 1] + ((2*t-2)/(2*t^2+5*t+2))*JackJ[2, 1] + (1/(2*t^2+3*t+1))*JackJ[3]
     688            ((-t^2+3*t-2)/(-6*t^2-18*t-12))*JackJ[1, 1, 1] + ((2*t-2)/(2*t^2+5*t+2))*JackJ[2, 1] + (1/(2*t^2+3*t+1))*JackJ[3]
    689689            sage: J(s([2,1]))
    690690            ((t-1)/(3*t+6))*JackJ[1, 1, 1] + (1/(t+2))*JackJ[2, 1]
    691691            sage: J(s([1,1,1]))
    class JackPolynomials_j(JackPolynomials_ 
    698698       
    699699            sage: J = JackPolynomialsJ(QQ)
    700700            sage: J([1])^2 #indirect doctest
    701             (t/(t+1))*JackJ[1, 1] + (1/(t+1))*JackJ[2]
     701            (-t/(-t-1))*JackJ[1, 1] + (1/(t+1))*JackJ[2]
    702702            sage: J([2])^2
    703             (2*t^2/(2*t^2+3*t+1))*JackJ[2, 2] + (4*t/(3*t^2+4*t+1))*JackJ[3, 1] + ((t+1)/(6*t^2+5*t+1))*JackJ[4]
     703            (-2*t^2/(-2*t^2-3*t-1))*JackJ[2, 2] + (-4*t/(-3*t^2-4*t-1))*JackJ[3, 1] + ((t+1)/(6*t^2+5*t+1))*JackJ[4]
    704704        """
    705705        return self( self._P(left) * self._P(right) )
    706706
  • sage/graphs/generic_graph.py

    diff --git a/sage/graphs/generic_graph.py b/sage/graphs/generic_graph.py
    a b class GenericGraph(GenericGraph_pyx): 
    65566556            sage: G.vertices(key = lambda x: (x[1], x[2], x[0]))
    65576557            [(0, 0, 0), (1, 0, 0), (2, 0, 0), (0, 0, 1), ... (1, 2, 2), (2, 2, 2)]
    65586558
    6559         ::
     6559        The discriminant of a polynomial is a function that returns an integer.
     6560        We build a graph whose vertices are polynomials, and use the discriminant
     6561        function to provide an ordering.  Note that since functions are first-class
     6562        objects in Python, we can specify precisely the function from the Sage library
     6563        that we wish to use as the key. ::
    65606564
    65616565            sage: t = polygen(QQ, 't')
    65626566            sage: K = Graph({5*t:[t^2], t^2:[t^2+2], t^2+2:[4*t^2-6], 4*t^2-6:[5*t]})
    6563             sage: dsc = sage.rings.polynomial.polynomial_element_generic.Polynomial_rational_dense.discriminant
    6564             sage: K.vertices(key=dsc)
     6567            sage: dsc = sage.rings.polynomial.polynomial_rational_flint.Polynomial_rational_flint.discriminant
     6568            sage: verts = K.vertices(key=dsc)
     6569            sage: verts
    65656570            [t^2 + 2, t^2, 5*t, 4*t^2 - 6]
     6571            sage: [x.discriminant() for x in verts]
     6572            [-8, 0, 1, 96]
    65666573
    65676574        If boundary vertices are requested first, then they are sorted
    65686575        separately from the remainder (which are also sorted). ::
  • new file sage/libs/flint/fmpq_poly.c

    diff --git a/sage/libs/flint/fmpq_poly.c b/sage/libs/flint/fmpq_poly.c
    new file mode 100644
    - +  
     1///////////////////////////////////////////////////////////////////////////////
     2// Copyright (C) 2010 Sebastian Pancratz                                     //
     3//                                                                           //
     4// Distributed under the terms of the GNU General Public License as          //
     5// published by the Free Software Foundation; either version 2 of the        //
     6// License, or (at your option) any later version.                           //
     7//                                                                           //
     8// http://www.gnu.org/licenses/                                              //
     9///////////////////////////////////////////////////////////////////////////////
     10
     11#ifdef __cplusplus
     12    extern "C" {
     13#endif
     14
     15#include "fmpq_poly.h"
     16
     17/**
     18 * \file     fmpq_poly.c
     19 * \brief    Fast implementation of the rational function field
     20 * \author   Sebastian Pancratz
     21 * \date     Mar 2010 -- July 2010
     22 * \version  0.1.8
     23 *
     24 * \mainpage
     25 *
     26 * <center>
     27 * A FLINT-based implementation of the rational polynomial ring.
     28 * </center>
     29 *
     30 * \section Overview
     31 *
     32 * The module <tt>fmpq_poly</tt> provides functions for performing
     33 * arithmetic on rational polynomials in \f$\mathbf{Q}[t]\f$, represented as
     34 * quotients of integer polynomials of type <tt>fmpz_poly_t</tt> and
     35 * denominators of type <tt>fmpz_t</tt>.  These functions start with the
     36 * prefix <tt>fmpq_poly_</tt>.
     37 *
     38 * Rational polynomials are stored in objects of type #fmpq_poly_t,
     39 * which is an array of #fmpq_poly_struct's of length one.  This
     40 * permits passing parameters of type #fmpq_poly_t by reference. 
     41 * We also define the type #fmpq_poly_ptr to be a pointer to
     42 * #fmpq_poly_struct's.
     43 *
     44 * The representation of a rational polynomial as the quotient of an integer
     45 * polynomial and an integer denominator can be made canonical by demanding
     46 * the numerator and denominator to be coprime and the denominator to be
     47 * positive.  As the only special case, we represent the zero function as
     48 * \f$0/1\f$.  All arithmetic functions assume that the operands are in this
     49 * canonical form, and canonicalize their result.  If the numerator or
     50 * denominator is modified individually, for example using the methods in the
     51 * group \ref NumDen, it is the user's responsibility to canonicalize the
     52 * rational function using the function #fmpq_poly_canonicalize() if necessary.
     53 *
     54 * All methods support aliasing of their inputs and outputs \e unless
     55 * explicitly stated otherwise, subject to the following caveat.  If
     56 * different rational polynomials (as objects in memory, not necessarily in
     57 * the mathematical sense) share some of the underlying integer polynomials 
     58 * or integers, the behaviour is undefined.
     59 *
     60 * \section Changelog  Version history
     61 * - 0.1.8
     62 *   - Extra checks for <tt>NULL</tt> returns of <tt>malloc</tt> calls
     63 * - 0.1.7
     64 *   - Explicit casts for return values of <tt>malloc</tt>
     65 *   - Changed a few calls to GMP and FLINT methods from <tt>_si</tt>
     66 *     to <tt>_ui</tt>
     67 * - 0.1.6
     68 *   - Added tons of testing code in the subdirectory <tt>test</tt>
     69 *   - Made a few changes throughout the code base indicated by the tests
     70 * - 0.1.5
     71 *   - Re-wrote the entire code to always initialise the denominator
     72 * - 0.1.4
     73 *   - Fixed a bug in #fmpq_poly_divrem()
     74 *   - Swapped calls to #fmpq_poly_degree to calls to #fmpq_poly_length in
     75 *     many places
     76 * - 0.1.3
     77 *   - Changed #fmpq_poly_inv()
     78 *   - Another sign check to <tt>fmpz_abs</tt> in #fmpq_poly_content()
     79 * - 0.1.2
     80 *   - Introduce a function #fmpq_poly_monic() and use this to simplify the
     81 *     code for the gcd and xgcd functions
     82 *   - Make further use of #fmpq_poly_is_zero()
     83 * - 0.1.1
     84 *   - Replaced a few sign checks and negations by <tt>fmpz_abs</tt>
     85 *   - Small changes to comments and the documentation
     86 *   - Moved some function bodies from #fmpq_poly.h to #fmpq_poly.c
     87 * - 0.1.0
     88 *   - First draft, based on the author's Sage code
     89 */
     90
     91/**
     92 * \defgroup Definitions        Type definitions
     93 * \defgroup MemoryManagement   Memory management
     94 * \defgroup NumDen             Accessing numerator and denominator
     95 * \defgroup Assignment         Assignment and basic manipulation
     96 * \defgroup Coefficients       Setting and retrieving individual coefficients
     97 * \defgroup PolyParameters     Polynomial parameters
     98 * \defgroup Comparison         Comparison
     99 * \defgroup Addition           Addition and subtraction
     100 * \defgroup ScalarMul          Scalar multiplication and division
     101 * \defgroup Multiplication     Multiplication
     102 * \defgroup Division           Euclidean division
     103 * \defgroup Powering           Powering
     104 * \defgroup GCD                Greatest common divisor
     105 * \defgroup Derivative         Derivative
     106 * \defgroup Evaluation         Evaluation
     107 * \defgroup Content            Gaussian content
     108 * \defgroup Resultant          Resultant
     109 * \defgroup Composition        Composition
     110 * \defgroup Squarefree         Square-free test
     111 * \defgroup Subpolynomials     Subpolynomials
     112 * \defgroup StringConversions  String conversions and I/O
     113 */
     114
     115///////////////////////////////////////////////////////////////////////////////
     116// Auxiliary functions                                                       //
     117///////////////////////////////////////////////////////////////////////////////
     118
     119/**
     120 * Returns the number of digits in the decimal representation of \c n.
     121 */
     122unsigned int fmpq_poly_places(ulong n)
     123{
     124    unsigned int count;
     125    if (n == 0)
     126        return 1u;
     127    count = 0;
     128    while (n > 0)
     129    {
     130        n /= 10ul;
     131        count++;
     132    }
     133    return count;
     134}
     135
     136///////////////////////////////////////////////////////////////////////////////
     137// Implementation                                                            //
     138///////////////////////////////////////////////////////////////////////////////
     139
     140/**
     141 * \brief    Puts <tt>f</tt> into canonical form.
     142 * \ingroup  Definitions
     143 *
     144 * This method ensures that the denominator is coprime to the content
     145 * of the numerator polynomial, and that the denominator is positive.
     146 *
     147 * The optional parameter <tt>temp</tt> is the temporary variable that this
     148 * method would otherwise need to allocate.  If <tt>temp</tt> is provided as
     149 * an initialized <tt>fmpz_t</tt>, it is assumed \e without further checks
     150 * that it is large enough.  For example,
     151 * <tt>max(f-\>num-\>limbs, fmpz_size(f-\>den))</tt> limbs will certainly
     152 * suffice.
     153 */
     154void fmpq_poly_canonicalize(fmpq_poly_ptr f, fmpz_t temp)
     155{
     156    if (fmpz_is_one(f->den))
     157        return;
     158   
     159    if (fmpq_poly_is_zero(f))
     160        fmpz_set_ui(f->den, 1ul);
     161   
     162    else if (fmpz_is_m1(f->den))
     163    {
     164        fmpz_poly_neg(f->num, f->num);
     165        fmpz_set_ui(f->den, 1ul);
     166    }
     167   
     168    else
     169    {
     170        int tcheck;  /* Whether temp == NULL */
     171        //if (temp == NULL)
     172        //{
     173            tcheck = 1;
     174            temp = fmpz_init(FLINT_MAX(f->num->limbs, fmpz_size(f->den)));
     175        //}
     176        //else
     177        //    tcheck = 0;
     178       
     179        fmpz_poly_content(temp, f->num);
     180        fmpz_abs(temp, temp);
     181       
     182        if (fmpz_is_one(temp))
     183        {
     184            if (fmpz_sgn(f->den) < 0)
     185            {
     186                fmpz_poly_neg(f->num, f->num);
     187                fmpz_neg(f->den, f->den);
     188            }
     189        }
     190        else
     191        {
     192            fmpz_t tempgcd;
     193            tempgcd = fmpz_init(FLINT_MAX(f->num->limbs, fmpz_size(f->den)));
     194           
     195            fmpz_gcd(tempgcd, temp, f->den);
     196            fmpz_abs(tempgcd, tempgcd);
     197           
     198            if (fmpz_is_one(tempgcd))
     199            {
     200                if (fmpz_sgn(f->den) < 0)
     201                {
     202                    fmpz_poly_neg(f->num, f->num);
     203                    fmpz_neg(f->den, f->den);
     204                }
     205            }
     206            else
     207            {
     208                if (fmpz_sgn(f->den) < 0)
     209                    fmpz_neg(tempgcd, tempgcd);
     210                fmpz_poly_scalar_div_fmpz(f->num, f->num, tempgcd);
     211                fmpz_tdiv(temp, f->den, tempgcd);
     212                fmpz_set(f->den, temp);
     213            }
     214           
     215            fmpz_clear(tempgcd);
     216        }
     217        if (tcheck)
     218            fmpz_clear(temp);
     219    }
     220}
     221
     222///////////////////////////////////////////////////////////////////////////////
     223// Assignment and basic manipulation
     224
     225/**
     226 * \ingroup  Assignment
     227 *
     228 * Sets the element \c rop to the additive inverse of \c op.
     229 */
     230void fmpq_poly_neg(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
     231{
     232    if (rop == op)
     233    {
     234        fmpz_poly_neg(rop->num, op->num);
     235    }
     236    else
     237    {
     238        fmpz_poly_neg(rop->num, op->num);
     239        rop->den = fmpz_realloc(rop->den, fmpz_size(op->den));
     240        fmpz_set(rop->den, op->den);
     241    }
     242}
     243
     244/**
     245 * \ingroup  Assignment
     246 *
     247 * Sets the element \c rop to the multiplicative inverse of \c op.
     248 *
     249 * Assumes that the element \c op is a unit.  Otherwise, an exception
     250 * is raised in the form of an <tt>abort</tt> statement.
     251 */
     252void fmpq_poly_inv(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
     253{
     254    /* Assertion! */
     255    if (fmpq_poly_length(op) != 1ul)
     256    {
     257        printf("ERROR (fmpq_poly_inv).  Element is not a unit.\n");
     258        abort();
     259    }
     260   
     261    if (rop == op)
     262    {
     263        fmpz_t t;
     264        t = fmpz_init(rop->num->limbs);
     265        fmpz_set(t, rop->num->coeffs);
     266        fmpz_poly_zero(rop->num);
     267        if (fmpz_sgn(t) > 0)
     268        {
     269            fmpz_poly_set_coeff_fmpz(rop->num, 0, rop->den);
     270            rop->den = fmpz_realloc(rop->den, fmpz_size(t));
     271            fmpz_set(rop->den, t);
     272        }
     273        else
     274        {
     275            fmpz_neg(rop->den, rop->den);
     276            fmpz_poly_set_coeff_fmpz(rop->num, 0, rop->den);
     277            rop->den = fmpz_realloc(rop->den, fmpz_size(t));
     278            fmpz_neg(rop->den, t);
     279        }
     280        fmpz_clear(t);
     281    }
     282    else
     283    {
     284        fmpz_poly_zero(rop->num);
     285        fmpz_poly_set_coeff_fmpz(rop->num, 0, op->den);
     286        rop->den = fmpz_realloc(rop->den, fmpz_size(op->num->coeffs));
     287        fmpz_set(rop->den, op->num->coeffs);
     288        if (fmpz_sgn(rop->den) < 0)
     289        {
     290            fmpz_poly_neg(rop->num, rop->num);
     291            fmpz_neg(rop->den, rop->den);
     292        }
     293    }
     294}
     295
     296///////////////////////////////////////////////////////////////////////////////
     297// Setting/ retrieving individual coefficients
     298
     299/**
     300 * \ingroup  Coefficients
     301 *
     302 * Obtains the <tt>i</tt>th coefficient from the polynomial <tt>op</tt>.
     303 */
     304void fmpq_poly_get_coeff_mpq(mpq_t rop, const fmpq_poly_ptr op, ulong i)
     305{
     306    fmpz_poly_get_coeff_mpz(mpq_numref(rop), op->num, i);
     307    fmpz_to_mpz(mpq_denref(rop), op->den);
     308    mpq_canonicalize(rop);
     309}
     310
     311/**
     312 * \ingroup  Coefficients
     313 *
     314 * Sets the <tt>i</tt>th coefficient of the polynomial <tt>op</tt> to \c x.
     315 */
     316void fmpq_poly_set_coeff_fmpz(fmpq_poly_ptr rop, ulong i, const fmpz_t x)
     317{
     318    fmpz_t t;
     319    int canonicalize;
     320   
     321    /* Is the denominator 1?  This includes the case when rop == 0. */
     322    if (fmpz_is_one(rop->den))
     323    {
     324        fmpz_poly_set_coeff_fmpz(rop->num, i, x);
     325        return;
     326    }
     327   
     328    t = fmpz_poly_get_coeff_ptr(rop->num, i);
     329    canonicalize = !(t == NULL || fmpz_is_zero(t));
     330   
     331    t = fmpz_init(fmpz_size(x) + fmpz_size(rop->den));
     332    fmpz_mul(t, x, rop->den);
     333    fmpz_poly_set_coeff_fmpz(rop->num, i, t);
     334    fmpz_clear(t);
     335    if (canonicalize)
     336        fmpq_poly_canonicalize(rop, NULL);
     337}
     338
     339/**
     340 * \ingroup  Coefficients
     341 *
     342 * Sets the <tt>i</tt>th coefficient of the polynomial <tt>op</tt> to \c x.
     343 */
     344void fmpq_poly_set_coeff_mpz(fmpq_poly_ptr rop, ulong i, const mpz_t x)
     345{
     346    mpz_t z;
     347    fmpz_t t;
     348    int canonicalize;
     349   
     350    /* Is the denominator 1?  This includes the case when rop == 0. */
     351    if (fmpz_is_one(rop->den))
     352    {
     353        fmpz_poly_set_coeff_mpz(rop->num, i, x);
     354        return;
     355    }
     356   
     357    t = fmpz_poly_get_coeff_ptr(rop->num, i);
     358    canonicalize = !(t == NULL || fmpz_is_zero(t));
     359   
     360    mpz_init(z);
     361    fmpz_to_mpz(z, rop->den);
     362    mpz_mul(z, z, x);
     363    fmpz_poly_set_coeff_mpz(rop->num, i, z);
     364    if (canonicalize)
     365        fmpq_poly_canonicalize(rop, NULL);
     366    mpz_clear(z);
     367}
     368
     369/**
     370 * \ingroup  Coefficients
     371 *
     372 * Sets the <tt>i</tt>th coefficient of the polynomial <tt>op</tt> to \c x.
     373 */
     374void fmpq_poly_set_coeff_mpq(fmpq_poly_ptr rop, ulong i, const mpq_t x)
     375{
     376    fmpz_t oldc;
     377    mpz_t den, gcd;
     378    int canonicalize;
     379   
     380    mpz_init(den);
     381    mpz_init(gcd);
     382   
     383    fmpz_to_mpz(den, rop->den);
     384    mpz_gcd(gcd, den, mpq_denref(x));
     385   
     386    oldc = fmpz_poly_get_coeff_ptr(rop->num, i);
     387    canonicalize = !(oldc == NULL || fmpz_is_zero(oldc));
     388   
     389    if (mpz_cmp(mpq_denref(x), gcd) == 0)
     390    {
     391        mpz_divexact(den, den, gcd);
     392        mpz_mul(gcd, den, mpq_numref(x));
     393        fmpz_poly_set_coeff_mpz(rop->num, i, gcd);
     394    }
     395    else
     396    {
     397        mpz_t t;
     398       
     399        mpz_init(t);
     400        mpz_divexact(t, mpq_denref(x), gcd);
     401        fmpz_poly_scalar_mul_mpz(rop->num, rop->num, t);
     402       
     403        mpz_divexact(gcd, den, gcd);
     404        mpz_mul(gcd, gcd, mpq_numref(x));
     405       
     406        fmpz_poly_set_coeff_mpz(rop->num, i, gcd);
     407       
     408        mpz_mul(den, den, t);
     409        rop->den = fmpz_realloc(rop->den, mpz_size(den));
     410        mpz_to_fmpz(rop->den, den);
     411       
     412        mpz_clear(t);
     413    }
     414   
     415    if (canonicalize)
     416        fmpq_poly_canonicalize(rop, NULL);
     417   
     418    mpz_clear(den);
     419    mpz_clear(gcd);
     420}
     421
     422/**
     423 * \ingroup  Coefficients
     424 *
     425 * Sets the <tt>i</tt>th coefficient of the polynomial <tt>op</tt> to \c x.
     426 */
     427void fmpq_poly_set_coeff_si(fmpq_poly_ptr rop, ulong i, long x)
     428{
     429    mpz_t z;
     430    fmpz_t t;
     431    int canonicalize;
     432   
     433    /* Is the denominator 1?  This includes the case when rop == 0. */
     434    if (fmpz_is_one(rop->den))
     435    {
     436        fmpz_poly_set_coeff_si(rop->num, i, x);
     437        return;
     438    }
     439   
     440    t = fmpz_poly_get_coeff_ptr(rop->num, i);
     441    canonicalize = !(t == NULL || fmpz_is_zero(t));
     442   
     443    mpz_init(z);
     444    fmpz_to_mpz(z, rop->den);
     445    mpz_mul_si(z, z, x);
     446    fmpz_poly_set_coeff_mpz(rop->num, i, z);
     447    if (canonicalize)
     448        fmpq_poly_canonicalize(rop, NULL);
     449    mpz_clear(z);
     450}
     451
     452///////////////////////////////////////////////////////////////////////////////
     453// Comparison
     454
     455/**
     456 * \brief    Returns whether \c op1 and \c op2 are equal.
     457 * \ingroup  Comparison
     458 *
     459 * Returns whether the two elements \c op1 and \c op2 are equal.
     460 */
     461int fmpq_poly_equal(const fmpq_poly_ptr op1, const fmpq_poly_ptr op2)
     462{
     463    return (op1->num->length == op2->num->length)
     464        && (fmpz_equal(op1->den, op2->den))
     465        && (fmpz_poly_equal(op1->num, op2->num));
     466}
     467
     468/**
     469 * \brief   Compares the two polynomials <tt>left</tt> and <tt>right</tt>.
     470 * \ingroup Comparison
     471 *
     472 * Compares the two polnomials <tt>left</tt> and <tt>right</tt>, returning
     473 * <tt>-1</tt>, <tt>0</tt>, or <tt>1</tt> as <tt>left</tt> is less than,
     474 * equal to, or greater than <tt>right</tt>.
     475 *
     476 * The comparison is first done by degree and then, in case of a tie, by
     477 * the individual coefficients, beginning with the highest.
     478 */
     479int fmpq_poly_cmp(const fmpq_poly_ptr left, const fmpq_poly_ptr right)
     480{
     481    long degdiff, i;
     482   
     483    /* Quick check whether left and right are the same object. */
     484    if (left == right)
     485        return 0;
     486   
     487    i = fmpz_poly_degree(left->num);
     488    degdiff = i - fmpz_poly_degree(right->num);
     489   
     490    if (degdiff < 0)
     491        return -1;
     492    else if (degdiff > 0)
     493        return 1;
     494    else
     495    {
     496        int ans;
     497        ulong limbs;
     498        fmpz_t lcoeff, rcoeff;
     499       
     500        if (fmpz_equal(left->den, right->den))
     501        {
     502            if (i == -1)  /* left and right are both zero */
     503                return 0;
     504            limbs = FLINT_MAX(left->num->limbs, right->num->limbs) + 1;
     505            lcoeff = fmpz_init(limbs);
     506            while (fmpz_equal(fmpz_poly_get_coeff_ptr(left->num, i),
     507                fmpz_poly_get_coeff_ptr(right->num, i)) && i > 0)
     508                i--;
     509            fmpz_sub(lcoeff, fmpz_poly_get_coeff_ptr(left->num, i),
     510                fmpz_poly_get_coeff_ptr(right->num, i));
     511            ans = fmpz_sgn(lcoeff);
     512            fmpz_clear(lcoeff);
     513            return ans;
     514        }
     515        else if (fmpz_is_one(left->den))
     516        {
     517            limbs = FLINT_MAX(left->num->limbs + fmpz_size(right->den),
     518                right->num->limbs) + 1;
     519            lcoeff = fmpz_init(limbs);
     520            fmpz_mul(lcoeff, fmpz_poly_get_coeff_ptr(left->num, i), right->den);
     521            while (fmpz_equal(lcoeff, fmpz_poly_get_coeff_ptr(right->num, i))
     522                && i > 0)
     523            {
     524                i--;
     525                fmpz_mul(lcoeff, fmpz_poly_get_coeff_ptr(left->num, i),
     526                    right->den);
     527            }
     528            fmpz_sub(lcoeff, lcoeff, fmpz_poly_get_coeff_ptr(right->num, i));
     529            ans = fmpz_sgn(lcoeff);
     530            fmpz_clear(lcoeff);
     531            return ans;
     532        }
     533        else if (fmpz_is_one(right->den))
     534        {
     535            limbs = FLINT_MAX(left->num->limbs,
     536                right->num->limbs + fmpz_size(left->den)) + 1;
     537            rcoeff = fmpz_init(limbs);
     538            fmpz_mul(rcoeff, fmpz_poly_get_coeff_ptr(right->num, i), left->den);
     539            while (fmpz_equal(fmpz_poly_get_coeff_ptr(left->num, i), rcoeff)
     540                && i > 0)
     541            {
     542                i--;
     543                fmpz_mul(rcoeff, fmpz_poly_get_coeff_ptr(right->num, i),
     544                    left->den);
     545            }
     546            fmpz_sub(rcoeff, fmpz_poly_get_coeff_ptr(left->num, i), rcoeff);
     547            ans = fmpz_sgn(rcoeff);
     548            fmpz_clear(rcoeff);
     549            return ans;
     550        }
     551        else
     552        {
     553            limbs = FLINT_MAX(left->num->limbs + fmpz_size(right->den),
     554                right->num->limbs + fmpz_size(left->den)) + 1;
     555            lcoeff = fmpz_init(limbs);
     556            rcoeff = fmpz_init(right->num->limbs + fmpz_size(left->den));
     557            fmpz_mul(lcoeff, fmpz_poly_get_coeff_ptr(left->num, i), right->den);
     558            fmpz_mul(rcoeff, fmpz_poly_get_coeff_ptr(right->num, i), left->den);
     559            while (fmpz_equal(lcoeff, rcoeff) && i > 0)
     560            {
     561                i--;
     562                fmpz_mul(lcoeff, fmpz_poly_get_coeff_ptr(left->num, i),
     563                    right->den);
     564                fmpz_mul(rcoeff, fmpz_poly_get_coeff_ptr(right->num, i),
     565                    left->den);
     566            }
     567            fmpz_sub(lcoeff, lcoeff, rcoeff);
     568            ans = fmpz_sgn(lcoeff);
     569            fmpz_clear(lcoeff);
     570            fmpz_clear(rcoeff);
     571            return ans;
     572        }
     573    }
     574}
     575
     576///////////////////////////////////////////////////////////////////////////////
     577// Addition/ subtraction
     578
     579/**
     580 * \ingroup  Addition
     581 *
     582 * Sets \c rop to the sum of \c rop and \c op.
     583 *
     584 * \todo  This is currently implemented by creating a copy!
     585 */
     586void _fmpq_poly_add_in_place(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
     587{
     588    if (rop == op)
     589    {
     590        fmpq_poly_scalar_mul_si(rop, rop, 2l);
     591        return;
     592    }
     593   
     594    if (fmpq_poly_is_zero(rop))
     595    {
     596        fmpq_poly_set(rop, op);
     597        return;
     598    }
     599    if (fmpq_poly_is_zero(op))
     600    {
     601        return;
     602    }
     603   
     604    /* Now we may assume that rop and op refer to distinct objects in        */
     605    /* memory and that both polynomials are non-zero.                        */
     606    fmpq_poly_t t;
     607    fmpq_poly_init(t);
     608    fmpq_poly_add(t, rop, op);
     609    fmpq_poly_swap(rop, t);
     610    fmpq_poly_clear(t);
     611}
     612
     613/**
     614 * \ingroup  Addition
     615 *
     616 * Sets \c rop to the sum of \c op1 and \c op2.
     617 */
     618void fmpq_poly_add(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2)
     619{
     620    if (op1 == op2)
     621    {
     622        fmpq_poly_scalar_mul_si(rop, op1, 2l);
     623        return;
     624    }
     625   
     626    if (rop == op1)
     627    {
     628        _fmpq_poly_add_in_place(rop, op2);
     629        return;
     630    }
     631    if (rop == op2)
     632    {
     633        _fmpq_poly_add_in_place(rop, op1);
     634        return;
     635    }
     636   
     637    /* From here on, we may assume that rop, op1 and op2 all refer to        */
     638    /* distinct objects in memory, although they may still be equal.         */
     639   
     640    if (fmpz_is_one(op1->den))
     641    {
     642        if (fmpz_is_one(op2->den))
     643        {
     644            fmpz_poly_add(rop->num, op1->num, op2->num);
     645            fmpz_set_ui(rop->den, 1ul);
     646        }
     647        else
     648        {
     649            fmpz_poly_scalar_mul_fmpz(rop->num, op1->num, op2->den);
     650            fmpz_poly_add(rop->num, rop->num, op2->num);
     651            rop->den = fmpz_realloc(rop->den, fmpz_size(op2->den));
     652            fmpz_set(rop->den, op2->den);
     653        }
     654    }
     655    else
     656    {
     657        if (fmpz_is_one(op2->den))
     658        {
     659            fmpz_poly_scalar_mul_fmpz(rop->num, op2->num, op1->den);
     660            fmpz_poly_add(rop->num, op1->num, rop->num);
     661            rop->den = fmpz_realloc(rop->den, fmpz_size(op1->den));
     662            fmpz_set(rop->den, op1->den);
     663        }
     664        else
     665        {
     666            fmpz_poly_t tpoly;
     667            fmpz_t tfmpz;
     668            ulong limbs;
     669           
     670            fmpz_poly_init(tpoly);
     671           
     672            limbs = fmpz_size(op1->den) + fmpz_size(op2->den);
     673            rop->den = fmpz_realloc(rop->den, limbs);
     674           
     675            limbs = FLINT_MAX(limbs, fmpz_size(op2->den) + op1->num->limbs);
     676            limbs = FLINT_MAX(limbs, fmpz_size(op1->den) + op2->num->limbs);
     677            tfmpz = fmpz_init(limbs);
     678           
     679            fmpz_gcd(rop->den, op1->den, op2->den);
     680            fmpz_tdiv(tfmpz, op2->den, rop->den);
     681            fmpz_poly_scalar_mul_fmpz(rop->num, op1->num, tfmpz);
     682            fmpz_tdiv(tfmpz, op1->den, rop->den);
     683            fmpz_poly_scalar_mul_fmpz(tpoly, op2->num, tfmpz);
     684            fmpz_poly_add(rop->num, rop->num, tpoly);
     685            fmpz_mul(rop->den, tfmpz, op2->den);
     686           
     687            fmpq_poly_canonicalize(rop, tfmpz);
     688           
     689            fmpz_poly_clear(tpoly);
     690            fmpz_clear(tfmpz);
     691        }
     692    }
     693}
     694
     695/**
     696 * \ingroup  Addition
     697 *
     698 * Sets \c rop to the difference of \c rop and \c op.
     699 *
     700 * \note  This is implemented using the methods #fmpq_poly_neg() and
     701 *        #_fmpq_poly_add_in_place().
     702 */
     703void _fmpq_poly_sub_in_place(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
     704{
     705    if (rop == op)
     706    {
     707        fmpq_poly_zero(rop);
     708        return;
     709    }
     710   
     711    fmpq_poly_neg(rop, rop);
     712    _fmpq_poly_add_in_place(rop, op);
     713    fmpq_poly_neg(rop, rop);
     714}
     715
     716/**
     717 * \ingroup  Addition
     718 *
     719 * Sets \c rop to the difference of \c op1 and \c op2.
     720 */
     721void fmpq_poly_sub(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2)
     722{
     723    if (op1 == op2)
     724    {
     725        fmpq_poly_zero(rop);
     726        return;
     727    }
     728    if (rop == op1)
     729    {
     730        _fmpq_poly_sub_in_place(rop, op2);
     731        return;
     732    }
     733    if (rop == op2)
     734    {
     735        _fmpq_poly_sub_in_place(rop, op1);
     736        fmpq_poly_neg(rop, rop);
     737        return;
     738    }
     739   
     740    /* From here on, we know that rop, op1 and op2 refer to distinct objects */
     741    /* in memory, although as rational functions they may still be equal     */
     742   
     743    if (fmpz_is_one(op1->den))
     744    {
     745        if (fmpz_is_one(op2->den))
     746        {
     747            fmpz_poly_sub(rop->num, op1->num, op2->num);
     748            fmpz_set_ui(rop->den, 1ul);
     749        }
     750        else
     751        {
     752            fmpz_poly_scalar_mul_fmpz(rop->num, op1->num, op2->den);
     753            fmpz_poly_sub(rop->num, rop->num, op2->num);
     754            rop->den = fmpz_realloc(rop->den, fmpz_size(op2->den));
     755            fmpz_set(rop->den, op2->den);
     756        }
     757    }
     758    else
     759    {
     760        if (fmpz_is_one(op2->den))
     761        {
     762            fmpz_poly_scalar_mul_fmpz(rop->num, op2->num, op1->den);
     763            fmpz_poly_sub(rop->num, op1->num, rop->num);
     764            rop->den = fmpz_realloc(rop->den, fmpz_size(op1->den));
     765            fmpz_set(rop->den, op1->den);
     766        }
     767        else
     768        {
     769            fmpz_poly_t tpoly;
     770            fmpz_t tfmpz;
     771            ulong limbs;
     772           
     773            fmpz_poly_init(tpoly);
     774           
     775            limbs = fmpz_size(op1->den) + fmpz_size(op2->den);
     776            rop->den = fmpz_realloc(rop->den, limbs);
     777           
     778            limbs = FLINT_MAX(limbs, fmpz_size(op2->den) + op1->num->limbs);
     779            limbs = FLINT_MAX(limbs, fmpz_size(op1->den) + op2->num->limbs);
     780            tfmpz = fmpz_init(limbs);
     781           
     782            fmpz_gcd(rop->den, op1->den, op2->den);
     783            fmpz_tdiv(tfmpz, op2->den, rop->den);
     784            fmpz_poly_scalar_mul_fmpz(rop->num, op1->num, tfmpz);
     785            fmpz_tdiv(tfmpz, op1->den, rop->den);
     786            fmpz_poly_scalar_mul_fmpz(tpoly, op2->num, tfmpz);
     787            fmpz_poly_sub(rop->num, rop->num, tpoly);
     788            fmpz_mul(rop->den, tfmpz, op2->den);
     789           
     790            fmpq_poly_canonicalize(rop, tfmpz);
     791           
     792            fmpz_poly_clear(tpoly);
     793            fmpz_clear(tfmpz);
     794        }
     795    }
     796}
     797
     798/**
     799 * \ingroup  Addition
     800 *
     801 * Sets \c rop to <tt>rop + op1 * op2</tt>.
     802 *
     803 * Currently, this method refers to the methods #fmpq_poly_mul() and
     804 * #fmpq_poly_add() to form the result in the naive way.
     805 *
     806 * \todo  Implement this method more efficiently.
     807 */
     808void fmpq_poly_addmul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2)
     809{
     810    fmpq_poly_t t;
     811    fmpq_poly_init(t);
     812    fmpq_poly_mul(t, op1, op2);
     813    fmpq_poly_add(rop, rop, t);
     814    fmpq_poly_clear(t);
     815}
     816
     817/**
     818 * \ingroup  Addition
     819 *
     820 * Sets \c rop to <tt>rop - op1 * op2</tt>.
     821 *
     822 * Currently, this method refers to the methods #fmpq_poly_mul() and
     823 * #fmpq_poly_sub() to form the result in the naive way.
     824 *
     825 * \todo  Implement this method more efficiently.
     826 */
     827void fmpq_poly_submul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2)
     828{
     829    fmpq_poly_t t;
     830    fmpq_poly_init(t);
     831    fmpq_poly_mul(t, op1, op2);
     832    fmpq_poly_sub(rop, rop, t);
     833    fmpq_poly_clear(t);
     834}
     835
     836///////////////////////////////////////////////////////////////////////////////
     837// Scalar multiplication and division
     838
     839/**
     840 * \ingroup  ScalarMul
     841 *
     842 * Sets \c rop to the scalar product of \c op and the integer \c x.
     843 */
     844void fmpq_poly_scalar_mul_si(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long x)
     845{
     846    if (fmpz_is_one(op->den))
     847    {
     848        fmpz_poly_scalar_mul_si(rop->num, op->num, x);
     849        fmpz_set_ui(rop->den, 1ul);
     850    }
     851    else
     852    {
     853        fmpz_t fx, g;
     854       
     855        g  = fmpz_init(fmpz_size(op->den));
     856        fx = fmpz_init(1);
     857        fmpz_set_si(fx, x);
     858        fmpz_gcd(g, op->den, fx);
     859        fmpz_abs(g, g);
     860        if (fmpz_is_one(g))
     861        {
     862            fmpz_poly_scalar_mul_si(rop->num, op->num, x);
     863            rop->den = fmpz_realloc(rop->den, fmpz_size(op->den));
     864            fmpz_set(rop->den, op->den);
     865        }
     866        else
     867        {
     868            if (rop == op)
     869            {
     870                fmpz_t t;
     871                t = fmpz_init(fmpz_size(op->den));
     872                fmpz_tdiv(t, fx, g);
     873                fmpz_poly_scalar_mul_fmpz(rop->num, op->num, t);
     874                fmpz_tdiv(t, op->den, g);
     875                fmpz_clear(rop->den);  // FIXME  Fixed
     876                rop->den = t;
     877            }
     878            else
     879            {
     880                rop->den = fmpz_realloc(rop->den, fmpz_size(op->den));
     881                fmpz_tdiv(rop->den, fx, g);
     882                fmpz_poly_scalar_mul_fmpz(rop->num, op->num, rop->den);
     883                fmpz_tdiv(rop->den, op->den, g);
     884            }
     885        }
     886        fmpz_clear(g);
     887        fmpz_clear(fx);
     888    }
     889}
     890
     891/**
     892 * \ingroup  ScalarMul
     893 *
     894 * Sets \c rop to the scalar multiple of \c op with the \c mpz_t \c x.
     895 */
     896void fmpq_poly_scalar_mul_mpz(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpz_t x)
     897{
     898    fmpz_t g, y;
     899
     900    if (fmpz_is_one(op->den))
     901    {
     902        fmpz_poly_scalar_mul_mpz(rop->num, op->num, x);
     903        fmpz_set_ui(rop->den, 1UL);
     904        return;
     905    }
     906   
     907    if (mpz_cmpabs_ui(x, 1) == 0)
     908    {
     909        if (mpz_sgn(x) > 0)
     910            fmpq_poly_set(rop, op);
     911        else
     912            fmpq_poly_neg(rop, op);
     913        return;
     914    }
     915    if (mpz_sgn(x) == 0)
     916    {
     917        fmpq_poly_zero(rop);
     918        return;
     919    }
     920   
     921    g = fmpz_init(FLINT_MAX(fmpz_size(op->den), mpz_size(x)));
     922    y = fmpz_init(mpz_size(x));
     923    mpz_to_fmpz(y, x);
     924   
     925    fmpz_gcd(g, op->den, y);
     926   
     927    if (fmpz_is_one(g))
     928    {
     929        fmpz_poly_scalar_mul_fmpz(rop->num, op->num, y);
     930        if (rop != op)
     931        {
     932            rop->den = fmpz_realloc(rop->den, fmpz_size(op->den));
     933            fmpz_set(rop->den, op->den);
     934        }
     935    }
     936    else
     937    {
     938        fmpz_t t;
     939        t = fmpz_init(FLINT_MAX(fmpz_size(y), fmpz_size(op->den)));
     940        fmpz_divides(t, y, g);
     941        fmpz_poly_scalar_mul_fmpz(rop->num, op->num, t);
     942        fmpz_divides(t, op->den, g);
     943        fmpz_clear(rop->den);
     944        rop->den = t;
     945    }
     946   
     947    fmpz_clear(g);
     948    fmpz_clear(y);
     949}
     950
     951/**
     952 * \ingroup  ScalarMul
     953 *
     954 * Sets \c rop to the scalar multiple of \c op with the \c mpq_t \c x.
     955 */
     956void fmpq_poly_scalar_mul_mpq(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x)
     957{
     958    fmpz_poly_scalar_mul_mpz(rop->num, op->num, mpq_numref(x));
     959    if (fmpz_is_one(op->den))
     960    {
     961        rop->den = fmpz_realloc(rop->den, mpz_size(mpq_denref(x)));
     962        mpz_to_fmpz(rop->den, mpq_denref(x));
     963    }
     964    else
     965    {
     966        fmpz_t s, t;
     967        s = fmpz_init(mpz_size(mpq_denref(x)));
     968        t = fmpz_init(fmpz_size(op->den));
     969        mpz_to_fmpz(s, mpq_denref(x));
     970        fmpz_set(t, op->den);
     971        rop->den = fmpz_realloc(rop->den, fmpz_size(s) + fmpz_size(t));
     972        fmpz_mul(rop->den, s, t);
     973        fmpz_clear(s);
     974        fmpz_clear(t);
     975    }
     976    fmpq_poly_canonicalize(rop, NULL);
     977}
     978
     979/**
     980 * /ingroup  ScalarMul
     981 *
     982 * Divides \c rop by the integer \c x. 
     983 *
     984 * Assumes that \c x is non-zero.  Otherwise, an exception is raised in the
     985 * form of an <tt>abort</tt> statement.
     986 */
     987void _fmpq_poly_scalar_div_si_in_place(fmpq_poly_ptr rop, long x)
     988{
     989    fmpz_t cont, fx, gcd;
     990   
     991    /* Assertion! */
     992    if (x == 0l)
     993    {
     994        printf("ERROR (_fmpq_poly_scalar_div_si_in_place).  Division by zero.\n");
     995        abort();
     996    }
     997   
     998    if (x == 1l)
     999    {
     1000        return;
     1001    }
     1002   
     1003    cont = fmpz_init(rop->num->limbs);
     1004    fmpz_poly_content(cont, rop->num);
     1005    fmpz_abs(cont, cont);
     1006   
     1007    if (fmpz_is_one(cont))
     1008    {
     1009        if (x > 0l)
     1010        {
     1011            if (fmpz_is_one(rop->den))
     1012            {
     1013                fmpz_set_si(rop->den, x);
     1014            }
     1015            else
     1016            {
     1017                fmpz_t t;
     1018                t = fmpz_init(fmpz_size(rop->den) + 1);
     1019                fmpz_mul_ui(t, rop->den, (ulong) x);
     1020                fmpz_clear(rop->den);  // FIXME  Fixed
     1021                rop->den = t;
     1022            }
     1023        }
     1024        else
     1025        {
     1026            fmpz_poly_neg(rop->num, rop->num);
     1027            if (fmpz_is_one(rop->den))
     1028            {
     1029                fmpz_set_si(rop->den, -x);
     1030            }
     1031            else
     1032            {
     1033                fmpz_t t;
     1034                t = fmpz_init(fmpz_size(rop->den) + 1);
     1035                fmpz_mul_ui(t, rop->den, (ulong) -x);
     1036                fmpz_clear(rop->den);  // FIXME  Fixed
     1037                rop->den = t;
     1038            }
     1039        }
     1040        fmpz_clear(cont);
     1041        return;
     1042    }
     1043   
     1044    fx = fmpz_init(1);
     1045    fmpz_set_si(fx, x);
     1046   
     1047    gcd = fmpz_init(FLINT_MAX(rop->num->limbs, fmpz_size(fx)));
     1048    fmpz_gcd(gcd, cont, fx);
     1049    fmpz_abs(gcd, gcd);
     1050   
     1051    if (fmpz_is_one(gcd))
     1052    {
     1053        if (x > 0l)
     1054        {
     1055            if (fmpz_is_one(rop->den))
     1056            {
     1057                fmpz_set_si(rop->den, x);
     1058            }
     1059            else
     1060            {
     1061                fmpz_t t;
     1062                t = fmpz_init(fmpz_size(rop->den) + 1);
     1063                fmpz_mul_ui(t, rop->den, (ulong) x);
     1064                fmpz_clear(rop->den);  // FIXME  Fixed
     1065                rop->den = t;
     1066            }
     1067        }
     1068        else
     1069        {
     1070            fmpz_poly_neg(rop->num, rop->num);
     1071            if (fmpz_is_one(rop->den))
     1072            {
     1073                fmpz_set_si(rop->den, -x);
     1074            }
     1075            else
     1076            {
     1077                fmpz_t t;
     1078                t = fmpz_init(fmpz_size(rop->den) + 1);
     1079                fmpz_mul_ui(t, rop->den, (ulong) -x);
     1080                fmpz_clear(rop->den);  // FIXME  Fixed
     1081                rop->den = t;
     1082            }
     1083        }
     1084    }
     1085    else
     1086    {
     1087        fmpz_poly_scalar_div_fmpz(rop->num, rop->num, gcd);
     1088        if (fmpz_is_one(rop->den))
     1089        {
     1090            fmpz_tdiv(rop->den, fx, gcd);
     1091        }
     1092        else
     1093        {
     1094            fmpz_tdiv(cont, fx, gcd);
     1095            fx = fmpz_realloc(fx, fmpz_size(rop->den));
     1096            fmpz_set(fx, rop->den);
     1097            rop->den = fmpz_realloc(rop->den, fmpz_size(rop->den) + 1);
     1098            fmpz_mul(rop->den, fx, cont);
     1099        }
     1100        if (x < 0l)
     1101        {
     1102            fmpz_poly_neg(rop->num, rop->num);
     1103            fmpz_neg(rop->den, rop->den);
     1104        }
     1105    }
     1106   
     1107    fmpz_clear(cont);
     1108    fmpz_clear(fx);
     1109    fmpz_clear(gcd);
     1110}
     1111
     1112
     1113/**
     1114 * \ingroup  ScalarMul
     1115 *
     1116 * Sets \c rop to the scalar multiple of \c op with the multiplicative inverse
     1117 * of the integer \c x.
     1118 *
     1119 * Assumes that \c x is non-zero.  Otherwise, an exception is raised in the
     1120 * form of an <tt>abort</tt> statement.
     1121 */
     1122void fmpq_poly_scalar_div_si(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long x)
     1123{
     1124    fmpz_t cont, fx, gcd;
     1125   
     1126    /* Assertion! */
     1127    if (x == 0l)
     1128    {
     1129        printf("ERROR (fmpq_poly_scalar_div_si).  Division by zero.\n");
     1130        abort();
     1131    }
     1132   
     1133    if (rop == op)
     1134    {
     1135        _fmpq_poly_scalar_div_si_in_place(rop, x);
     1136        return;
     1137    }
     1138   
     1139    /* From here on, we may assume that rop and op denote two different      */
     1140    /* rational polynomials (as objects in memory).                          */
     1141   
     1142    if (x == 1l)
     1143    {
     1144        fmpq_poly_set(rop, op);
     1145        return;
     1146    }
     1147   
     1148    cont = fmpz_init(op->num->limbs);
     1149    fmpz_poly_content(cont, op->num);
     1150    fmpz_abs(cont, cont);
     1151   
     1152    if (fmpz_is_one(cont))
     1153    {
     1154        if (x > 0l)
     1155        {
     1156            fmpz_poly_set(rop->num, op->num);
     1157            if (fmpz_is_one(op->den))
     1158            {
     1159                fmpz_set_si(rop->den, x);
     1160            }
     1161            else
     1162            {
     1163                rop->den = fmpz_realloc(rop->den, fmpz_size(op->den) + 1);
     1164                fmpz_mul_ui(rop->den, op->den, (ulong) x);
     1165            }
     1166        }
     1167        else
     1168        {
     1169            fmpz_poly_neg(rop->num, op->num);
     1170            if (fmpz_is_one(op->den))
     1171            {
     1172                fmpz_set_si(rop->den, -x);
     1173            }
     1174            else
     1175            {
     1176                rop->den = fmpz_realloc(rop->den, fmpz_size(op->den) + 1);
     1177                fmpz_mul_ui(rop->den, op->den, (ulong) -x);
     1178            }
     1179        }
     1180        fmpz_clear(cont);
     1181        return;
     1182    }
     1183   
     1184    fx = fmpz_init(1);
     1185    fmpz_set_si(fx, x);
     1186   
     1187    gcd = fmpz_init(FLINT_MAX(op->num->limbs, fmpz_size(fx)));
     1188    fmpz_gcd(gcd, cont, fx);
     1189    fmpz_abs(gcd, gcd);
     1190   
     1191    if (fmpz_is_one(gcd))
     1192    {
     1193        if (x > 0l)
     1194        {
     1195            fmpz_poly_set(rop->num, op->num);
     1196            if (fmpz_is_one(op->den))
     1197            {
     1198                fmpz_set_si(rop->den, x);
     1199            }
     1200            else
     1201            {
     1202                rop->den = fmpz_realloc(rop->den, fmpz_size(op->den) + 1);
     1203                fmpz_mul_ui(rop->den, op->den, (ulong) x);
     1204            }
     1205        }
     1206        else
     1207        {
     1208            fmpz_poly_neg(rop->num, op->num);
     1209            if (fmpz_is_one(op->den))
     1210            {
     1211                fmpz_set_si(rop->den, -x);
     1212            }
     1213            else
     1214            {
     1215                rop->den = fmpz_realloc(rop->den, fmpz_size(op->den) + 1);
     1216                fmpz_mul_ui(rop->den, op->den, (ulong) -x);
     1217            }
     1218        }
     1219    }
     1220    else
     1221    {
     1222        fmpz_poly_scalar_div_fmpz(rop->num, op->num, gcd);
     1223        if (fmpz_is_one(op->den))
     1224        {
     1225            fmpz_tdiv(rop->den, fx, gcd);
     1226        }
     1227        else
     1228        {
     1229            rop->den = fmpz_realloc(rop->den, fmpz_size(op->den) + 1);
     1230            fmpz_tdiv(cont, fx, gcd);  /* fx and gcd are word-sized */
     1231            fmpz_mul(rop->den, op->den, cont);
     1232        }
     1233        if (x < 0l)
     1234        {
     1235            fmpz_poly_neg(rop->num, rop->num);
     1236            fmpz_neg(rop->den, rop->den);
     1237        }
     1238    }
     1239   
     1240    fmpz_clear(cont);
     1241    fmpz_clear(fx);
     1242    fmpz_clear(gcd);
     1243}
     1244
     1245/**
     1246 * \ingroup  ScalarMul
     1247 *
     1248 * Sets \c rop to the scalar multiple of \c op with the multiplicative inverse
     1249 * of the integer \c x.
     1250 *
     1251 * Assumes that \c x is non-zero.  Otherwise, an exception is raised in the
     1252 * form of an <tt>abort</tt> statement.
     1253 */
     1254void _fmpq_poly_scalar_div_mpz_in_place(fmpq_poly_ptr rop, const mpz_t x)
     1255{
     1256    fmpz_t cont, fx, gcd;
     1257   
     1258    /* Assertion! */
     1259    if (mpz_sgn(x) == 0)
     1260    {
     1261        printf("ERROR (_fmpq_poly_scalar_div_mpz_in_place).  Division by zero.\n");
     1262        abort();
     1263    }
     1264   
     1265    if (mpz_cmp_ui(x, 1) == 0)
     1266        return;
     1267   
     1268    cont = fmpz_init(rop->num->limbs);
     1269    fmpz_poly_content(cont, rop->num);
     1270    fmpz_abs(cont, cont);
     1271   
     1272    if (fmpz_is_one(cont))
     1273    {
     1274        if (fmpz_is_one(rop->den))
     1275        {
     1276            rop->den = fmpz_realloc(rop->den, mpz_size(x));
     1277            mpz_to_fmpz(rop->den, x);
     1278        }
     1279        else
     1280        {
     1281            fmpz_t t;
     1282            fx = fmpz_init(mpz_size(x));
     1283            mpz_to_fmpz(fx, x);
     1284            t = fmpz_init(fmpz_size(rop->den) + mpz_size(x));
     1285            fmpz_mul(t, rop->den, fx);
     1286            fmpz_clear(rop->den);  // FIXME  Fixed
     1287            rop->den = t;
     1288            fmpz_clear(fx);
     1289        }
     1290        if (mpz_sgn(x) < 0)
     1291        {
     1292            fmpz_poly_neg(rop->num, rop->num);
     1293            fmpz_neg(rop->den, rop->den);
     1294        }
     1295        fmpz_clear(cont);
     1296        return;
     1297    }
     1298   
     1299    fx = fmpz_init(mpz_size(x));
     1300    mpz_to_fmpz(fx, x);
     1301   
     1302    gcd = fmpz_init(FLINT_MAX(rop->num->limbs, fmpz_size(fx)));
     1303    fmpz_gcd(gcd, cont, fx);
     1304    fmpz_abs(gcd, gcd);
     1305   
     1306    if (fmpz_is_one(gcd))
     1307    {
     1308        if (fmpz_is_one(rop->den))
     1309        {
     1310            rop->den = fmpz_realloc(rop->den, mpz_size(x));
     1311            mpz_to_fmpz(rop->den, x);
     1312        }
     1313        else
     1314        {
     1315            fmpz_t t;
     1316            t = fmpz_init(fmpz_size(rop->den) + mpz_size(x));
     1317            fmpz_mul(t, rop->den, fx);
     1318            fmpz_clear(rop->den);  // FIXME  Fixed
     1319            rop->den = t;
     1320        }
     1321        if (mpz_sgn(x) < 0)
     1322        {
     1323            fmpz_poly_neg(rop->num, rop->num);
     1324            fmpz_neg(rop->den, rop->den);
     1325        }
     1326    }
     1327    else
     1328    {
     1329        fmpz_poly_scalar_div_fmpz(rop->num, rop->num, gcd);
     1330        if (fmpz_is_one(rop->den))
     1331        {
     1332            rop->den = fmpz_realloc(rop->den, fmpz_size(fx));
     1333            fmpz_tdiv(rop->den, fx, gcd);
     1334        }
     1335        else
     1336        {
     1337            cont = fmpz_realloc(cont, fmpz_size(fx));
     1338            fmpz_tdiv(cont, fx, gcd);
     1339            fx = fmpz_realloc(fx, fmpz_size(rop->den));
     1340            fmpz_set(fx, rop->den);
     1341            rop->den = fmpz_realloc(rop->den, fmpz_size(fx) + fmpz_size(cont));
     1342            fmpz_mul(rop->den, fx, cont);
     1343        }
     1344        if (mpz_sgn(x) < 0)
     1345        {
     1346            fmpz_poly_neg(rop->num, rop->num);
     1347            fmpz_neg(rop->den, rop->den);
     1348        }
     1349    }
     1350   
     1351    fmpz_clear(cont);
     1352    fmpz_clear(fx);
     1353    fmpz_clear(gcd);
     1354}
     1355
     1356/**
     1357 * \ingroup  ScalarMul
     1358 *
     1359 * Sets \c rop to the scalar multiple of \c op with the multiplicative inverse
     1360 * of the integer \c x.
     1361 *
     1362 * Assumes that \c x is non-zero.  Otherwise, an exception is raised in the
     1363 * form of an <tt>abort</tt> statement.
     1364 */
     1365void fmpq_poly_scalar_div_mpz(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpz_t x)
     1366{
     1367    fmpz_t cont, fx, gcd;
     1368   
     1369    /* Assertion! */
     1370    if (mpz_sgn(x) == 0)
     1371    {
     1372        printf("ERROR (fmpq_poly_scalar_div_mpz).  Division by zero.\n");
     1373        abort();
     1374    }
     1375   
     1376    if (rop == op)
     1377    {
     1378        _fmpq_poly_scalar_div_mpz_in_place(rop, x);
     1379        return;
     1380    }
     1381   
     1382    /* From here on, we may assume that rop and op denote two different      */
     1383    /* rational polynomials (as objects in memory).                          */
     1384   
     1385    if (mpz_cmp_ui(x, 1) == 0)
     1386    {
     1387        fmpq_poly_set(rop, op);
     1388        return;
     1389    }
     1390   
     1391    cont = fmpz_init(op->num->limbs);
     1392    fmpz_poly_content(cont, op->num);
     1393    fmpz_abs(cont, cont);
     1394   
     1395    if (fmpz_is_one(cont))
     1396    {
     1397        if (fmpz_is_one(op->den))
     1398        {
     1399            rop->den = fmpz_realloc(rop->den, mpz_size(x));
     1400            mpz_to_fmpz(rop->den, x);
     1401        }
     1402        else
     1403        {
     1404            fmpz_t t;
     1405            t = fmpz_init(mpz_size(x));
     1406            mpz_to_fmpz(t, x);
     1407            rop->den = fmpz_realloc(rop->den, fmpz_size(op->den) + fmpz_size(t));
     1408            fmpz_mul(rop->den, op->den, t);
     1409            fmpz_clear(t);
     1410        }
     1411        if (mpz_sgn(x) > 0)
     1412        {
     1413            fmpz_poly_set(rop->num, op->num);
     1414        }
     1415        else
     1416        {
     1417            fmpz_poly_neg(rop->num, op->num);
     1418            fmpz_neg(rop->den, rop->den);
     1419        }
     1420        fmpz_clear(cont);
     1421        return;
     1422    }
     1423   
     1424    fx = fmpz_init(mpz_size(x));
     1425    mpz_to_fmpz(fx, x);
     1426   
     1427    gcd = fmpz_init(FLINT_MAX(op->num->limbs, fmpz_size(fx)));
     1428    fmpz_gcd(gcd, cont, fx);
     1429    fmpz_abs(gcd, gcd);
     1430   
     1431    if (fmpz_is_one(gcd))
     1432    {
     1433        if (fmpz_is_one(op->den))
     1434        {
     1435            rop->den = fmpz_realloc(rop->den, mpz_size(x));
     1436            mpz_to_fmpz(rop->den, x);
     1437        }
     1438        else
     1439        {
     1440            fmpz_t t;
     1441            t = fmpz_init(mpz_size(x));
     1442            mpz_to_fmpz(t, x);
     1443            rop->den = fmpz_realloc(rop->den, fmpz_size(op->den) + fmpz_size(t));
     1444            fmpz_mul(rop->den, op->den, t);
     1445            fmpz_clear(t);
     1446        }
     1447        if (mpz_sgn(x) > 0)
     1448        {
     1449            fmpz_poly_set(rop->num, op->num);
     1450        }
     1451        else
     1452        {
     1453            fmpz_poly_neg(rop->num, op->num);
     1454            fmpz_neg(rop->den, rop->den);
     1455        }
     1456    }
     1457    else
     1458    {
     1459        fmpz_poly_scalar_div_fmpz(rop->num, op->num, gcd);
     1460        if (fmpz_is_one(op->den))
     1461        {
     1462            rop->den = fmpz_realloc(rop->den, fmpz_size(fx));
     1463            fmpz_tdiv(rop->den, fx, gcd);
     1464        }
     1465        else
     1466        {
     1467            cont = fmpz_realloc(cont, fmpz_size(fx));
     1468            fmpz_tdiv(cont, fx, gcd);
     1469            rop->den = fmpz_realloc(rop->den, fmpz_size(op->den) + fmpz_size(cont));
     1470            fmpz_mul(rop->den, op->den, cont);
     1471        }
     1472        if (mpz_sgn(x) < 0)
     1473        {
     1474            fmpz_poly_neg(rop->num, rop->num);
     1475            fmpz_neg(rop->den, rop->den);
     1476        }
     1477    }
     1478   
     1479    fmpz_clear(cont);
     1480    fmpz_clear(fx);
     1481    fmpz_clear(gcd);
     1482}
     1483
     1484/**
     1485 * \ingroup  ScalarMul
     1486 *
     1487 * Sets \c rop to the scalar multiple of \c op with the multiplicative inverse
     1488 * of the rational \c x.
     1489 *
     1490 * Assumes that the rational \c x is in lowest terms and non-zero.  If the
     1491 * rational is not in lowest terms, the resulting value of <tt>rop</tt> is
     1492 * undefined.  If <tt>x</tt> is zero, an exception is raised in the form
     1493 * of an <tt>abort</tt> statement.
     1494 */
     1495void fmpq_poly_scalar_div_mpq(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x)
     1496{
     1497    /* Assertion! */
     1498    if (mpz_sgn(mpq_numref(x)) == 0)
     1499    {
     1500        printf("ERROR (fmpq_poly_scalar_div_mpq).  Division by zero.\n");
     1501        abort();
     1502    }
     1503   
     1504    fmpz_poly_scalar_mul_mpz(rop->num, op->num, mpq_denref(x));
     1505    if (fmpz_is_one(op->den))
     1506    {
     1507        rop->den = fmpz_realloc(rop->den, mpz_size(mpq_numref(x)));
     1508        mpz_to_fmpz(rop->den, mpq_numref(x));
     1509    }
     1510    else
     1511    {
     1512        fmpz_t s, t;
     1513        s = fmpz_init(mpz_size(mpq_numref(x)));
     1514        t = fmpz_init(fmpz_size(op->den));
     1515        mpz_to_fmpz(s, mpq_numref(x));
     1516        fmpz_set(t, op->den);
     1517        rop->den = fmpz_realloc(rop->den, fmpz_size(s) + fmpz_size(t));
     1518        fmpz_mul(rop->den, s, t);
     1519        fmpz_clear(s);
     1520        fmpz_clear(t);
     1521    }
     1522    fmpq_poly_canonicalize(rop, NULL);
     1523}
     1524
     1525///////////////////////////////////////////////////////////////////////////////
     1526// Multiplication
     1527
     1528/**
     1529 * \ingroup  Multiplication
     1530 *
     1531 * Multiplies <tt>rop</tt> by <tt>op</tt>.
     1532 */
     1533void _fmpq_poly_mul_in_place(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
     1534{
     1535    fmpq_poly_t t;
     1536   
     1537    if (rop == op)
     1538    {
     1539        fmpz_poly_power(rop->num, op->num, 2ul);
     1540        if (!fmpz_is_one(rop->den))
     1541        {
     1542            rop->den = fmpz_realloc(rop->den, 2 * fmpz_size(rop->den));
     1543            fmpz_pow_ui(rop->den, op->den, 2ul);
     1544        }
     1545        return;
     1546    }
     1547   
     1548    if (fmpq_poly_is_zero(rop) || fmpq_poly_is_zero(op))
     1549    {
     1550        fmpq_poly_zero(rop);
     1551        return;
     1552    }
     1553   
     1554    /* From here on, rop and op point to two different objects in memory,    */
     1555    /* and these are both non-zero rational polynomials.                     */
     1556   
     1557    fmpq_poly_init(t);
     1558    fmpq_poly_mul(t, rop, op);
     1559    fmpq_poly_swap(rop, t);
     1560    fmpq_poly_clear(t);
     1561}
     1562
     1563/**
     1564 * \ingroup  Multiplication
     1565 *
     1566 * Sets \c rop to the product of \c op1 and \c op2.
     1567 */
     1568void fmpq_poly_mul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2)
     1569{
     1570    fmpz_t gcd1, gcd2;
     1571   
     1572    if (op1 == op2)
     1573    {
     1574        fmpz_poly_power(rop->num, op1->num, 2ul);
     1575        if (fmpz_is_one(op1->den))
     1576        {
     1577            fmpz_set_ui(rop->den, 1);
     1578        }
     1579        else
     1580        {
     1581            if (rop == op1)
     1582            {
     1583                fmpz_t t;
     1584                t = fmpz_init(2 * fmpz_size(op1->den));
     1585                fmpz_pow_ui(t, op1->den, 2ul);
     1586                fmpz_clear(rop->den);  // FIXME  Fixed
     1587                rop->den = t;
     1588            }
     1589            else
     1590            {
     1591                rop->den = fmpz_realloc(rop->den, 2 * fmpz_size(op1->den));
     1592                fmpz_pow_ui(rop->den, op1->den, 2ul);
     1593            }
     1594        }
     1595        return;
     1596    }
     1597   
     1598    if (rop == op1)
     1599    {
     1600        _fmpq_poly_mul_in_place(rop, op2);
     1601        return;
     1602    }
     1603    if (rop == op2)
     1604    {
     1605        _fmpq_poly_mul_in_place(rop, op1);
     1606        return;
     1607    }
     1608   
     1609    /* From here on, we may assume that rop, op1 and op2 all refer to        */
     1610    /* distinct objects in memory, although they may still be equal.         */
     1611
     1612    gcd1 = NULL;
     1613    gcd2 = NULL;
     1614
     1615    if (!fmpz_is_one(op2->den))
     1616    {
     1617        gcd1 = fmpz_init(FLINT_MAX(op1->num->limbs, fmpz_size(op2->den)));
     1618        fmpz_poly_content(gcd1, op1->num);
     1619        if (!fmpz_is_one(gcd1))
     1620        {
     1621            fmpz_t t;
     1622            t = fmpz_init(FLINT_MAX(fmpz_size(gcd1), fmpz_size(op2->den)));
     1623            fmpz_gcd(t, gcd1, op2->den);
     1624            fmpz_clear(gcd1);
     1625            gcd1 = t;
     1626        }
     1627    }
     1628    if (!fmpz_is_one(op1->den))
     1629    {
     1630        gcd2 = fmpz_init(FLINT_MAX(op2->num->limbs, fmpz_size(op1->den)));
     1631        fmpz_poly_content(gcd2, op2->num);
     1632        if (!fmpz_is_one(gcd2))
     1633        {
     1634            fmpz_t t;
     1635            t = fmpz_init(FLINT_MAX(fmpz_size(gcd2), fmpz_size(op1->den)));
     1636            fmpz_gcd(t, gcd2, op1->den);
     1637            fmpz_clear(gcd2);
     1638            gcd2 = t;
     1639        }
     1640    }
     1641   
     1642    /* TODO:  If gcd1 and gcd2 are very large compared to the degrees of  */
     1643    /* poly1 and poly2, we might want to create copies of the polynomials */
     1644    /* and divide out the common factors *before* the multiplication.     */
     1645   
     1646    fmpz_poly_mul(rop->num, op1->num, op2->num);
     1647    rop->den = fmpz_realloc(rop->den, fmpz_size(op1->den)
     1648             + fmpz_size(op2->den));
     1649    fmpz_mul(rop->den, op1->den, op2->den);
     1650   
     1651    if (gcd1 == NULL || fmpz_is_one(gcd1))
     1652    {
     1653        if (gcd2 != NULL && !fmpz_is_one(gcd2))
     1654        {
     1655            fmpz_t h;
     1656            h = fmpz_init(fmpz_size(rop->den));
     1657            fmpz_poly_scalar_div_fmpz(rop->num, rop->num, gcd2);
     1658            fmpz_divides(h, rop->den, gcd2);
     1659            fmpz_clear(rop->den);
     1660            rop->den = h;
     1661        }
     1662    }
     1663    else
     1664    {
     1665        if (gcd2 == NULL || fmpz_is_one(gcd2))
     1666        {
     1667            fmpz_t h;
     1668            h = fmpz_init(fmpz_size(rop->den));
     1669            fmpz_poly_scalar_div_fmpz(rop->num, rop->num, gcd1);
     1670            fmpz_divides(h, rop->den, gcd1);
     1671            fmpz_clear(rop->den);
     1672            rop->den = h;
     1673        }
     1674        else
     1675        {
     1676            fmpz_t g, h;
     1677            g = fmpz_init(fmpz_size(gcd1) + fmpz_size(gcd2));
     1678            h = fmpz_init(fmpz_size(rop->den));
     1679            fmpz_mul(g, gcd1, gcd2);
     1680            fmpz_poly_scalar_div_fmpz(rop->num, rop->num, g);
     1681            fmpz_divides(h, rop->den, g);
     1682            fmpz_clear(rop->den);
     1683            rop->den = h;
     1684            fmpz_clear(g);
     1685        }
     1686    }
     1687   
     1688    if (gcd1 != NULL) fmpz_clear(gcd1);
     1689    if (gcd2 != NULL) fmpz_clear(gcd2);
     1690}
     1691
     1692///////////////////////////////////////////////////////////////////////////////
     1693// Division
     1694
     1695/**
     1696 * \ingroup  Division
     1697 *
     1698 * Returns the quotient of the Euclidean division of \c a by \c b.
     1699 *
     1700 * Assumes that \c b is non-zero.  Otherwise, an exception is raised in the
     1701 * form of an <tt>abort</tt> statement.
     1702 */
     1703void fmpq_poly_floordiv(fmpq_poly_ptr q, const fmpq_poly_ptr a, const fmpq_poly_ptr b)
     1704{
     1705    ulong m;
     1706    fmpz_t lead;
     1707   
     1708    /* Assertion! */
     1709    if (fmpq_poly_is_zero(b))
     1710    {
     1711        printf("ERROR (fmpq_poly_floordiv).  Division by zero.\n");
     1712        abort();
     1713    }
     1714   
     1715    /* Catch the case when a and b are the same objects in memory. */
     1716    /* As of FLINT version 1.5.1, there is a bug in this case.     */
     1717    if (a == b)
     1718    {
     1719        fmpq_poly_set_si(q, 1);
     1720        return;
     1721    }
     1722   
     1723    /* Deal with the various other cases of aliasing. */
     1724    if (q == a | q == b)
     1725    {
     1726        fmpq_poly_t tpoly;
     1727        fmpq_poly_init(tpoly);
     1728        fmpq_poly_floordiv(tpoly, a, b);
     1729        fmpq_poly_swap(q, tpoly);
     1730        fmpq_poly_clear(tpoly);
     1731        return;
     1732    }
     1733   
     1734    /* Deal separately with the case deg(b) = 0. */
     1735    if (fmpq_poly_length(b) == 1ul)
     1736    {
     1737        lead = fmpz_poly_get_coeff_ptr(b->num, 0);
     1738        if (fmpz_is_one(a->den))
     1739        {
     1740            if (fmpz_is_one(b->den))  /* a->den == b->den == 1 */
     1741                fmpz_poly_set(q->num, a->num);
     1742            else  /* a->den == 1, b->den > 1 */
     1743                fmpz_poly_scalar_mul_fmpz(q->num, a->num, b->den);
     1744           
     1745            q->den = fmpz_realloc(q->den, fmpz_size(lead));
     1746            fmpz_set(q->den, lead);
     1747            fmpq_poly_canonicalize(q, NULL);
     1748        }
     1749        else
     1750        {
     1751            if (fmpz_is_one(b->den))  /* a->den > 1, b->den == 1 */
     1752                fmpz_poly_set(q->num, a->num);
     1753            else  /* a->den, b->den > 1 */
     1754                fmpz_poly_scalar_mul_fmpz(q->num, a->num, b->den);
     1755           
     1756            q->den = fmpz_realloc(q->den, fmpz_size(a->den) + fmpz_size(lead));
     1757            fmpz_mul(q->den, a->den, lead);
     1758            fmpq_poly_canonicalize(q, NULL);
     1759        }
     1760        return;
     1761    }
     1762   
     1763    /* General case..                            */
     1764    /* Set q to b->den q->num / (a->den lead^m). */
     1765   
     1766    /* Since the fmpz_poly_pseudo_div method is broken as of FLINT 1.5.0, we */
     1767    /* use the fmpz_poly_pseudo_divrem method instead.                       */
     1768    /* fmpz_poly_pseudo_div(q->num, &m, a->num, b->num);                     */
     1769    {
     1770        fmpz_poly_t r;
     1771        fmpz_poly_init(r);
     1772        fmpz_poly_pseudo_divrem(q->num, r, &m, a->num, b->num);
     1773        fmpz_poly_clear(r);
     1774    }
     1775   
     1776    lead = fmpz_poly_lead(b->num);
     1777   
     1778    if (!fmpz_is_one(b->den))
     1779        fmpz_poly_scalar_mul_fmpz(q->num, q->num, b->den);
     1780   
     1781    /* Case 1:  lead^m is 1 */
     1782    if (fmpz_is_one(lead) || m == 0ul || (fmpz_is_m1(lead) & m % 2 == 0))
     1783    {
     1784        q->den = fmpz_realloc(q->den, fmpz_size(a->den));
     1785        fmpz_set(q->den, a->den);
     1786        fmpq_poly_canonicalize(q, NULL);
     1787    }
     1788    /* Case 2:  lead^m is -1 */
     1789    else if (fmpz_is_m1(lead) & m % 2)
     1790    {
     1791        fmpz_poly_neg(q->num, q->num);
     1792        q->den = fmpz_realloc(q->den, fmpz_size(a->den));
     1793        fmpz_set(q->den, a->den);
     1794        fmpq_poly_canonicalize(q, NULL);
     1795    }
     1796    /* Case 3:  lead^m is not +-1 */
     1797    else
     1798    {
     1799        ulong limbs;
     1800        limbs = m * fmpz_size(lead);
     1801       
     1802        if (fmpz_is_one(a->den))
     1803        {
     1804            q->den = fmpz_realloc(q->den, limbs);
     1805            fmpz_pow_ui(q->den, lead, m);
     1806        }
     1807        else
     1808        {
     1809            fmpz_t t;
     1810            t = fmpz_init(limbs);
     1811            q->den = fmpz_realloc(q->den, limbs + fmpz_size(a->den));
     1812            fmpz_pow_ui(t, lead, m);
     1813            fmpz_mul(q->den, t, a->den);
     1814            fmpz_clear(t);
     1815        }
     1816        fmpq_poly_canonicalize(q, NULL);
     1817    }
     1818}
     1819
     1820/**
     1821 * \ingroup  Division
     1822 *
     1823 * Sets \c r to the remainder of the Euclidean division of \c a by \c b.
     1824 *
     1825 * Assumes that \c b is non-zero.  Otherwise, an exception is raised in the
     1826 * form of an <tt>abort</tt> statement.
     1827 */
     1828void fmpq_poly_mod(fmpq_poly_ptr r, const fmpq_poly_ptr a, const fmpq_poly_ptr b)
     1829{
     1830    ulong m;
     1831    fmpz_t lead;
     1832   
     1833    /* Assertion! */
     1834    if (fmpq_poly_is_zero(b))
     1835    {
     1836        printf("ERROR (fmpq_poly_mod).  Division by zero.\n");
     1837        abort();
     1838    }
     1839   
     1840    /* Catch the case when a and b are the same objects in memory. */
     1841    /* As of FLINT version 1.5.1, there is a bug in this case.     */
     1842    if (a == b)
     1843    {
     1844        fmpq_poly_set_si(r, 0);
     1845        return;
     1846    }
     1847   
     1848    /* Deal with the various other cases of aliasing. */
     1849    if (r == a | r == b)
     1850    {
     1851        fmpq_poly_t tpoly;
     1852        fmpq_poly_init(tpoly);
     1853        fmpq_poly_mod(tpoly, a, b);
     1854        fmpq_poly_swap(r, tpoly);
     1855        fmpq_poly_clear(tpoly);
     1856        return;
     1857    }
     1858   
     1859    /* Deal separately with the case deg(b) = 0. */
     1860    if (fmpq_poly_length(b) == 1ul)
     1861    {
     1862        fmpz_poly_zero(r->num);
     1863        fmpz_set_ui(r->den, 1);
     1864        return;
     1865    }
     1866   
     1867    /* General case..                     */
     1868    /* Set r to r->num / (a->den lead^m). */
     1869   
     1870
     1871    /* As of FLINT 1.5.0, the fmpz_poly_pseudo_mod method is not             */
     1872    /* asymptotically fast and hence we swap to the fmpz_poly_pseudo_divrem  */
     1873    /* method if one of the operands' lengths is at least 32.                */
     1874    if (fmpq_poly_length(a) < 32 && fmpq_poly_length(b) < 32)
     1875        fmpz_poly_pseudo_rem(r->num, &m, a->num, b->num);
     1876    else
     1877    {
     1878        fmpz_poly_t q;
     1879        fmpz_poly_init(q);
     1880        fmpz_poly_pseudo_divrem(q, r->num, &m, a->num, b->num);
     1881        fmpz_poly_clear(q);
     1882    }
     1883   
     1884    lead = fmpz_poly_lead(b->num);
     1885   
     1886    /* Case 1:  lead^m is 1 */
     1887    if (fmpz_is_one(lead) || m == 0ul || (fmpz_is_m1(lead) & m % 2 == 0))
     1888    {
     1889        r->den = fmpz_realloc(r->den, fmpz_size(a->den));
     1890        fmpz_set(r->den, a->den);
     1891    }
     1892    /* Case 2:  lead^m is -1 */
     1893    else if (fmpz_is_m1(lead) & m % 2)
     1894    {
     1895        r->den = fmpz_realloc(r->den, fmpz_size(a->den));
     1896        fmpz_neg(r->den, a->den);
     1897    }
     1898    /* Case 3:  lead^m is not +-1 */
     1899    else
     1900    {
     1901        ulong limbs;
     1902        limbs = m * fmpz_size(lead);
     1903       
     1904        if (fmpz_is_one(a->den))
     1905        {
     1906            r->den = fmpz_realloc(r->den, limbs);
     1907            fmpz_pow_ui(r->den, lead, m);
     1908        }
     1909        else
     1910        {
     1911            fmpz_t t;
     1912            t = fmpz_init(limbs);
     1913            r->den = fmpz_realloc(r->den, limbs + fmpz_size(a->den));
     1914            fmpz_pow_ui(t, lead, m);
     1915            fmpz_mul(r->den, t, a->den);
     1916            fmpz_clear(t);
     1917        }
     1918    }
     1919    fmpq_poly_canonicalize(r, NULL);
     1920}
     1921
     1922/**
     1923 * \ingroup  Division
     1924 *
     1925 * Sets \c q and \c r to the quotient and remainder of the Euclidean
     1926 * division of \c a by \c b.
     1927 *
     1928 * Assumes that \c b is non-zero, and that \c q and \c r refer to distinct
     1929 * objects in memory.
     1930 */
     1931void fmpq_poly_divrem(fmpq_poly_ptr q, fmpq_poly_ptr r, const fmpq_poly_ptr a, const fmpq_poly_ptr b)
     1932{
     1933    ulong m;
     1934    fmpz_t lead;
     1935   
     1936    /* Assertion! */
     1937    if (fmpq_poly_is_zero(b))
     1938    {
     1939        printf("ERROR (fmpq_poly_divrem).  Division by zero.\n");
     1940        abort();
     1941    }
     1942    if (q == r)
     1943    {
     1944        printf("ERROR (fmpq_poly_divrem).  Output arguments aliased.\n");
     1945        abort();
     1946    }
     1947   
     1948    /* Catch the case when a and b are the same objects in memory.           */
     1949    /* As of FLINT version 1.5.1, there is a bug in this case.               */
     1950    if (a == b)
     1951    {
     1952        fmpq_poly_set_si(q, 1);
     1953        fmpq_poly_set_si(r, 0);
     1954        return;
     1955    }
     1956   
     1957    /* Deal with the various other cases of aliasing.                        */
     1958    if (r == a | r == b)
     1959    {
     1960        if (q == a | q == b)
     1961        {
     1962            fmpq_poly_t tempq, tempr;
     1963            fmpq_poly_init(tempq);
     1964            fmpq_poly_init(tempr);
     1965            fmpq_poly_divrem(tempq, tempr, a, b);
     1966            fmpq_poly_swap(q, tempq);
     1967            fmpq_poly_swap(r, tempr);
     1968            fmpq_poly_clear(tempq);
     1969            fmpq_poly_clear(tempr);
     1970            return;
     1971        }
     1972        else
     1973        {
     1974            fmpq_poly_t tempr;
     1975            fmpq_poly_init(tempr);
     1976            fmpq_poly_divrem(q, tempr, a, b);
     1977            fmpq_poly_swap(r, tempr);
     1978            fmpq_poly_clear(tempr);
     1979            return;
     1980        }
     1981    }
     1982    else
     1983    {
     1984        if (q == a | q == b)
     1985        {
     1986            fmpq_poly_t tempq;
     1987            fmpq_poly_init(tempq);
     1988            fmpq_poly_divrem(tempq, r, a, b);
     1989            fmpq_poly_swap(q, tempq);
     1990            fmpq_poly_clear(tempq);
     1991            return;
     1992        }
     1993    }
     1994   
     1995    // TODO: Implement the case `\deg(b) = 0` more efficiently!
     1996   
     1997    /* General case..                                                        */
     1998    /* Set q to b->den q->num / (a->den lead^m)                              */
     1999    /* and r to r->num / (a->den lead^m).                                    */
     2000   
     2001    fmpz_poly_pseudo_divrem(q->num, r->num, &m, a->num, b->num);
     2002   
     2003    lead = fmpz_poly_lead(b->num);
     2004   
     2005    if (!fmpz_is_one(b->den))
     2006        fmpz_poly_scalar_mul_fmpz(q->num, q->num, b->den);
     2007   
     2008    /* Case 1.  lead^m is 1 */
     2009    if (fmpz_is_one(lead) || m == 0ul || (fmpz_is_m1(lead) & m % 2 == 0))
     2010    {
     2011        q->den = fmpz_realloc(q->den, fmpz_size(a->den));
     2012        fmpz_set(q->den, a->den);
     2013    }
     2014    /* Case 2.  lead^m is -1 */
     2015    else if (fmpz_is_m1(lead) & m % 2)
     2016    {
     2017        fmpz_poly_neg(q->num, q->num);
     2018        q->den = fmpz_realloc(q->den, fmpz_size(a->den));
     2019        fmpz_set(q->den, a->den);
     2020       
     2021        fmpz_poly_neg(r->num, r->num);
     2022    }
     2023    /* Case 3.  lead^m is not +-1 */
     2024    else
     2025    {
     2026        ulong limbs;
     2027        limbs = m * fmpz_size(lead);
     2028       
     2029        if (fmpz_is_one(a->den))
     2030        {
     2031            q->den = fmpz_realloc(q->den, limbs);
     2032            fmpz_pow_ui(q->den, lead, m);
     2033        }
     2034        else
     2035        {
     2036            fmpz_t t;
     2037            t = fmpz_init(limbs);
     2038            q->den = fmpz_realloc(q->den, limbs + fmpz_size(a->den));
     2039            fmpz_pow_ui(t, lead, m);
     2040            fmpz_mul(q->den, t, a->den);
     2041            fmpz_clear(t);
     2042        }
     2043    }
     2044    r->den = fmpz_realloc(r->den, fmpz_size(q->den));
     2045    fmpz_set(r->den, q->den);
     2046   
     2047    fmpq_poly_canonicalize(q, NULL);
     2048    fmpq_poly_canonicalize(r, NULL);
     2049}
     2050
     2051///////////////////////////////////////////////////////////////////////////////
     2052// Powering
     2053
     2054/**
     2055 * \ingroup  Powering
     2056 *
     2057 * Sets \c rop to the <tt>exp</tt>th power of \c op.
     2058 *
     2059 * The corner case of <tt>exp == 0</tt> is handled by setting \c rop to the
     2060 * constant function \f$1\f$.  Note that this includes the case \f$0^0 = 1\f$.
     2061 */
     2062void fmpq_poly_power(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong exp)
     2063{
     2064    if (exp == 0ul)
     2065    {
     2066        fmpq_poly_one(rop);
     2067    }
     2068    else
     2069    {
     2070        fmpz_poly_power(rop->num, op->num, exp);
     2071        if (fmpz_is_one(op->den))
     2072        {
     2073            fmpz_set_ui(rop->den, 1);
     2074        }
     2075        else
     2076        {
     2077            if (rop == op)
     2078            {
     2079                fmpz_t t;
     2080                t = fmpz_init(exp * fmpz_size(op->den));
     2081                fmpz_pow_ui(t, op->den, exp);
     2082                fmpz_clear(rop->den);  // FIXME  Fixed
     2083                rop->den = t;
     2084            }
     2085            else
     2086            {
     2087                rop->den = fmpz_realloc(rop->den, exp * fmpz_size(op->den));
     2088                fmpz_pow_ui(rop->den, op->den, exp);
     2089            }
     2090        }
     2091    }
     2092}
     2093
     2094///////////////////////////////////////////////////////////////////////////////
     2095// Greatest common divisor
     2096
     2097/**
     2098 * \ingroup  GCD
     2099 *
     2100 * Returns the (monic) greatest common divisor \c res of \c a and \c b.
     2101 *
     2102 * Corner cases:  If \c a and \c b are both zero, returns zero.  If only
     2103 * one of them is zero, returns the other polynomial, up to normalisation.
     2104 */
     2105void fmpq_poly_gcd(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b)
     2106{
     2107    fmpz_t lead;
     2108    fmpz_poly_t num;
     2109   
     2110    /* Deal with aliasing */
     2111    if (rop == a | rop == b)
     2112    {
     2113        fmpq_poly_t tpoly;
     2114        fmpq_poly_init(tpoly);
     2115        fmpq_poly_gcd(tpoly, a, b);
     2116        fmpq_poly_swap(rop, tpoly);
     2117        fmpq_poly_clear(tpoly);
     2118        return;
     2119    }
     2120   
     2121    /* Deal with corner cases */
     2122    if (fmpq_poly_is_zero(a))
     2123    {
     2124        if (fmpq_poly_is_zero(b))  /* a == b == 0 */
     2125            fmpq_poly_zero(rop);
     2126        else                       /* a == 0, b != 0 */
     2127            fmpq_poly_monic(rop, b);
     2128        return;
     2129    }
     2130    else
     2131    {
     2132        if (fmpq_poly_is_zero(b))  /* a != 0, b == 0 */
     2133        {
     2134            fmpq_poly_monic(rop, a);
     2135            return;
     2136        }
     2137    }
     2138   
     2139    /* General case.. */
     2140    fmpz_poly_init(num);
     2141    fmpz_poly_primitive_part(rop->num, a->num);
     2142    fmpz_poly_primitive_part(num, b->num);
     2143   
     2144    /* Since rop.num is the greatest common divisor of the primitive parts   */
     2145    /* of a.num and b.num, it is also primitive.  But as of FLINT 1.4.0, the */
     2146    /* leading term *might* be negative.                                     */
     2147    fmpz_poly_gcd(rop->num, rop->num, num);
     2148   
     2149    lead = fmpz_poly_lead(rop->num);
     2150    rop->den = fmpz_realloc(rop->den, fmpz_size(lead));
     2151    if (fmpz_sgn(lead) < 0)
     2152    {
     2153        fmpz_poly_neg(rop->num, rop->num);
     2154        fmpz_neg(rop->den, lead);
     2155    }
     2156    else
     2157        fmpz_set(rop->den, lead);
     2158   
     2159    fmpz_poly_clear(num);
     2160}
     2161
     2162/**
     2163 * \ingroup  GCD
     2164 *
     2165 * Returns polynomials \c s, \c t and \c rop such that \c rop is
     2166 * (monic) greatest common divisor of \c a and \c b, and such that
     2167 * <tt>rop = s a + t b</tt>.
     2168 *
     2169 * Corner cases:  If \c a and \c b are zero, returns zero polynomials.
     2170 * Otherwise, if only \c a is zero, returns <tt>(res, s, t) = (b, 0, 1)</tt>
     2171 * up to normalisation, and similarly if only \c b is zero.
     2172 *
     2173 * Assumes that the output parameters \c rop, \c s, and \c t refer to
     2174 * distinct objects in memory.  Otherwise, an exception is raised in the
     2175 * form of an <tt>abort</tt> statement.
     2176 */
     2177void 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)
     2178{
     2179    fmpz_t lead, temp;
     2180    ulong bound, limbs;
     2181   
     2182    /* Assertion! */
     2183    if (rop == s | rop == t | s == t)
     2184    {
     2185        printf("ERROR (fmpq_poly_xgcd).  Output arguments aliased.\n");
     2186        abort();
     2187    }
     2188   
     2189    /* Deal with aliasing */
     2190    if (rop == a | rop == b)
     2191    {
     2192        if (s == a | s == b)
     2193        {
     2194            /* We know that t does not coincide with a or b, since otherwise */
     2195            /* two of rop, s, and t coincide, too.                           */
     2196            fmpq_poly_t tempg, temps;
     2197            fmpq_poly_init(tempg);
     2198            fmpq_poly_init(temps);
     2199            fmpq_poly_xgcd(tempg, temps, t, a, b);
     2200            fmpq_poly_swap(rop, tempg);
     2201            fmpq_poly_swap(s, temps);
     2202            fmpq_poly_clear(tempg);
     2203            fmpq_poly_clear(temps);
     2204            return;
     2205        }
     2206        else
     2207        {
     2208            if (t == a | t == b)
     2209            {
     2210                fmpq_poly_t tempg, tempt;
     2211                fmpq_poly_init(tempg);
     2212                fmpq_poly_init(tempt);
     2213                fmpq_poly_xgcd(tempg, s, tempt, a, b);
     2214                fmpq_poly_swap(rop, tempg);
     2215                fmpq_poly_swap(t, tempt);
     2216                fmpq_poly_clear(tempg);
     2217                fmpq_poly_clear(tempt);
     2218                return;
     2219            }
     2220            else
     2221            {
     2222                fmpq_poly_t tempg;
     2223                fmpq_poly_init(tempg);
     2224                fmpq_poly_xgcd(tempg, s, t, a, b);
     2225                fmpq_poly_swap(rop, tempg);
     2226                fmpq_poly_clear(tempg);
     2227                return;
     2228            }
     2229        }
     2230    }
     2231    else
     2232    {
     2233        if (s == a | s == b)
     2234        {
     2235            if (t == a | t == b)
     2236            {
     2237                fmpq_poly_t temps, tempt;
     2238                fmpq_poly_init(temps);
     2239                fmpq_poly_init(tempt);
     2240                fmpq_poly_xgcd(rop, temps, tempt, a, b);
     2241                fmpq_poly_swap(s, temps);
     2242                fmpq_poly_swap(t, tempt);
     2243                fmpq_poly_clear(temps);
     2244                fmpq_poly_clear(tempt);
     2245                return;
     2246            }
     2247            else
     2248            {
     2249                fmpq_poly_t temps;
     2250                fmpq_poly_init(temps);
     2251                fmpq_poly_xgcd(rop, temps, t, a, b);
     2252                fmpq_poly_swap(s, temps);
     2253                fmpq_poly_clear(temps);
     2254                return;
     2255            }
     2256        }
     2257        else
     2258        {
     2259            if (t == a | t == b)
     2260            {
     2261                fmpq_poly_t tempt;
     2262                fmpq_poly_init(tempt);
     2263                fmpq_poly_xgcd(rop, s, tempt, a, b);
     2264                fmpq_poly_swap(t, tempt);
     2265                fmpq_poly_clear(tempt);
     2266                return;
     2267            }
     2268        }
     2269    }
     2270   
     2271    /* From here on, we may assume that none of the output variables are */
     2272    /* aliases for the input variables.                                  */
     2273   
     2274    /* Deal with the following three corner cases: */
     2275    /*   a == 0, b == 0                            */
     2276    /*   a == 0, b =! 0                            */
     2277    /*   a != 0, b == 0                            */
     2278    if (fmpq_poly_is_zero(a))
     2279    {
     2280        if (fmpq_poly_is_zero(b))  /* Case 1.  a == b == 0 */
     2281        {
     2282            fmpq_poly_zero(rop);
     2283            fmpq_poly_zero(s);
     2284            fmpq_poly_zero(t);
     2285            return;
     2286        }
     2287        else  /* Case 2.  a == 0, b != 0 */
     2288        {
     2289            fmpz_t blead = fmpz_poly_lead(b->num);
     2290           
     2291            fmpq_poly_monic(rop, b);
     2292            fmpq_poly_zero(s);
     2293           
     2294            fmpz_poly_zero(t->num);
     2295            fmpz_poly_set_coeff_fmpz(t->num, 0, b->den);
     2296            if (fmpz_is_one(blead))
     2297                fmpz_set_ui(t->den, 1);
     2298            else
     2299            {
     2300                fmpz_t g;
     2301                g = fmpz_init(FLINT_MAX(b->num->limbs, fmpz_size(b->den)));
     2302                fmpz_gcd(g, blead, b->den);
     2303                if (fmpz_sgn(g) != fmpz_sgn(blead))
     2304                    fmpz_neg(g, g);
     2305                fmpz_poly_scalar_div_fmpz(t->num, t->num, g);
     2306                t->den = fmpz_realloc(t->den, fmpz_size(blead));
     2307                fmpz_divides(t->den, blead, g);
     2308                fmpz_clear(g);
     2309            }
     2310            return;
     2311        }
     2312    }
     2313    else
     2314    {
     2315        if (fmpq_poly_is_zero(b))  /* Case 3.  a != 0, b == 0 */
     2316        {
     2317            fmpq_poly_xgcd(rop, t, s, b, a);
     2318            return;
     2319        }
     2320    }
     2321   
     2322    /* We are now in the general case where a and b are non-zero. */
     2323   
     2324    s->den = fmpz_realloc(s->den, a->num->limbs);
     2325    t->den = fmpz_realloc(t->den, b->num->limbs);
     2326    fmpz_poly_content(s->den, a->num);
     2327    fmpz_poly_content(t->den, b->num);
     2328    fmpz_poly_scalar_div_fmpz(s->num, a->num, s->den);
     2329    fmpz_poly_scalar_div_fmpz(t->num, b->num, t->den);
     2330   
     2331    /* Note that, since s->num and t->num are primitive, rop->num is         */
     2332    /* primitive, too.  In fact, it is the rational greatest common divisor  */
     2333    /* of a and b.  As of FLINT 1.4.0, the leading coefficient of res.num    */
     2334    /* *might* be negative.                                                  */
     2335   
     2336    fmpz_poly_gcd(rop->num, s->num, t->num);
     2337    if (fmpz_sgn(fmpz_poly_lead(rop->num)) < 0)
     2338        fmpz_poly_neg(rop->num, rop->num);
     2339    lead = fmpz_poly_lead(rop->num);
     2340   
     2341    /* Now rop->num is a (primitive) rational greatest common divisor of */
     2342    /* a and b.                                                          */
     2343   
     2344    if (fmpz_poly_degree(rop->num) > 0)
     2345    {
     2346        fmpz_poly_div(s->num, s->num, rop->num);
     2347        fmpz_poly_div(t->num, t->num, rop->num);
     2348    }
     2349   
     2350    bound = fmpz_poly_resultant_bound(s->num, t->num);
     2351    if (fmpz_is_one(lead))
     2352        bound = bound / FLINT_BITS + 2;
     2353    else
     2354        bound = FLINT_MAX(bound / FLINT_BITS + 2, fmpz_size(lead));
     2355    rop->den = fmpz_realloc(rop->den, bound);
     2356   
     2357    fmpz_poly_xgcd(rop->den, s->num, t->num, s->num, t->num);
     2358   
     2359    /* Now the following equation holds:                                     */
     2360    /*   rop->den rop->num ==                                                */
     2361    /*             (s->num a->den / s->den) a +  (t->num b->den / t->den) b. */
     2362   
     2363    limbs = FLINT_MAX(s->num->limbs, t->num->limbs);
     2364    limbs = FLINT_MAX(limbs, fmpz_size(s->den));
     2365    limbs = FLINT_MAX(limbs, fmpz_size(t->den) + fmpz_size(rop->den) + fmpz_size(lead));
     2366    temp = fmpz_init(limbs);
     2367   
     2368    s->den = fmpz_realloc(s->den, fmpz_size(s->den) + fmpz_size(rop->den)
     2369                                                    + fmpz_size(lead));
     2370    if (!fmpz_is_one(a->den))
     2371        fmpz_poly_scalar_mul_fmpz(s->num, s->num, a->den);
     2372    fmpz_mul(temp, s->den, rop->den);
     2373    fmpz_mul(s->den, temp, lead);
     2374   
     2375    t->den = fmpz_realloc(t->den, fmpz_size(t->den) + fmpz_size(rop->den)
     2376                                                    + fmpz_size(lead));
     2377    if (!fmpz_is_one(b->den))
     2378        fmpz_poly_scalar_mul_fmpz(t->num, t->num, b->den);
     2379    fmpz_mul(temp, t->den, rop->den);
     2380    fmpz_mul(t->den, temp, lead);
     2381   
     2382    fmpq_poly_canonicalize(s, temp);
     2383    fmpq_poly_canonicalize(t, temp);
     2384   
     2385    fmpz_set(rop->den, lead);
     2386   
     2387    fmpz_clear(temp);
     2388}
     2389
     2390/**
     2391 * \ingroup  GCD
     2392 *
     2393 * Computes the monic (or zero) least common multiple of \c a and \c b.
     2394 *
     2395 * If either of \c a and \c b is zero, returns zero.  This behaviour ensures
     2396 * that the relation
     2397 * \f[
     2398 * \text{lcm}(a,b) \gcd(a,b) \sim a b
     2399 * \f]
     2400 * holds, where \f$\sim\f$ denotes equality up to units.
     2401 */
     2402void fmpq_poly_lcm(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b)
     2403{
     2404    /* Handle aliasing */
     2405    if (rop == a | rop == b)
     2406    {
     2407        fmpq_poly_t tpoly;
     2408        fmpq_poly_init(tpoly);
     2409        fmpq_poly_lcm(tpoly, a, b);
     2410        fmpq_poly_swap(rop, tpoly);
     2411        fmpq_poly_clear(tpoly);
     2412        return;
     2413    }
     2414   
     2415    if (fmpq_poly_is_zero(a))
     2416        fmpq_poly_zero(rop);
     2417    else if (fmpq_poly_is_zero(b))
     2418        fmpq_poly_zero(rop);
     2419    else
     2420    {
     2421        fmpz_poly_t prod;
     2422        fmpz_t lead;
     2423       
     2424        fmpz_poly_init(prod);
     2425        fmpq_poly_gcd(rop, a, b);
     2426        fmpz_poly_mul(prod, a->num, b->num);
     2427        fmpz_poly_primitive_part(prod, prod);
     2428        fmpz_poly_div(rop->num, prod, rop->num);
     2429       
     2430        /* Note that GCD returns a monic rop and so a primitive rop.num.   */
     2431        /* Dividing the primitive prod by this yields a primitive quotient */
     2432        /* rop->num.                                                       */
     2433       
     2434        lead = fmpz_poly_lead(rop->num);
     2435        rop->den = fmpz_realloc(rop->den, fmpz_size(lead));
     2436        if (fmpz_sgn(lead) < 0)
     2437            fmpz_poly_neg(rop->num, rop->num);
     2438        fmpz_set(rop->den, fmpz_poly_lead(rop->num));
     2439       
     2440        fmpz_poly_clear(prod);
     2441    }
     2442}
     2443
     2444///////////////////////////////////////////////////////////////////////////////
     2445// Derivative
     2446
     2447/**
     2448 * \ingroup  Derivative
     2449 *
     2450 * Sets \c rop to the derivative of \c op.
     2451 *
     2452 * \todo  The second argument should be declared \c const, but as of
     2453 *        FLINT 1.5.0 this generates a compile-time warning.
     2454 */
     2455void fmpq_poly_derivative(fmpq_poly_ptr rop, fmpq_poly_ptr op)
     2456{
     2457    if (fmpq_poly_length(op) < 2ul)
     2458        fmpq_poly_zero(rop);
     2459    else
     2460    {
     2461        fmpz_poly_derivative(rop->num, op->num);
     2462        rop->den = fmpz_realloc(rop->den, fmpz_size(op->den));
     2463        fmpz_set(rop->den, op->den);
     2464        fmpq_poly_canonicalize(rop, NULL);
     2465    }
     2466}
     2467
     2468///////////////////////////////////////////////////////////////////////////////
     2469// Evaluation
     2470
     2471/**
     2472 * \ingroup  Evaluation
     2473 *
     2474 * Evaluates the integer polynomial \c f at the rational \c a using Horner's
     2475 * method.
     2476 */
     2477void _fmpz_poly_evaluate_mpq_horner(mpq_t rop, const fmpz_poly_t f, const mpq_t a)
     2478{
     2479    mpq_t temp;
     2480    ulong n;
     2481   
     2482    /* Handle aliasing */
     2483    if (rop == a)
     2484    {
     2485        mpq_t tempr;
     2486        mpq_init(tempr);
     2487        _fmpz_poly_evaluate_mpq_horner(tempr, f, a);
     2488        mpq_swap(rop, tempr);
     2489        mpq_clear(tempr);
     2490        return;
     2491    }
     2492   
     2493    n = fmpz_poly_length(f);
     2494   
     2495    if (n == 0ul)
     2496    {
     2497        mpq_set_ui(rop, 0, 1);
     2498    }
     2499    else if (n == 1ul)
     2500    {
     2501        fmpz_poly_get_coeff_mpz(mpq_numref(rop), f, 0);
     2502        mpz_set_ui(mpq_denref(rop), 1);
     2503    }
     2504    else
     2505    {
     2506        n--;
     2507        fmpz_poly_get_coeff_mpz(mpq_numref(rop), f, n);
     2508        mpz_set_ui(mpq_denref(rop), 1);
     2509        mpq_init(temp);
     2510        do {
     2511            n--;
     2512            mpq_mul(temp, rop, a);
     2513            fmpz_poly_get_coeff_mpz(mpq_numref(rop), f, n);
     2514            mpz_set_ui(mpq_denref(rop), 1);
     2515            mpq_add(rop, rop, temp);
     2516        } while (n);
     2517        mpq_clear(temp);
     2518    }
     2519}
     2520
     2521/**
     2522 * \ingroup  Evaluation
     2523 *
     2524 * Evaluates the rational polynomial \c f at the integer \c a.
     2525 *
     2526 * Assumes that the numerator and denominator of the <tt>mpq_t</tt>
     2527 * \c rop are distinct (as objects in memory) from the <tt>mpz_t</tt>
     2528 * \c a.
     2529 */
     2530void fmpq_poly_evaluate_mpz(mpq_t rop, fmpq_poly_ptr f, const mpz_t a)
     2531{
     2532    fmpz_t num, t;
     2533    ulong limbs, max;
     2534   
     2535    if (fmpq_poly_is_zero(f))
     2536    {
     2537        mpq_set_ui(rop, 0, 1);
     2538        return;
     2539    }
     2540   
     2541    /* Establish a bound on the size of f->num evaluated at a */
     2542    max = (f->num->length) * (f->num->limbs + f->num->length * mpz_size(a));
     2543   
     2544    /* Compute the result */
     2545    num = fmpz_init(max);
     2546    t   = fmpz_init(mpz_size(a));
     2547    mpz_to_fmpz(t, a);
     2548    fmpz_poly_evaluate(num, f->num, t);
     2549    fmpz_to_mpz(mpq_numref(rop), num);
     2550    if (fmpz_is_one(f->den))
     2551    {
     2552        mpz_set_ui(mpq_denref(rop), 1);
     2553    }
     2554    else
     2555    {
     2556        fmpz_to_mpz(mpq_denref(rop), f->den);
     2557        mpq_canonicalize(rop);
     2558    }
     2559   
     2560    fmpz_clear(num);
     2561    fmpz_clear(t);
     2562}
     2563
     2564/**
     2565 * \ingroup  Evaluation
     2566 *
     2567 * Evaluates the rational polynomial \c f at the rational \c a.
     2568 */
     2569void fmpq_poly_evaluate_mpq(mpq_t rop, fmpq_poly_ptr f, const mpq_t a)
     2570{
     2571    if (rop == a)
     2572    {
     2573        mpq_t tempr;
     2574        mpq_init(tempr);
     2575        fmpq_poly_evaluate_mpq(tempr, f, a);
     2576        mpq_swap(rop, tempr);
     2577        mpq_clear(tempr);
     2578        return;
     2579    }
     2580   
     2581    _fmpz_poly_evaluate_mpq_horner(rop, f->num, a);
     2582    if (!fmpz_is_one(f->den))
     2583    {
     2584        mpz_t den;
     2585        mpz_init(den);
     2586        fmpz_to_mpz(den, f->den);
     2587        mpz_mul(mpq_denref(rop), mpq_denref(rop), den);
     2588        mpq_canonicalize(rop);
     2589        mpz_clear(den);
     2590    }
     2591}
     2592
     2593///////////////////////////////////////////////////////////////////////////////
     2594// Gaussian content
     2595
     2596/**
     2597 * \ingroup  Content
     2598 *
     2599 * Returns the non-negative content of \c op.
     2600 *
     2601 * The content of \f$0\f$ is defined to be \f$0\f$.
     2602 */
     2603void fmpq_poly_content(mpq_t rop, const fmpq_poly_ptr op)
     2604{
     2605    fmpz_t numc;
     2606   
     2607    numc = fmpz_init(op->num->limbs);
     2608    fmpz_poly_content(numc, op->num);
     2609    fmpz_abs(numc, numc);
     2610    fmpz_to_mpz(mpq_numref(rop), numc);
     2611    if (fmpz_is_one(op->den))
     2612        mpz_set_ui(mpq_denref(rop), 1);
     2613    else
     2614        fmpz_to_mpz(mpq_denref(rop), op->den);
     2615    fmpz_clear(numc);
     2616}
     2617
     2618/**
     2619 * \ingroup  Content
     2620 *
     2621 * Returns the primitive part (with non-negative leading coefficient) of
     2622 * \c op as an element of type #fmpq_poly_t.
     2623 */
     2624void fmpq_poly_primitive_part(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
     2625{
     2626    if (fmpq_poly_is_zero(op))
     2627        fmpq_poly_zero(rop);
     2628    else
     2629    {
     2630        fmpz_poly_primitive_part(rop->num, op->num);
     2631        if (fmpz_sgn(fmpz_poly_lead(rop->num)) < 0)
     2632            fmpz_poly_neg(rop->num, rop->num);
     2633        fmpz_set_ui(rop->den, 1);
     2634    }
     2635}
     2636
     2637/**
     2638 * \brief    Returns whether \c op is monic.
     2639 * \ingroup  Content
     2640 *
     2641 * Returns whether \c op is monic.
     2642 *
     2643 * By definition, the zero polynomial is \e not monic.
     2644 */
     2645int fmpq_poly_is_monic(const fmpq_poly_ptr op)
     2646{
     2647    fmpz_t lead;
     2648   
     2649    if (fmpq_poly_is_zero(op))
     2650        return 0;
     2651    else
     2652    {
     2653        lead = fmpz_poly_lead(op->num);
     2654        if (fmpz_is_one(op->den))
     2655            return fmpz_is_one(lead);
     2656        else
     2657            return (fmpz_sgn(lead) > 0) && (fmpz_cmpabs(lead, op->den) == 0);
     2658    }
     2659}
     2660
     2661/**
     2662 * Sets \c rop to the unique monic scalar multiple of \c op.
     2663 *
     2664 * As the only special case, if \c op is the zero polynomial, \c rop is set
     2665 * to zero, too.
     2666 */
     2667void fmpq_poly_monic(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
     2668{
     2669    fmpz_t lead;
     2670   
     2671    if (fmpq_poly_is_zero(op))
     2672    {
     2673        fmpq_poly_zero(rop);
     2674        return;
     2675    }
     2676   
     2677    fmpz_poly_primitive_part(rop->num, op->num);
     2678    lead = fmpz_poly_lead(rop->num);
     2679    rop->den = fmpz_realloc(rop->den, fmpz_size(lead));
     2680    if (fmpz_sgn(lead) < 0)
     2681        fmpz_poly_neg(rop->num, rop->num);
     2682    fmpz_set(rop->den, fmpz_poly_lead(rop->num));
     2683}
     2684
     2685///////////////////////////////////////////////////////////////////////////////
     2686// Resultant
     2687
     2688/**
     2689 * \brief    Returns the resultant of \c a and \c b.
     2690 * \ingroup  Resultant
     2691 *
     2692 * Returns the resultant of \c a and \c b.
     2693 *
     2694 * Enumerating the roots of \c a and \c b over \f$\bar{\mathbf{Q}}\f$ as
     2695 * \f$r_1, \dotsc, r_m\f$ and \f$s_1, \dotsc, s_n\f$, respectively, and
     2696 * letting \f$x\f$ and \f$y\f$ denote the leading coefficients, the resultant
     2697 * is defined as
     2698 * \f[
     2699 * x^{\deg(b)} y^{\deg(a)} \prod_{1 \leq i, j \leq n} (r_i - s_j).
     2700 * \f]
     2701 *
     2702 * We handle special cases as follows:  if one of the polynomials is zero,
     2703 * the resultant is zero.  Note that otherwise if one of the polynomials is
     2704 * constant, the last term in the above expression is the empty product.
     2705 */
     2706void fmpq_poly_resultant(mpq_t rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b)
     2707{
     2708    fmpz_t rest, t1, t2;
     2709    fmpz_poly_t anum, bnum;
     2710    fmpz_t acont, bcont;
     2711    long d1, d2;
     2712    ulong bound, denbound, numbound;
     2713    fmpz_poly_t g;
     2714   
     2715    d1 = fmpq_poly_degree(a);
     2716    d2 = fmpq_poly_degree(b);
     2717   
     2718    /* We first handle special cases. */
     2719    if (d1 < 0l | d2 < 0l)
     2720    {
     2721        mpq_set_ui(rop, 0, 1);
     2722        return;
     2723    }
     2724    if (d1 == 0l)
     2725    {
     2726        if (d2 == 0l)
     2727            mpq_set_ui(rop, 1, 1);
     2728        else if (d2 == 1l)
     2729        {
     2730            fmpz_to_mpz(mpq_numref(rop), fmpz_poly_lead(a->num));
     2731            fmpz_to_mpz(mpq_denref(rop), a->den);
     2732        }
     2733        else
     2734        {
     2735            if (fmpz_is_one(a->den))
     2736                bound = a->num->limbs;
     2737            else
     2738                bound = FLINT_MAX(a->num->limbs, fmpz_size(a->den));
     2739            t1 = fmpz_init(d2 * bound);
     2740            fmpz_pow_ui(t1, fmpz_poly_lead(a->num), d2);
     2741            fmpz_to_mpz(mpq_numref(rop), t1);
     2742            if (fmpz_is_one(a->den))
     2743                mpz_set_ui(mpq_denref(rop), 1);
     2744            else
     2745            {
     2746                fmpz_pow_ui(t1, a->den, d2);
     2747                fmpz_to_mpz(mpq_denref(rop), t1);
     2748            }
     2749            fmpz_clear(t1);
     2750        }
     2751        return;
     2752    }
     2753    if (d2 == 0l)
     2754    {
     2755        fmpq_poly_resultant(rop, b, a);
     2756        return;
     2757    }
     2758   
     2759    /* We are now in the general case, with both polynomials of degree at */
     2760    /* least 1.                                                           */
     2761   
     2762    /* We set a->num =: acont anum with acont > 0 and anum primitive. */
     2763    acont = fmpz_init(a->num->limbs);
     2764    fmpz_poly_content(acont, a->num);
     2765    fmpz_abs(acont, acont);
     2766    fmpz_poly_init(anum);
     2767    fmpz_poly_scalar_div_fmpz(anum, a->num, acont);
     2768   
     2769    /* We set b->num =: bcont bnum with bcont > 0 and bnum primitive. */
     2770    bcont = fmpz_init(b->num->limbs);
     2771    fmpz_poly_content(bcont, b->num);
     2772    fmpz_abs(bcont, bcont);
     2773    fmpz_poly_init(bnum);
     2774    fmpz_poly_scalar_div_fmpz(bnum, b->num, bcont);
     2775   
     2776    fmpz_poly_init(g);
     2777    fmpz_poly_gcd(g, anum, bnum);
     2778   
     2779    /* If the gcd has positive degree, the resultant is zero. */
     2780    if (fmpz_poly_degree(g) > 0)
     2781    {
     2782        mpq_set_ui(rop, 0, 1);
     2783       
     2784        /* Clean up */
     2785        fmpz_clear(acont);
     2786        fmpz_clear(bcont);
     2787        fmpz_poly_clear(anum);
     2788        fmpz_poly_clear(bnum);
     2789        fmpz_poly_clear(g);
     2790        return;
     2791    }
     2792   
     2793    /* Set some bounds */
     2794    if (fmpz_is_one(a->den))
     2795    {
     2796        if (fmpz_is_one(b->den))
     2797        {
     2798            numbound = FLINT_MAX(d2 * fmpz_size(acont), d1 * fmpz_size(bcont));
     2799            denbound = 1;
     2800        }
     2801        else
     2802        {
     2803            numbound = FLINT_MAX(d2 * fmpz_size(acont), d1 * fmpz_size(bcont));
     2804            denbound = d1 * fmpz_size(b->den);
     2805        }
     2806    }
     2807    else
     2808    {
     2809        if (fmpz_is_one(b->den))
     2810        {
     2811            numbound = FLINT_MAX(d2 * fmpz_size(acont), d1 * fmpz_size(bcont));
     2812            denbound = d2 * fmpz_size(a->den);
     2813        }
     2814        else
     2815        {
     2816            numbound = FLINT_MAX(d2 * fmpz_size(acont), d1 * fmpz_size(bcont));
     2817            denbound = FLINT_MAX(d2 * fmpz_size(a->den), d1 * fmpz_size(b->den));
     2818        }
     2819    }
     2820   
     2821    /* Now anum and bnum are coprime and we compute their resultant using */
     2822    /* the method from the fmpz_poly module.                              */
     2823    bound = fmpz_poly_resultant_bound(anum, bnum);
     2824    bound = bound/FLINT_BITS + 2 + d2 * fmpz_size(acont)
     2825        + d1 * fmpz_size(bcont);
     2826    rest = fmpz_init(FLINT_MAX(bound, denbound));
     2827    fmpz_poly_resultant(rest, anum, bnum);
     2828   
     2829    /* Finally, w take into account the factors acont/a.den and bcont/b.den. */
     2830    t1 = fmpz_init(FLINT_MAX(bound, denbound));
     2831    t2 = fmpz_init(FLINT_MAX(numbound, denbound));
     2832   
     2833    if (!fmpz_is_one(acont))
     2834    {
     2835        fmpz_pow_ui(t2, acont, d2);
     2836        fmpz_set(t1, rest);
     2837        fmpz_mul(rest, t1, t2);
     2838    }
     2839    if (!fmpz_is_one(bcont))
     2840    {
     2841        fmpz_pow_ui(t2, bcont, d1);
     2842        fmpz_set(t1, rest);
     2843        fmpz_mul(rest, t1, t2);
     2844    }
     2845   
     2846    fmpz_to_mpz(mpq_numref(rop), rest);
     2847   
     2848    if (fmpz_is_one(a->den))
     2849    {
     2850        if (fmpz_is_one(b->den))
     2851            fmpz_set_ui(rest, 1);
     2852        else
     2853            fmpz_pow_ui(rest, b->den, d1);
     2854    }
     2855    else
     2856    {
     2857        if (fmpz_is_one(b->den))
     2858            fmpz_pow_ui(rest, a->den, d2);
     2859        else
     2860        {
     2861            fmpz_pow_ui(t1, a->den, d2);
     2862            fmpz_pow_ui(t2, b->den, d1);
     2863            fmpz_mul(rest, t1, t2);
     2864        }
     2865    }
     2866   
     2867    fmpz_to_mpz(mpq_denref(rop), rest);
     2868    mpq_canonicalize(rop);
     2869   
     2870    /* Clean up */
     2871    fmpz_clear(rest);
     2872    fmpz_clear(t1);
     2873    fmpz_clear(t2);
     2874    fmpz_poly_clear(anum);
     2875    fmpz_poly_clear(bnum);
     2876    fmpz_clear(acont);
     2877    fmpz_clear(bcont);
     2878    fmpz_poly_clear(g);
     2879}
     2880
     2881/**
     2882 * \brief    Computes the discriminant of \c a.
     2883 * \ingroup  Discriminant
     2884 *
     2885 * Computes the discriminant of the polynomial \f$a\f$.
     2886 *
     2887 * The discriminant \f$R_n\f$ is defined as
     2888 * \f[
     2889 * R_n = a_n^{2 n-2} \prod_{1 \le i < j \le n} (r_i - r_j)^2,
     2890 * \f]
     2891 * where \f$n\f$ is the degree of this polynomial, \f$a_n\f$ is the leading
     2892 * coefficient and \f$r_1, \ldots, r_n\f$ are the roots over
     2893 * \f$\bar{\mathbf{Q}}\f$ are.
     2894 *
     2895 * The discriminant of constant polynomials is defined to be \f$0\f$.
     2896 *
     2897 * This implementation uses the identity
     2898 * \f[
     2899 * R_n(f) := (-1)^(n (n-1)/2) R(f,f') / a_n,
     2900 * \f]
     2901 * where \f$n\f$ is the degree of this polynomial, \f$a_n\f$ is the leading
     2902 * coefficient and \f$f'\f$ is the derivative of \f$f\f$.
     2903 *
     2904 * \see  #fmpq_poly_resultant()
     2905 */
     2906void fmpq_poly_discriminant(mpq_t disc, fmpq_poly_t a)
     2907{
     2908    fmpq_poly_t der;
     2909    mpq_t t;
     2910    long n;
     2911   
     2912    n = fmpq_poly_degree(a);
     2913    if (n < 1L)
     2914    {
     2915        mpq_set_ui(disc, 0, 1);
     2916        return;
     2917    }
     2918   
     2919    fmpq_poly_init(der);
     2920    fmpq_poly_derivative(der, a);
     2921    fmpq_poly_resultant(disc, a, der);
     2922    mpq_init(t);
     2923   
     2924    fmpz_to_mpz(mpq_numref(t), a->den);
     2925    fmpz_to_mpz(mpq_denref(t), fmpz_poly_lead(a->num));
     2926    mpq_canonicalize(t);
     2927    mpq_mul(disc, disc, t);
     2928   
     2929    if (n % 4 > 1)
     2930        mpz_neg(mpq_numref(disc), mpq_numref(disc));
     2931   
     2932    fmpq_poly_clear(der);
     2933    mpq_clear(t);
     2934}
     2935
     2936///////////////////////////////////////////////////////////////////////////////
     2937// Composition
     2938
     2939/**
     2940 * \brief    Returns the composition of \c a and \c b.
     2941 * \ingroup  Composition
     2942 *
     2943 * Returns the composition of \c a and \c b.
     2944 *
     2945 * To be clear about the order of the composition, denoting the polynomials
     2946 * \f$a(T)\f$ and \f$b(T)\f$, returns the polynomial \f$a(b(T))\f$.
     2947 */
     2948void fmpq_poly_compose(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b)
     2949{
     2950    mpq_t x;            /* Rational to hold the inverse of b->den */
     2951   
     2952    if (fmpz_is_one(b->den))
     2953    {
     2954        fmpz_poly_compose(rop->num, a->num, b->num);
     2955        rop->den = fmpz_realloc(rop->den, fmpz_size(a->den));
     2956        fmpz_set(rop->den, a->den);
     2957        fmpq_poly_canonicalize(rop, NULL);
     2958        return;
     2959    }
     2960   
     2961    /* Aliasing.                                                            */
     2962    /*                                                                      */
     2963    /* Note that rop and a, as well as a and b may be aliased, but rop and  */
     2964    /* b may not be aliased.                                                */
     2965    if (rop == b)
     2966    {
     2967        fmpq_poly_t tempr;
     2968        fmpq_poly_init(tempr);
     2969        fmpq_poly_compose(tempr, a, b);
     2970        fmpq_poly_swap(rop, tempr);
     2971        fmpq_poly_clear(tempr);
     2972        return;
     2973    }
     2974   
     2975    /* Set x = 1/b.den, and note this is already in canonical form. */
     2976    mpq_init(x);
     2977    mpz_set_ui(mpq_numref(x), 1);
     2978    fmpz_to_mpz(mpq_denref(x), b->den);
     2979   
     2980    /* First set rop = a(T / b.den) and then use FLINT's composition to */
     2981    /* set rop->num = rop->num(b->num).                                 */
     2982    fmpq_poly_rescale(rop, a, x);
     2983    fmpz_poly_compose(rop->num, rop->num, b->num);
     2984    fmpq_poly_canonicalize(rop, NULL);
     2985   
     2986    mpq_clear(x);
     2987}
     2988
     2989/**
     2990 * \brief    Rescales the co-ordinate.
     2991 * \ingroup  Composition
     2992 *
     2993 * Denoting this polynomial \f$a(T)\f$, computes the polynomial \f$a(x T)\f$.
     2994 */
     2995void fmpq_poly_rescale(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x)
     2996{
     2997    ulong numsize, densize;
     2998    ulong limbs;
     2999    fmpz_t num, den;
     3000    fmpz_t coeff, power, t;
     3001    long i, n;
     3002   
     3003    fmpq_poly_set(rop, op);
     3004   
     3005    if (fmpq_poly_length(rop) < 2ul)
     3006        return;
     3007   
     3008    num = fmpz_init(mpz_size(mpq_numref(x)));
     3009    den = fmpz_init(mpz_size(mpq_denref(x)));
     3010    mpz_to_fmpz(num, mpq_numref(x));
     3011    mpz_to_fmpz(den, mpq_denref(x));
     3012    numsize = fmpz_size(num);
     3013    densize = fmpz_size(den);
     3014   
     3015    n = fmpz_poly_degree(rop->num);
     3016   
     3017    if (fmpz_is_one(den))
     3018    {
     3019        coeff = fmpz_init(rop->num->limbs + n * numsize);
     3020        power = fmpz_init(n * numsize);
     3021        t = fmpz_init(rop->num->limbs + n * numsize);
     3022       
     3023        fmpz_set(power, num);
     3024       
     3025        fmpz_poly_get_coeff_fmpz(t, rop->num, 1);
     3026        fmpz_mul(coeff, t, power);
     3027        fmpz_poly_set_coeff_fmpz(rop->num, 1, coeff);
     3028       
     3029        for (i = 2; i <= n; i++)
     3030        {
     3031            fmpz_set(t, power);
     3032            fmpz_mul(power, t, num);
     3033            fmpz_poly_get_coeff_fmpz(t, rop->num, i);
     3034            fmpz_mul(coeff, t, power);
     3035            fmpz_poly_set_coeff_fmpz(rop->num, i, coeff);
     3036        }
     3037    }
     3038    else
     3039    {
     3040        coeff = fmpz_init(rop->num->limbs + n * (numsize + densize));
     3041        power = fmpz_init(n * (numsize + densize));
     3042        limbs = rop->num->limbs + n * (numsize + densize);
     3043        limbs = FLINT_MAX(limbs, fmpz_size(rop->den));
     3044        t = fmpz_init(limbs);
     3045       
     3046        fmpz_pow_ui(power, den, n);
     3047       
     3048        if (fmpz_is_one(rop->den))
     3049        {
     3050            rop->den = fmpz_realloc(rop->den, n * densize);
     3051            fmpz_set(rop->den, power);
     3052        }
     3053        else
     3054        {
     3055            fmpz_set(t, rop->den);
     3056            limbs = n * densize + fmpz_size(rop->den);
     3057            rop->den = fmpz_realloc(rop->den, limbs);
     3058            fmpz_mul(rop->den, power, t);
     3059        }
     3060       
     3061        fmpz_set_ui(power, 1);
     3062        for (i = n - 1; i >= 0; i--)
     3063        {
     3064            fmpz_set(t, power);
     3065            fmpz_mul(power, t, den);
     3066            fmpz_poly_get_coeff_fmpz(t, rop->num, i);
     3067            fmpz_mul(coeff, t, power);
     3068            fmpz_poly_set_coeff_fmpz(rop->num, i, coeff);
     3069        }
     3070       
     3071        fmpz_set_ui(power, 1);
     3072        for (i = 1; i <= n; i++)
     3073        {
     3074            fmpz_set(t, power);
     3075            fmpz_mul(power, t, num);
     3076            fmpz_poly_get_coeff_fmpz(t, rop->num, i);
     3077            fmpz_mul(coeff, t, power);
     3078            fmpz_poly_set_coeff_fmpz(rop->num, i, coeff);
     3079        }
     3080    }
     3081    fmpq_poly_canonicalize(rop, NULL);
     3082    fmpz_clear(num);
     3083    fmpz_clear(den);
     3084    fmpz_clear(coeff);
     3085    fmpz_clear(power);
     3086    fmpz_clear(t);
     3087}
     3088
     3089///////////////////////////////////////////////////////////////////////////////
     3090// Square-free
     3091
     3092/**
     3093 * \brief    Returns whether \c op is squarefree.
     3094 * \ingroup  Squarefree
     3095 *
     3096 * Returns whether \c op is squarefree.
     3097 *
     3098 * By definition, a polynomial is square-free if it is not a multiple of a
     3099 * non-unit factor.  In particular, polynomials up to degree 1 (including)
     3100 * are square-free.
     3101 */
     3102int fmpq_poly_is_squarefree(const fmpq_poly_ptr op)
     3103{
     3104    if (fmpq_poly_length(op) < 3ul)
     3105        return 1;
     3106    else
     3107    {
     3108        int ans;
     3109        fmpz_poly_t prim;
     3110       
     3111        fmpz_poly_init(prim);
     3112        fmpz_poly_primitive_part(prim, op->num);
     3113        ans = fmpz_poly_is_squarefree(prim);
     3114        fmpz_poly_clear(prim);
     3115        return ans;
     3116    }
     3117}
     3118
     3119///////////////////////////////////////////////////////////////////////////////
     3120// Subpolynomials
     3121
     3122/**
     3123 * \brief    Returns a contiguous subpolynomial.
     3124 * \ingroup  Subpolynomials
     3125 *
     3126 * Returns the slice with coefficients from \f$x^i\f$ (including) to
     3127 * \f$x^j\f$ (excluding).
     3128 */
     3129void fmpq_poly_getslice(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong i, ulong j)
     3130{
     3131    ulong k;
     3132   
     3133    j = FLINT_MIN(fmpq_poly_length(op), j);
     3134   
     3135    /* Aliasing */
     3136    if (rop == op)
     3137    {
     3138        for (k = 0; k < i; k++)
     3139            fmpz_poly_set_coeff_ui(rop->num, k, 0);
     3140        for (k = j; k < fmpz_poly_length(rop->num); k++)
     3141            fmpz_poly_set_coeff_ui(rop->num, k, 0);
     3142        fmpq_poly_canonicalize(rop, NULL);
     3143        return;
     3144    }
     3145   
     3146    fmpz_poly_zero(rop->num);
     3147    if (i < j)
     3148    {
     3149        for (k = j; k != i; )
     3150        {
     3151            k--;
     3152            fmpz_poly_set_coeff_fmpz(rop->num, k,
     3153                                     fmpz_poly_get_coeff_ptr(op->num, k));
     3154        }
     3155        rop->den = fmpz_realloc(rop->den, fmpz_size(op->den));
     3156        fmpz_set(rop->den, op->den);
     3157        fmpq_poly_canonicalize(rop, NULL);
     3158    }
     3159    else
     3160    {
     3161        fmpz_set_ui(rop->den, 1);
     3162    }
     3163}
     3164
     3165/**
     3166 * \brief    Shifts this polynomial to the left by \f$n\f$ coefficients.
     3167 * \ingroup  Subpolynomials
     3168 *
     3169 * Notionally multiplies \c op by \f$t^n\f$ and stores the result in \c rop.
     3170 */
     3171void fmpq_poly_left_shift(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong n)
     3172{
     3173    /* XXX:  As a workaround for a bug in FLINT 1.5.1, we need to handle */
     3174    /* the zero polynomial separately.                                   */
     3175    if (fmpq_poly_is_zero(op))
     3176    {
     3177        fmpq_poly_zero(rop);
     3178        return;
     3179    }
     3180   
     3181    if (n == 0ul)
     3182    {
     3183        fmpq_poly_set(rop, op);
     3184        return;
     3185    }
     3186   
     3187    fmpz_poly_left_shift(rop->num, op->num, n);
     3188    if (rop != op)
     3189    {
     3190        rop->den = fmpz_realloc(rop->den, fmpz_size(op->den));
     3191        fmpz_set(rop->den, op->den);
     3192    }
     3193}
     3194
     3195/**
     3196 * \brief    Shifts this polynomial to the right by \f$n\f$ coefficients.
     3197 * \ingroup  Subpolynomials
     3198 *
     3199 * Notationally returns the quotient of floor division of \c rop by \c op.
     3200 */
     3201void fmpq_poly_right_shift(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong n)
     3202{
     3203    /* XXX:  As a workaround for a bug in FLINT 1.5.1, we need to handle */
     3204    /* the zero polynomial separately.                                   */
     3205    if (fmpq_poly_is_zero(op))
     3206    {
     3207        fmpq_poly_zero(rop);
     3208        return;
     3209    }
     3210   
     3211    if (n == 0ul)
     3212    {
     3213        fmpq_poly_set(rop, op);
     3214        return;
     3215    }
     3216   
     3217    fmpz_poly_right_shift(rop->num, op->num, n);
     3218    if (rop != op)
     3219    {
     3220        rop->den = fmpz_realloc(rop->den, fmpz_size(op->den));
     3221        fmpz_set(rop->den, op->den);
     3222    }
     3223    fmpq_poly_canonicalize(rop, NULL);
     3224}
     3225
     3226/**
     3227 * \brief    Truncates this polynomials.
     3228 * \ingroup  Subpolynomials
     3229 *
     3230 * Truncates <tt>op</tt>  modulo \f$x^n\f$ whenever \f$n\f$ is positive. 
     3231 * Returns zero otherwise.
     3232 */
     3233void fmpq_poly_truncate(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong n)
     3234{
     3235    fmpq_poly_set(rop, op);
     3236    if (fmpq_poly_length(rop) > n)
     3237    {
     3238        fmpz_poly_truncate(rop->num, n);
     3239        fmpq_poly_canonicalize(rop, NULL);
     3240    }
     3241}
     3242
     3243/**
     3244 * \brief    Reverses this polynomial.
     3245 * \ingroup  Subpolynomials
     3246 *
     3247 * Reverses the coefficients of \c op - thought of as a polynomial of
     3248 * length \c n - and places the result in <tt>rop</tt>.
     3249 */
     3250void fmpq_poly_reverse(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong n)
     3251{
     3252    ulong len;
     3253    len = fmpq_poly_length(op);
     3254   
     3255    if (n == 0ul || len == 0ul)
     3256    {
     3257        fmpq_poly_zero(rop);
     3258        return;
     3259    }
     3260   
     3261    fmpz_poly_reverse(rop->num, op->num, n);
     3262    if (rop != op)
     3263    {
     3264        rop->den = fmpz_realloc(rop->den, fmpz_size(op->den));
     3265        fmpz_set(rop->den, op->den);
     3266    }
     3267   
     3268    if (n < len)
     3269        fmpq_poly_canonicalize(rop, NULL);
     3270}
     3271
     3272///////////////////////////////////////////////////////////////////////////////
     3273// String conversion
     3274
     3275/**
     3276 * \addtogroup StringConversions
     3277 *
     3278 * The following three methods enable users to construct elements of type
     3279 * \c fmpq_poly_ptr from strings or to obtain string representations of such
     3280 * elements.
     3281 *
     3282 * The format used is based on the FLINT format for integer polynomials of
     3283 * type <tt>fmpz_poly_t</tt>, which we recall first:
     3284 *
     3285 * A non-zero polynomial \f$a_0 + a_1 X + \dotsb + a_n X^n\f$ of length
     3286 * \f$n + 1\f$ is represented by the string <tt>n+1  a_0 a_1 ... a_n</tt>,
     3287 * where there are two space characters following the length and single space
     3288 * characters separating the individual coefficients.  There is no leading or
     3289 * trailing white-space.  In contrast, the zero polynomial is represented by
     3290 * <tt>0</tt>.
     3291 *
     3292 * We adapt this notation for rational polynomials by using the <tt>mpq_t</tt>
     3293 * notation for the coefficients without any additional white-space.
     3294 *
     3295 * There is also a <tt>_pretty</tt> variant available.
     3296 *
     3297 * Note that currently these functions are not optimized for performance and
     3298 * are intended to be used only for debugging purposes or one-off input and
     3299 * output, rather than as a low-level parser.
     3300 */
     3301
     3302/**
     3303 * \ingroup StringConversions
     3304 * \brief Constructs a polynomial from a list of rationals.
     3305 *
     3306 * Given a list of length <tt>n</tt> containing rationals of type
     3307 * <tt>mpq_t</tt>, this method constructs a polynomial with these
     3308 * coefficients, beginning with the constant term.
     3309 */
     3310void _fmpq_poly_from_list(fmpq_poly_ptr rop, mpq_t * a, ulong n)
     3311{
     3312    mpz_t den, t;
     3313    ulong i;
     3314   
     3315    mpz_init(t);
     3316    mpz_init_set_ui(den, 1);
     3317    for (i = 0; i < n; i++)
     3318        mpz_lcm(den, den, mpq_denref(a[i]));
     3319   
     3320    for (i = 0; i < n; i++)
     3321    {
     3322        mpz_divexact(t, den, mpq_denref(a[i]));
     3323        mpz_mul(mpq_numref(a[i]), mpq_numref(a[i]), t);
     3324    }
     3325   
     3326    fmpz_poly_realloc(rop->num, n);
     3327    for (i = n - 1ul; i != -1ul; i--)
     3328        fmpz_poly_set_coeff_mpz(rop->num, i, mpq_numref(a[i]));
     3329   
     3330    if (mpz_cmp_ui(den, 1) == 0)
     3331        fmpz_set_ui(rop->den, 1);
     3332    else
     3333    {
     3334        rop->den = fmpz_realloc(rop->den, mpz_size(den));
     3335        mpz_to_fmpz(rop->den, den);
     3336    }
     3337    mpz_clear(den);
     3338    mpz_clear(t);
     3339}
     3340
     3341/**
     3342 * \ingroup  StringConversions
     3343 *
     3344 * Sets the rational polynomial \c rop to the value specified by the
     3345 * null-terminated string \c str.
     3346 *
     3347 * The behaviour is undefined if the format of the string \c str does not
     3348 * conform to the specification.
     3349 *
     3350 * \todo  In a future version, it would be nice to have this function return
     3351 *        <tt>0</tt> or <tt>1</tt> depending on whether the input format
     3352 *        was correct.  Currently, the method always returns <tt>1</tt>.
     3353 */
     3354int fmpq_poly_from_string(fmpq_poly_ptr rop, const char * str)
     3355{
     3356    mpq_t * a;
     3357    char * strcopy;
     3358    ulong strcopy_len;
     3359    ulong i, j, k, n;
     3360   
     3361    n = atoi(str);
     3362   
     3363    /* Handle the case of the zero polynomial */
     3364    if (n == 0ul)
     3365    {
     3366        fmpq_poly_zero(rop);
     3367        return 1;
     3368    }
     3369   
     3370    /* Compute the maximum length that the copy buffer has to be */
     3371    strcopy_len = 0;
     3372    for (j = 0; str[j] != ' '; j++);
     3373    j += 2;
     3374    for (i = 0; i < n; i++)
     3375    {
     3376        for (k = j; !(str[k] == ' ' || str[k] == '\0'); k++);
     3377        strcopy_len = FLINT_MAX(strcopy_len, k - j + 1);
     3378        j = k + 1;
     3379    }
     3380   
     3381    strcopy = (char *) malloc(strcopy_len * sizeof(char));
     3382    if (!strcopy)
     3383    {
     3384        printf("ERROR (fmpq_poly_from_string).  Memory allocation failed.\n");
     3385        abort();
     3386    }
     3387   
     3388    /* Read the data into the array a of mpq_t's */
     3389    a = (mpq_t *) malloc(n * sizeof(mpq_t));
     3390    for (j = 0; str[j] != ' '; j++);
     3391    j += 2;
     3392    for (i = 0; i < n; i++)
     3393    {
     3394        for (k = j; !(str[k] == ' ' || str[k] == '\0'); k++);
     3395        memcpy(strcopy, str + j, k - j + 1);
     3396        strcopy[k - j] = '\0';
     3397       
     3398        mpq_init(a[i]);
     3399        mpq_set_str(a[i], strcopy, 10);
     3400        mpq_canonicalize(a[i]);
     3401        j = k + 1;
     3402    }
     3403   
     3404    _fmpq_poly_from_list(rop, a, n);
     3405   
     3406    /* Clean-up */
     3407    free(strcopy);
     3408    for (i = 0; i < n; i++)
     3409        mpq_clear(a[i]);
     3410    free(a);
     3411   
     3412    return 1;
     3413}
     3414
     3415/**
     3416 * \ingroup  StringConversions
     3417 *
     3418 * Returns the string representation of the rational polynomial \c op.
     3419 */
     3420char * fmpq_poly_to_string(const fmpq_poly_ptr op)
     3421{
     3422    ulong i, j;
     3423    ulong len;     /* Upper bound on the length          */
     3424    ulong denlen;  /* Size of the denominator in base 10 */
     3425    mpz_t z;
     3426    mpq_t q;
     3427   
     3428    char * str;
     3429   
     3430    if (fmpq_poly_is_zero(op))
     3431    {
     3432        str = (char *) malloc(2 * sizeof(char));
     3433        if (!str)
     3434        {
     3435            printf("ERROR (fmpq_poly_to_string).  Memory allocation failed.\n");
     3436            abort();
     3437        }
     3438        str[0] = '0';
     3439        str[1] = '\0';
     3440        return str;
     3441    }
     3442   
     3443    mpz_init(z);
     3444    if (fmpz_is_one(op->den))
     3445    {
     3446        denlen = 0;
     3447    }
     3448    else
     3449    {
     3450        fmpz_to_mpz(z, op->den);
     3451        denlen = mpz_sizeinbase(z, 10);
     3452    }
     3453    len = fmpq_poly_places(fmpq_poly_length(op)) + 2;
     3454    for (i = 0; i < fmpq_poly_length(op); i++)
     3455    {
     3456        fmpz_poly_get_coeff_mpz(z, op->num, i);
     3457        len += mpz_sizeinbase(z, 10) + 1;
     3458        if (mpz_sgn(z) != 0)
     3459            len += 2 + denlen;
     3460    }
     3461   
     3462    mpq_init(q);
     3463    str = (char *) malloc(len * sizeof(char));
     3464    if (!str)
     3465    {
     3466        printf("ERROR (fmpq_poly_to_string).  Memory allocation failed.\n");
     3467        abort();
     3468    }
     3469    sprintf(str, "%lu", fmpq_poly_length(op));
     3470    for (j = 0; str[j] != '\0'; j++);
     3471    str[j++] = ' ';
     3472    for (i = 0; i < fmpq_poly_length(op); i++)
     3473    {
     3474        str[j++] = ' ';
     3475        fmpz_poly_get_coeff_mpz(mpq_numref(q), op->num, i);
     3476        if (op->den == NULL)
     3477            mpz_set_ui(mpq_denref(q), 1);
     3478        else
     3479            fmpz_to_mpz(mpq_denref(q), op->den);
     3480        mpq_canonicalize(q);
     3481        mpq_get_str(str + j, 10, q);
     3482        for ( ; str[j] != '\0'; j++);
     3483    }
     3484   
     3485    mpq_clear(q);
     3486    mpz_clear(z);
     3487   
     3488    return str;
     3489}
     3490
     3491/**
     3492 * \ingroup  StringConversions
     3493 *
     3494 * Returns the pretty string representation of \c op.
     3495 *
     3496 * Returns the pretty string representation of the rational polynomial \c op,
     3497 * using the string \c var as the variable name.
     3498 */
     3499char * fmpq_poly_to_string_pretty(const fmpq_poly_ptr op, const char * var)
     3500{
     3501    long i;
     3502    ulong j;
     3503    ulong len;     /* Upper bound on the length          */
     3504    ulong denlen;  /* Size of the denominator in base 10 */
     3505    ulong varlen;  /* Length of the variable name        */
     3506    mpz_t z;       /* op->den (if this is not 1)         */
     3507    mpq_t q;
     3508    char * str;
     3509   
     3510    if (fmpq_poly_is_zero(op))
     3511    {
     3512        str = (char *) malloc(2 * sizeof(char));
     3513        if (!str)
     3514        {
     3515            printf("ERROR (fmpq_poly_to_string_pretty).  Memory allocation failed.\n");
     3516            abort();
     3517        }
     3518        str[0] = '0';
     3519        str[1] = '\0';
     3520        return str;
     3521    }
     3522   
     3523    if (fmpq_poly_length(op) == 1ul)
     3524    {
     3525        mpq_init(q);
     3526        fmpz_to_mpz(mpq_numref(q), fmpz_poly_lead(op->num));
     3527        if (fmpz_is_one(op->den))
     3528            mpz_set_ui(mpq_denref(q), 1);
     3529        else
     3530        {
     3531            fmpz_to_mpz(mpq_denref(q), op->den);
     3532            mpq_canonicalize(q);
     3533        }
     3534        str = mpq_get_str(NULL, 10, q);
     3535        mpq_clear(q);
     3536        return str;
     3537    }
     3538   
     3539    varlen = strlen(var);
     3540   
     3541    /* Copy the denominator into an mpz_t */
     3542    mpz_init(z);
     3543    if (fmpz_is_one(op->den))
     3544    {
     3545        denlen = 0;
     3546    }
     3547    else
     3548    {
     3549        fmpz_to_mpz(z, op->den);
     3550        denlen = mpz_sizeinbase(z, 10);
     3551    }
     3552   
     3553    /* Estimate the length */
     3554    len = 0;
     3555    for (i = 0; i < fmpq_poly_length(op); i++)
     3556    {
     3557        fmpz_poly_get_coeff_mpz(z, op->num, i);
     3558        len += mpz_sizeinbase(z, 10);            /* Numerator                */
     3559        if (mpz_sgn(z) != 0)
     3560            len += 1 + denlen;                   /* Denominator and /        */
     3561        len += 3;                                /* Operator and whitespace  */
     3562        len += 1 + varlen + 1;                   /* *, x and ^               */
     3563        len += fmpq_poly_places(i);              /* Exponent                 */
     3564    }
     3565   
     3566    mpq_init(q);
     3567    str = (char *) malloc(len * sizeof(char));
     3568    if (!str)
     3569    {
     3570        printf("ERROR (fmpq_poly_to_string_pretty).  Memory allocation failed.\n");
     3571        abort();
     3572    }
     3573    j = 0;
     3574   
     3575    /* Print the leading term */
     3576    fmpz_to_mpz(mpq_numref(q), fmpz_poly_lead(op->num));
     3577    fmpz_to_mpz(mpq_denref(q), op->den);
     3578    mpq_canonicalize(q);
     3579   
     3580    if (mpq_cmp_ui(q, 1, 1) != 0)
     3581    {
     3582        if (mpq_cmp_si(q, -1, 1) == 0)
     3583            str[j++] = '-';
     3584        else
     3585        {
     3586            mpq_get_str(str, 10, q);
     3587            for ( ; str[j] != '\0'; j++);
     3588            str[j++] = '*';
     3589        }
     3590    }
     3591    sprintf(str + j, "%s", var);
     3592    j += varlen;
     3593    if (fmpz_poly_degree(op->num) > 1)
     3594    {
     3595        str[j++] = '^';
     3596        sprintf(str + j, "%li", fmpz_poly_degree(op->num));
     3597        for ( ; str[j] != '\0'; j++);
     3598    }
     3599   
     3600    for (i = fmpq_poly_degree(op) - 1; i >= 0; i--)
     3601    {
     3602        if (fmpz_is_zero(fmpz_poly_get_coeff_ptr(op->num, i)))
     3603            continue;
     3604       
     3605        fmpz_poly_get_coeff_mpz(mpq_numref(q), op->num, i);
     3606        fmpz_to_mpz(mpq_denref(q), op->den);
     3607        mpq_canonicalize(q);
     3608       
     3609        str[j++] = ' ';
     3610        if (mpq_sgn(q) < 0)
     3611        {
     3612            mpq_abs(q, q);
     3613            str[j++] = '-';
     3614        }
     3615        else
     3616            str[j++] = '+';
     3617        str[j++] = ' ';
     3618       
     3619        mpq_get_str(str + j, 10, q);
     3620        for ( ; str[j] != '\0'; j++);
     3621       
     3622        if (i > 0)
     3623        {
     3624            str[j++] = '*';
     3625            sprintf(str + j, "%s", var);
     3626            j += varlen;
     3627            if (i > 1)
     3628            {
     3629                str[j++] = '^';
     3630                sprintf(str + j, "%li", i);
     3631                for ( ; str[j] != '\0'; j++);
     3632            }
     3633        }
     3634    }
     3635   
     3636    mpq_clear(q);
     3637    mpz_clear(z);
     3638    return str;
     3639}
     3640
     3641#ifdef __cplusplus
     3642    }
     3643#endif
     3644
  • new file sage/libs/flint/fmpq_poly.h

    diff --git a/sage/libs/flint/fmpq_poly.h b/sage/libs/flint/fmpq_poly.h
    new file mode 100644
    - +  
     1///////////////////////////////////////////////////////////////////////////////
     2// Copyright (C) 2010 Sebastian Pancratz                                     //
     3//                                                                           //
     4// Distributed under the terms of the GNU General Public License as          //
     5// published by the Free Software Foundation; either version 2 of the        //
     6// License, or (at your option) any later version.                           //
     7//                                                                           //
     8// http://www.gnu.org/licenses/                                              //
     9///////////////////////////////////////////////////////////////////////////////
     10
     11#ifndef _FMPQ_POLY_H_
     12#define _FMPQ_POLY_H_
     13
     14#ifdef __cplusplus
     15    extern "C" {
     16#endif
     17
     18#include <stdio.h>
     19#include <gmp.h>
     20
     21#include "fmpz.h"
     22#include "fmpz_poly.h"
     23
     24/**
     25 * \file     fmpq_poly.h
     26 * \brief    Fast implementation of the rational polynomial ring
     27 * \author   Sebastian Pancratz
     28 * \date     Mar 2010 -- July 2010
     29 * \version  0.1.8
     30 */
     31
     32/**
     33 * \ingroup  Definitions
     34 * \brief    Data type for a rational polynomial.
     35 * \details  We represent a rational polynomial as the quotient of an integer
     36 *           polynomial of type <tt>fmpz_poly_t</tt> and a denominator of type
     37 *           <tt>fmpz_t</tt>, enforcing coprimality and that the denominator
     38 *           is positive.  The zero polynomial is represented as \f$0/1\f$.
     39 */
     40typedef struct
     41{
     42    fmpz_poly_t num;  //!< Numerator
     43    fmpz_t den;       //!< Denominator
     44} fmpq_poly_struct;
     45
     46/**
     47 * \ingroup  Definitions
     48 * \brief    Array of #fmpq_poly_struct of length one.
     49 * \details  This allows passing <em>by reference</em> without having to
     50 *           explicitly allocate memory for the structure, as one would have
     51 *           to when using pointers.
     52 */
     53typedef fmpq_poly_struct fmpq_poly_t[1];
     54
     55/**
     56 * \ingroup  Definitions
     57 * \brief    Pointer to #fmpq_poly_struct.
     58 */
     59typedef fmpq_poly_struct * fmpq_poly_ptr;
     60
     61void fmpq_poly_canonicalize(fmpq_poly_ptr f, fmpz_t temp);
     62
     63///////////////////////////////////////////////////////////////////////////////
     64// fmpq_poly_*                                                               //
     65///////////////////////////////////////////////////////////////////////////////
     66
     67///////////////////////////////////////////////////////////////////////////////
     68// Accessing numerator and denominator
     69
     70/**
     71 * \def      fmpq_poly_numref(op)
     72 * \brief    Returns a reference to the numerator of \c op.
     73 * \ingroup  NumDen
     74 * \details  The <tt>fmpz_poly</tt> functions can be used on the result of
     75 *           this macro.
     76 */
     77#define fmpq_poly_numref(op)  ((op)->num)
     78
     79/**
     80 * \def      fmpq_poly_denref(op)
     81 * \brief    Returns a reference to the denominator of \c op.
     82 * \ingroup  NumDen
     83 * \details  The <tt>fmpz</tt> functions can be used on the result of
     84 *           this macro.
     85 */
     86#define fmpq_poly_denref(op)  ((op)->den)
     87
     88///////////////////////////////////////////////////////////////////////////////
     89// Memory management
     90
     91/**
     92 * \ingroup  MemoryManagement
     93 *
     94 * Initializes the element \c rop to zero.
     95 *
     96 * This function should be called once for the #fmpq_poly_ptr \c rop, before
     97 * using it with other <tt>fmpq_poly</tt> functions, or following a
     98 * preceeding call to #fmpq_poly_clear().
     99 */
     100static inline
     101void fmpq_poly_init(fmpq_poly_ptr rop)
     102{
     103    fmpz_poly_init(rop->num);
     104    rop->den = fmpz_init(1ul);
     105    fmpz_set_ui(rop->den, 1ul);
     106}
     107
     108/**
     109 * \ingroup  MemoryManagement
     110 *
     111 * Clears the element \c rop.
     112 *
     113 * This function should only be called on an element which has been
     114 * initialised.
     115 */
     116static inline
     117void fmpq_poly_clear(fmpq_poly_ptr rop)
     118{
     119    fmpz_poly_clear(rop->num);
     120    fmpz_clear(rop->den);
     121}
     122
     123///////////////////////////////////////////////////////////////////////////////
     124// Assignment and basic manipulation
     125
     126/**
     127 * \ingroup  Assignment
     128 *
     129 * Sets the element \c rop to the same value as the element \c op.
     130 */
     131static inline
     132void fmpq_poly_set(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
     133{
     134    if (rop != op)
     135    {
     136        fmpz_poly_set(rop->num, op->num);
     137       
     138        rop->den = fmpz_realloc(rop->den, fmpz_size(op->den));
     139        fmpz_set(rop->den, op->den);
     140    }
     141}
     142
     143/**
     144 * \ingroup  Assignment
     145 *
     146 * Sets the element \c rop to the value given by the \c int \c op.
     147 */
     148static inline 
     149void fmpq_poly_set_si(fmpq_poly_ptr rop, long op)
     150{
     151    fmpz_poly_zero(rop->num);
     152    fmpz_poly_set_coeff_si(rop->num, 0, op);
     153    fmpz_set_ui(rop->den, 1ul);
     154}
     155
     156/**
     157 * \ingroup  Assignment
     158 *
     159 * Sets the element \c rop to the integer \c x.
     160 */
     161static inline
     162void fmpq_poly_set_fmpz(fmpq_poly_ptr rop, const fmpz_t x)
     163{
     164    fmpz_poly_zero(rop->num);
     165    fmpz_poly_set_coeff_fmpz(rop->num, 0, x);
     166    fmpz_set_ui(rop->den, 1ul);
     167}
     168
     169/**
     170 * \ingroup  Assignment
     171 *
     172 * Sets the element \c rop to the integer \c x.
     173 */
     174static inline
     175void fmpq_poly_set_mpz(fmpq_poly_ptr rop, const mpz_t x)
     176{
     177    fmpz_poly_zero(rop->num);
     178    fmpz_poly_set_coeff_mpz(rop->num, 0, x);
     179    fmpz_set_ui(rop->den, 1ul);
     180}
     181
     182/**
     183 * \ingroup  Assignment
     184 *
     185 * Sets the element \c rop to the rational \c x, assumed to be given
     186 * in lowest terms.
     187 */
     188static inline
     189void fmpq_poly_set_mpq(fmpq_poly_ptr rop, const mpq_t x)
     190{
     191    fmpz_poly_zero(rop->num);
     192    fmpz_poly_set_coeff_mpz(rop->num, 0, mpq_numref(x));
     193    rop->den = fmpz_realloc(rop->den, mpz_size(mpq_denref(x)));
     194    mpz_to_fmpz(rop->den, mpq_denref(x));
     195}
     196
     197/**
     198 * \ingroup  Assignment
     199 *
     200 * Swaps the elements \c op1 and \c op2.
     201 *
     202 * This is done efficiently by swapping pointers.
     203 */
     204static inline
     205void fmpq_poly_swap(fmpq_poly_ptr op1, fmpq_poly_ptr op2)
     206{
     207    if (op1 != op2)
     208    {
     209        fmpz_t t;
     210        fmpz_poly_swap(op1->num, op2->num);
     211        t        = op1->den;
     212        op1->den = op2->den;
     213        op2->den = t;
     214    }
     215}
     216
     217/**
     218 * \ingroup  Assignment
     219 *
     220 * Sets the element \c rop to zero.
     221 */
     222static inline
     223void fmpq_poly_zero(fmpq_poly_ptr rop)
     224{
     225    fmpz_poly_zero(rop->num);
     226    fmpz_set_ui(rop->den, 1ul);
     227}
     228
     229/**
     230 * \ingroup  Assignment
     231 *
     232 * Sets the element \c rop to one.
     233 */
     234static inline
     235void fmpq_poly_one(fmpq_poly_ptr rop)
     236{
     237    fmpz_poly_zero(rop->num);
     238    fmpz_poly_set_coeff_si(rop->num, 0, 1);
     239    fmpz_set_ui(rop->den, 1ul);
     240}
     241
     242void fmpq_poly_neg(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
     243void fmpq_poly_inv(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
     244
     245///////////////////////////////////////////////////////////////////////////////
     246// Setting/ retrieving individual coefficients
     247
     248void fmpq_poly_get_coeff_mpq(mpq_t rop, const fmpq_poly_ptr op, ulong i);
     249
     250void fmpq_poly_set_coeff_fmpz(fmpq_poly_ptr rop, ulong i, const fmpz_t x);
     251void fmpq_poly_set_coeff_mpz(fmpq_poly_ptr rop, ulong i, const mpz_t x);
     252void fmpq_poly_set_coeff_mpq(fmpq_poly_ptr rop, ulong i, const mpq_t x);
     253void fmpq_poly_set_coeff_si(fmpq_poly_ptr rop, ulong i, long x);
     254
     255///////////////////////////////////////////////////////////////////////////////
     256// Comparison
     257
     258/**
     259 * \brief    Returns whether \c op is zero.
     260 * \ingroup  Comparison
     261 *
     262 * Returns whether the element \c op is zero.
     263 */
     264static inline
     265int fmpq_poly_is_zero(const fmpq_poly_ptr op)
     266{
     267    return op->num->length == 0ul;
     268}
     269
     270/**
     271 * \brief    Returns whether \c op is equal to \f$1\f$.
     272 * \ingroup  Comparison
     273 *
     274 * Returns whether the element \c op is equal to the constant polynomial
     275 * \f$1\f$.
     276 */
     277static inline
     278int fmpq_poly_is_one(const fmpq_poly_ptr op)
     279{
     280    return (op->num->length == 1ul) && fmpz_is_one(op->num->coeffs)
     281                                    && fmpz_is_one(op->den);
     282}
     283
     284int fmpq_poly_equal(const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
     285int fmpq_poly_cmp(const fmpq_poly_ptr left, const fmpq_poly_ptr right);
     286
     287///////////////////////////////////////////////////////////////////////////////
     288// Polynomial parameters
     289
     290/**
     291 * \brief    Returns the length of \c op.
     292 * \ingroup  PolyParameters
     293 * \details  Returns the length of the polynomial \c op, which is one greater
     294 *           than its degree.
     295 */
     296static inline
     297ulong fmpq_poly_length(const fmpq_poly_ptr op)
     298{
     299    return op->num->length;
     300}
     301
     302/**
     303 * \brief    Returns the degree of \c op.
     304 * \ingroup  Parameters
     305 * \details  Returns the degree of the polynomial \c op.
     306 */
     307static inline
     308long fmpq_poly_degree(const fmpq_poly_ptr op)
     309{
     310    return (long) (op->num->length - 1ul);
     311}
     312
     313///////////////////////////////////////////////////////////////////////////////
     314// Addition/ subtraction
     315
     316void fmpq_poly_add(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
     317void fmpq_poly_sub(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
     318
     319void fmpq_poly_addmul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
     320void fmpq_poly_submul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
     321
     322///////////////////////////////////////////////////////////////////////////////
     323// Scalar multiplication and division
     324
     325void fmpq_poly_scalar_mul_si(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long x);
     326void fmpq_poly_scalar_mul_mpz(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpz_t x);
     327void fmpq_poly_scalar_mul_mpq(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x);
     328
     329void fmpq_poly_scalar_div_si(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long x);
     330void fmpq_poly_scalar_div_mpz(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpz_t x);
     331void fmpq_poly_scalar_div_mpq(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x);
     332
     333///////////////////////////////////////////////////////////////////////////////
     334// Multiplication
     335
     336void fmpq_poly_mul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
     337
     338///////////////////////////////////////////////////////////////////////////////
     339// Division
     340
     341void fmpq_poly_floordiv(fmpq_poly_ptr q, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
     342void fmpq_poly_mod(fmpq_poly_ptr r, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
     343void fmpq_poly_divrem(fmpq_poly_ptr q, fmpq_poly_ptr r, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
     344
     345///////////////////////////////////////////////////////////////////////////////
     346// Powering
     347
     348void fmpq_poly_power(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong exp);
     349
     350///////////////////////////////////////////////////////////////////////////////
     351// Greatest common divisor
     352
     353void fmpq_poly_gcd(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
     354void 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);
     355void fmpq_poly_lcm(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
     356
     357///////////////////////////////////////////////////////////////////////////////
     358// Derivative
     359
     360void fmpq_poly_derivative(fmpq_poly_ptr rop, fmpq_poly_ptr op);
     361
     362///////////////////////////////////////////////////////////////////////////////
     363// Evaluation
     364
     365void fmpq_poly_evaluate_mpz(mpq_t rop, fmpq_poly_ptr f, const mpz_t a);
     366void fmpq_poly_evaluate_mpq(mpq_t rop, fmpq_poly_ptr f, const mpq_t a);
     367
     368///////////////////////////////////////////////////////////////////////////////
     369// Gaussian content
     370
     371void fmpq_poly_content(mpq_t rop, const fmpq_poly_ptr op);
     372void fmpq_poly_primitive_part(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
     373int fmpq_poly_is_monic(const fmpq_poly_ptr op);
     374void fmpq_poly_monic(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
     375
     376///////////////////////////////////////////////////////////////////////////////
     377// Resultant
     378
     379void fmpq_poly_resultant(mpq_t rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
     380void fmpq_poly_discriminant(mpq_t disc, fmpq_poly_t a);
     381
     382///////////////////////////////////////////////////////////////////////////////
     383// Composition
     384
     385void fmpq_poly_compose(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
     386void fmpq_poly_rescale(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x);
     387
     388///////////////////////////////////////////////////////////////////////////////
     389// Square-free
     390
     391int fmpq_poly_is_squarefree(const fmpq_poly_ptr op);
     392
     393///////////////////////////////////////////////////////////////////////////////
     394// Subpolynomials
     395
     396void fmpq_poly_getslice(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong i, ulong j);
     397void fmpq_poly_left_shift(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong n);
     398void fmpq_poly_right_shift(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong n);
     399void fmpq_poly_truncate(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong n);
     400void fmpq_poly_reverse(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong n);
     401
     402///////////////////////////////////////////////////////////////////////////////
     403// String conversion
     404
     405void _fmpq_poly_from_list(fmpq_poly_ptr rop, mpq_t * a, ulong n);
     406int fmpq_poly_from_string(fmpq_poly_ptr rop, const char * s);
     407char * fmpq_poly_to_string(const fmpq_poly_ptr op);
     408char * fmpq_poly_to_string_pretty(const fmpq_poly_ptr op, const char * x);
     409
     410#ifdef __cplusplus
     411    }
     412#endif
     413
     414#endif
     415
  • new file sage/libs/flint/fmpq_poly.pxd

    diff --git a/sage/libs/flint/fmpq_poly.pxd b/sage/libs/flint/fmpq_poly.pxd
    new file mode 100644
    - +  
     1###############################################################################
     2#          Copyright (C) 2010 Sebastian Pancratz <sfp@pancratz.org>           #
     3#                                                                             #
     4#     Distributed under the terms of the GNU General Public License (GPL)     #
     5#                                                                             #
     6#                        http://www.gnu.org/licenses/                         #
     7###############################################################################
     8
     9#include "fmpz_poly.pxi"
     10#include "fmpz.pxi"
     11
     12cdef extern from "gmp.h":
     13    ctypedef void * mpz_t
     14    ctypedef void * mpq_t
     15
     16cdef extern from "fmpq_poly.h":
     17    ctypedef void * fmpz_t
     18    ctypedef void * fmpz_poly_p
     19    struct fmpq_poly:
     20        fmpz_poly_p num
     21        fmpz_t den
     22       
     23    ctypedef fmpq_poly fmpq_poly_struct
     24    ctypedef fmpq_poly_struct fmpq_poly_t[1]
     25
     26    void * fmpq_poly_canonicalize(fmpq_poly_t, fmpz_t)
     27
     28    void * fmpq_poly_numref(fmpq_poly_t)
     29    void * fmpq_poly_denref(fmpq_poly_t)
     30
     31    void fmpq_poly_init(fmpq_poly_t)
     32    void fmpq_poly_clear(fmpq_poly_t)
     33
     34    void fmpq_poly_set(fmpq_poly_t, fmpq_poly_t)
     35    void fmpq_poly_set_si(fmpq_poly_t, long)
     36    void fmpq_poly_set_mpz(fmpq_poly_t, mpz_t)
     37    void fmpq_poly_set_mpq(fmpq_poly_t, mpq_t)
     38    void fmpq_poly_swap(fmpq_poly_t, fmpq_poly_t)
     39    void fmpq_poly_zero(fmpq_poly_t)
     40    void fmpq_poly_neg(fmpq_poly_t, fmpq_poly_t)
     41
     42    void fmpq_poly_get_coeff_mpq(mpq_t, fmpq_poly_t, unsigned long)
     43    void fmpq_poly_set_coeff_si(fmpq_poly_t, unsigned long, long)
     44    void fmpq_poly_set_coeff_mpq(fmpq_poly_t, unsigned long, mpq_t)
     45    void fmpq_poly_set_coeff_mpz(fmpq_poly_t, unsigned long, mpz_t)
     46   
     47    int fmpq_poly_equal(fmpq_poly_t, fmpq_poly_t)
     48    int fmpq_poly_cmp(fmpq_poly_t, fmpq_poly_t)
     49    int fmpq_poly_is_zero(fmpq_poly_t)
     50
     51    long fmpq_poly_degree(fmpq_poly_t)
     52    unsigned long fmpq_poly_length(fmpq_poly_t)
     53   
     54    void fmpq_poly_add(fmpq_poly_t, fmpq_poly_t, fmpq_poly_t)
     55    void fmpq_poly_sub(fmpq_poly_t, fmpq_poly_t, fmpq_poly_t)
     56   
     57    void fmpq_poly_scalar_mul_mpq(fmpq_poly_t, fmpq_poly_t, mpq_t)
     58    void fmpq_poly_scalar_div_mpq(fmpq_poly_t, fmpq_poly_t, mpq_t)
     59   
     60    void fmpq_poly_mul(fmpq_poly_t, fmpq_poly_t, fmpq_poly_t)
     61
     62    void fmpq_poly_divrem(fmpq_poly_t, fmpq_poly_t, fmpq_poly_t, fmpq_poly_t)
     63    void fmpq_poly_floordiv(fmpq_poly_t, fmpq_poly_t, fmpq_poly_t)
     64    void fmpq_poly_mod(fmpq_poly_t, fmpq_poly_t, fmpq_poly_t)
     65
     66    void fmpq_poly_power(fmpq_poly_t, fmpq_poly_t, unsigned long)
     67
     68    void fmpq_poly_left_shift(fmpq_poly_t, fmpq_poly_t, unsigned long)
     69    void fmpq_poly_right_shift(fmpq_poly_t, fmpq_poly_t, unsigned long)
     70   
     71    void fmpq_poly_gcd(fmpq_poly_t, fmpq_poly_t, fmpq_poly_t)
     72    void fmpq_poly_xgcd(fmpq_poly_t, fmpq_poly_t, fmpq_poly_t, fmpq_poly_t, fmpq_poly_t)
     73    void fmpq_poly_lcm(fmpq_poly_t, fmpq_poly_t, fmpq_poly_t)
     74
     75    void fmpq_poly_derivative(fmpq_poly_t, fmpq_poly_t)
     76
     77    void fmpq_poly_evaluate_mpz(mpq_t, fmpq_poly_t, mpz_t)
     78    void fmpq_poly_evaluate_mpq(mpq_t, fmpq_poly_t, mpq_t)
     79
     80    void fmpq_poly_content(mpq_t, fmpq_poly_t)
     81    void fmpq_poly_primitive_part(fmpq_poly_t, fmpq_poly_t)
     82
     83    void fmpq_poly_resultant(mpq_t, fmpq_poly_t, fmpq_poly_t)
     84
     85    void fmpq_poly_compose(fmpq_poly_t, fmpq_poly_t, fmpq_poly_t)
     86
     87    void fmpq_poly_getslice(fmpq_poly_t, fmpq_poly_t, unsigned long, unsigned long)
     88    void fmpq_poly_truncate(fmpq_poly_t, fmpq_poly_t, unsigned long)
     89    void fmpq_poly_reverse(fmpq_poly_t, fmpq_poly_t, unsigned long)
     90
     91    void _fmpq_poly_from_list(fmpq_poly_t, mpq_t *, unsigned long)
     92    void fmpq_poly_from_string(fmpq_poly_t, char *)
     93    char * fmpq_poly_to_string(fmpq_poly_t, char *)
     94    char * fmpq_poly_to_string_pretty(fmpq_poly_t, char *)
     95
  • sage/libs/flint/fmpz.pxi

    diff --git a/sage/libs/flint/fmpz.pxi b/sage/libs/flint/fmpz.pxi
    a b include "../ntl/decl.pxi" 
    22
    33cdef extern from "FLINT/fmpz.h":
    44   
    5     ctypedef void* fmpz_t
     5    ctypedef void * fmpz_t
     6    ctypedef void * mpz_t
    67
    78    fmpz_t fmpz_init(unsigned long limbs)
    89    void fmpz_clear(fmpz_t f)
    910    void fmpz_print(fmpz_t f)
    1011    int fmpz_is_one(fmpz_t f)
     12   
     13    void fmpz_to_mpz(mpz_t rop, fmpz_t op)
  • sage/misc/explain_pickle.py

    diff --git a/sage/misc/explain_pickle.py b/sage/misc/explain_pickle.py
    a b EXAMPLES:: 
    1717    pg_make_integer = unpickle_global('sage.rings.integer', 'make_integer')
    1818    pg_make_integer('c1p')
    1919    sage: explain_pickle(dumps(polygen(QQ)))
    20     pg_Polynomial_rational_dense = unpickle_global('sage.rings.polynomial.polynomial_element_generic', 'Polynomial_rational_dense')
     20    pg_Polynomial_rational_flint = unpickle_global('sage.rings.polynomial.polynomial_rational_flint', 'Polynomial_rational_flint')
    2121    pg_PolynomialRing = unpickle_global('sage.rings.polynomial.polynomial_ring_constructor', 'PolynomialRing')
    2222    pg_RationalField = unpickle_global('sage.rings.rational_field', 'RationalField')
    2323    pg = unpickle_instantiate(pg_RationalField, ())
    2424    pg_make_rational = unpickle_global('sage.rings.rational', 'make_rational')
    25     pg_Polynomial_rational_dense(pg_PolynomialRing(pg, 'x', None, False), [pg_make_rational('0'), pg_make_rational('1')], False, True)
     25    pg_Polynomial_rational_flint(pg_PolynomialRing(pg, 'x', None, False), [pg_make_rational('0'), pg_make_rational('1')], False, True)
    2626    sage: sage_eval(explain_pickle(dumps(polygen(QQ)))) == polygen(QQ)
    2727    True
    2828
    version of Sage; here are the above two  
    3939    from sage.rings.integer import make_integer
    4040    make_integer('c1p')
    4141    sage: explain_pickle(dumps(polygen(QQ)), in_current_sage=True)
    42     from sage.rings.polynomial.polynomial_element_generic import Polynomial_rational_dense
     42    from sage.rings.polynomial.polynomial_rational_flint import Polynomial_rational_flint
    4343    from sage.rings.rational import make_rational
    44     Polynomial_rational_dense(PolynomialRing(RationalField(), 'x', None, False), [make_rational('0'), make_rational('1')], False, True)
     44    Polynomial_rational_flint(PolynomialRing(RationalField(), 'x', None, False), [make_rational('0'), make_rational('1')], False, True)
    4545
    4646The explain_pickle function has several use cases.
    4747
  • sage/rings/arith.py

    diff --git a/sage/rings/arith.py b/sage/rings/arith.py
    a b def __LCM_sequence(v): 
    15731573        2*X^3 + 4*X^2
    15741574        sage: R.<X>=QQ[]
    15751575        sage: __LCM_sequence(Sequence((2*X+4,2*X^2,2)))
    1576         4*X^3 + 8*X^2
     1576        X^3 + 2*X^2
    15771577    """
    15781578    if len(v) == 0:
    15791579        return ZZ(1)
  • sage/rings/complex_field.py

    diff --git a/sage/rings/complex_field.py b/sage/rings/complex_field.py
    a b class ComplexField_class(field.Field): 
    126126        sage: C(S.gen())
    127127        Traceback (most recent call last):
    128128        ...
    129         TypeError: unable to coerce to a ComplexNumber: <class 'sage.rings.polynomial.polynomial_element_generic.Polynomial_rational_dense'>
     129        TypeError: unable to coerce to a ComplexNumber: <type 'sage.rings.polynomial.polynomial_rational_flint.Polynomial_rational_flint'>
    130130   
    131131    This illustrates precision.
    132132   
  • sage/rings/number_field/number_field_ideal.py

    diff --git a/sage/rings/number_field/number_field_ideal.py b/sage/rings/number_field/number_field_ideal.py
    a b class NumberFieldIdeal(Ideal_generic): 
    192192        EXAMPLES::
    193193
    194194            sage: NumberField(x^2 + 1, 'a').ideal(7).__hash__()
    195             -9223372036854775779                 # 64-bit
    196             -2147483619                          # 32-bit
     195            -288230376151711715                 # 64-bit
     196            -67108835                           # 32-bit
    197197        """
    198198        try: return self._hash
    199199        # At some point in the future (e.g., for relative extensions), we'll likely
  • sage/rings/polynomial/polynomial_element.pyx

    diff --git a/sage/rings/polynomial/polynomial_element.pyx b/sage/rings/polynomial/polynomial_element.pyx
    a b cdef class Polynomial(CommutativeAlgebra 
    470470       
    471471        ::
    472472       
    473             sage: R.<x> = QQ[]
     473            sage: R.<x> = PolynomialRing(QQ, sparse=True)
    474474            sage: f = x^3-2*x
    475475            sage: f(x) is f
    476476            True
    cdef class Polynomial(CommutativeAlgebra 
    14841484       
    14851485        EXAMPLES::
    14861486       
    1487             sage: x = polygen(QQ)
     1487            sage: R.<x> = PolynomialRing(QQ, implementation="FLINT")
    14881488            sage: f = x^3+2/3*x^2 - 5/3
    14891489            sage: f._repr_()
    14901490            'x^3 + 2/3*x^2 - 5/3'
    14911491            sage: f.rename('vaughn')
    1492             sage: f
    1493             vaughn
     1492            Traceback (most recent call last):
     1493            ...
     1494            NotImplementedError: object does not support renaming: x^3 + 2/3*x^2 - 5/3
    14941495        """
    14951496        return self._repr()
    14961497
    cdef class Polynomial(CommutativeAlgebra 
    20632064
    20642065        Note that some subclasses may implement their own
    20652066        denominator function. For example, see
    2066         :class:`sage.rings.polynomial.polynomial_element_generic.Polynomial_rational_dense`
     2067        :class:`sage.rings.polynomial.polynomial_rational_flint.Polynomial_rational_flint`
    20672068
    20682069        .. warning::
    20692070
    cdef class Polynomial(CommutativeAlgebra 
    21272128
    21282129        Note that some subclases may implement its own numerator
    21292130        function. For example, see
    2130         :class:`sage.rings.polynomial.polynomial_element_generic.Polynomial_rational_dense`
     2131        :class:`sage.rings.polynomial.polynomial_rational_flint.Polynomial_rational_flint`
    21312132
    21322133        .. warning::
    21332134
    cdef class Polynomial(CommutativeAlgebra 
    49874988    def is_irreducible(self):
    49884989        """
    49894990        Return True precisely if this polynomial is irreducible over its
    4990         base ring. Testing irreducibility over
    4991         `\ZZ/n\ZZ` for composite `n` is not
     4991        base ring.
     4992
     4993        Testing irreducibility over `\ZZ/n\ZZ` for composite `n` is not
    49924994        implemented.
    4993        
    4994         The function returns False for polynomials which are units,
    4995         and raises an exception for the zero polynomial.
    49964995
    49974996        EXAMPLES::
    49984997       
    cdef class Polynomial(CommutativeAlgebra 
    50045003            sage: (x^3 + 2).is_irreducible()
    50055004            True
    50065005            sage: R(0).is_irreducible()
    5007             Traceback (most recent call last):
    5008             ...
    5009             ValueError: self must be nonzero
    5010        
    5011         See \#5140:
     5006            False
     5007       
     5008        See \#5140,
     5009       
     5010        ::
     5011       
    50125012            sage: R(1).is_irreducible()
    50135013            False
    50145014            sage: R(4).is_irreducible()
    cdef class Polynomial(CommutativeAlgebra 
    50165016            sage: R(5).is_irreducible()
    50175017            True
    50185018
    5019         The base ring does matter: for example, 2x is irreducible as a
    5020         polynomial in QQ[x], but not in ZZ[x]:
     5019        The base ring does matter:  for example, 2x is irreducible as a
     5020        polynomial in QQ[x], but not in ZZ[x],
     5021       
     5022        ::
    50215023       
    50225024            sage: R.<x> = ZZ[]
    50235025            sage: R(2*x).is_irreducible()
    cdef class Polynomial(CommutativeAlgebra 
    50385040            True
    50395041        """
    50405042        if self.is_zero():
    5041             raise ValueError, "self must be nonzero"
     5043            return False
    50425044        if self.is_unit():
    50435045            return False
    50445046        if self.degree() == 0:
  • sage/rings/polynomial/polynomial_element_generic.py

    diff --git a/sage/rings/polynomial/polynomial_element_generic.py b/sage/rings/polynomial/polynomial_element_generic.py
    a b class Polynomial_generic_dense_field(Pol 
    604604    def __init__(self, parent, x=None, check=True, is_gen = False, construct=False):
    605605        Polynomial_generic_dense.__init__(self, parent, x, check, is_gen)
    606606
     607############################################################################
     608# XXX:  Ensures that the generic polynomials implemented in SAGE via PARI  #
     609# until at least until 4.5.0 unpickle correctly as polynomials implemented #
     610# via FLINT.                                                               #
     611from sage.structure.sage_object import register_unpickle_override
     612from sage.rings.polynomial.polynomial_rational_flint import Polynomial_rational_flint
    607613
    608 class Polynomial_rational_dense(Polynomial_generic_field):
    609     """
    610     A dense polynomial over the rational numbers.
    611     """
    612     def __init__(self, parent, x=None, check=True, is_gen=False, construct=False):
    613         Polynomial.__init__(self, parent, is_gen=is_gen)
    614 
    615         if construct:
    616             self.__poly = x
    617             return
    618 
    619         self.__poly = pari([]).Polrev()
    620 
    621         if x is None:
    622             return         # leave initialized to 0 polynomial.
    623 
    624 
    625         if fraction_field_element.is_FractionFieldElement(x):
    626             if not x.denominator().is_unit():
    627                 raise TypeError, "denominator must be a unit"
    628             elif x.denominator() != 1:
    629                 x = x.numerator() * x.denominator().inverse_of_unit()
    630             else:
    631                 x = x.numerator()
    632        
    633         if isinstance(x, Polynomial):
    634             if x.parent() == self.parent():
    635                 from copy import copy
    636                 self.__poly = copy(x.__poly)
    637                 return
    638             else:
    639                 x = [QQ(a) for a in x.list()]
    640                 check = False
    641 
    642         elif isinstance(x, dict):
    643             x = self._dict_to_list(x, QQ(0))
    644            
    645            
    646         elif isinstance(x, pari_gen):
    647             f = x.Polrev()
    648             self.__poly = f
    649             assert self.__poly.type() == "t_POL"           
    650             return
    651            
    652         elif not isinstance(x, list):
    653             x = [x]   # constant polynomials
    654 
    655         if check:
    656             x = [QQ(z) for z in x]
    657 
    658         self.__list = list(x)
    659         while len(self.__list) > 0 and not self.__list[-1]:
    660             del self.__list[-1]
    661        
    662         # NOTE: It is much faster to convert to string and let pari's parser at it,
    663         # which is why we pass str(x) in.
    664         self.__poly = pari(str(x)).Polrev()
    665         assert self.__poly.type() == "t_POL"
    666 
    667     def _repr(self, name=None):
    668         if name is None:
    669             name = self.parent().variable_name()
    670         return str(self.__poly).replace("x", name)
    671 
    672     def _repr_(self):
    673         return self._repr()
    674 
    675     def __reduce__(self):
    676         return Polynomial_rational_dense, \
    677                (self.parent(), self.list(), False, self.is_gen())
    678        
    679     def __getitem__(self, n):
    680         return QQ(self.__poly[n])
    681 
    682     def __getslice__(self, i, j):
    683         if i < 0:
    684             i = 0
    685         v = [QQ(x) for x in self.__poly[i:j]]
    686         P = self.parent()
    687         return P([0]*int(i) + v)
    688 
    689     def _pow(self, n):
    690         if self.degree() <= 0:
    691             return self.parent()(self[0]**n)
    692         if n < 0:
    693             return (~self)**(-n)
    694         return Polynomial_rational_dense(self.parent(), self.__poly**n, construct=True)
    695      
    696     def _add_(self, right):
    697         return Polynomial_rational_dense(self.parent(),
    698                                          self.__poly + right.__poly, construct=True)
    699 
    700     def is_irreducible(self):
    701         """
    702         EXAMPLES::
    703 
    704             sage: R.<x> = QQ[]
    705             sage: (x^2 + 2).is_irreducible()
    706             True
    707             sage: (x^2 - 1).is_irreducible()
    708             False
    709 
    710         See trac #5140::
    711 
    712             sage: (2*x).is_irreducible()
    713             True
    714 
    715         TESTS::
    716             sage: R(0).is_irreducible()
    717             Traceback (most recent call last):
    718             ...
    719             ValueError: must be nonzero
    720         """
    721         if self.is_zero():
    722             raise ValueError, "must be nonzero"
    723         S = PolynomialRing(ZZ, self.variable_name())
    724         f = S(self.denominator()*self)
    725         f //= f.content()
    726         return f.is_irreducible()
    727 
    728     def galois_group(self, pari_group=False, algorithm='pari'):
    729         r"""
    730         Return the Galois group of f as a permutation group.
    731 
    732         INPUT:
    733 
    734         - ``self`` -- an irreducible polynomial
    735            
    736         - ``pari_group`` -- bool (default: False); if True instead return the
    737           Galois group as a PARI group.  This has a useful label in it, and may
    738           be slightly faster since it doesn't require looking up a group in
    739           Gap.  To get a permutation group from a PARI group P, type
    740           ``PermutationGroup(P)``.
    741 
    742         - ``algorithm`` -- ``'pari'``, ``'kash'``, ``'magma'`` (default: ``'pari'``,
    743           except when the degree is `\ge 12` when ``'kash'`` is tried). NOTE:
    744           ``'magma'`` also does not return a proven correct result.  Please see
    745           the Magma docs for how to get a proven result.
    746 
    747         ALGORITHM: The Galois group is computed using PARI in C library mode,
    748         or possibly kash or magma.
    749 
    750         .. note::
    751        
    752             The PARI documentation contains the following warning: The method
    753             used is that of resolvent polynomials and is sensitive to the
    754             current precision. The precision is updated internally but, in very
    755             rare cases, a wrong result may be returned if the initial precision
    756             was not sufficient.
    757 
    758         EXAMPLES::
    759 
    760             sage: R.<x> = PolynomialRing(QQ)
    761             sage: f = x^4 - 17*x^3 - 2*x + 1
    762             sage: G = f.galois_group(); G            # optional - database_gap
    763             Transitive group number 5 of degree 4
    764             sage: G.gens()                           # optional - database_gap
    765             [(1,2,3,4), (1,2)]
    766             sage: G.order()                          # optional - database_gap
    767             24
    768 
    769         It is potentially useful to instead obtain the corresponding
    770         PARI group, which is little more than a 4-tuple.  See the
    771         PARI manual for the exact details.  (Note that the third
    772         entry in the tuple is in the new standard ordering.)
    773        
    774         ::
    775 
    776             sage: f = x^4 - 17*x^3 - 2*x + 1
    777             sage: G = f.galois_group(pari_group=True); G
    778             PARI group [24, -1, 5, "S4"] of degree 4
    779             sage: PermutationGroup(G)                # optional - database_gap
    780             Transitive group number 5 of degree 4
    781 
    782         You can use KASH to compute Galois groups as well.  The
    783         advantage is that KASH can compute Galois groups of fields up
    784         to degree 23, whereas PARI only goes to degree 11.  (In my
    785         not-so-thorough experiments PARI is faster than KASH.)
    786 
    787         ::
    788 
    789             sage: f = x^4 - 17*x^3 - 2*x + 1
    790             sage: f.galois_group(algorithm='kash')      # optional - kash
    791             Transitive group number 5 of degree 4
    792        
    793             sage: f = x^4 - 17*x^3 - 2*x + 1
    794             sage: f.galois_group(algorithm='magma')     # optional - magma, database_gap
    795             Transitive group number 5 of degree 4
    796         """
    797         from sage.groups.all import PariGroup, PermutationGroup, TransitiveGroup
    798 
    799         if not self.is_irreducible():
    800             raise ValueError, "polynomial must be irreducible"
    801        
    802         if self.degree() > 11 and algorithm=='pari':
    803             algorithm = 'kash'
    804 
    805         if algorithm == 'kash':
    806             try:
    807                 from sage.interfaces.all import kash
    808                 kash.eval('X := PolynomialRing(RationalField()).1')
    809                 s = self._repr(name='X')
    810                 G = kash('Galois(%s)'%s)
    811                 d = int(kash.eval('%s.ext1'%G.name()))
    812                 n = int(kash.eval('%s.ext2'%G.name()))
    813                 return TransitiveGroup(d, n)
    814             except RuntimeError, msg:
    815                 raise NotImplementedError, str(msg) + "\nSorry, computation of Galois groups of fields of degree bigger than 11 is not yet implemented.  Try installing the optional free (closed source) KASH package, which supports larger degrees, or use algorithm='magma' if you have magma."
    816         elif algorithm == 'magma':
    817             from sage.interfaces.all import magma
    818             X = magma(self).GaloisGroup()
    819             try:
    820                 n, d = X.TransitiveGroupIdentification(nvals=2)
    821                 d = int(d)
    822                 n = int(n)
    823             except RuntimeError, msg:
    824                 raise RuntimeError, str(msg) + "\nUnable to lookup description of Galois group as a transitive group.\n%s"%X
    825             return TransitiveGroup(d, n)
    826            
    827         G = self.__poly.polgalois()
    828         H = PariGroup(G, self.degree())
    829         if pari_group:
    830             return H
    831         else:
    832             return PermutationGroup(H)
    833 
    834     @coerce_binop
    835     def quo_rem(self, right):
    836         """
    837         Returns a tuple (quotient, remainder) where
    838         self = quotient*right + remainder.
    839 
    840         EXAMPLES::
    841 
    842             sage: R.<x> = QQ[]
    843             sage: f = x^5 + 17*x + 3
    844             sage: g = x^3 - 19
    845             sage: q,r = f.quo_rem(g)
    846             sage: q*g + r
    847             x^5 + 17*x + 3       
    848         """
    849         v = self.__poly.divrem(right.__poly)
    850         return Polynomial_rational_dense(self.parent(), v[0], construct=True), \
    851                Polynomial_rational_dense(self.parent(), v[1], construct=True)
    852        
    853                                    
    854     def _mul_(self, right):
    855         """
    856         EXAMPLES::
    857 
    858             sage: R.<x> = QQ[]
    859             sage: (x - QQ('2/3'))*(x^2 - 8*x + 16)
    860             x^3 - 26/3*x^2 + 64/3*x - 32/3
    861         """
    862         return self.parent()(self.__poly * right.__poly, construct=True)       
    863 
    864     def _sub_(self, right):
    865         """
    866         EXAMPLES::
    867 
    868             sage: R.<x> = QQ[]
    869             sage: x^5 + 17*x^3 + x+ 3 - (x^3 - 19)
    870             x^5 + 16*x^3 + x + 22
    871         """
    872         return self.parent()(self.__poly - right.__poly, construct=True)
    873                                
    874     def _unsafe_mutate(self, n, value):
    875         try:
    876             del self.__list
    877         except AttributeError:
    878             pass
    879         n = int(n)
    880         if n < 0:
    881             raise IndexError, "n must be >= 0"
    882         if n <= self.__poly.poldegree():
    883             self.__poly[n] = QQ(value)
    884         else:
    885             self.__poly = self.__poly + pari('(%s)*x^%s'%(QQ(value),n))
    886         if hasattr(self, "__list"):
    887             del self.__list
    888      
    889 
    890     def real_root_intervals(self):
    891         """
    892         Returns isolating intervals for the real roots of this polynomial.
    893 
    894         EXAMPLE::
    895 
    896             sage: R.<x> = PolynomialRing(QQ)
    897             sage: f = (x - 1/2) * (x - 3/4) * (x - 3/2)
    898             sage: f.real_root_intervals()
    899             [((243/512, 1215/2048), 1), ((729/1024, 1701/2048), 1), ((243/256, 1011/512), 1)]
    900         """
    901 
    902         from sage.rings.polynomial.real_roots import real_roots
    903 
    904         return real_roots(self)
    905 
    906     def __copy__(self):
    907         """
    908         Return a copy of this polynomial.
    909         """
    910         f = Polynomial_rational_dense(self.parent())
    911         from copy import copy
    912         f.__poly = copy(self.__poly)
    913         return f
    914        
    915     def degree(self, gen=None):
    916         """
    917         Return the degree of this polynomial. The zero polynomial
    918         has degree -1.
    919 
    920         EXAMPLES::
    921 
    922             sage: R.<x> = QQ[]
    923             sage: (x^5 + 17*x^3 + x+ 3).degree()
    924             5
    925             sage: R(0).degree()
    926             -1
    927             sage: type(x.degree())
    928             <type 'sage.rings.integer.Integer'>
    929         """
    930         return ZZ(max(self.__poly.poldegree(), -1))
    931 
    932     def discriminant(self):
    933         """
    934         EXAMPLES::
    935 
    936             sage: _.<x> = PolynomialRing(QQ)
    937             sage: f = x^3 + 3*x - 17
    938             sage: f.discriminant()
    939             -7911
    940         """
    941         return QQ(self.__poly.poldisc())
    942 
    943     def disc(self):
    944         """
    945         Same as :meth:`.discriminant`.
    946 
    947         EXAMPLES::
    948 
    949             sage: _.<x> = PolynomialRing(QQ)
    950             sage: f = x^3 + 3*x - 17
    951             sage: f.disc()
    952             -7911
    953         """
    954         return self.discriminant()
    955        
    956     def numerator(self):
    957         """
    958         Returns the numerator of self as a polynomial in `\ZZ[x]`.
    959        
    960         EXAMPLES::
    961 
    962             sage: R.<x> = QQ[]
    963             sage: (x/2).numerator()
    964             x
    965             sage: (x + 1/2).numerator()
    966             2*x + 1
    967             sage: (x^2/12 - 1/15).numerator()
    968             5*x^2 - 4
    969             sage: f = R.random_element(60)
    970             sage: f.numerator() in ZZ['x']
    971             True
    972             sage: f.numerator() / f.denominator() == f
    973             True
    974         """
    975         return ZZ[self.variable_name()](self.denominator() * self)
    976        
    977     def denominator(self):
    978         """
    979         Returns the denominator of self as an element of `\ZZ`.
    980        
    981         EXAMPLES::
    982 
    983             sage: R.<x> = QQ[]
    984             sage: (x/2).denominator()
    985             2
    986             sage: (x/2 + 1/3).denominator()
    987             6
    988             sage: R.<x> = QQ[]
    989             sage: f = R.random_element(50)
    990             sage: f * f.denominator() in ZZ['x']
    991             True
    992         """
    993         return integer.LCM_list([a.denominator() for a in self])
    994      
    995     def factor_mod(self, p):
    996         """
    997         Return the factorization of self modulo the prime p.
    998 
    999         INPUT:
    1000 
    1001         - ``p`` -- prime
    1002 
    1003         OUTPUT:
    1004        
    1005         factorization of self reduced modulo p.
    1006 
    1007         EXAMPLES::
    1008 
    1009             sage: R.<x> = QQ[]
    1010             sage: (x^5 + 17*x^3 + x+ 3).factor_mod(3)
    1011             x * (x^2 + 1)^2
    1012             sage: (x^5 + 2).factor_mod(5)
    1013             (x + 2)^5       
    1014         """
    1015         import sage.rings.finite_rings.constructor as finite_field
    1016         p = integer.Integer(p)
    1017         if not p.is_prime():
    1018             raise ValueError, "p must be prime"
    1019         G = self._pari_().factormod(p)
    1020         K = finite_field.FiniteField(p)
    1021         R = K[self.parent().variable_name()]
    1022         return R(1)._factor_pari_helper(G, unit=R(self).leading_coefficient())
    1023 
    1024     def factor_padic(self, p, prec=10):
    1025         """
    1026         Return p-adic factorization of self to given precision.
    1027 
    1028         INPUT:
    1029 
    1030         - p -- prime
    1031         - prec -- integer; the precision
    1032 
    1033         OUTPUT:
    1034 
    1035         factorization of self viewed as a polynomial over the p-adics
    1036 
    1037         EXAMPLES::
    1038 
    1039             sage: R.<x> = QQ[]
    1040             sage: f = x^3 - 2
    1041             sage: f.factor_padic(2)
    1042             (1 + O(2^10))*x^3 + (O(2^10))*x^2 + (O(2^10))*x + (2 + 2^2 + 2^3 + 2^4 + 2^5 + 2^6 + 2^7 + 2^8 + 2^9 + O(2^10))
    1043             sage: f.factor_padic(3)
    1044             (1 + O(3^10))*x^3 + (O(3^10))*x^2 + (O(3^10))*x + (1 + 2*3 + 2*3^2 + 2*3^3 + 2*3^4 + 2*3^5 + 2*3^6 + 2*3^7 + 2*3^8 + 2*3^9 + O(3^10))
    1045             sage: f.factor_padic(5)
    1046             ((1 + O(5^10))*x + (2 + 4*5 + 2*5^2 + 2*5^3 + 5^4 + 3*5^5 + 4*5^7 + 2*5^8 + 5^9 + O(5^10))) * ((1 + O(5^10))*x^2 + (3 + 2*5^2 + 2*5^3 + 3*5^4 + 5^5 + 4*5^6 + 2*5^8 + 3*5^9 + O(5^10))*x + (4 + 5 + 2*5^2 + 4*5^3 + 4*5^4 + 3*5^5 + 3*5^6 + 4*5^7 + 4*5^9 + O(5^10)))
    1047         """
    1048         from sage.rings.padics.factory import Qp
    1049         p = integer.Integer(p)
    1050         if not p.is_prime():
    1051             raise ValueError, "p must be prime"
    1052         prec = integer.Integer(prec)
    1053         if prec <= 0:
    1054             raise ValueError, "prec must be positive"
    1055         K = Qp(p, prec, type='capped-rel')
    1056         R = K[self.parent().variable_name()]
    1057         return R(self).factor() # absprec = prec)
    1058 
    1059     def list(self):
    1060         """
    1061         Return a new copy of the list of the underlying
    1062         elements of self.
    1063 
    1064         EXAMPLES::
    1065 
    1066             sage: _.<x> = PolynomialRing(QQ)
    1067             sage: f = x^3 + 3*x - 17/13; f
    1068             x^3 + 3*x - 17/13
    1069             sage: v = f.list(); v
    1070             [-17/13, 3, 0, 1]
    1071             sage: v[0] = 0
    1072             sage: f
    1073             x^3 + 3*x - 17/13
    1074             sage: f.list()
    1075             [-17/13, 3, 0, 1]           
    1076         """
    1077         return [QQ(x) for x in self.__poly.Vecrev()]
    1078 
    1079    
    1080 
    1081 ##     def partial_fraction(self, g):
    1082 ##         """
    1083 ##         Return partial fraction decomposition of self/g, where g
    1084 ##         has the same parent as self.
    1085 ##         """
    1086 ##         g = self.parent()(g)
    1087 ##         from sage.interfaces.maxima import maxima
    1088 ##         h = maxima(self)/maxima(g)
    1089 ##         k = h.partfrac(self.parent().variable())
    1090 
    1091     def rescale(self, a):
    1092         """
    1093         Return `f(a*X)`.
    1094         """
    1095         b = 1
    1096         c = []
    1097         for i in range(self.degree()+1):
    1098             c.append(b*self[i])
    1099             b *= a
    1100         return self.parent()(c)
    1101 
    1102     def hensel_lift(self, p, e):
    1103         """
    1104         Assuming that self factors modulo `p` into distinct factors,
    1105         computes the Hensel lifts of these factors modulo `p^e`.  We
    1106         assume that `p` has integer coefficients.
    1107         """
    1108         p = integer.Integer(p)
    1109         if not p.is_prime():
    1110             raise ValueError, "p must be prime"
    1111         e = integer.Integer(e)
    1112         if e < 1:
    1113             raise ValueError, "e must be at least 1"
    1114         F = self.factor_mod(p)
    1115         y = []
    1116         for g, n in F:
    1117             if n > 1:
    1118                 raise ArithmeticError, "The polynomial must be square free modulo p."
    1119             y.append(g)
    1120         H = self._pari_().polhensellift(y, p, e)
    1121         from sage.rings.finite_rings.integer_mod_ring import IntegerModRing
    1122         R = IntegerModRing(p**e)
    1123         S = R[self.parent().variable_name()]
    1124         return [S(eval(str(m.Vec().Polrev().Vec()))) for m in H]
     614register_unpickle_override( \
     615    'sage.rings.polynomial.polynomial_element_generic', \
     616    'Polynomial_rational_dense', Polynomial_rational_flint)
    1125617
    1126618class Polynomial_padic_generic_dense(Polynomial_generic_dense, Polynomial_generic_domain):
    1127619    def __init__(self, parent, x=None, check=True, is_gen = False, construct=False, absprec=None):
  • new file sage/rings/polynomial/polynomial_rational_flint.pxd

    diff --git a/sage/rings/polynomial/polynomial_rational_flint.pxd b/sage/rings/polynomial/polynomial_rational_flint.pxd
    new file mode 100644
    - +  
     1###############################################################################
     2#          Copyright (C) 2010 Sebastian Pancratz <sfp@pancratz.org>           #
     3#                                                                             #
     4#     Distributed under the terms of the GNU General Public License (GPL)     #
     5#                                                                             #
     6#                        http://www.gnu.org/licenses/                         #
     7###############################################################################
     8
     9include "../../ext/cdefs.pxi"
     10include "../../libs/flint/fmpz.pxi"
     11include "../../libs/flint/fmpz_poly.pxi"
     12include "../../libs/flint/fmpq_poly.pxd"
     13
     14from sage.rings.polynomial.polynomial_element cimport Polynomial
     15
     16cdef class Polynomial_rational_flint(Polynomial):
     17    cdef fmpq_poly_t __poly
     18
     19    cdef Polynomial_rational_flint _new(self)
     20    cpdef _unsafe_mutate(self, unsigned long n, value)
     21    cpdef Polynomial truncate(self, long n)
     22
  • new file sage/rings/polynomial/polynomial_rational_flint.pyx

    diff --git a/sage/rings/polynomial/polynomial_rational_flint.pyx b/sage/rings/polynomial/polynomial_rational_flint.pyx
    new file mode 100644
    - +  
     1r"""
     2Univariate polynomials over `\QQ` implemented via FLINT
     3
     4AUTHOR:
     5
     6- Sebastian Pancratz
     7"""
     8
     9###############################################################################
     10#          Copyright (C) 2010 Sebastian Pancratz <sfp@pancratz.org>           #
     11#                                                                             #
     12#     Distributed under the terms of the GNU General Public License (GPL)     #
     13#                                                                             #
     14#                        http://www.gnu.org/licenses/                         #
     15###############################################################################
     16
     17include "../../ext/stdsage.pxi"
     18include "../../ext/interrupt.pxi"
     19include "../../ext/gmp.pxi"
     20include "../../libs/ntl/decl.pxi"
     21
     22include "../../ext/cdefs.pxi"
     23include "../../libs/flint/fmpz.pxi"
     24include "../../libs/flint/fmpz_poly.pxi"
     25include "../../libs/flint/fmpq_poly.pxd"
     26
     27from sage.interfaces.all import singular as singular_default
     28
     29from sage.libs.all import pari, pari_gen
     30from sage.libs.flint.ntl_interface cimport *
     31from sage.libs.flint.fmpz_poly cimport fmpz_poly_set
     32
     33from sage.rings.integer cimport Integer
     34from sage.rings.integer_ring import ZZ
     35from sage.rings.fraction_field_element import FractionFieldElement
     36from sage.rings.rational cimport Rational
     37from sage.rings.rational_field import QQ
     38from sage.rings.polynomial.polynomial_element cimport Polynomial
     39from sage.rings.polynomial.polynomial_element import is_Polynomial
     40from sage.rings.polynomial.polynomial_integer_dense_flint cimport Polynomial_integer_dense_flint
     41
     42from sage.structure.element cimport Element, ModuleElement, RingElement
     43from sage.structure.element import coerce_binop
     44from sage.structure.factorization import Factorization
     45
     46
     47cdef inline bint _do_sig(fmpq_poly_t op):
     48    """
     49    Returns 1 when signal handling should be carried out for an operation
     50    on this polynomial and 0 otherwise.
     51
     52    Strictly speaking, whether or not signal handling should be carried
     53    ought to depend on the operation as well as the operands in question.
     54    For simplicity we carry out signal handling for all but the simplest
     55    of operands regardless of the operation.
     56
     57    TESTS::
     58
     59        sage: R.<t> = QQ[]
     60        sage: f = 1 + t/2
     61        sage: g = 2/3 + t^2
     62        sage: _ = f * g      # indirect doctest
     63    """
     64    return fmpz_poly_max_limbs(fmpq_poly_numref(op)) > 1 \
     65                                                 or fmpq_poly_degree(op) > 1000
     66
     67cdef class Polynomial_rational_flint(Polynomial):
     68    """
     69    Univariate polynomials over the rationals, implemented via FLINT.
     70   
     71    Internally, we represent rational polynomial as the quotient of an integer
     72    polynomial and a positive denominator which is coprime to the content of
     73    the numerator.
     74   
     75    The implementation is based on the C module fmpq_poly written on top of
     76    FLINT 1.
     77    """
     78
     79    ###########################################################################
     80    # Allocation & initialisation                                             #
     81    ###########################################################################
     82   
     83    cdef Polynomial_rational_flint _new(self):
     84        """
     85        Quickly creates a new polynomial object in this class.
     86       
     87        OUTPUT:
     88       
     89        - Polynomial of type Polynomial_rational_flint
     90       
     91        TESTS::
     92       
     93            sage: R.<t> = QQ[]
     94            sage: f = 2/3*t^2
     95            sage: g = -1/2*t + 2
     96            sage: f + g           # indirect doctest
     97            2/3*t^2 - 1/2*t + 2
     98        """
     99        cdef Polynomial_rational_flint res = PY_NEW(Polynomial_rational_flint)
     100        res._parent = self._parent
     101        res._is_gen = 0
     102        return res
     103
     104    def __cinit__(self):
     105        """
     106        Initialises the underlying data structure.
     107       
     108        TESTS::
     109       
     110            sage: R.<t> = QQ[]
     111            sage: f = 2/3 * t - 7  #indirect doctest
     112        """
     113        fmpq_poly_init(self.__poly)
     114
     115    def __dealloc__(self):
     116        """
     117        Deallocates the underlying data structure.
     118       
     119        TESTS::
     120       
     121            sage: R.<t> = QQ[]
     122            sage: f = 1/3 * t
     123            sage: del f        # untested
     124        """
     125        fmpq_poly_clear(self.__poly)
     126
     127    def __init__(self, parent, x=None, check=True, is_gen=False, construct=False):
     128        """
     129        Initialises the associated data for the polynomial self.
     130       
     131        INPUT:
     132       
     133        - ``parent`` - Polynomial ring, the parent of self
     134        - ``x`` - Data for the new polynomial self, e.g. a polynomial, an
     135          integer, a rational, a list of rationals, a dictionary with keys
     136          the degrees and the rational coefficients, etc (default: ``None``)
     137        - `check`` - Whether the integrity of the data needs to be verified,
     138          largely ignored by this method (default: ``True``)
     139        - ``is_gen`` - Whether self shall be initialised as the generator of
     140          the parent polynomial ring
     141        - ``construct`` - Whether the element shall always be constructed
     142          as an independent copy of any input data (default: ``False``)
     143       
     144        TESTS::
     145       
     146            sage: R.<t> = QQ[]
     147            sage: f = -4 * t^2 + 1/3 * t - 1/7  # indirect doctest
     148        """
     149        cdef long deg
     150        cdef unsigned long n
     151        cdef Rational c
     152        cdef list L1
     153        cdef mpq_t * L2
     154       
     155        Polynomial.__init__(self, parent, is_gen=is_gen)
     156       
     157        if is_gen:
     158            fmpq_poly_set_coeff_si(self.__poly, 1, 1)
     159       
     160        elif PY_TYPE_CHECK(x, Polynomial_rational_flint):
     161            fmpq_poly_set(self.__poly, (<Polynomial_rational_flint> x).__poly)
     162
     163        elif PY_TYPE_CHECK(x, int):
     164            fmpq_poly_set_si(self.__poly, <int> x)
     165       
     166        elif PY_TYPE_CHECK(x, Integer):
     167            fmpq_poly_set_mpz(self.__poly, (<Integer> x).value)
     168
     169        elif PY_TYPE_CHECK(x, Rational):
     170            fmpq_poly_set_mpq(self.__poly, (<Rational> x).value)
     171
     172        elif PY_TYPE_CHECK(x, list) or PY_TYPE_CHECK(x, tuple):
     173           
     174            if len(x) == 0:
     175                return
     176            elif len(x) == 1:
     177                Polynomial_rational_flint.__init__(self, parent, x[0], \
     178                                check=check, is_gen=False, construct=construct)
     179                return
     180           
     181            L1 = [e if isinstance(e, Rational) else Rational(e) for e in x]
     182            n  = <unsigned long> len(x)
     183            _sig_on
     184            L2 = <mpq_t *> malloc(n * sizeof(mpq_t))
     185            for deg from 0 <= deg < n:
     186                mpq_init(L2[deg])
     187                mpq_set(L2[deg], (<Rational> L1[deg]).value)
     188            _fmpq_poly_from_list(self.__poly, L2, n)
     189            for deg from 0 <= deg < n:
     190                mpq_clear(L2[deg])
     191            free(L2)
     192            _sig_off
     193           
     194#           deg = 0
     195#           for e in x:
     196#               c = Rational(e)
     197#               fmpq_poly_set_coeff_mpq(self.__poly, deg, c.value)
     198#               deg += 1
     199           
     200        elif PY_TYPE_CHECK(x, dict):
     201            for deg, e in x.iteritems():
     202                c = Rational(e)
     203                fmpq_poly_set_coeff_mpq(self.__poly, deg, c.value)
     204
     205        elif PY_TYPE_CHECK(x, pari_gen):
     206            k = self._parent.base_ring()
     207            x = [k(w) for w in x.Vecrev()]
     208            Polynomial_rational_flint.__init__(self, parent, x, check=True, \
     209                                             is_gen=False, construct=construct)
     210
     211        elif PY_TYPE_CHECK(x, Polynomial):
     212            k = self._parent.base_ring()
     213            x = [k(w) for w in list(x)]
     214            Polynomial_rational_flint.__init__(self, parent, x, check=True, \
     215                                             is_gen=False, construct=construct)
     216
     217        elif PY_TYPE_CHECK(x, FractionFieldElement) and (x.parent().base() is parent or x.parent().base() == parent) and x.denominator() == 1:
     218            x = x.numerator()
     219            Polynomial_rational_flint.__init__(self, parent, x, check=check, \
     220                                            is_gen=is_gen, construct=construct)
     221
     222        else:
     223            x = parent.base_ring()(x)
     224            Polynomial_rational_flint.__init__(self, parent, x, check=check, \
     225                                            is_gen=is_gen, construct=construct)
     226   
     227    def __reduce__(self):
     228        """
     229        This is used when pickling polynomials.
     230       
     231        TESTS::
     232       
     233            sage: R.<t> = QQ[]
     234            sage: f = 2/3 * t^2 + 1
     235            sage: r = f.__reduce__(); r
     236            (<type 'sage.rings.polynomial.polynomial_rational_flint.Polynomial_rational_flint'>, (Univariate Polynomial Ring in t over Rational Field, [1, 0, 2/3], False, False))
     237            sage: r[0](*r[1])
     238            2/3*t^2 + 1
     239            sage: loads(dumps(f)) == f
     240            True
     241        """
     242        return Polynomial_rational_flint, \
     243               (self.parent(), self.list(), False, self.is_gen())
     244
     245    def __copy__(self):
     246        """
     247        Returns a copy of self.
     248
     249        TESTS::
     250       
     251            sage: R.<t> = QQ[]
     252            sage: f = 4/5 * t^3 - 1/17
     253            sage: copy(f) == f
     254            True
     255        """
     256        cdef Polynomial_rational_flint res = self._new()
     257        fmpq_poly_set(res.__poly, self.__poly)
     258        return res
     259   
     260    def _singular_(self, singular=singular_default, have_ring=False):
     261        """
     262        Returns a Singular representation of self.
     263
     264        INPUT:
     265
     266        - ``singular`` - Singular interpreter (default: default interpreter)
     267        - ``have_ring`` - set to True if the ring was already set in Singular
     268
     269        EXAMPLES::
     270
     271            sage: P.<x> = PolynomialRing(QQ)
     272            sage: f = 3*x^2 + 2*x + 5
     273            sage: singular(f)
     274            3*x^2+2*x+5
     275        """
     276        if not have_ring:
     277            self._parent._singular_(singular).set_ring()  # Expensive!
     278        return singular(self._singular_init_())
     279
     280    def list(self):
     281        """
     282        Returns a list with the coefficients of self.
     283
     284        EXAMPLES::
     285
     286            sage: R.<t> = QQ[]
     287            sage: f = 1 + t + t^2/2 + t^3/3 + t^4/4
     288            sage: f.list()
     289            [1, 1, 1/2, 1/3, 1/4]
     290            sage: g = R(0)
     291            sage: g.list()
     292            []
     293        """
     294        cdef unsigned long length = fmpq_poly_length(self.__poly)
     295        return [self[n] for n in range(length)]
     296
     297    ###########################################################################
     298    # Basis access                                                            #
     299    ###########################################################################
     300
     301    def degree(self):
     302        """
     303        Returns the degree of self.
     304       
     305        By convention, the degree of the zero polynomial is -1.
     306       
     307        EXAMPLES::
     308       
     309            sage: R.<t> = QQ[]
     310            sage: f = 1 + t + t^2/2 + t^3/3 + t^4/4
     311            sage: f.degree()
     312            4
     313            sage: g = R(0)
     314            sage: g.degree()
     315            -1
     316        """
     317        cdef Integer deg = PY_NEW(Integer)
     318        mpz_set_si(deg.value, fmpq_poly_degree(self.__poly))
     319        return deg
     320
     321    def __getitem__(self, long n):
     322        """
     323        Returns the `n`th coefficient of self.
     324       
     325        INPUT:
     326       
     327        - ``n`` - Degree of the monomial whose coefficient is to be returned
     328       
     329        EXAMPLES::
     330           
     331            sage: R.<t> = QQ[]
     332            sage: f = 1 + t + t^2/2 + t^3/3 + t^4/4
     333            sage: f[-1], f[0], f[3], f[5]            # indirect doctest
     334            (0, 1, 1/3, 0)
     335        """
     336        cdef Rational z = PY_NEW(Rational)
     337        if 0 <= n and n < fmpq_poly_length(self.__poly):
     338            fmpq_poly_get_coeff_mpq(z.value, self.__poly, n)
     339        return z
     340
     341    def __getslice__(self, long i, long j):
     342        """
     343        Returns the subpolynomial of self from the `i`th to the `j`th
     344        coefficient, where the lower bound is inclusive and the upper bound
     345        is exclusive.
     346       
     347        INPUT:
     348       
     349        - ``i`` - Lower index for the slice
     350        - ``j`` - Upper index for the slice
     351       
     352        EXAMPLES::
     353       
     354            sage: R.<t> = QQ[]
     355            sage: f = 1 + t + t^2/2 + t^3/3 + t^4/4
     356            sage: f[1:3]                             # indirect doctest
     357            1/2*t^2 + t
     358        """
     359        cdef Polynomial_rational_flint res = self._new()
     360        cdef bint do_sig = _do_sig(self.__poly)
     361       
     362        if do_sig: _sig_on
     363        fmpq_poly_getslice(res.__poly, self.__poly, i, j)
     364        if do_sig: _sig_off
     365        return res
     366
     367    cpdef _unsafe_mutate(self, unsigned long n, value):
     368        """
     369        Sets the `n`th coefficient of self to value.
     370       
     371        TESTS::
     372       
     373            sage: R.<t> = QQ[]
     374            sage: f = 1 + t + t^2/2 + t^3/3 + t^4/4
     375            sage: f._unsafe_mutate(4, 1/5)
     376            sage: f
     377            1/5*t^4 + 1/3*t^3 + 1/2*t^2 + t + 1
     378       
     379        WARNING:
     380       
     381        Polynomials in Sage are meant to be immutable, and some methods may
     382        rely on this convention.  This method should be used only with the
     383        utmost care.
     384        """
     385        cdef bint do_sig = _do_sig(self.__poly)
     386       
     387        if PY_TYPE_CHECK(value, int):
     388            if do_sig: _sig_on
     389            fmpq_poly_set_coeff_si(self.__poly, n, value)
     390            if do_sig: _sig_off
     391        elif PY_TYPE_CHECK(value, Integer):
     392            if do_sig: _sig_on
     393            fmpq_poly_set_coeff_mpz(self.__poly, n, (<Integer> value).value)
     394            if do_sig: _sig_off
     395        elif PY_TYPE_CHECK(value, Rational):
     396            if do_sig: _sig_on
     397            fmpq_poly_set_coeff_mpq(self.__poly, n, (<Rational> value).value)
     398            if do_sig: _sig_off
     399        else:
     400            value = Rational(value)
     401            if do_sig: _sig_on
     402            fmpq_poly_set_coeff_mpq(self.__poly, n, (<Rational> value).value)
     403            if do_sig: _sig_off
     404
     405    def __call__(self, *x, **kwds):
     406        """
     407        Calls this polynomial with the given parameters, which can be
     408        interpreted as polynomial composition or evaluation by this
     409        method.
     410
     411        If the argument is not simply an integer, a rational, or a
     412        polynomial, the call is passed on to the generic implementation
     413        in the Polynomial class.
     414       
     415        EXAMPLES:
     416
     417        The first example illustrates polynomial composition::
     418       
     419            sage: R.<t> = QQ[]
     420            sage: f = t^2 - 1
     421            sage: g = t + 1
     422            sage: f(g)          # indirect doctest
     423            t^2 + 2*t
     424               
     425        Now we illustrate how a polynomial can be evaluated at a rational
     426        number::
     427       
     428            sage: f(-2/3)       # indirect doctest
     429            -5/9
     430        """
     431        cdef Polynomial_rational_flint f
     432        cdef Rational r
     433       
     434        if len(x) == 1:
     435            a = x[0]
     436            if isinstance(a, Polynomial_rational_flint):
     437                f = (<Polynomial_rational_flint> a)._new()
     438                _sig_on
     439                fmpq_poly_compose(f.__poly, self.__poly, \
     440                    (<Polynomial_rational_flint> a).__poly)
     441                _sig_off
     442                return f
     443            if isinstance(a, Rational):
     444                r = PY_NEW(Rational)
     445                _sig_on
     446                fmpq_poly_evaluate_mpq(r.value, self.__poly, (<Rational> a).value)
     447                _sig_off
     448                return r
     449            if isinstance(a, Integer):
     450                r = PY_NEW(Rational)
     451                _sig_on
     452                fmpq_poly_evaluate_mpz(r.value, self.__poly, (<Integer> a).value)
     453                _sig_off
     454                return r
     455
     456        return Polynomial.__call__(self, *x, **kwds)
     457
     458    cpdef Polynomial truncate(self, long n):
     459        """
     460        Returns self truncated modulo `t^n`.
     461       
     462        INPUT:
     463       
     464        - ``n`` - The power of `t` modulo which self is truncated
     465       
     466        EXAMPLES::
     467       
     468            sage: R.<t> = QQ[]
     469            sage: f = 1 - t + 1/2*t^2 - 1/3*t^3
     470            sage: f.truncate(0)
     471            0
     472            sage: f.truncate(2)
     473            -t + 1
     474        """
     475        cdef Polynomial_rational_flint res
     476        cdef bint do_sig
     477       
     478        if (n >= fmpq_poly_length(self.__poly)):
     479            return self
     480        else:
     481            res = self._new()
     482            if n > 0:
     483                do_sig = _do_sig(self.__poly)
     484                if do_sig: _sig_on
     485                fmpq_poly_truncate(res.__poly, self.__poly, n)
     486                if do_sig: _sig_off
     487            return res
     488
     489    def reverse(self, n = None):
     490        """
     491        Reverses the coefficients of self - thought of as a polynomial of
     492        length `n`.
     493       
     494        Assumes that whenever `n` is not ``None`` it is an integral value
     495        that fits into an unsigned long.  Otherwise, a ValueError is raised.
     496       
     497        INPUT:
     498       
     499        - ``n`` - Integral value that fits in an unsigned long:  the power
     500          of `t` modulo which we consider self before reversing it
     501          (default:  ``None``, interpreted as the length of self)
     502       
     503        OUTPUT:
     504       
     505        - The reversed polynomial as a Polynomial_rational_flint
     506       
     507        EXAMPLES:
     508       
     509        We first consider the simplest case, where we reverse all coefficients
     510        of a polynomial and obtain a polynomial of the same degree::
     511       
     512            sage: R.<t> = QQ[]
     513            sage: f = 1 + t + t^2 / 2 + t^3 / 3 + t^4 / 4
     514            sage: f.reverse()
     515            t^4 + t^3 + 1/2*t^2 + 1/3*t + 1/4
     516       
     517        Next, an example we the returned polynomial has lower degree because
     518        the original polynomial has low coefficients equal to zero::
     519       
     520            sage: R.<t> = QQ[]
     521            sage: f = 3/4*t^2 + 6*t^7
     522            sage: f.reverse()
     523            3/4*t^5 + 6
     524       
     525        The next example illustrates the passing of a value for `n` less than
     526        the length of self, notationally resulting in truncation prior to
     527        reversing::
     528       
     529            sage: R.<t> = QQ[]
     530            sage: f = 1 + t + t^2 / 2 + t^3 / 3 + t^4 / 4
     531            sage: f.reverse(3)
     532            t^2 + t + 1/2
     533       
     534        Now we illustrate the passing of a value for `n` greater than the
     535        length of self, notationally resulting in zero padding at the top
     536        end prior to reversing::
     537       
     538            sage: R.<t> = QQ[]
     539            sage: f = 1 + t + t^2 / 2 + t^3 / 3
     540            sage: f.reverse(5)
     541            t^4 + t^3 + 1/2*t^2 + 1/3*t
     542       
     543        TESTS::
     544       
     545        We illustrate two ways in which the interpretation of `n` as an
     546        unsigned long int may fail.  Firstly, an integral value which is
     547        too large, yielding an OverflowError::
     548       
     549            sage: R.<t> = QQ[]
     550            sage: f = 1 + t/2
     551            sage: f.reverse(2**64)
     552            Traceback (most recent call last):
     553            ...
     554            OverflowError: long int too large to convert
     555       
     556        Secondly, a value which cannot be converted to an integral value,
     557        resulting in a TypeError::
     558           
     559            sage: R.<t> = QQ[]
     560            sage: f = 1 + t/2
     561            sage: f.reverse(I)
     562            Traceback (most recent call last):
     563            ...
     564            TypeError: can't convert complex to int; use int(abs(z))
     565        """
     566        cdef unsigned long len
     567        cdef Polynomial_rational_flint res
     568        cdef bint do_sig
     569       
     570        if n is None:
     571            len = fmpq_poly_length(self.__poly)
     572        else:
     573            len = <unsigned long> n
     574       
     575        res = self._new()
     576        do_sig = _do_sig(self.__poly)
     577        if do_sig: _sig_on
     578        fmpq_poly_reverse(res.__poly, self.__poly, len)
     579        if do_sig: _sig_off
     580        return res
     581   
     582    ###########################################################################
     583    # Comparisons                                                             #
     584    ###########################################################################
     585
     586    def is_zero(self):
     587        """
     588        Returns whether or not self is the zero polynomial.
     589       
     590        EXAMPLES::
     591       
     592            sage: R.<t> = QQ[]
     593            sage: f = 1 - t + 1/2*t^2 - 1/3*t^3
     594            sage: f.is_zero()
     595            False
     596            sage: R(0).is_zero()
     597            True
     598        """
     599        return bool(fmpq_poly_is_zero(self.__poly))
     600
     601    def __nonzero__(self):
     602        """
     603        Returns whether or not self is non-zero.
     604       
     605        EXAMPLES::
     606       
     607            sage: R.<t> = QQ[]
     608            sage: f = 1 - t + 1/2*t^2 - 1/3*t^3
     609            sage: bool(f)
     610            True
     611            sage: bool(R(0))
     612            False
     613        """
     614        return not fmpq_poly_is_zero(self.__poly)
     615
     616    ###########################################################################
     617    # Shifting                                                                #
     618    ###########################################################################
     619
     620    def __lshift__(self, n):
     621        """
     622        Notationally multiplies self by `t^n`.
     623       
     624        EXAMPLES::
     625       
     626            sage: R.<t> = QQ[]
     627            sage: t << 10                     # indirect doctest
     628            t^11
     629       
     630        TESTS::
     631       
     632            sage: R.<t> = QQ[]
     633            sage: f = R.random_element(1000)
     634            sage: (f << 23) >> 23 == f        # indirect doctest
     635            True
     636        """
     637        cdef unsigned long k = <unsigned long> n
     638        cdef Polynomial_rational_flint f = <Polynomial_rational_flint> self
     639        cdef Polynomial_rational_flint res
     640        cdef bint do_sig
     641       
     642        if k == 0 or fmpq_poly_is_zero(f.__poly):
     643            return self
     644        else:
     645            res = f._new()
     646            do_sig = fmpq_poly_length(f.__poly) > 5000 or n > 5000
     647           
     648            if do_sig: _sig_on
     649            fmpq_poly_left_shift(res.__poly, f.__poly, k)
     650            if do_sig: _sig_off
     651            return res
     652
     653    def __rshift__(self, n):
     654        """
     655        Notationally returns the quotient of Euclidean division of self
     656        by `t^n`.
     657       
     658        EXAMPLES::
     659       
     660            sage: R.<t> = QQ[]
     661            sage: f = 1 + t + t^2/2 + t^3/3 + t^4/4
     662            sage: f >> 2
     663            1/4*t^2 + 1/3*t + 1/2
     664        """
     665        cdef unsigned long k = <unsigned long> n
     666        cdef Polynomial_rational_flint f = <Polynomial_rational_flint> self
     667        cdef Polynomial_rational_flint res
     668        cdef bint do_sig
     669       
     670        if k == 0 or fmpq_poly_is_zero(f.__poly):
     671            return self
     672        else:
     673            res = f._new()
     674            do_sig = _do_sig(f.__poly)
     675           
     676            if do_sig: _sig_on
     677            fmpq_poly_right_shift(res.__poly, f.__poly, k)
     678            if do_sig: _sig_off
     679            return res
     680
     681    ###########################################################################
     682    # Arithmetic                                                              #
     683    ###########################################################################
     684
     685    cpdef ModuleElement _add_(self, ModuleElement right):
     686        """
     687        Returns the sum of two rational polynomials.
     688       
     689        EXAMPLES::
     690       
     691            sage: R.<t> = QQ[]
     692            sage: f = 2/3 + t + 2*t^3
     693            sage: g = -1 + t/3 - 10/11*t^4
     694            sage: f + g
     695            -10/11*t^4 + 2*t^3 + 4/3*t - 1/3
     696       
     697        TESTS::
     698       
     699            sage: R.<t> = QQ[]
     700            sage: f = R.random_element(2000)
     701            sage: f + f == 2 * f              # indirect doctest
     702            True
     703        """
     704        cdef Polynomial_rational_flint op2 = <Polynomial_rational_flint> right
     705        cdef Polynomial_rational_flint res = self._new()
     706        cdef bint do_sig = _do_sig(self.__poly) or _do_sig(op2.__poly)
     707       
     708        if do_sig: _sig_on
     709        fmpq_poly_add(res.__poly, self.__poly, op2.__poly)
     710        if do_sig: _sig_off
     711        return res
     712
     713    cpdef ModuleElement _sub_(self, ModuleElement right):
     714        """
     715        Returns the difference of two rational polynomials.
     716       
     717        EXAMPLES::
     718
     719            sage: R.<t> = QQ[]
     720            sage: f = -10/11*t^4 + 2*t^3 + 4/3*t - 1/3
     721            sage: g = 2*t^3
     722            sage: f - g                                 # indirect doctest
     723            -10/11*t^4 + 4/3*t - 1/3
     724       
     725        TESTS::
     726       
     727            sage: R.<t> = QQ[]
     728            sage: f = R.random_element(2000)
     729            sage: f - f/2 == 1/2 * f          # indirect doctest
     730            True
     731            sage: f[:1000] == f - f[1000:]    # indirect doctest
     732            True
     733        """
     734        cdef Polynomial_rational_flint op2 = <Polynomial_rational_flint> right
     735        cdef Polynomial_rational_flint res = self._new()
     736        cdef bint do_sig = _do_sig(self.__poly) or _do_sig(op2.__poly)
     737       
     738        if do_sig: _sig_on
     739        fmpq_poly_sub(res.__poly, self.__poly, op2.__poly)
     740        if do_sig: _sig_off
     741        return res
     742
     743    cpdef ModuleElement _neg_(self):
     744        """
     745        Returns the difference of two rational polynomials.
     746       
     747        EXAMPLES::
     748       
     749            sage: R.<t> = QQ[]
     750            sage: f = 3*t/2
     751            sage: -f            # indirect doctest
     752            -3/2*t
     753       
     754        TESTS::
     755       
     756            sage: R.<t> = QQ[]
     757            sage: f = R.random_element(2000)
     758            sage: f + (-f) == 0               # indirect doctest
     759            True
     760        """
     761        cdef Polynomial_rational_flint res = self._new()
     762        cdef bint do_sig = _do_sig(self.__poly)
     763       
     764        if do_sig: _sig_on
     765        fmpq_poly_neg(res.__poly, self.__poly)
     766        if do_sig: _sig_off
     767        return res
     768
     769    @coerce_binop
     770    def quo_rem(self, right):
     771        """
     772        Returns the quotient and remainder of the Euclidean division of
     773        self and right.
     774       
     775        Raises a ZerodivisionError if right is zero.
     776       
     777        EXAMPLES::
     778       
     779            sage: R.<t> = QQ[]
     780            sage: f = R.random_element(2000)
     781            sage: g = R.random_element(1000)
     782            sage: q, r = f.quo_rem(g)
     783            sage: f == q*g + r
     784            True
     785        """
     786        if right.is_zero():
     787            raise ZeroDivisionError, "division by zero polynomial"
     788        if self.is_zero():
     789            return self, self
     790
     791        cdef Polynomial_rational_flint qq = self._new()
     792        cdef Polynomial_rational_flint rr = self._new()
     793
     794        _sig_on
     795        fmpq_poly_divrem(qq.__poly, rr.__poly, self.__poly, \
     796                         (<Polynomial_rational_flint> right).__poly)
     797        _sig_off
     798        return qq, rr
     799
     800    @coerce_binop
     801    def gcd(self, right):
     802        """
     803        Returns the (monic) greatest common divisor of self and right.
     804       
     805        Corner cases:  if self and right are both zero, returns zero.  If
     806        only one of them is zero, returns the other polynomial, up to
     807        normalisation.
     808       
     809        EXAMPLES::
     810       
     811            sage: R.<t> = QQ[]
     812            sage: f = -2 + 3*t/2 + 4*t^2/7 - t^3
     813            sage: g = 1/2 + 4*t + 2*t^4/3
     814            sage: f.gcd(g)
     815            1
     816            sage: f = (-3*t + 1/2) * f
     817            sage: g = (-3*t + 1/2) * (4*t^2/3 - 1) * g
     818            sage: f.gcd(g)
     819            t - 1/6
     820        """
     821        cdef Polynomial_rational_flint res = self._new()
     822       
     823        _sig_on
     824        fmpq_poly_gcd(res.__poly, self.__poly, \
     825                (<Polynomial_rational_flint> right).__poly)
     826        _sig_off
     827        return res
     828
     829    @coerce_binop
     830    def lcm(self, right):
     831        """
     832        Returns the monic (or zero) least common multiple of self and right.
     833       
     834        Corner cases:  if either of self and right are zero, returns zero.
     835        This behaviour is ensures that the relation lcm(a,b) gcd(a,b) == a b
     836        holds up to multiplication by rationals.
     837
     838        EXAMPLES::
     839
     840            sage: R.<t> = QQ[]
     841            sage: f = -2 + 3*t/2 + 4*t^2/7 - t^3
     842            sage: g = 1/2 + 4*t + 2*t^4/3
     843            sage: f.lcm(g)
     844            t^7 - 4/7*t^6 - 3/2*t^5 + 8*t^4 - 75/28*t^3 - 66/7*t^2 + 87/8*t + 3/2
     845            sage: f.lcm(g) * f.gcd(g) // (f * g)
     846            -3/2
     847        """
     848        cdef Polynomial_rational_flint res = self._new()
     849       
     850        _sig_on
     851        fmpq_poly_lcm(res.__poly, self.__poly, \
     852                      (<Polynomial_rational_flint> right).__poly)
     853        _sig_off
     854        return res
     855
     856    @coerce_binop
     857    def xgcd(self, right):
     858        """
     859        Returns polynomials d, s, and t such that d == s * self + t * right,
     860        where d is the (monic) greatest common divisor of self and right.
     861        The choice of s and t is not specified any further.
     862       
     863        Corner cases:  if self and right are zero, returns zero polynomials.
     864        Otherwise, if only self is zero, returns (d, s, t) = (right, 0, 1) up
     865        to normalisation, and similarly if only right is zero.
     866       
     867        EXAMPLES::
     868       
     869            sage: R.<t> = QQ[]
     870            sage: f = 2/3 + 3/4 * t - t^2
     871            sage: g = -3 + 1/7 * t
     872            sage: f.xgcd(g)
     873            (1, -12/5095, -84/5095*t - 1701/5095)
     874        """
     875        cdef Polynomial_rational_flint d = self._new()
     876        cdef Polynomial_rational_flint s = self._new()
     877        cdef Polynomial_rational_flint t = self._new()
     878       
     879        _sig_on
     880        fmpq_poly_xgcd(d.__poly, s.__poly, t.__poly, self.__poly, (<Polynomial_rational_flint>right).__poly)
     881        _sig_off
     882        return d, s, t
     883
     884    cpdef RingElement _mul_(self, RingElement right):
     885        """
     886        Returns the product of self and right.
     887       
     888        EXAMPLES::
     889
     890            sage: R.<t> = QQ[]
     891            sage: f = -1 + 3*t/2 - t^3
     892            sage: g = 2/3 + 7/3*t + 3*t^2
     893            sage: f * g                           # indirect doctest
     894            -3*t^5 - 7/3*t^4 + 23/6*t^3 + 1/2*t^2 - 4/3*t - 2/3
     895       
     896        TESTS::
     897       
     898            sage: R.<t> = QQ[]
     899            sage: f = R.random_element(2000)
     900            sage: g = R.random_element(2000)
     901            sage: (f + g) * (f - g) == f^2 - g^2  # indirect doctest
     902            True
     903        """
     904        cdef Polynomial_rational_flint op2 = <Polynomial_rational_flint> right
     905        cdef Polynomial_rational_flint res = self._new()
     906        cdef bint do_sig = _do_sig(self.__poly) or _do_sig(op2.__poly)
     907       
     908        if do_sig: _sig_on
     909        fmpq_poly_mul(res.__poly, self.__poly, op2.__poly)
     910        if do_sig: _sig_off
     911        return res
     912
     913    cpdef ModuleElement _lmul_(self, RingElement right):
     914        r"""
     915        Returns right times self, where right is a rational number.
     916       
     917        EXAMPLES::
     918       
     919            sage: R.<t> = QQ[]
     920            sage: f = 3/2*t^3 - t + 1/3
     921            sage: 6 * f                  # indirect doctest
     922            9*t^3 - 6*t + 2
     923        """
     924        cdef Polynomial_rational_flint res = self._new()
     925        cdef bint do_sig = _do_sig(self.__poly)
     926       
     927        if do_sig: _sig_on
     928        fmpq_poly_scalar_mul_mpq(res.__poly, self.__poly, \
     929                                 (<Rational> right).value)
     930        if do_sig: _sig_off
     931        return res
     932
     933    cpdef ModuleElement _rmul_(self, RingElement right):
     934        r"""
     935        Returns self * right, where right is a rational number.
     936       
     937        EXAMPLES::
     938       
     939            sage: R.<t> = QQ[]
     940            sage: f = 3/2*t^3 - t + 1/3
     941            sage: f * 6                   # indirect doctest
     942            9*t^3 - 6*t + 2
     943        """
     944        cdef Polynomial_rational_flint res = self._new()
     945        cdef bint do_sig = _do_sig(self.__poly)
     946       
     947        if do_sig: _sig_on
     948        fmpq_poly_scalar_mul_mpq(res.__poly, self.__poly, \
     949                                 (<Rational> right).value)
     950        if do_sig: _sig_off
     951        return res
     952
     953    def __pow__(Polynomial_rational_flint self, exp, ignored):
     954        """
     955        Returns self raised to the power of exp.
     956       
     957        The corner case of exp == 0 is handled by returning the constant
     958        polynomial 1.  Note that this includes the case 0^0 == 1.
     959       
     960        This method only supports integral values for exp that fit into
     961        a signed long int.
     962       
     963        INPUT:
     964       
     965        - exp - Exponent
     966       
     967        OUTPUT:
     968       
     969        - Polynomial; self raised to the power of exp
     970       
     971        EXAMPLES::
     972
     973            sage: R.<t> = QQ[]
     974            sage: f = 1/2 + 2*t - t^2/3
     975            sage: f^0
     976            1
     977            sage: f^3
     978            -1/27*t^6 + 2/3*t^5 - 23/6*t^4 + 6*t^3 + 23/4*t^2 + 3/2*t + 1/8
     979
     980        TESTS::
     981       
     982            sage: R.<t> = QQ[]
     983            sage: f = 1/2 + t
     984            sage: f^0
     985            1
     986            sage: R(0)^0
     987            1
     988           
     989        We verify the size condition on the exponent::
     990       
     991            sage: R.<t> = QQ[]
     992            sage: (1 + t)^(2^63)
     993            Traceback (most recent call last):
     994            ...
     995            OverflowError: long int too large to convert to int
     996        """
     997        cdef long n
     998        cdef Polynomial_rational_flint res
     999       
     1000        if not PY_TYPE_CHECK(exp, Integer):
     1001            try:
     1002                exp = Integer(exp)
     1003            except TypeError:
     1004                raise TypeError, "non-integral exponents not supported"
     1005       
     1006        n = <long> exp
     1007       
     1008        if n < 0:
     1009            if fmpq_poly_is_zero(self.__poly):
     1010                raise ZeroDivisionError, "negative exponent in power of zero"
     1011            res = self._new()
     1012            _sig_on
     1013            fmpq_poly_power(res.__poly, self.__poly, -n)
     1014            _sig_off
     1015            return ~res
     1016        else:
     1017            res = self._new()
     1018            if self._is_gen:
     1019                fmpq_poly_set_coeff_si(res.__poly, n, 1)
     1020            else:
     1021                _sig_on
     1022                fmpq_poly_power(res.__poly, self.__poly, n)
     1023                _sig_off
     1024            return res
     1025
     1026    def __floordiv__(Polynomial_rational_flint self, right):
     1027        """
     1028        Returns the quotient of self and right obtain by Euclidean division.
     1029       
     1030        EXAMPLES::
     1031       
     1032            sage: R.<t> = QQ[]
     1033            sage: f = t^3 - t/2 + 1/5
     1034            sage: g = 2/3*t - 1
     1035            sage: f // g                       # indirect doctest
     1036            3/2*t^2 + 9/4*t + 21/8
     1037       
     1038        TESTS::
     1039       
     1040            sage: R.<t> = QQ[]
     1041            sage: f  = R.random_element(1000)
     1042            sage: g  = R.random_element(500)
     1043            sage: if g == 0: g = R(1)
     1044            sage: qr = f.quo_rem(g)
     1045            sage: q  = f // g                  # indirect doctest
     1046            sage: qr[0] == q
     1047            True
     1048        """
     1049        cdef Polynomial_rational_flint res
     1050        cdef bint do_sig
     1051       
     1052        if right == 0:
     1053            raise ZeroDivisionError, "division by zero polynomial"
     1054
     1055        if not PY_TYPE_CHECK(right, Polynomial_rational_flint):
     1056            if right in QQ:
     1057                res = self._new()
     1058                do_sig = _do_sig(self.__poly)
     1059               
     1060                if do_sig: _sig_on
     1061                fmpq_poly_scalar_div_mpq(res.__poly, self.__poly,
     1062                                                  (<Rational> QQ(right)).value)
     1063                if do_sig: _sig_off
     1064                return res
     1065           
     1066            right = self._parent(right)
     1067       
     1068        res = self._new()
     1069        _sig_on
     1070        fmpq_poly_floordiv(res.__poly, self.__poly,
     1071                                     (<Polynomial_rational_flint>right).__poly)
     1072        _sig_off
     1073        return res
     1074
     1075    def __mod__(Polynomial_rational_flint self, right):
     1076        """
     1077        Returns the remainder of self and right obtain by Euclidean division.
     1078       
     1079        EXAMPLES::
     1080       
     1081            sage: R.<t> = QQ[]
     1082            sage: f = t^3 - t/2 + 1/5
     1083            sage: g = 2/3*t - 1
     1084            sage: f % g                        # indirect doctest
     1085            113/40
     1086       
     1087        TESTS::
     1088       
     1089            sage: R.<t> = QQ[]
     1090            sage: f  = R.random_element(1000)
     1091            sage: g  = R.random_element(500)
     1092            sage: if g == 0: g = R(1)
     1093            sage: qr = f.quo_rem(g)
     1094            sage: r  = f % g                   # indirect doctest
     1095            sage: qr[1] == r
     1096            True
     1097        """
     1098        cdef Polynomial_rational_flint res
     1099       
     1100        if right == 0:
     1101            raise ZeroDivisionError, "division by zero polynomial"
     1102
     1103        if not PY_TYPE_CHECK(right, Polynomial_rational_flint):
     1104            right = self._parent(right)
     1105
     1106        res = self._new()
     1107        _sig_on
     1108        fmpq_poly_mod(res.__poly, self.__poly,
     1109                                     (<Polynomial_rational_flint>right).__poly)
     1110        _sig_off
     1111        return res
     1112
     1113    ###########################################################################
     1114    # Further methods                                                         #
     1115    ###########################################################################
     1116
     1117    def numerator(self):
     1118        """
     1119        Returns the numerator of self.
     1120       
     1121        Representing self as the quotient of an integer polynomial and
     1122        a positive integer denominator (coprime to the content of the
     1123        polynomial), returns the integer polynomial.
     1124       
     1125        EXAMPLE::
     1126       
     1127            sage: R.<t> = QQ[]
     1128            sage: f = (3 * t^3 + 1) / -3
     1129            sage: f.numerator()
     1130            -3*t^3 - 1
     1131        """
     1132        cdef Polynomial_integer_dense_flint num = \
     1133                                         PY_NEW(Polynomial_integer_dense_flint)
     1134        parent = ZZ[self.variable_name()]
     1135        _sig_on
     1136        Polynomial_integer_dense_flint.__init__(num, parent, x=None, \
     1137                                    check=False, is_gen=False, construct=False)
     1138        fmpz_poly_set(num.__poly, fmpq_poly_numref(self.__poly))
     1139        _sig_off
     1140        return num
     1141
     1142    def denominator(self):
     1143        """
     1144        Returns the denominator of self.
     1145
     1146        EXAMPLE::
     1147       
     1148            sage: R.<t> = QQ[]
     1149            sage: f = (3 * t^3 + 1) / -3
     1150            sage: f.denominator()
     1151            3
     1152        """
     1153        cdef Integer den = PY_NEW(Integer)
     1154        if fmpq_poly_denref(self.__poly) is NULL:
     1155            mpz_set_ui(den.value, 1)
     1156        else:
     1157            fmpz_to_mpz(den.value, fmpq_poly_denref(self.__poly))
     1158        return den
     1159   
     1160    def _derivative(self, var = None):
     1161        """
     1162        Returns the derivative of self with respect to ``var``.
     1163       
     1164        INPUT:
     1165       
     1166        -  ``var`` - Must be either (equal to) the generator of the polynomial
     1167           ring to which this polynomial belongs, or ``None``; either way the
     1168           behaviour is the same.
     1169       
     1170        OUTPUT:
     1171       
     1172        -  Derivative as a ``Polynomial_rational_flint``
     1173       
     1174        .. seealso:: :meth:`.derivative`
     1175       
     1176        EXAMPLES::
     1177           
     1178            sage: R.<x> = QQ[]
     1179            sage: f = x^4 - x - 1
     1180            sage: f._derivative()
     1181            4*x^3 - 1
     1182            sage: f._derivative(None)
     1183            4*x^3 - 1
     1184            sage: f._derivative(2*x)
     1185            Traceback (most recent call last):
     1186            ...
     1187            ValueError: Cannot differentiate with respect to 2*x
     1188            sage: y = var("y")
     1189            sage: f._derivative(y)
     1190            Traceback (most recent call last):
     1191            ...
     1192            ValueError: Cannot differentiate with respect to y
     1193        """
     1194        cdef Polynomial_rational_flint der
     1195        cdef bint do_sig
     1196       
     1197        if var is not None and var != self._parent.gen():
     1198            raise ValueError, "Cannot differentiate with respect to %s" %var
     1199       
     1200        der = self._new()
     1201        do_sig = _do_sig(self.__poly)
     1202       
     1203        if do_sig: _sig_on
     1204        fmpq_poly_derivative(der.__poly, self.__poly)
     1205        if do_sig: _sig_off
     1206        return der
     1207   
     1208    def real_root_intervals(self):
     1209        """
     1210        Returns isolating intervals for the real roots of self.
     1211
     1212        EXAMPLES:
     1213       
     1214        We compute the roots of the characteristic polynomial of some
     1215        Salem numbers::
     1216
     1217            sage: R.<t> = QQ[]
     1218            sage: f = 1 - t^2 - t^3 - t^4 + t^6
     1219            sage: f.real_root_intervals()
     1220            [((1/2, 3/4), 1), ((1, 3/2), 1)]
     1221        """
     1222        from sage.rings.polynomial.real_roots import real_roots
     1223        return real_roots(self)
     1224
     1225    @coerce_binop
     1226    def resultant(Polynomial_rational_flint self, right):
     1227        r"""
     1228        Returns the resultant of self and right.
     1229       
     1230        Enumerating the roots over `\QQ` as `r_1, \cdots, r_m` and
     1231        `s_1, \cdots, s_n` and letting `x` and `y` denote the leading
     1232        coefficients of `f` and `g`, the resultant of the two polynomials
     1233        is defined by
     1234       
     1235        .. math::
     1236       
     1237            x^{\deg g} y^{\deg f} \prod_{i,j} (r_i - s_j).
     1238       
     1239        Corner cases:  if one of the polynomials is zero, the resultant
     1240        is zero.  Note that otherwise if one of the polynomials is constant,
     1241        the last term in the above is the empty product.
     1242       
     1243        EXAMPLES::
     1244       
     1245            sage: R.<t> = QQ[]
     1246            sage: f = (t - 2/3) * (t + 4/5) * (t - 1)
     1247            sage: g = (t - 1/3) * (t + 1/2) * (t + 1)
     1248            sage: f.resultant(g)
     1249            119/1350
     1250            sage: h = (t - 1/3) * (t + 1/2) * (t - 1)
     1251            sage: f.resultant(h)
     1252            0
     1253        """
     1254        cdef Rational res = PY_NEW(Rational)
     1255        _sig_on
     1256        fmpq_poly_resultant(res.value, self.__poly, \
     1257                            (<Polynomial_rational_flint>right).__poly)
     1258        _sig_off
     1259        return res
     1260   
     1261    def is_irreducible(self):
     1262        r"""
     1263        Returns whether self is irreducible.
     1264
     1265        This method computes the primitive part as an element of `\ZZ[t]` and
     1266        calls the method ``is_irreducible`` for elements of that polynomial
     1267        ring.
     1268
     1269        By definition, over any integral domain, an element `r` is irreducible
     1270        if and only if it is non-zero, not a unit and whenever `r = ab` then 
     1271        `a` or `b` is a unit.
     1272
     1273        OUTPUT:
     1274
     1275        -  ``bool`` - Whether this polynomial is irreducible
     1276
     1277        EXAMPLES::
     1278
     1279            sage: R.<t> = QQ[]
     1280            sage: (t^2 + 2).is_irreducible()
     1281            True
     1282            sage: (t^2 - 1).is_irreducible()
     1283            False
     1284
     1285        TESTS::
     1286
     1287            sage: R.<t> = QQ[]
     1288            sage: R(0).is_irreducible()
     1289            False
     1290            sage: R(-1/2).is_irreducible()
     1291            False
     1292            sage: (t+1).is_irreducible()
     1293            True
     1294        """
     1295        cdef Polynomial_integer_dense_flint primitive
     1296        cdef unsigned long length = fmpq_poly_length(self.__poly)
     1297       
     1298        if length < 2:
     1299            return False
     1300        elif length == 2:
     1301            return True
     1302        else:
     1303            primitive = PY_NEW(Polynomial_integer_dense_flint)
     1304            parent = ZZ[self.variable_name()]
     1305            _sig_on
     1306            Polynomial_integer_dense_flint.__init__(primitive, parent, \
     1307                             x=None, check=True, is_gen=False, construct=False)
     1308            fmpz_poly_primitive_part(primitive.__poly, \
     1309                                     fmpq_poly_numref(self.__poly))
     1310            _sig_off
     1311            return primitive.is_irreducible()
     1312
     1313    ###########################################################################
     1314    # Methods using PARI                                                      #
     1315    ###########################################################################
     1316
     1317    def galois_group(self, pari_group = False, algorithm = 'pari'):
     1318        """
     1319        Returns the Galois group of self as a permutation group.
     1320       
     1321        INPUT:
     1322       
     1323        -  ``self`` - Irreducible polynomial
     1324       
     1325        -  ``pari_group`` - bool (default: ``False``); if ``True`` instead
     1326           return the Galois group as a PARI group.  This has a useful label
     1327           in it, and may be slightly faster since it doesn't require looking
     1328           up a group in Gap.  To get a permutation group from a PARI
     1329           group ``P``, type ``PermutationGroup(P)``.
     1330       
     1331        -  ``algorithm`` - ``'pari'``, ``'kash'``, ``'magma'`` (default:
     1332           ``'pari'``, except when the degree is at least 12 in which case
     1333           ``'kash'`` is tried).
     1334       
     1335        OUTPUT:
     1336       
     1337        -  Galois group
     1338       
     1339        ALGORITHM:
     1340       
     1341        The Galois group is computed using PARI in C library mode, or possibly
     1342        KASH or MAGMA.
     1343       
     1344        .. note::
     1345       
     1346            The PARI documentation contains the following warning: The method
     1347            used is that of resolvent polynomials and is sensitive to the
     1348            current precision. The precision is updated internally but, in very
     1349            rare cases, a wrong result may be returned if the initial precision
     1350            was not sufficient.
     1351           
     1352            MAGMA does not return a provably correct result.  Please see the
     1353            MAGMA documentation for how to obtain a provably correct result.
     1354       
     1355        EXAMPLES::
     1356       
     1357            sage: R.<x> = QQ[]
     1358            sage: f = x^4 - 17*x^3 - 2*x + 1
     1359            sage: G = f.galois_group(); G            # optional - database_gap
     1360            Transitive group number 5 of degree 4
     1361            sage: G.gens()                           # optional - database_gap
     1362            [(1,2,3,4), (1,2)]
     1363            sage: G.order()                          # optional - database_gap
     1364            24
     1365       
     1366        It is potentially useful to instead obtain the corresponding PARI
     1367        group, which is little more than a 4-tuple.  See the PARI manual for
     1368        the exact details.  (Note that the third entry in the tuple is in the
     1369        new standard ordering.)
     1370       
     1371        ::
     1372       
     1373            sage: f = x^4 - 17*x^3 - 2*x + 1
     1374            sage: G = f.galois_group(pari_group=True); G
     1375            PARI group [24, -1, 5, "S4"] of degree 4
     1376            sage: PermutationGroup(G)                # optional - database_gap
     1377            Transitive group number 5 of degree 4
     1378       
     1379        You can use KASH to compute Galois groups as well.  The advantage is
     1380        that KASH can compute Galois groups of fields up to degree 23, whereas
     1381        PARI only goes to degree 11.  (In my not-so-thorough experiments PARI
     1382        is faster than KASH.)
     1383       
     1384        ::
     1385       
     1386            sage: f = x^4 - 17*x^3 - 2*x + 1
     1387            sage: f.galois_group(algorithm='kash')   # optional - kash
     1388            Transitive group number 5 of degree 4
     1389       
     1390            sage: f = x^4 - 17*x^3 - 2*x + 1
     1391            sage: f.galois_group(algorithm='magma')  # optional - magma, database_gap
     1392            Transitive group number 5 of degree 4
     1393       
     1394        TESTS:
     1395
     1396        We illustrate the behaviour in the case of reducible polynomials::
     1397       
     1398            sage: R.<t> = QQ[]
     1399            sage: f = (1 + t)^2
     1400            sage: f.galois_group()
     1401            Traceback (most recent call last):
     1402            ...
     1403            ValueError: The polynomial must be irreducible
     1404        """
     1405        from sage.groups.all import PariGroup, PermutationGroup, TransitiveGroup
     1406       
     1407        if not self.is_irreducible():
     1408            raise ValueError, "The polynomial must be irreducible"
     1409       
     1410        if self.degree() > 11 and algorithm=='pari':
     1411            algorithm = 'kash'
     1412       
     1413        if algorithm == 'pari':
     1414            G = self._pari_().Polrev().polgalois()
     1415            H = PariGroup(G, self.degree())
     1416            if pari_group:
     1417                return H
     1418            else:
     1419                return PermutationGroup(H)
     1420       
     1421        elif algorithm == 'kash':
     1422            try:
     1423                from sage.interfaces.all import kash
     1424                kash.eval('X := PolynomialRing(RationalField()).1')
     1425                s = self._repr(name='X')
     1426                G = kash('Galois(%s)'%s)
     1427                d = int(kash.eval('%s.ext1'%G.name()))
     1428                n = int(kash.eval('%s.ext2'%G.name()))
     1429                return TransitiveGroup(d, n)
     1430            except RuntimeError, msg:
     1431                raise NotImplementedError, (str(msg) + "\nSorry, " +
     1432                    "computation of Galois groups of fields of degree " +
     1433                    "bigger than 11 is not yet implemented.  Try installing " +
     1434                    "the optional free (closed source) KASH package, which " +
     1435                    "supports larger degrees, or use algorithm='magma' if " +
     1436                    "you have magma.")
     1437       
     1438        elif algorithm == 'magma':
     1439            from sage.interfaces.all import magma
     1440            X = magma(self).GaloisGroup()
     1441            try:
     1442                n, d = X.TransitiveGroupIdentification(nvals=2)
     1443                d = int(d)
     1444                n = int(n)
     1445