Ticket #4000: trac4000_fmpq_poly_c.patch

File trac4000_fmpq_poly_c.patch, 117.9 KB (added by spancratz, 8 years ago)

fmpq_poly.c and header

  • new file sage/libs/flint/fmpq_poly.c

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

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