Ticket #9541: trac_9541-add_fmpq.patch

File trac_9541-add_fmpq.patch, 135.2 KB (added by was, 7 years ago)
  • module_list.py

    # HG changeset patch
    # User William Stein <wstein@gmail.com>
    # Date 1279483889 -7200
    # Node ID 17f4c5d3dd4b40cfc0168648d1f2d5c844a84ac5
    # Parent  97f2f6b0bc622cf1ded811a229240c1376024c69
    [mq]: trac_9541-add_fmpq.patch
    
    diff -r 97f2f6b0bc62 -r 17f4c5d3dd4b module_list.py
    a b  
    468468              extra_compile_args=["-std=c99", "-D_XPG6"],
    469469              depends = [SAGE_ROOT + "/local/include/FLINT/flint.h"]),
    470470
     471    Extension('sage.libs.flint.fmpq_poly_example',
     472              sources = ["sage/libs/flint/fmpq_poly_example.pyx", "sage/libs/flint/fmpq_poly.c"],
     473              libraries = ["csage", "flint", "gmp", "gmpxx", "m", "stdc++"],
     474              include_dirs = [SAGE_ROOT+'/local/include/FLINT/'],
     475              extra_compile_args=["-std=c99", "-D_XPG6"],
     476              depends = [SAGE_ROOT + "/local/include/FLINT/flint.h"]),
     477
    471478    Extension('sage.libs.flint.fmpz_poly',
    472479              sources = ["sage/libs/flint/fmpz_poly.pyx"],
    473480              libraries = ["csage", "flint", "gmp", "gmpxx", "m", "stdc++"],
  • new file sage/libs/flint/fmpq_poly.c

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

    diff -r 97f2f6b0bc62 -r 17f4c5d3dd4b sage/libs/flint/fmpq_poly.h
    - +  
     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#include <stdio.h>
     15#include <gmp.h>
     16
     17#include "fmpz.h"
     18#include "fmpz_poly.h"
     19
     20/**
     21 * \file     fmpq_poly.h
     22 * \brief    Fast implementation of the rational polynomial ring
     23 * \author   Sebastian Pancratz
     24 * \date     Mar 2010 -- July 2010
     25 * \version  0.1.2
     26 */
     27
     28/**
     29 * \ingroup  Definitions
     30 * \brief    Data type for a rational polynomial.
     31 * \details  We represent a rational polynomial as the quotient of an integer
     32 *           polynomial of type <tt>fmpz_poly_t</tt> and a denominator of type
     33 *           <tt>fmpz_t</tt>, enforcing coprimality and that the denominator
     34 *           is positive.  Moreover, the denominator may be <tt>NULL</tt>,
     35 *           which is interpreted as \f$1\f$.  The zero polynomial is
     36 *           represented as \f$0/1\f$.
     37 */
     38typedef struct
     39{
     40    fmpz_poly_p num; //!< Numerator
     41    fmpz_t den;      //!< Denominator
     42} fmpq_poly_struct;
     43
     44/**
     45 * \ingroup  Definitions
     46 * \brief    Array of #fmpq_poly_struct of length one.
     47 * \details  This allows passing <em>by reference</em> without having to
     48 *           explicitly allocate memory for the structure, as one would have
     49 *           to when using pointers.
     50 */
     51typedef fmpq_poly_struct fmpq_poly_t[1];
     52
     53/**
     54 * \ingroup  Definitions
     55 * \brief    Pointer to #fmpq_poly_struct.
     56 */
     57typedef fmpq_poly_struct * fmpq_poly_ptr;
     58
     59void fmpq_poly_canonicalize(fmpq_poly_ptr f, fmpz_t temp);
     60
     61///////////////////////////////////////////////////////////////////////////////
     62// _fmpq_poly_*                                                              //
     63///////////////////////////////////////////////////////////////////////////////
     64
     65/**
     66 * \ingroup  Helper
     67 *
     68 * Ensures that <tt>op-\>den</tt> can hold <tt>limbs</tt> limbs (plus a
     69 * sign/ size limb).
     70 *
     71 * This method never shrinks <tt>op-\>den</tt>.  Moreover, after a call to
     72 * this method, <tt>op-\>den</tt> is guaranteed to be initialised.
     73 */
     74static inline
     75void _fmpq_poly_den_fit_limbs(fmpq_poly_ptr op, unsigned long limbs)
     76{
     77    if (op->den == NULL)
     78    {
     79        op->den = fmpz_init(limbs > 1 ? limbs : 1);
     80        fmpz_set_si(op->den, 1);  /* XXX: FLINT, Valgrind problem */
     81    }
     82    else if (fmpz_size(op->den) < limbs)
     83        op->den = fmpz_realloc(op->den, limbs);
     84}
     85
     86/**
     87 * \ingroup  Helper
     88 *
     89 * Compares the denominators of <tt>op1</tt> and <tt>op2</tt>.
     90 *
     91 * This method handles <tt>NULL</tt> values correctly.  This method only deals
     92 * with the <tt>den</tt> attributes of the arguments.  In particular, it does
     93 * not ensure that the contents of the numerators are coprime to the
     94 * denominators.
     95 *
     96 * All arguments may be aliased.
     97 */
     98static inline
     99int _fmpq_poly_den_equal(fmpq_poly_ptr op1, fmpq_poly_ptr op2)
     100{
     101    if (op1->den == NULL)
     102        return op2->den == NULL || fmpz_is_one(op2->den);
     103    else
     104        return op2->den == NULL ? fmpz_is_one(op1->den)
     105                                : fmpz_equal(op1->den, op2->den);
     106}
     107
     108/**
     109 * \ingroup  Helper
     110 *
     111 * Returns whether the denominator of <tt>op</tt> is one.
     112 *
     113 * This method does not check that the denominator is coprime to the
     114 * numerator.
     115 */
     116static inline
     117int _fmpq_poly_den_is_one(fmpq_poly_ptr op)
     118{
     119    return op->den == NULL || fmpz_is_one(op->den);
     120}
     121
     122///////////////////////////////////////////////////////////////////////////////
     123// fmpq_poly_*                                                               //
     124///////////////////////////////////////////////////////////////////////////////
     125
     126///////////////////////////////////////////////////////////////////////////////
     127// Accessing numerator and denominator
     128
     129/**
     130 * \def      fmpq_poly_numref(op)
     131 * \brief    Returns a reference to the numerator of \c op.
     132 * \ingroup  NumDen
     133 * \details  The <tt>fmpz_poly</tt> functions can be used on the result of
     134 *           this macro.
     135 */
     136#define fmpq_poly_numref(op)  ((op)->num)
     137
     138/**
     139 * \def      fmpq_poly_denref(op)
     140 * \brief    Returns a reference to the denominator of \c op.
     141 * \ingroup  NumDen
     142 * \details  The <tt>fmpz</tt> functions can be used on the result of
     143 *           this macro \e if the result is not <tt>NULL</tt>.
     144 */
     145#define fmpq_poly_denref(op)  ((op)->den)
     146
     147/**
     148 * \brief    Sets \c rop to the numerator of \c op.
     149 * \ingroup  NumDen
     150 *
     151 * \note  The arguments <tt>rop</tt> and <tt>op-\>num</tt> may be aliased.
     152 *
     153 * \param[out] rop  Integer polynomial
     154 * \param[in] op    Rational polynomial
     155 */
     156static inline
     157void fmpq_poly_get_num(fmpz_poly_p rop, const fmpq_poly_ptr op)
     158{
     159    fmpz_poly_set(rop, op->num);
     160}
     161
     162/**
     163 * \brief    Sets \c rop to the denominator of \c op.
     164 * \ingroup  NumDen
     165 *
     166 * This method handles <tt>NULL</tt> values correctly.
     167 *
     168 * \note  The arguments <tt>rop</tt> and <tt>op-\>den</tt> may be aliased.
     169 *
     170 * \param[out] rop  Integer
     171 * \param[in] op    Rational polynomial
     172 */
     173static inline
     174void fmpq_poly_get_den(fmpz_t rop, const fmpq_poly_ptr op)
     175{
     176    if (rop == op->den)  /* Aliasing */
     177        return;
     178   
     179    if (op->den == NULL)
     180        fmpz_set_si(rop, 1);
     181    else
     182    {
     183        fmpz_set_si(rop, 1);
     184        if (fmpz_size(rop) < fmpz_size(op->den))
     185            rop = fmpz_realloc(rop, fmpz_size(op->den));
     186        fmpz_set(rop, op->den);
     187    }
     188}
     189
     190/**
     191 * \brief    Sets the numerator of \c rop to the polynomial \c op.
     192 * \ingroup  NumDen
     193 *
     194 * \note  The arguments <tt>rop-\>num</tt> and <tt>op</tt> may be aliased.
     195 *
     196 * \param[out] rop  Rational polynomial
     197 * \param[in] op    Integer polynomial
     198 */
     199static inline
     200void fmpq_poly_set_num(fmpq_poly_ptr rop, const fmpz_poly_p op)
     201{
     202    fmpz_poly_set(rop->num, op);
     203}
     204
     205/**
     206 * \brief    Sets the denominator of \c rop to the integer \c op.
     207 * \ingroup  NumDen
     208 *
     209 * This method handles <tt>NULL</tt> values correctly.  It does \e not ensure
     210 * that \c rop ends up in canonical form.
     211 *
     212 * The arguments <tt>rop-\>den</tt> and <tt>op</tt> may be aliased.
     213 *
     214 * \param[out] rop  Rational polynomial
     215 * \param[in] op    Integer
     216 */
     217static inline
     218void fmpq_poly_set_den(fmpq_poly_ptr rop, const fmpz_t op)
     219{
     220    if (op == NULL)
     221    {
     222        if (rop->den != NULL)
     223            fmpz_set_si(rop->den, 1);
     224    }
     225    else
     226    {
     227        if (rop->den == NULL)
     228        {
     229            if (!fmpz_is_one(op))
     230            {
     231                rop->den = fmpz_init(fmpz_size(op));
     232                fmpz_set(rop->den, op);
     233            }
     234        }
     235        else
     236        {
     237            if (fmpz_size(rop->den) < fmpz_size(op))
     238                rop->den = fmpz_realloc(rop->den, fmpz_size(op));
     239            fmpz_set(rop->den, op);
     240        }
     241    }
     242}
     243
     244///////////////////////////////////////////////////////////////////////////////
     245// Memory management
     246
     247/**
     248 * \ingroup  MemoryManagement
     249 *
     250 * Initializes the element \c rop to zero.
     251 *
     252 * This function should be called once for the #fmpq_poly_ptr \c rop, before
     253 * using it with other <tt>fmpq_poly</tt> functions, or following a
     254 * preceeding call to #fmpq_poly_clear().
     255 */
     256static inline
     257void fmpq_poly_init(fmpq_poly_ptr rop)
     258{
     259    rop->num = (fmpz_poly_p) malloc(sizeof(fmpz_poly_struct));
     260    fmpz_poly_init(rop->num);
     261    rop->den = NULL;
     262}
     263
     264/**
     265 * \ingroup  MemoryManagement
     266 *
     267 * Clears the element \c rop.
     268 *
     269 * This function should only be called on an element which has been
     270 * initialised.
     271 */
     272static inline
     273void fmpq_poly_clear(fmpq_poly_ptr rop)
     274{
     275    if (rop->num != NULL)
     276    {
     277        fmpz_poly_clear(rop->num);
     278        free(rop->num);
     279        rop->num = NULL;
     280    }
     281    if (rop->den != NULL);
     282    {
     283        fmpz_clear(rop->den);
     284        rop->den = NULL;
     285    }
     286}
     287
     288///////////////////////////////////////////////////////////////////////////////
     289// Assignment and basic manipulation
     290
     291/**
     292 * \ingroup  Assignment
     293 *
     294 * Sets the element \c rop to the same value as the element \c op.
     295 */
     296static inline
     297void fmpq_poly_set(fmpq_poly_ptr rop, const fmpq_poly_ptr op)
     298{
     299    if (rop != op)
     300    {
     301        fmpz_poly_set(rop->num, op->num);
     302        if (op->den == NULL)
     303        {
     304            if (rop->den != NULL)
     305            {
     306                fmpz_clear(rop->den);
     307                rop->den = NULL;
     308            }
     309        }
     310        else
     311        {
     312            _fmpq_poly_den_fit_limbs(rop, fmpz_size(op->den));
     313            fmpz_set(rop->den, op->den);
     314        }
     315    }
     316}
     317
     318/**
     319 * \ingroup  Assignment
     320 *
     321 * Sets the element \c rop to the value given by the \c int \c op.
     322 */
     323static inline 
     324void fmpq_poly_set_si(fmpq_poly_ptr rop, const long op)
     325{
     326    fmpz_poly_zero(rop->num);
     327    fmpz_poly_set_coeff_si(rop->num, 0, op);
     328    if (rop->den != NULL)
     329        fmpz_set_si(rop->den, 1);
     330}
     331
     332/**
     333 * \ingroup  Assignment
     334 *
     335 * Sets the element \c rop to the integer \c x.
     336 */
     337static inline
     338void fmpq_poly_set_fmpz(fmpq_poly_ptr rop, const fmpz_t x)
     339{
     340    fmpz_poly_zero(rop->num);
     341    fmpz_poly_set_coeff_fmpz(rop->num, 0, x);
     342    if (rop->den != NULL)
     343        fmpz_set_si(rop->den, 1);
     344}
     345
     346/**
     347 * \ingroup  Assignment
     348 *
     349 * Sets the element \c rop to the integer \c x.
     350 */
     351static inline
     352void fmpq_poly_set_mpz(fmpq_poly_ptr rop, const mpz_t x)
     353{
     354    fmpz_poly_zero(rop->num);
     355    fmpz_poly_set_coeff_mpz(rop->num, 0, x);
     356    if (rop->den != NULL)
     357        fmpz_set_si(rop->den, 1);
     358}
     359
     360/**
     361 * \ingroup  Assignment
     362 *
     363 * Sets the element \c rop to the rational \c x.
     364 */
     365static inline
     366void fmpq_poly_set_mpq(fmpq_poly_ptr rop, const mpq_t x)
     367{
     368    fmpz_poly_zero(rop->num);
     369    fmpz_poly_set_coeff_mpz(rop->num, 0, mpq_numref(x));
     370    if (mpz_cmp_si(mpq_denref(x), 1))  /* x->den > 1 */
     371    {
     372        _fmpq_poly_den_fit_limbs(rop, mpz_size(mpq_denref(x)));
     373        mpz_to_fmpz(rop->den, mpq_denref(x));
     374    }
     375    else                               /* x.den == 1 */
     376        if (rop->den != NULL)
     377            fmpz_set_si(rop->den, 1);
     378}
     379
     380/**
     381 * \ingroup  Assignment
     382 *
     383 * Swaps the elements \c op1 and \c op2.
     384 *
     385 * This is done efficiently by swapping pointers.
     386 */
     387static inline
     388void fmpq_poly_swap(fmpq_poly_ptr op1, fmpq_poly_ptr op2)
     389{
     390    fmpz_t t;
     391    if (op1 != op2)
     392    {
     393        fmpz_poly_swap(op1->num, op2->num);
     394        t        = op1->den;
     395        op1->den = op2->den;
     396        op2->den = t;
     397    }
     398}
     399
     400/**
     401 * \ingroup  Assignment
     402 *
     403 * Sets the element \c rop to zero.
     404 */
     405static inline
     406void fmpq_poly_zero(fmpq_poly_ptr rop)
     407{
     408    fmpz_poly_zero(rop->num);
     409    if (rop->den != NULL)
     410        fmpz_set_si(rop->den, 1);
     411}
     412
     413/**
     414 * \ingroup  Assignment
     415 *
     416 * Sets the element \c rop to one.
     417 */
     418static inline
     419void fmpq_poly_one(fmpq_poly_ptr rop)
     420{
     421    fmpz_poly_zero(rop->num);
     422    fmpz_poly_set_coeff_si(rop->num, 0, 1);
     423    if (rop->den != NULL)
     424        fmpz_set_si(rop->den, 1);
     425}
     426
     427void fmpq_poly_neg(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
     428void fmpq_poly_inv(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
     429
     430///////////////////////////////////////////////////////////////////////////////
     431// Setting/ retrieving individual coefficients
     432
     433void fmpq_poly_get_coeff_mpq(mpq_t rop, const fmpq_poly_ptr op, unsigned long i);
     434
     435void fmpq_poly_set_coeff_fmpz(fmpq_poly_ptr rop, const unsigned long i, const fmpz_t x);
     436void fmpq_poly_set_coeff_mpz(fmpq_poly_ptr rop, const unsigned long i, const mpz_t x);
     437void fmpq_poly_set_coeff_mpq(fmpq_poly_ptr rop, const unsigned long i, const mpq_t x);
     438void fmpq_poly_set_coeff_si(fmpq_poly_ptr rop, const unsigned long i, long x);
     439
     440///////////////////////////////////////////////////////////////////////////////
     441// Comparison
     442
     443/**
     444 * \brief    Returns whether \c op is zero.
     445 * \ingroup  Comparison
     446 *
     447 * Returns whether the element \c op is zero.
     448 */
     449static inline
     450int fmpq_poly_is_zero(const fmpq_poly_ptr op)
     451{
     452    return fmpz_poly_degree(op->num) < 0;
     453}
     454
     455/**
     456 * \brief    Returns whether \c op is equal to \f$1\f$.
     457 * \ingroup  Comparison
     458 *
     459 * Returns whether the element \c op is equal to the constant polynomial
     460 * \f$1\f$.
     461 */
     462static inline
     463int fmpq_poly_is_one(const fmpq_poly_ptr op)
     464{
     465    return (op->num->length == 1) && fmpz_is_one(op->num->coeffs)
     466        && _fmpq_poly_den_is_one(op);
     467}
     468
     469int fmpq_poly_equal(const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
     470int fmpq_poly_cmp(const fmpq_poly_ptr left, const fmpq_poly_ptr right);
     471
     472///////////////////////////////////////////////////////////////////////////////
     473// Polynomial parameters
     474
     475/**
     476 * \brief    Returns the length of \c op.
     477 * \ingroup  PolyParameters
     478 * \details  Returns the length of the polynomial \c op, which is one greater
     479 *           than its degree.
     480 */
     481static inline
     482unsigned long fmpq_poly_length(const fmpq_poly_ptr op)
     483{
     484    return fmpz_poly_length(op->num);
     485}
     486
     487/**
     488 * \brief    Returns the degree of \c op.
     489 * \ingroup  Parameters
     490 * \details  Returns the degree of the polynomial \c op.
     491 */
     492static inline
     493long fmpq_poly_degree(const fmpq_poly_ptr op)
     494{
     495    return fmpz_poly_degree(op->num);
     496}
     497
     498///////////////////////////////////////////////////////////////////////////////
     499// Addition/ subtraction
     500
     501void fmpq_poly_add(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
     502void fmpq_poly_sub(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
     503
     504void fmpq_poly_addmul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
     505void fmpq_poly_submul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
     506
     507///////////////////////////////////////////////////////////////////////////////
     508// Scalar multiplication and division
     509
     510void fmpq_poly_scalar_mul_si(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long x);
     511void fmpq_poly_scalar_mul_mpz(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpz_t x);
     512void fmpq_poly_scalar_mul_mpq(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x);
     513
     514void fmpq_poly_scalar_div_si(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long x);
     515void fmpq_poly_scalar_div_mpz(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpz_t x);
     516void fmpq_poly_scalar_div_mpq(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x);
     517
     518///////////////////////////////////////////////////////////////////////////////
     519// Multiplication
     520
     521void fmpq_poly_mul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
     522
     523///////////////////////////////////////////////////////////////////////////////
     524// Division
     525
     526void fmpq_poly_floordiv(fmpq_poly_ptr q, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
     527void fmpq_poly_mod(fmpq_poly_ptr r, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
     528void fmpq_poly_divrem(fmpq_poly_ptr q, fmpq_poly_ptr r, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
     529
     530///////////////////////////////////////////////////////////////////////////////
     531// Powering
     532
     533void fmpq_poly_power(fmpq_poly_ptr rop, const fmpq_poly_ptr op, unsigned long exp);
     534
     535///////////////////////////////////////////////////////////////////////////////
     536// Greatest common divisor
     537
     538void fmpq_poly_gcd(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
     539void 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);
     540void fmpq_poly_lcm(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
     541
     542///////////////////////////////////////////////////////////////////////////////
     543// Derivative
     544
     545void fmpq_poly_derivative(fmpq_poly_ptr rop, fmpq_poly_ptr op);
     546
     547///////////////////////////////////////////////////////////////////////////////
     548// Evaluation
     549
     550void fmpq_poly_evaluate_mpz(mpq_t rop, fmpq_poly_ptr f, const mpz_t a);
     551void fmpq_poly_evaluate_mpq(mpq_t rop, fmpq_poly_ptr f, const mpq_t a);
     552
     553///////////////////////////////////////////////////////////////////////////////
     554// Gaussian content
     555
     556void fmpq_poly_content(mpq_t rop, const fmpq_poly_ptr op);
     557void fmpq_poly_primitive_part(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
     558int fmpq_poly_is_monic(const fmpq_poly_ptr op);
     559void fmpq_poly_monic(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
     560
     561///////////////////////////////////////////////////////////////////////////////
     562// Resultant
     563
     564void fmpq_poly_resultant(mpq_t rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
     565void fmpq_poly_discriminant(mpq_t disc, fmpq_poly_t a);
     566
     567///////////////////////////////////////////////////////////////////////////////
     568// Composition
     569
     570void fmpq_poly_compose(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
     571void fmpq_poly_rescale(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x);
     572
     573///////////////////////////////////////////////////////////////////////////////
     574// Square-free
     575
     576int fmpq_poly_is_squarefree(const fmpq_poly_ptr op);
     577
     578///////////////////////////////////////////////////////////////////////////////
     579// Subpolynomials
     580
     581void fmpq_poly_getslice(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long i, long j);
     582void fmpq_poly_shift(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long n);
     583void fmpq_poly_truncate(fmpq_poly_ptr rop, const fmpq_poly_ptr op, unsigned long n);
     584void fmpq_poly_reverse(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
     585
     586///////////////////////////////////////////////////////////////////////////////
     587// String conversion
     588
     589int fmpq_poly_from_string(fmpq_poly_ptr rop, const char * s);
     590char * fmpq_poly_to_string(const fmpq_poly_ptr op);
     591char * fmpq_poly_to_string_pretty(const fmpq_poly_ptr op, const char * x);
     592
     593#endif
     594
  • new file sage/libs/flint/fmpq_poly.pxd

    diff -r 97f2f6b0bc62 -r 17f4c5d3dd4b sage/libs/flint/fmpq_poly.pxd
    - +  
     1#include "fmpz_poly.pxi"
     2#include "fmpz.pxi"
     3
     4cdef extern from "fmpq_poly.h":
     5    ctypedef void* fmpz_t
     6    ctypedef void* fmpz_poly_p
     7    struct fmpq_poly:
     8        fmpz_poly_p num
     9        fmpz_t den
     10       
     11    ctypedef fmpq_poly fmpq_poly_struct
     12    ctypedef fmpq_poly_struct fmpq_poly_t[1]
     13
     14    void fmpq_poly_init(fmpq_poly_t)
     15    void fmpq_poly_from_string(fmpq_poly_t, char*)
     16    char* fmpq_poly_to_string_pretty(fmpq_poly_t, char*)
     17    void fmpq_poly_mul(fmpq_poly_t, fmpq_poly_t, fmpq_poly_t)
     18    void fmpq_poly_clear(fmpq_poly_t)
  • new file sage/libs/flint/fmpq_poly_example.pyx

    diff -r 97f2f6b0bc62 -r 17f4c5d3dd4b sage/libs/flint/fmpq_poly_example.pyx
    - +  
     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# This file contains a simple example for how to use fmpq_poly package.
     12
     13include "fmpq_poly.pxd"
     14
     15cdef extern from "stdlib.h":
     16    void free(void*)
     17
     18def example():
     19    """
     20    This is an example of how to use fmpq from Cython.
     21   
     22    EXAMPLES::
     23   
     24        sage: import sage.libs.flint.fmpq_poly_example as e
     25        sage: e.example()
     26        f     = 3/5*t + 1/2
     27        g     = 2/7*t^2 - 3/2*t - 1
     28        f * g = 6/35*t^3 - 53/70*t^2 - 27/20*t - 1/2
     29    """
     30    cdef char *str, *strf, *strg
     31    cdef fmpq_poly_t f, g
     32    fmpq_poly_init(f)
     33    fmpq_poly_init(g)
     34    fmpq_poly_from_string(f, "2  1/2 3/5")
     35    fmpq_poly_from_string(g, "3  -1 -3/2 2/7")
     36    strf = fmpq_poly_to_string_pretty(f, "t")
     37    strg = fmpq_poly_to_string_pretty(g, "t")
     38    print("f     = %s"%strf)
     39    print("g     = %s"%strg)
     40    fmpq_poly_mul(f, f, g)
     41    str  = fmpq_poly_to_string_pretty(f, "t")
     42    print("f * g = %s"%str);
     43    free(str)
     44    free(strf)
     45    free(strg)
     46    fmpq_poly_clear(f)
     47    fmpq_poly_clear(g)
     48
     49