Ticket #12173: flint-2.3-sage-5.1.beta2.patch

File flint-2.3-sage-5.1.beta2.patch, 248.5 KB (added by jpflori, 10 years ago)

Fredrik's patch; rebase on 5.1.beta2.

  • module_list.py

    # HG changeset patch
    # User Jean-Pierre Flori <jean-pierre.flor@ssi.gouv.fr>
    # Date 1338464095 -7200
    # Node ID 3d9a05cac4204fbf2008f1125a9577bae5da57bf
    # Parent  6b81b33d7cd4b8a338117932159da75570fb5da5
    [mq]: flint-2.3-sage-5.0.diff
    
    diff --git a/module_list.py b/module_list.py
    a b  
    16371637              include_dirs = ['sage/libs/ntl/']),
    16381638
    16391639    Extension('sage.rings.polynomial.polynomial_rational_flint',
    1640               sources = ['sage/rings/polynomial/polynomial_rational_flint.pyx', 'sage/libs/flint/fmpq_poly.c'],
     1640              sources = ['sage/rings/polynomial/polynomial_rational_flint.pyx'],
    16411641              language = 'c++',
    16421642              extra_compile_args=["-std=c99"] + uname_specific('SunOS', [], ['-D_XPG6']),
    16431643              libraries = ["csage", "flint", "ntl", "gmpxx", "gmp"],
  • sage/algebras/quatalg/quaternion_algebra_element.pyx

    diff --git a/sage/algebras/quatalg/quaternion_algebra_element.pyx b/sage/algebras/quatalg/quaternion_algebra_element.pyx
    a b  
    16131613        x, y, z, w = to_quaternion(parent._base, v)
    16141614        cdef NumberFieldElement a = <NumberFieldElement>(parent._base(parent._a))
    16151615        cdef NumberFieldElement b = <NumberFieldElement>(parent._base(parent._b))
    1616         ZZX_to_fmpz_poly(self.x, (<NumberFieldElement>x).__numerator)
    1617         ZZX_to_fmpz_poly(self.y, (<NumberFieldElement>y).__numerator)
    1618         ZZX_to_fmpz_poly(self.z, (<NumberFieldElement>z).__numerator)
    1619         ZZX_to_fmpz_poly(self.w, (<NumberFieldElement>w).__numerator)
     1616        fmpz_poly_set_ZZX(self.x, (<NumberFieldElement>x).__numerator)
     1617        fmpz_poly_set_ZZX(self.y, (<NumberFieldElement>y).__numerator)
     1618        fmpz_poly_set_ZZX(self.z, (<NumberFieldElement>z).__numerator)
     1619        fmpz_poly_set_ZZX(self.w, (<NumberFieldElement>w).__numerator)
    16201620       
    16211621        ZZ_to_mpz(&T1, &(<NumberFieldElement>x).__denominator)
    16221622        ZZ_to_mpz(&T2, &(<NumberFieldElement>y).__denominator)
     
    16371637        fmpz_poly_scalar_mul_mpz(self.z, self.z, t3)
    16381638        fmpz_poly_scalar_mul_mpz(self.w, self.w, t4)
    16391639
    1640         ZZX_to_fmpz_poly(self.a, a.__numerator)     # we will assume that the denominator of a and b are 1
    1641         ZZX_to_fmpz_poly(self.b, b.__numerator)
     1640        fmpz_poly_set_ZZX(self.a, a.__numerator)     # we will assume that the denominator of a and b are 1
     1641        fmpz_poly_set_ZZX(self.b, b.__numerator)
    16421642
    1643         ZZX_to_fmpz_poly(self.modulus, (<NumberFieldElement>x).__fld_numerator.x) # and same for the modulus
     1643        fmpz_poly_set_ZZX(self.modulus, (<NumberFieldElement>x).__fld_numerator.x) # and same for the modulus
    16441644
    16451645    def __getitem__(self, int i):
    16461646        """
     
    16691669        cdef NumberFieldElement item = el._new()
    16701670       
    16711671        if i == 0:
    1672             fmpz_poly_to_ZZX(item.__numerator, self.x)
     1672            fmpz_poly_get_ZZX(item.__numerator, self.x)
    16731673        elif i == 1:
    1674             fmpz_poly_to_ZZX(item.__numerator, self.y)
     1674            fmpz_poly_get_ZZX(item.__numerator, self.y)
    16751675        elif i == 2:
    1676             fmpz_poly_to_ZZX(item.__numerator, self.z)
     1676            fmpz_poly_get_ZZX(item.__numerator, self.z)
    16771677        elif i == 3:
    1678             fmpz_poly_to_ZZX(item.__numerator, self.w)
     1678            fmpz_poly_get_ZZX(item.__numerator, self.w)
    16791679        else:
    16801680            raise IndexError, "quaternion element index out of range"
    16811681
     
    19881988        # Implementation-wise, we compute the GCD's one at a time,
    19891989        # and quit if it ever becomes one
    19901990 
    1991         cdef fmpz_t content = fmpz_init(fmpz_poly_max_limbs(self.x)) # TODO: think about how big this should be (probably the size of d)
    1992                                                                      # Note that we have to allocate this here, and not
    1993                                                                      # as a global variable, because fmpz_t's do not
    1994                                                                      # self allocate memory
     1991        cdef fmpz_t content
     1992        fmpz_init(content)
    19951993        fmpz_poly_content(content, self.x)
    1996         fmpz_to_mpz(U1, content)
     1994        fmpz_get_mpz(U1, content)
    19971995        mpz_gcd(U1, self.d, U1)
    19981996        if mpz_cmp_ui(U1, 1) != 0:
    19991997            fmpz_poly_content(content, self.y)
    2000             fmpz_to_mpz(U2, content)
     1998            fmpz_get_mpz(U2, content)
    20011999            mpz_gcd(U1, U1, U2)
    20022000            if mpz_cmp_ui(U1, 1) != 0:
    20032001                fmpz_poly_content(content, self.z)
    2004                 fmpz_to_mpz(U2, content)
     2002                fmpz_get_mpz(U2, content)
    20052003                mpz_gcd(U1, U1, U2)
    20062004                if mpz_cmp_ui(U1, 1) != 0:
    20072005                    fmpz_poly_content(content, self.w)
    2008                     fmpz_to_mpz(U2, content)
     2006                    fmpz_get_mpz(U2, content)
    20092007                    mpz_gcd(U1, U1, U2)
    20102008                    if mpz_cmp_ui(U1, 1) != 0:
    2011                         fmpz_poly_scalar_div_mpz(self.x, self.x, U1)
    2012                         fmpz_poly_scalar_div_mpz(self.y, self.y, U1)
    2013                         fmpz_poly_scalar_div_mpz(self.z, self.z, U1)
    2014                         fmpz_poly_scalar_div_mpz(self.w, self.w, U1)
     2009                        fmpz_poly_scalar_divexact_mpz(self.x, self.x, U1)
     2010                        fmpz_poly_scalar_divexact_mpz(self.y, self.y, U1)
     2011                        fmpz_poly_scalar_divexact_mpz(self.z, self.z, U1)
     2012                        fmpz_poly_scalar_divexact_mpz(self.w, self.w, U1)
    20152013                        mpz_divexact(self.d, self.d, U1)
    20162014       
    20172015        fmpz_clear(content)
  • sage/graphs/matchpoly.pyx

    diff --git a/sage/graphs/matchpoly.pyx b/sage/graphs/matchpoly.pyx
    a b  
    364364    """
    365365    cdef int i, j, k, edge1, edge2, new_edge1, new_edge2, new_nedges
    366366    cdef int *edges1, *edges2, *new_edges1, *new_edges2
    367     cdef fmpz_t coeff
     367    cdef fmpz * coeff
    368368
    369369    if nverts == 3:
    370370        coeff = fmpz_poly_get_coeff_ptr(pol, 3)
    371371        if coeff is NULL:
    372372            fmpz_poly_set_coeff_ui(pol, 3, 1)
    373373        else:
    374             fmpz_add_ui_inplace(coeff, 1)
     374            fmpz_add_ui(coeff, coeff, 1)
    375375        coeff = fmpz_poly_get_coeff_ptr(pol, 1)
    376376        if coeff is NULL:
    377377            fmpz_poly_set_coeff_ui(pol, 1, nedges)
    378378        else:
    379             fmpz_add_ui_inplace(coeff, nedges)
     379            fmpz_add_ui(coeff, coeff, nedges)
    380380        return
    381381
    382382    if nedges == 0:
     
    384384        if coeff is NULL:
    385385            fmpz_poly_set_coeff_ui(pol, nverts, 1)
    386386        else:
    387             fmpz_add_ui_inplace(coeff, 1)
     387            fmpz_add_ui(coeff, coeff, 1)
    388388        return
    389389
    390390    edges1 = edges[2*depth]
  • sage/groups/perm_gps/partn_ref/data_structures_pxd.pxi

    diff --git a/sage/groups/perm_gps/partn_ref/data_structures_pxd.pxi b/sage/groups/perm_gps/partn_ref/data_structures_pxd.pxi
    a b  
    1313cdef extern from "stdlib.h":
    1414    int rand()
    1515
    16 cdef extern from "FLINT/long_extras.h":
    17     int z_isprime(unsigned long n)
     16cdef extern from "FLINT/ulong_extras.h":
     17    int n_is_prime(unsigned long n)
    1818
    1919cdef struct OrbitPartition:
    2020    # Disjoint set data structure for representing the orbits of the generators
  • sage/groups/perm_gps/partn_ref/data_structures_pyx.pxi

    diff --git a/sage/groups/perm_gps/partn_ref/data_structures_pyx.pxi b/sage/groups/perm_gps/partn_ref/data_structures_pyx.pxi
    a b  
    14321432            if OP.parent[i] == i:
    14331433                q = OP.size[i]
    14341434                if m < 2*q and q < m-2:
    1435                     if z_isprime(q):
     1435                    if n_is_prime(q):
    14361436                        sage_free(perm)
    14371437                        OP_dealloc(OP)
    14381438                        return True
  • deleted file sage/libs/flint/fmpq_poly.c

    diff --git a/sage/libs/flint/fmpq_poly.c b/sage/libs/flint/fmpq_poly.c
    deleted file mode 100644
    + -  
    1 ///////////////////////////////////////////////////////////////////////////////
    2 // Copyright (C) 2010 Sebastian Pancratz                                     //
    3 //                                                                           //
    4 // Distributed under the terms of the GNU General Public License as          //
    5 // published by the Free Software Foundation; either version 2 of the        //
    6 // License, or (at your option) any later version.                           //
    7 //                                                                           //
    8 // http://www.gnu.org/licenses/                                              //
    9 ///////////////////////////////////////////////////////////////////////////////
    10 
    11 #ifdef __cplusplus
    12     extern "C" {
    13 #endif
    14 
    15 #include "fmpq_poly.h"
    16 
    17 /**
    18  * \file     fmpq_poly.c
    19  * \brief    Fast implementation of the rational function field
    20  * \author   Sebastian Pancratz
    21  * \date     Mar 2010 -- July 2010
    22  * \version  0.1.8
    23  *
    24  * \mainpage
    25  *
    26  * <center>
    27  * A FLINT-based implementation of the rational polynomial ring.
    28  * </center>
    29  *
    30  * \section Overview
    31  *
    32  * The module <tt>fmpq_poly</tt> provides functions for performing
    33  * arithmetic on rational polynomials in \f$\mathbf{Q}[t]\f$, represented as
    34  * quotients of integer polynomials of type <tt>fmpz_poly_t</tt> and
    35  * denominators of type <tt>fmpz_t</tt>.  These functions start with the
    36  * prefix <tt>fmpq_poly_</tt>.
    37  *
    38  * Rational polynomials are stored in objects of type #fmpq_poly_t,
    39  * which is an array of #fmpq_poly_struct's of length one.  This
    40  * permits passing parameters of type #fmpq_poly_t by reference. 
    41  * We also define the type #fmpq_poly_ptr to be a pointer to
    42  * #fmpq_poly_struct's.
    43  *
    44  * The representation of a rational polynomial as the quotient of an integer
    45  * polynomial and an integer denominator can be made canonical by demanding
    46  * the numerator and denominator to be coprime and the denominator to be
    47  * positive.  As the only special case, we represent the zero function as
    48  * \f$0/1\f$.  All arithmetic functions assume that the operands are in this
    49  * canonical form, and canonicalize their result.  If the numerator or
    50  * denominator is modified individually, for example using the methods in the
    51  * group \ref NumDen, it is the user's responsibility to canonicalize the
    52  * rational function using the function #fmpq_poly_canonicalize() if necessary.
    53  *
    54  * All methods support aliasing of their inputs and outputs \e unless
    55  * explicitly stated otherwise, subject to the following caveat.  If
    56  * different rational polynomials (as objects in memory, not necessarily in
    57  * the mathematical sense) share some of the underlying integer polynomials 
    58  * or integers, the behaviour is undefined.
    59  *
    60  * \section Changelog  Version history
    61  * - 0.1.8
    62  *   - Extra checks for <tt>NULL</tt> returns of <tt>malloc</tt> calls
    63  * - 0.1.7
    64  *   - Explicit casts for return values of <tt>malloc</tt>
    65  *   - Changed a few calls to GMP and FLINT methods from <tt>_si</tt>
    66  *     to <tt>_ui</tt>
    67  * - 0.1.6
    68  *   - Added tons of testing code in the subdirectory <tt>test</tt>
    69  *   - Made a few changes throughout the code base indicated by the tests
    70  * - 0.1.5
    71  *   - Re-wrote the entire code to always initialise the denominator
    72  * - 0.1.4
    73  *   - Fixed a bug in #fmpq_poly_divrem()
    74  *   - Swapped calls to #fmpq_poly_degree to calls to #fmpq_poly_length in
    75  *     many places
    76  * - 0.1.3
    77  *   - Changed #fmpq_poly_inv()
    78  *   - Another sign check to <tt>fmpz_abs</tt> in #fmpq_poly_content()
    79  * - 0.1.2
    80  *   - Introduce a function #fmpq_poly_monic() and use this to simplify the
    81  *     code for the gcd and xgcd functions
    82  *   - Make further use of #fmpq_poly_is_zero()
    83  * - 0.1.1
    84  *   - Replaced a few sign checks and negations by <tt>fmpz_abs</tt>
    85  *   - Small changes to comments and the documentation
    86  *   - Moved some function bodies from #fmpq_poly.h to #fmpq_poly.c
    87  * - 0.1.0
    88  *   - First draft, based on the author's Sage code
    89  */
    90 
    91 /**
    92  * \defgroup Definitions        Type definitions
    93  * \defgroup MemoryManagement   Memory management
    94  * \defgroup NumDen             Accessing numerator and denominator
    95  * \defgroup Assignment         Assignment and basic manipulation
    96  * \defgroup Coefficients       Setting and retrieving individual coefficients
    97  * \defgroup PolyParameters     Polynomial parameters
    98  * \defgroup Comparison         Comparison
    99  * \defgroup Addition           Addition and subtraction
    100  * \defgroup ScalarMul          Scalar multiplication and division
    101  * \defgroup Multiplication     Multiplication
    102  * \defgroup Division           Euclidean division
    103  * \defgroup Powering           Powering
    104  * \defgroup GCD                Greatest common divisor
    105  * \defgroup Derivative         Derivative
    106  * \defgroup Evaluation         Evaluation
    107  * \defgroup Content            Gaussian content
    108  * \defgroup Resultant          Resultant
    109  * \defgroup Composition        Composition
    110  * \defgroup Squarefree         Square-free test
    111  * \defgroup Subpolynomials     Subpolynomials
    112  * \defgroup StringConversions  String conversions and I/O
    113  */
    114 
    115 ///////////////////////////////////////////////////////////////////////////////
    116 // Auxiliary functions                                                       //
    117 ///////////////////////////////////////////////////////////////////////////////
    118 
    119 /**
    120  * Returns the number of digits in the decimal representation of \c n.
    121  */
    122 unsigned 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  */
    154 void 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  */
    230 void 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  */
    252 void 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  */
    304 void 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  */
    316 void 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  */
    344 void 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  */
    374 void 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  */
    427 void 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  */
    461 int 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  */
    479 int 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  */
    586 void _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  */
    618 void 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  */
    703 void _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  */
    721 void 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  */
    808 void 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  */
    827 void 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  */
    844 void 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  */
    896 void 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  */
    956 void 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  */
    987 void _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  */
    1122 void 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  */
    1254 void _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  */
    1365 void 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  */
    1495 void 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  */
    1533 void _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  */
    1568 void 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  */
    1703 void 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  */
    1828 void 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  */
    1931 void 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  */
    2062 void 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  */
    2105 void 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  */
    2177 void fmpq_poly_xgcd(fmpq_poly_ptr rop, fmpq_poly_ptr s, fmpq_poly_ptr t, const fmpq_poly_ptr a, const fmpq_poly_ptr b)
    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  */
    2402 void 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  */
    2455 void 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  */
    2477 void _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  */
    2530 void 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  */
    2569 void 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  */
    2603 void 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  */
    2624 void 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  */
    2645 int 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  */
    2667 void 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  */
    2706 void 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  */
    2906 void 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  */
    2948 void 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  */
    2995 void 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  */
    3102 int 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  */
    3129 void 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  */
    3171 void 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  */
    3201 void 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  */
    3233 void 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  */
    3250 void 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  */
    3310 void _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  */
    3354 int 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  */
    3420 char * 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  */
    3499 char * 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 
  • deleted file sage/libs/flint/fmpq_poly.h

    diff --git a/sage/libs/flint/fmpq_poly.h b/sage/libs/flint/fmpq_poly.h
    deleted file mode 100644
    + -  
    1 ///////////////////////////////////////////////////////////////////////////////
    2 // Copyright (C) 2010 Sebastian Pancratz                                     //
    3 //                                                                           //
    4 // Distributed under the terms of the GNU General Public License as          //
    5 // published by the Free Software Foundation; either version 2 of the        //
    6 // License, or (at your option) any later version.                           //
    7 //                                                                           //
    8 // http://www.gnu.org/licenses/                                              //
    9 ///////////////////////////////////////////////////////////////////////////////
    10 
    11 #ifndef _FMPQ_POLY_H_
    12 #define _FMPQ_POLY_H_
    13 
    14 #ifdef __cplusplus
    15     extern "C" {
    16 #endif
    17 
    18 #include <stdio.h>
    19 #include <gmp.h>
    20 
    21 #include "fmpz.h"
    22 #include "fmpz_poly.h"
    23 
    24 /**
    25  * \file     fmpq_poly.h
    26  * \brief    Fast implementation of the rational polynomial ring
    27  * \author   Sebastian Pancratz
    28  * \date     Mar 2010 -- July 2010
    29  * \version  0.1.8
    30  */
    31 
    32 /**
    33  * \ingroup  Definitions
    34  * \brief    Data type for a rational polynomial.
    35  * \details  We represent a rational polynomial as the quotient of an integer
    36  *           polynomial of type <tt>fmpz_poly_t</tt> and a denominator of type
    37  *           <tt>fmpz_t</tt>, enforcing coprimality and that the denominator
    38  *           is positive.  The zero polynomial is represented as \f$0/1\f$.
    39  */
    40 typedef 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  */
    53 typedef fmpq_poly_struct fmpq_poly_t[1];
    54 
    55 /**
    56  * \ingroup  Definitions
    57  * \brief    Pointer to #fmpq_poly_struct.
    58  */
    59 typedef fmpq_poly_struct * fmpq_poly_ptr;
    60 
    61 void 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  */
    100 static inline
    101 void 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  */
    116 static inline
    117 void 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  */
    131 static inline
    132 void 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 same value as the element \c op.
    147  */
    148 static inline
    149 void fmpq_poly_set_fmpz_poly(fmpq_poly_ptr rop, const fmpz_poly_t op)
    150 {
    151     fmpz_poly_set(rop->num, op);
    152     fmpz_set_ui(rop->den, 1);
    153 }
    154 
    155 /**
    156  * \ingroup  Assignment
    157  *
    158  * Sets the element \c rop to the value given by the \c int \c op.
    159  */
    160 static inline 
    161 void fmpq_poly_set_si(fmpq_poly_ptr rop, long op)
    162 {
    163     fmpz_poly_zero(rop->num);
    164     fmpz_poly_set_coeff_si(rop->num, 0, op);
    165     fmpz_set_ui(rop->den, 1ul);
    166 }
    167 
    168 /**
    169  * \ingroup  Assignment
    170  *
    171  * Sets the element \c rop to the integer \c x.
    172  */
    173 static inline
    174 void fmpq_poly_set_fmpz(fmpq_poly_ptr rop, const fmpz_t x)
    175 {
    176     fmpz_poly_zero(rop->num);
    177     fmpz_poly_set_coeff_fmpz(rop->num, 0, x);
    178     fmpz_set_ui(rop->den, 1ul);
    179 }
    180 
    181 /**
    182  * \ingroup  Assignment
    183  *
    184  * Sets the element \c rop to the integer \c x.
    185  */
    186 static inline
    187 void fmpq_poly_set_mpz(fmpq_poly_ptr rop, const mpz_t x)
    188 {
    189     fmpz_poly_zero(rop->num);
    190     fmpz_poly_set_coeff_mpz(rop->num, 0, x);
    191     fmpz_set_ui(rop->den, 1ul);
    192 }
    193 
    194 /**
    195  * \ingroup  Assignment
    196  *
    197  * Sets the element \c rop to the rational \c x, assumed to be given
    198  * in lowest terms.
    199  */
    200 static inline
    201 void fmpq_poly_set_mpq(fmpq_poly_ptr rop, const mpq_t x)
    202 {
    203     fmpz_poly_zero(rop->num);
    204     fmpz_poly_set_coeff_mpz(rop->num, 0, mpq_numref(x));
    205     rop->den = fmpz_realloc(rop->den, mpz_size(mpq_denref(x)));
    206     mpz_to_fmpz(rop->den, mpq_denref(x));
    207 }
    208 
    209 /**
    210  * \ingroup  Assignment
    211  *
    212  * Swaps the elements \c op1 and \c op2.
    213  *
    214  * This is done efficiently by swapping pointers.
    215  */
    216 static inline
    217 void fmpq_poly_swap(fmpq_poly_ptr op1, fmpq_poly_ptr op2)
    218 {
    219     if (op1 != op2)
    220     {
    221         fmpz_t t;
    222         fmpz_poly_swap(op1->num, op2->num);
    223         t        = op1->den;
    224         op1->den = op2->den;
    225         op2->den = t;
    226     }
    227 }
    228 
    229 /**
    230  * \ingroup  Assignment
    231  *
    232  * Sets the element \c rop to zero.
    233  */
    234 static inline
    235 void fmpq_poly_zero(fmpq_poly_ptr rop)
    236 {
    237     fmpz_poly_zero(rop->num);
    238     fmpz_set_ui(rop->den, 1ul);
    239 }
    240 
    241 /**
    242  * \ingroup  Assignment
    243  *
    244  * Sets the element \c rop to one.
    245  */
    246 static inline
    247 void fmpq_poly_one(fmpq_poly_ptr rop)
    248 {
    249     fmpz_poly_zero(rop->num);
    250     fmpz_poly_set_coeff_si(rop->num, 0, 1);
    251     fmpz_set_ui(rop->den, 1ul);
    252 }
    253 
    254 void fmpq_poly_neg(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
    255 void fmpq_poly_inv(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
    256 
    257 ///////////////////////////////////////////////////////////////////////////////
    258 // Setting/ retrieving individual coefficients
    259 
    260 void fmpq_poly_get_coeff_mpq(mpq_t rop, const fmpq_poly_ptr op, ulong i);
    261 
    262 void fmpq_poly_set_coeff_fmpz(fmpq_poly_ptr rop, ulong i, const fmpz_t x);
    263 void fmpq_poly_set_coeff_mpz(fmpq_poly_ptr rop, ulong i, const mpz_t x);
    264 void fmpq_poly_set_coeff_mpq(fmpq_poly_ptr rop, ulong i, const mpq_t x);
    265 void fmpq_poly_set_coeff_si(fmpq_poly_ptr rop, ulong i, long x);
    266 
    267 ///////////////////////////////////////////////////////////////////////////////
    268 // Comparison
    269 
    270 /**
    271  * \brief    Returns whether \c op is zero.
    272  * \ingroup  Comparison
    273  *
    274  * Returns whether the element \c op is zero.
    275  */
    276 static inline
    277 int fmpq_poly_is_zero(const fmpq_poly_ptr op)
    278 {
    279     return op->num->length == 0ul;
    280 }
    281 
    282 /**
    283  * \brief    Returns whether \c op is equal to \f$1\f$.
    284  * \ingroup  Comparison
    285  *
    286  * Returns whether the element \c op is equal to the constant polynomial
    287  * \f$1\f$.
    288  */
    289 static inline
    290 int fmpq_poly_is_one(const fmpq_poly_ptr op)
    291 {
    292     return (op->num->length == 1ul) && fmpz_is_one(op->num->coeffs)
    293                                     && fmpz_is_one(op->den);
    294 }
    295 
    296 int fmpq_poly_equal(const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
    297 int fmpq_poly_cmp(const fmpq_poly_ptr left, const fmpq_poly_ptr right);
    298 
    299 ///////////////////////////////////////////////////////////////////////////////
    300 // Polynomial parameters
    301 
    302 /**
    303  * \brief    Returns the length of \c op.
    304  * \ingroup  PolyParameters
    305  * \details  Returns the length of the polynomial \c op, which is one greater
    306  *           than its degree.
    307  */
    308 static inline
    309 ulong fmpq_poly_length(const fmpq_poly_ptr op)
    310 {
    311     return op->num->length;
    312 }
    313 
    314 /**
    315  * \brief    Returns the degree of \c op.
    316  * \ingroup  Parameters
    317  * \details  Returns the degree of the polynomial \c op.
    318  */
    319 static inline
    320 long fmpq_poly_degree(const fmpq_poly_ptr op)
    321 {
    322     return (long) (op->num->length - 1ul);
    323 }
    324 
    325 ///////////////////////////////////////////////////////////////////////////////
    326 // Addition/ subtraction
    327 
    328 void fmpq_poly_add(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
    329 void fmpq_poly_sub(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
    330 
    331 void fmpq_poly_addmul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
    332 void fmpq_poly_submul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
    333 
    334 ///////////////////////////////////////////////////////////////////////////////
    335 // Scalar multiplication and division
    336 
    337 void fmpq_poly_scalar_mul_si(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long x);
    338 void fmpq_poly_scalar_mul_mpz(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpz_t x);
    339 void fmpq_poly_scalar_mul_mpq(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x);
    340 
    341 void fmpq_poly_scalar_div_si(fmpq_poly_ptr rop, const fmpq_poly_ptr op, long x);
    342 void fmpq_poly_scalar_div_mpz(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpz_t x);
    343 void fmpq_poly_scalar_div_mpq(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x);
    344 
    345 ///////////////////////////////////////////////////////////////////////////////
    346 // Multiplication
    347 
    348 void fmpq_poly_mul(fmpq_poly_ptr rop, const fmpq_poly_ptr op1, const fmpq_poly_ptr op2);
    349 
    350 ///////////////////////////////////////////////////////////////////////////////
    351 // Division
    352 
    353 void fmpq_poly_floordiv(fmpq_poly_ptr q, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
    354 void fmpq_poly_mod(fmpq_poly_ptr r, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
    355 void fmpq_poly_divrem(fmpq_poly_ptr q, fmpq_poly_ptr r, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
    356 
    357 ///////////////////////////////////////////////////////////////////////////////
    358 // Powering
    359 
    360 void fmpq_poly_power(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong exp);
    361 
    362 ///////////////////////////////////////////////////////////////////////////////
    363 // Greatest common divisor
    364 
    365 void fmpq_poly_gcd(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
    366 void fmpq_poly_xgcd(fmpq_poly_ptr rop, fmpq_poly_ptr s, fmpq_poly_ptr t, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
    367 void fmpq_poly_lcm(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
    368 
    369 ///////////////////////////////////////////////////////////////////////////////
    370 // Derivative
    371 
    372 void fmpq_poly_derivative(fmpq_poly_ptr rop, fmpq_poly_ptr op);
    373 
    374 ///////////////////////////////////////////////////////////////////////////////
    375 // Evaluation
    376 
    377 void fmpq_poly_evaluate_mpz(mpq_t rop, fmpq_poly_ptr f, const mpz_t a);
    378 void fmpq_poly_evaluate_mpq(mpq_t rop, fmpq_poly_ptr f, const mpq_t a);
    379 
    380 ///////////////////////////////////////////////////////////////////////////////
    381 // Gaussian content
    382 
    383 void fmpq_poly_content(mpq_t rop, const fmpq_poly_ptr op);
    384 void fmpq_poly_primitive_part(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
    385 int fmpq_poly_is_monic(const fmpq_poly_ptr op);
    386 void fmpq_poly_monic(fmpq_poly_ptr rop, const fmpq_poly_ptr op);
    387 
    388 ///////////////////////////////////////////////////////////////////////////////
    389 // Resultant
    390 
    391 void fmpq_poly_resultant(mpq_t rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
    392 void fmpq_poly_discriminant(mpq_t disc, fmpq_poly_t a);
    393 
    394 ///////////////////////////////////////////////////////////////////////////////
    395 // Composition
    396 
    397 void fmpq_poly_compose(fmpq_poly_ptr rop, const fmpq_poly_ptr a, const fmpq_poly_ptr b);
    398 void fmpq_poly_rescale(fmpq_poly_ptr rop, const fmpq_poly_ptr op, const mpq_t x);
    399 
    400 ///////////////////////////////////////////////////////////////////////////////
    401 // Square-free
    402 
    403 int fmpq_poly_is_squarefree(const fmpq_poly_ptr op);
    404 
    405 ///////////////////////////////////////////////////////////////////////////////
    406 // Subpolynomials
    407 
    408 void fmpq_poly_getslice(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong i, ulong j);
    409 void fmpq_poly_left_shift(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong n);
    410 void fmpq_poly_right_shift(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong n);
    411 void fmpq_poly_truncate(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong n);
    412 void fmpq_poly_reverse(fmpq_poly_ptr rop, const fmpq_poly_ptr op, ulong n);
    413 
    414 ///////////////////////////////////////////////////////////////////////////////
    415 // String conversion
    416 
    417 void _fmpq_poly_from_list(fmpq_poly_ptr rop, mpq_t * a, ulong n);
    418 int fmpq_poly_from_string(fmpq_poly_ptr rop, const char * s);
    419 char * fmpq_poly_to_string(const fmpq_poly_ptr op);
    420 char * fmpq_poly_to_string_pretty(const fmpq_poly_ptr op, const char * x);
    421 
    422 #ifdef __cplusplus
    423     }
    424 #endif
    425 
    426 #endif
    427 
  • sage/libs/flint/fmpq_poly.pxd

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

    diff --git a/sage/libs/flint/fmpz.pxi b/sage/libs/flint/fmpz.pxi
    a b  
    11include "../ntl/decl.pxi"
    22
    33cdef extern from "FLINT/fmpz.h":
    4    
    5     ctypedef void * fmpz_t
     4
     5    ctypedef long fmpz   
     6    ctypedef long * fmpz_t
    67    ctypedef void * mpz_t
    78
    8     fmpz_t fmpz_init(unsigned long limbs)
     9    void fmpz_init(fmpz_t x)
    910
    1011    void fmpz_set_ui(fmpz_t res, unsigned long x)
    1112    void fmpz_set_si(fmpz_t res, long x)
     
    1314    void fmpz_clear(fmpz_t f)
    1415    void fmpz_print(fmpz_t f)
    1516    int fmpz_is_one(fmpz_t f)
     17   
     18    void fmpz_get_mpz(mpz_t rop, fmpz_t op)
     19    void fmpz_set_mpz(fmpz_t rop, mpz_t op)
    1620
    17     void fmpz_add_ui_inplace(fmpz_t output, unsigned long x)
    18     void fmpz_sub_ui_inplace(fmpz_t output, unsigned long x)
    19    
    20     void fmpz_to_mpz(mpz_t rop, fmpz_t op)
     21    void fmpz_add_ui(fmpz_t f, fmpz_t g, unsigned long c)
     22
  • sage/libs/flint/fmpz_poly.pxi

    diff --git a/sage/libs/flint/fmpz_poly.pxi b/sage/libs/flint/fmpz_poly.pxi
    a b  
    66    ctypedef void* fmpz_poly_t
    77   
    88    void fmpz_poly_init(fmpz_poly_t poly)
    9     void fmpz_poly_init2(fmpz_poly_t poly, unsigned long alloc, \
    10             unsigned long limbs)
     9    void fmpz_poly_init2(fmpz_poly_t poly, unsigned long alloc)
    1110    void fmpz_poly_realloc(fmpz_poly_t poly, unsigned long alloc)
    1211   
    1312    void fmpz_poly_fit_length(fmpz_poly_t poly, unsigned long alloc)
    14     void fmpz_poly_resize_limbs(fmpz_poly_t poly, unsigned long limbs)
    15     void fmpz_poly_fit_limbs(fmpz_poly_t poly, unsigned long limbs)
    16     unsigned long fmpz_poly_limbs(fmpz_poly_t poly)
    1713
    1814    void fmpz_poly_clear(fmpz_poly_t poly)
    1915   
     
    5450            unsigned long x)
    5551    void fmpz_poly_scalar_mul_si(fmpz_poly_t output, fmpz_poly_t input, long x)
    5652
    57     void fmpz_poly_scalar_mul_mpz(fmpz_poly_t output, fmpz_poly_t poly,
    58             mpz_t x)
     53    void fmpz_poly_scalar_mul_fmpz(fmpz_poly_t output, fmpz_poly_t poly,
     54            fmpz_t x)
     55    void fmpz_poly_scalar_mul_mpz(fmpz_poly_t output, fmpz_poly_t poly, mpz_t x)
    5956
    60     void fmpz_poly_scalar_div_exact_ui(fmpz_poly_t output, fmpz_poly_t poly, \
     57    void fmpz_poly_scalar_divexact_ui(fmpz_poly_t output, fmpz_poly_t poly, \
    6158            unsigned long x)
    62     void fmpz_poly_scalar_div_exact_si(fmpz_poly_t output, fmpz_poly_t poly, \
     59    void fmpz_poly_scalar_divexact_si(fmpz_poly_t output, fmpz_poly_t poly, \
    6360            long x)
    64     void fmpz_poly_scalar_div_exact_fmpz(fmpz_poly_t output, fmpz_poly_t poly, \
     61    void fmpz_poly_scalar_divexact_fmpz(fmpz_poly_t output, fmpz_poly_t poly, \
    6562            fmpz_t x)
    6663   
    67     void fmpz_poly_scalar_div_mpz( fmpz_poly_t output, fmpz_poly_t poly, mpz_t x)
     64    void fmpz_poly_scalar_divexact_mpz( fmpz_poly_t output, fmpz_poly_t poly, mpz_t x)
     65
     66    void fmpz_poly_scalar_fdiv_ui(fmpz_poly_t output, fmpz_poly_t poly, unsigned long x)
     67    void fmpz_poly_scalar_fdiv_mpz(fmpz_poly_t output, fmpz_poly_t poly, mpz_t x)
    6868
    6969    void fmpz_poly_div(fmpz_poly_t Q, fmpz_poly_t A, fmpz_poly_t B)
    7070    void fmpz_poly_divrem(fmpz_poly_t Q, fmpz_poly_t R, fmpz_poly_t A, \
     
    7777   
    7878    int fmpz_poly_equal(fmpz_poly_t poly1, fmpz_poly_t poly2)
    7979
    80     bint fmpz_poly_from_string(fmpz_poly_t poly, char* s)
    81     char* fmpz_poly_to_string(fmpz_poly_t poly)
     80    int fmpz_poly_set_str(fmpz_poly_t poly, char* s)
     81    char* fmpz_poly_get_str(fmpz_poly_t poly)
    8282    void fmpz_poly_print(fmpz_poly_t poly)
    8383    bint fmpz_poly_read(fmpz_poly_t poly)
    8484
    85     void fmpz_poly_power(fmpz_poly_t output, fmpz_poly_t poly, \
     85    void fmpz_poly_pow(fmpz_poly_t output, fmpz_poly_t poly, \
    8686            unsigned long exp)
    87     void fmpz_poly_power_trunc_n(fmpz_poly_t output, fmpz_poly_t poly, \
     87    void fmpz_poly_pow_trunc(fmpz_poly_t output, fmpz_poly_t poly, \
    8888            unsigned long exp, unsigned long n)
    8989
    90     void fmpz_poly_left_shift ( fmpz_poly_t output ,
     90    void fmpz_poly_shift_left ( fmpz_poly_t output ,
    9191                                fmpz_poly_t poly , unsigned long n )
    92     void fmpz_poly_right_shift ( fmpz_poly_t output ,
     92    void fmpz_poly_shift_right ( fmpz_poly_t output ,
    9393                                 fmpz_poly_t poly , unsigned long n )
    9494
    9595    void fmpz_poly_derivative ( fmpz_poly_t der , fmpz_poly_t poly )
  • sage/libs/flint/fmpz_poly.pyx

    diff --git a/sage/libs/flint/fmpz_poly.pyx b/sage/libs/flint/fmpz_poly.pyx
    a b  
    5353        cdef long c
    5454        cdef Integer w
    5555        if PY_TYPE_CHECK(v, str):
    56             if fmpz_poly_from_string(self.poly, v):
     56            if not fmpz_poly_set_str(self.poly, v):
    5757                return
    5858            else:
    5959                raise ValueError, "Unable to create Fmpz_poly from that string."
     
    115115            sage: f = Fmpz_poly([0,1]); f^7
    116116            8  0 0 0 0 0 0 0 1
    117117        """
    118         cdef char* ss = fmpz_poly_to_string(self.poly)
     118        cdef char* ss = fmpz_poly_get_str(self.poly)
    119119        cdef object s = ss
    120120        sage_free(ss)
    121121        return s
     
    246246        if not PY_TYPE_CHECK(self, Fmpz_poly):
    247247            raise TypeError
    248248        cdef Fmpz_poly res = <Fmpz_poly>PY_NEW(Fmpz_poly)
    249         fmpz_poly_power(res.poly, (<Fmpz_poly>self).poly, nn)
     249        fmpz_poly_pow(res.poly, (<Fmpz_poly>self).poly, nn)
    250250        return res
    251251       
    252252    def pow_truncate(self, exp, n):
     
    267267            raise ValueError, "Exponent must be at least 0"
    268268        cdef long exp_c = exp, nn = n
    269269        cdef Fmpz_poly res = <Fmpz_poly>PY_NEW(Fmpz_poly)
    270         fmpz_poly_power_trunc_n(res.poly, (<Fmpz_poly>self).poly, exp_c, nn)
     270        fmpz_poly_pow_trunc(res.poly, (<Fmpz_poly>self).poly, exp_c, nn)
    271271        return res
    272272       
    273273    def __floordiv__(left, right):
     
    329329        """
    330330        cdef Fmpz_poly res = <Fmpz_poly>PY_NEW(Fmpz_poly)
    331331
    332         fmpz_poly_left_shift(res.poly, self.poly, n)
     332        fmpz_poly_shift_left(res.poly, self.poly, n)
    333333
    334334        return res
    335335
     
    345345        """
    346346        cdef Fmpz_poly res = <Fmpz_poly>PY_NEW(Fmpz_poly)
    347347
    348         fmpz_poly_right_shift(res.poly, self.poly, n)
     348        fmpz_poly_shift_right(res.poly, self.poly, n)
    349349
    350350        return res
    351351
  • deleted file sage/libs/flint/long_extras.pxd

    diff --git a/sage/libs/flint/long_extras.pxd b/sage/libs/flint/long_extras.pxd
    deleted file mode 100644
    + -  
    1 include "../../ext/stdsage.pxi"
    2 include "../../ext/cdefs.pxi"
    3 
    4 from sage.libs.flint.flint cimport *
    5 
    6 
    7 cdef extern from "FLINT/long_extras.h":
    8 
    9     ctypedef struct factor_t:
    10         int num
    11         unsigned long p[15]
    12         unsigned long exp[15]
    13 
    14     cdef unsigned long z_randint(unsigned long limit)
    15 
    16     cdef unsigned long z_randbits(unsigned long bits)
    17 
    18     cdef unsigned long z_randprime(unsigned long bits, int proved)
    19 
    20     cdef double z_precompute_inverse(unsigned long n)
    21 
    22     cdef double z_precompute_inverse2(unsigned long n)
    23 
    24     cdef double z_ll_precompute_inverse(unsigned long n)
    25 
    26     cdef unsigned long z_addmod(unsigned long a, unsigned long b, unsigned long p)
    27 
    28     cdef unsigned long z_submod(unsigned long a, unsigned long b, unsigned long p)
    29 
    30     cdef unsigned long z_negmod(unsigned long a, unsigned long p)
    31 
    32     cdef unsigned long z_mod_precomp(unsigned long a, unsigned long n, double ninv)
    33 
    34     cdef unsigned long z_div2_precomp(unsigned long a, unsigned long n, double ninv)
    35 
    36     cdef unsigned long z_mod2_precomp(unsigned long a, unsigned long n, double ninv)
    37 
    38     cdef unsigned long z_ll_mod_precomp(unsigned long a_hi, unsigned long a_lo, unsigned long n, double ninv)
    39 
    40     cdef unsigned long z_mulmod_precomp(unsigned long a, unsigned long b, unsigned long n, double ninv)
    41 
    42     cdef unsigned long z_mulmod2_precomp(unsigned long a, unsigned long b, unsigned long n, double ninv)
    43 
    44     cdef unsigned long z_powmod(unsigned long a, long exp, unsigned long n)
    45 
    46     cdef unsigned long z_powmod2(unsigned long a, long exp, unsigned long n)
    47 
    48     cdef unsigned long z_powmod_precomp(unsigned long a, long exp, unsigned long n, double ninv)
    49 
    50     cdef unsigned long z_powmod2_precomp(unsigned long a, long exp, unsigned long n, double ninv)
    51 
    52     cdef int z_legendre_precomp(unsigned long a, unsigned long p, double pinv)
    53 
    54     cdef int z_jacobi(long x, unsigned long y)
    55 
    56     cdef int z_ispseudoprime_fermat(unsigned long n, unsigned long i)
    57 
    58     cdef int z_isprime(unsigned long n)
    59 
    60     cdef int z_isprime_precomp(unsigned long n, double ninv)
    61 
    62     cdef unsigned long z_nextprime(unsigned long n, int proved)
    63 
    64     cdef int z_isprime_pocklington(unsigned long n, unsigned long iterations)
    65 
    66     cdef int z_ispseudoprime_lucas_ab(unsigned long n, int a, int b)
    67 
    68     cdef int z_ispseudoprime_lucas(unsigned long n)
    69 
    70     cdef unsigned long z_pow(unsigned long a, unsigned long exp)
    71 
    72     cdef unsigned long z_sqrtmod(unsigned long a, unsigned long p)
    73 
    74     cdef unsigned long z_cuberootmod(unsigned long * cuberoot1, unsigned long a, unsigned long p)
    75 
    76     cdef unsigned long z_invert(unsigned long a, unsigned long p)
    77    
    78     cdef long z_gcd_invert(long* a, long x, long y)
    79 
    80     cdef long z_extgcd(long* a, long* b, long x, long y)
    81 
    82     cdef unsigned long z_gcd(long x, long y)
    83 
    84     cdef unsigned long z_intsqrt(unsigned long r)
    85 
    86     cdef int z_issquare(unsigned long x)
    87 
    88     cdef unsigned long z_CRT(unsigned long x1, unsigned long n1, unsigned long x2, unsigned long n2)
    89 
    90     cdef int z_issquarefree(unsigned long n, int proved)
    91 
    92     cdef int z_remove_precomp(unsigned long * n, unsigned long p, double pinv)
    93 
    94     cdef int z_remove(unsigned long * n, unsigned long p)
    95 
    96     cdef unsigned long z_factor_SQUFOF(unsigned long n)
    97 
    98     cdef unsigned long z_primitive_root(unsigned long p)
    99 
    100     cdef unsigned long z_primitive_root_precomp(unsigned long p, double p_inv)
    101 
    102     cdef void z_factor(factor_t * factors, unsigned long n, int proved)
    103 
    104 
  • sage/libs/flint/ntl_interface.pxd

    diff --git a/sage/libs/flint/ntl_interface.pxd b/sage/libs/flint/ntl_interface.pxd
    a b  
    66from sage.libs.ntl.ntl_ZZX_decl cimport ZZX_c
    77
    88cdef extern from "FLINT/NTL-interface.h":
    9     unsigned long ZZ_limbs(ZZ_c z)
    109
    11     void fmpz_poly_to_ZZX(ZZX_c output, fmpz_poly_t poly)
    12     void ZZX_to_fmpz_poly(fmpz_poly_t output, ZZX_c poly)
     10    void fmpz_poly_get_ZZX(ZZX_c output, fmpz_poly_t poly)
     11    void fmpz_poly_set_ZZX(fmpz_poly_t output, ZZX_c poly)
    1312
    14     void fmpz_to_mpz(mpz_t res, fmpz_t f)
    15     void ZZ_to_fmpz(fmpz_t output, ZZ_c z)
    16    
     13    void fmpz_get_ZZ(ZZ_c output, fmpz_t f)
     14    void fmpz_set_ZZ(fmpz_t output, ZZ_c z)
     15
  • new file sage/libs/flint/ulong_extras.pxd

    diff --git a/sage/libs/flint/ulong_extras.pxd b/sage/libs/flint/ulong_extras.pxd
    new file mode 100644
    - +  
     1include "../../ext/stdsage.pxi"
     2include "../../ext/cdefs.pxi"
     3
     4from sage.libs.flint.flint cimport *
     5
     6
     7cdef extern from "FLINT/ulong_extras.h":
     8
     9    ctypedef struct n_factor_t:
     10        int num
     11        unsigned long p[15]
     12        unsigned long exp[15]
     13
     14    cdef int n_jacobi(long x, unsigned long y)
     15
     16    cdef int n_is_prime(unsigned long n)
     17
     18    cdef unsigned long n_gcd(long x, long y)
     19
     20    cdef void n_factor(n_factor_t * factors, unsigned long n, int proved)
  • sage/libs/flint/zmod_poly.pxd

    diff --git a/sage/libs/flint/zmod_poly.pxd b/sage/libs/flint/zmod_poly.pxd
    a b  
    55
    66from flint import *
    77
    8 cdef extern from "FLINT/zmod_poly.h":
    9     ctypedef struct zmod_poly_struct:
    10         unsigned long *coeffs
    11         unsigned long alloc
    12         unsigned long length
    13         unsigned long p
    14         double p_inv
     8cdef extern from "FLINT/nmod_poly.h":
    159
    16     ctypedef zmod_poly_struct* zmod_poly_t
     10    ctypedef unsigned long mp_bitcnt_t
     11    ctypedef void * mp_srcptr
    1712
    18     ctypedef struct zmod_poly_factor_struct:
    19         zmod_poly_t *factors
    20         unsigned long *exponents
    21         unsigned long alloc
    22         unsigned long num_factors
     13    ctypedef struct nmod_t:
     14        mp_limb_t n
     15        mp_limb_t ninv
     16        mp_bitcnt_t norm
    2317
    24     ctypedef zmod_poly_factor_struct* zmod_poly_factor_t
     18    ctypedef struct nmod_poly_struct:
     19        mp_limb_t *coeffs
     20        long alloc
     21        long length
     22        nmod_t mod
    2523
    26     cdef void zmod_poly_init(zmod_poly_t poly, unsigned long p)
    27     cdef void zmod_poly_init_precomp(zmod_poly_t poly, unsigned long p, double p_inv)
    28     cdef void zmod_poly_init2(zmod_poly_t poly, unsigned long p, unsigned long alloc)
    29     cdef void zmod_poly_init2_precomp(zmod_poly_t poly, unsigned long p, double p_inv, unsigned long alloc)
    30     cdef void zmod_poly_clear(zmod_poly_t poly)
     24    ctypedef nmod_poly_struct* nmod_poly_t
    3125
    32     cdef void zmod_poly_realloc(zmod_poly_t poly, unsigned long alloc)
    33     # _bits_ only applies to newly allocated coefficients, not existing ones...
     26    ctypedef struct nmod_poly_factor_struct:
     27        nmod_poly_struct *p
     28        long *exp
     29        long num
     30        long alloc
    3431
    35     # this non-inlined version REQUIRES that alloc > poly->alloc
    36     void __zmod_poly_fit_length(zmod_poly_t poly, unsigned long alloc)
     32    ctypedef nmod_poly_factor_struct* nmod_poly_factor_t
    3733
    38     # this is arranged so that the initial comparison (very frequent) is inlined,
    39     # but the actual allocation (infrequent) is not
    40     cdef void zmod_poly_fit_length(zmod_poly_t poly, unsigned long alloc)
     34    # Memory management
    4135
    42     # ------------------------------------------------------
    43     # Setting/retrieving coefficients
     36    cdef void nmod_poly_init(nmod_poly_t poly, mp_limb_t n)
     37    cdef void nmod_poly_init_preinv(nmod_poly_t poly, mp_limb_t n, mp_limb_t ninv)
     38    cdef void nmod_poly_init2(nmod_poly_t poly, mp_limb_t n, long alloc)
     39    cdef void nmod_poly_init2_preinv(nmod_poly_t poly, mp_limb_t n, mp_limb_t ninv, long alloc)
     40    cdef void nmod_poly_realloc(nmod_poly_t poly, long alloc)
     41    cdef void nmod_poly_clear(nmod_poly_t poly)
     42    cdef void nmod_poly_fit_length(nmod_poly_t poly, long alloc)
     43    cdef void _nmod_poly_normalise(nmod_poly_t poly)
    4444
    45     cdef unsigned long zmod_poly_get_coeff_ui(zmod_poly_t poly, unsigned long n)
     45    # Polynomial parameters
    4646
    47     cdef unsigned long _zmod_poly_get_coeff_ui(zmod_poly_t poly, unsigned long n)
     47    cdef long nmod_poly_length(nmod_poly_t poly)
     48    cdef long nmod_poly_degree(nmod_poly_t poly)
     49    cdef mp_limb_t nmod_poly_modulus(nmod_poly_t poly)
     50    cdef mp_bitcnt_t nmod_poly_max_bits(nmod_poly_t poly)
    4851
    49     cdef void zmod_poly_set_coeff_ui(zmod_poly_t poly, unsigned long n, unsigned long c)
     52    # Assignment and basic manipulation
    5053
    51     cdef void _zmod_poly_set_coeff_ui(zmod_poly_t poly, unsigned long n, unsigned long c)
     54    cdef void nmod_poly_set(nmod_poly_t a, nmod_poly_t b)
     55    cdef void nmod_poly_swap(nmod_poly_t poly1, nmod_poly_t poly2)
     56    cdef void nmod_poly_zero(nmod_poly_t res)
     57    cdef void nmod_poly_one(nmod_poly_t res)
     58    cdef void nmod_poly_truncate(nmod_poly_t poly, long len)
     59    cdef void nmod_poly_reverse(nmod_poly_t output, nmod_poly_t input, long m)
    5260
    53     # ------------------------------------------------------
    54     # String conversions and I/O
     61    # Comparison
    5562
    56     cdef int zmod_poly_from_string(zmod_poly_t poly, char* s)
    57     cdef char* zmod_poly_to_string(zmod_poly_t poly)
    58     cdef void zmod_poly_print(zmod_poly_t poly)
    59     cdef void zmod_poly_fprint(zmod_poly_t poly, FILE* f)
    60     cdef int zmod_poly_read(zmod_poly_t poly)
    61     cdef int zmod_poly_fread(zmod_poly_t poly, FILE* f)
     63    cdef int nmod_poly_equal(nmod_poly_t a, nmod_poly_t b)
     64    cdef int nmod_poly_is_zero(nmod_poly_t poly)
     65    cdef int nmod_poly_is_one(nmod_poly_t poly)
    6266
    63     # ------------------------------------------------------
    64     # Length and degree
     67    # Getting and setting coefficients
    6568
    66     cdef void __zmod_poly_normalise(zmod_poly_t poly)
    67     cdef int __zmod_poly_normalised(zmod_poly_t poly)
    68     cdef void zmod_poly_truncate(zmod_poly_t poly, unsigned long length)
     69    cdef unsigned long nmod_poly_get_coeff_ui(nmod_poly_t poly, unsigned long j)
     70    cdef void nmod_poly_set_coeff_ui(nmod_poly_t poly, unsigned long j, unsigned long c)
    6971
    70     cdef unsigned long zmod_poly_length(zmod_poly_t poly)
     72    # Input and output
    7173
    72     cdef long zmod_poly_degree(zmod_poly_t poly)
     74    cdef char * nmod_poly_get_str(nmod_poly_t poly)
     75    cdef int nmod_poly_set_str(char * s, nmod_poly_t poly)
     76    cdef int nmod_poly_print(nmod_poly_t a)
     77    cdef int nmod_poly_fread(FILE * f, nmod_poly_t poly)
     78    cdef int nmod_poly_fprint(FILE * f, nmod_poly_t poly)
     79    cdef int nmod_poly_read(nmod_poly_t poly)
    7380
    74     cdef unsigned long zmod_poly_modulus(zmod_poly_t poly)
     81    # Shifting
    7582
    76     cdef double zmod_poly_precomputed_inverse(zmod_poly_t poly)
     83    cdef void nmod_poly_shift_left(nmod_poly_t res, nmod_poly_t poly, long k)
     84    cdef void nmod_poly_shift_right(nmod_poly_t res, nmod_poly_t poly, long k)
    7785
    78     # ------------------------------------------------------
    79     # Assignment
     86    # Addition and subtraction
    8087
    81     cdef void _zmod_poly_set(zmod_poly_t res, zmod_poly_t poly)
    82     cdef void zmod_poly_set(zmod_poly_t res, zmod_poly_t poly)
     88    cdef void nmod_poly_add(nmod_poly_t res, nmod_poly_t poly1, nmod_poly_t poly2)
     89    cdef void nmod_poly_sub(nmod_poly_t res, nmod_poly_t poly1, nmod_poly_t poly2)
     90    cdef void nmod_poly_neg(nmod_poly_t res, nmod_poly_t poly1)
    8391
    84     cdef void zmod_poly_zero(zmod_poly_t poly)
     92    # Scalar multiplication and division
    8593
    86     cdef void zmod_poly_swap(zmod_poly_t poly1, zmod_poly_t poly2)
     94    cdef void nmod_poly_scalar_mul_nmod(nmod_poly_t res, nmod_poly_t poly1, mp_limb_t c)
     95    cdef void nmod_poly_make_monic(nmod_poly_t output, nmod_poly_t input)
    8796
    88     #
    89     # Subpolynomials
    90     #
     97    # Multiplication
    9198
    92     cdef void _zmod_poly_attach(zmod_poly_t output, zmod_poly_t input)
     99    cdef void nmod_poly_mul(nmod_poly_t res, nmod_poly_t poly1, nmod_poly_t poly2)
     100    cdef void nmod_poly_mullow(nmod_poly_t res, nmod_poly_t poly1, nmod_poly_t poly2, long trunc)
     101    cdef void nmod_poly_mulhigh(nmod_poly_t res, nmod_poly_t poly1, nmod_poly_t poly2, long n)
     102    cdef void nmod_poly_mulmod(nmod_poly_t res, nmod_poly_t poly1, nmod_poly_t poly2, nmod_poly_t f)
    93103
    94     cdef void zmod_poly_attach(zmod_poly_t output, zmod_poly_t input)
     104    # Powering
    95105
    96     #
    97     # Attach input shifted right by n to output
    98     #
     106    cdef void nmod_poly_pow(nmod_poly_t res, nmod_poly_t poly, unsigned long e)
     107    cdef void nmod_poly_pow_trunc(nmod_poly_t res, nmod_poly_t poly, unsigned long e, long trunc)
    99108
    100     cdef void _zmod_poly_attach_shift(zmod_poly_t output, zmod_poly_t input, unsigned long n)
     109    # Division
    101110
    102     cdef void zmod_poly_attach_shift(zmod_poly_t output, zmod_poly_t input, unsigned long n)
     111    cdef void nmod_poly_divrem(nmod_poly_t Q, nmod_poly_t R, nmod_poly_t A, nmod_poly_t B)
     112    cdef void nmod_poly_div(nmod_poly_t Q, nmod_poly_t A, nmod_poly_t B)
     113    cdef void nmod_poly_inv_series(nmod_poly_t Qinv, nmod_poly_t Q, long n)
     114    cdef void nmod_poly_div_series(nmod_poly_t Q, nmod_poly_t A, nmod_poly_t B, long n)
    103115
    104     #
    105     # Attach input to first n coefficients of input
    106     #
     116    # Derivative
    107117
    108     cdef void _zmod_poly_attach_truncate(zmod_poly_t output,  zmod_poly_t input, unsigned long n)
     118    cdef void nmod_poly_derivative(nmod_poly_t x_prime, nmod_poly_t x)
     119    cdef void nmod_poly_integral(nmod_poly_t x_int, nmod_poly_t x)
    109120
    110     cdef void zmod_poly_attach_truncate(zmod_poly_t output, zmod_poly_t input, unsigned long n)
     121    # Evaluation
    111122
    112     #
    113     # Comparison functions
    114     #
     123    cdef mp_limb_t nmod_poly_evaluate_nmod(nmod_poly_t poly, mp_limb_t c)
     124    cdef void nmod_poly_evaluate_nmod_vec(mp_ptr ys, nmod_poly_t poly, mp_srcptr xs, long n)
    115125
    116     cdef int zmod_poly_equal(zmod_poly_t poly1, zmod_poly_t poly2)
     126    # Interpolation
    117127
    118     cdef int zmod_poly_is_one(zmod_poly_t poly1)
     128    cdef void nmod_poly_interpolate_nmod_vec(nmod_poly_t poly, mp_srcptr xs, mp_srcptr ys, long n)
    119129
    120     #
    121     # Reversal
    122     #
     130    # Composition
    123131
    124     cdef void _zmod_poly_reverse(zmod_poly_t output, zmod_poly_t input, unsigned long length)
    125     cdef void zmod_poly_reverse(zmod_poly_t output, zmod_poly_t input, unsigned long length)
     132    cdef void nmod_poly_compose(nmod_poly_t res, nmod_poly_t poly1, nmod_poly_t poly2)
    126133
    127     #
    128     # Monic polys
    129     #
     134    # Power series composition and reversion
    130135
    131     cdef void zmod_poly_make_monic(zmod_poly_t output, zmod_poly_t pol)
     136    cdef void nmod_poly_compose_series(nmod_poly_t res, nmod_poly_t poly1, nmod_poly_t poly2, long n)
     137    cdef void nmod_poly_reverse_series(nmod_poly_t Qinv, nmod_poly_t Q, long n)
    132138
    133     #
    134     # Addition and subtraction
    135     #
     139    # GCD
    136140
    137     cdef void zmod_poly_add(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2)
    138     cdef void zmod_poly_add_without_mod(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2)
    139     cdef void zmod_poly_sub(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2)
    140     cdef void _zmod_poly_sub(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2)
    141     cdef void zmod_poly_neg(zmod_poly_t res, zmod_poly_t poly)
     141    cdef void nmod_poly_gcd(nmod_poly_t G, nmod_poly_t A, nmod_poly_t B)
     142    cdef void nmod_poly_xgcd(nmod_poly_t G, nmod_poly_t S, nmod_poly_t T, nmod_poly_t A, nmod_poly_t B)
     143    mp_limb_t nmod_poly_resultant(nmod_poly_t f, nmod_poly_t g)
    142144
    143     #
    144     # Shifting functions
    145     #
     145    # Square roots
    146146
    147     cdef void zmod_poly_left_shift(zmod_poly_t res, zmod_poly_t poly, unsigned long k)
    148     cdef void zmod_poly_right_shift(zmod_poly_t res, zmod_poly_t poly, unsigned long k)
     147    cdef void nmod_poly_invsqrt_series(nmod_poly_t g, nmod_poly_t h, long n)
     148    cdef void nmod_poly_sqrt_series(nmod_poly_t g, nmod_poly_t h, long n)
     149    int nmod_poly_sqrt(nmod_poly_t b, nmod_poly_t a)
    149150
    150     #
    151     # Polynomial multiplication
    152     #
    153     # All multiplication functions require that the modulus be no more than FLINT_BITS-1 bits
    154     #
     151    # Transcendental functions
    155152
    156     cdef void zmod_poly_mul(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2)
    157     cdef void zmod_poly_sqr(zmod_poly_t res, zmod_poly_t poly)
     153    cdef void nmod_poly_atan_series(nmod_poly_t g, nmod_poly_t h, long n)
     154    cdef void nmod_poly_tan_series(nmod_poly_t g, nmod_poly_t h, long n)
     155    cdef void nmod_poly_asin_series(nmod_poly_t g, nmod_poly_t h, long n)
     156    cdef void nmod_poly_sin_series(nmod_poly_t g, nmod_poly_t h, long n)
     157    cdef void nmod_poly_cos_series(nmod_poly_t g, nmod_poly_t h, long n)
     158    cdef void nmod_poly_asinh_series(nmod_poly_t g, nmod_poly_t h, long n)
     159    cdef void nmod_poly_atanh_series(nmod_poly_t g, nmod_poly_t h, long n)
     160    cdef void nmod_poly_sinh_series(nmod_poly_t g, nmod_poly_t h, long n)
     161    cdef void nmod_poly_cosh_series(nmod_poly_t g, nmod_poly_t h, long n)
     162    cdef void nmod_poly_tanh_series(nmod_poly_t g, nmod_poly_t h, long n)
     163    cdef void nmod_poly_log_series(nmod_poly_t res, nmod_poly_t f, long n)
     164    cdef void nmod_poly_exp_series(nmod_poly_t f, nmod_poly_t h, long)
    158165
    159     # Requires that poly1 bits + poly2 bits + log_length is not greater than 2*FLINT_BITS
     166    # Inflation and deflation
    160167
    161     cdef void zmod_poly_mul_KS(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, unsigned long bits_input)
    162     cdef void _zmod_poly_mul_KS(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, unsigned long bits_input)
     168    cdef unsigned long nmod_poly_deflation(nmod_poly_t input)
     169    cdef void nmod_poly_deflate(nmod_poly_t result, nmod_poly_t input, unsigned long deflation)
     170    cdef void nmod_poly_inflate(nmod_poly_t result, nmod_poly_t input, unsigned long inflation)
    163171
    164     cdef void zmod_poly_mul_KS_trunc(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, unsigned long bits_input, unsigned long trunc)
    165     cdef void _zmod_poly_mul_KS_trunc(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, unsigned long bits_input, unsigned long trunc)
     172    # Factoring
    166173
    167     cdef void _zmod_poly_mul_classical(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2)
    168     cdef void __zmod_poly_mul_classical_mod_last(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, unsigned long bits)
    169     cdef void __zmod_poly_mul_classical_mod_throughout(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, unsigned long bits)
    170     cdef void zmod_poly_mul_classical(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2)
    171     cdef void _zmod_poly_sqr_classical(zmod_poly_t res, zmod_poly_t poly)
    172     cdef void zmod_poly_sqr_classical(zmod_poly_t res, zmod_poly_t poly)
    173 
    174     cdef void _zmod_poly_mul_classical_trunc(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, unsigned long trunc)
    175     cdef void __zmod_poly_mul_classical_trunc_mod_last(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, unsigned long bits, unsigned long trunc)
    176     cdef void __zmod_poly_mul_classical_trunc_mod_throughout(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, unsigned long bits, unsigned long trunc)
    177     cdef void zmod_poly_mul_classical_trunc(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, unsigned long trunc)
    178 
    179     cdef void _zmod_poly_mul_classical_trunc_left(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, unsigned long trunc)
    180     cdef void __zmod_poly_mul_classical_trunc_left_mod_last(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, unsigned long bits, unsigned long trunc)
    181     cdef void __zmod_poly_mul_classical_trunc_left_mod_throughout(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, unsigned long bits, unsigned long trunc)
    182     cdef void zmod_poly_mul_classical_trunc_left(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, unsigned long trunc)
    183 
    184     cdef void zmod_poly_mul_trunc_n(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, unsigned long trunc)
    185     cdef void zmod_poly_mul_trunc_left_n(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, unsigned long trunc)
    186 
    187     #
    188     # Bit packing functions
    189     #
    190 
    191     cdef unsigned long zmod_poly_bits(zmod_poly_t poly)
    192     cdef void _zmod_poly_bit_pack_mpn(mp_limb_t * res, zmod_poly_t poly, unsigned long bits, unsigned long length)
    193     cdef void _zmod_poly_bit_unpack_mpn(zmod_poly_t poly, mp_limb_t *mpn, unsigned long length, unsigned long bits)
    194 
    195     cdef void print_binary(unsigned long n, unsigned long len)
    196     cdef void print_binary2(unsigned long n, unsigned long len, unsigned long space_bit)
    197 
    198     #
    199     # Scalar multiplication
    200     #
    201 
    202     cdef void _zmod_poly_scalar_mul(zmod_poly_t res, zmod_poly_t poly, unsigned long scalar)
    203     cdef void zmod_poly_scalar_mul(zmod_poly_t res, zmod_poly_t poly, unsigned long scalar)
    204     cdef void __zmod_poly_scalar_mul_without_mod(zmod_poly_t res, zmod_poly_t poly, unsigned long scalar)
    205 
    206     #
    207     # Division
    208     #
    209 
    210     cdef void zmod_poly_divrem_classical(zmod_poly_t Q, zmod_poly_t R, zmod_poly_t A, zmod_poly_t B)
    211     cdef void __zmod_poly_divrem_classical_mod_last(zmod_poly_t Q, zmod_poly_t R, zmod_poly_t A, zmod_poly_t B)
    212     cdef void zmod_poly_div_classical(zmod_poly_t Q, zmod_poly_t A, zmod_poly_t B)
    213     cdef void __zmod_poly_div_classical_mod_last(zmod_poly_t Q, zmod_poly_t A, zmod_poly_t B)
    214     cdef void zmod_poly_div_divconquer_recursive(zmod_poly_t Q, zmod_poly_t BQ, zmod_poly_t A, zmod_poly_t B)
    215     cdef void zmod_poly_divrem_divconquer(zmod_poly_t Q, zmod_poly_t R, zmod_poly_t A, zmod_poly_t B)
    216     cdef void zmod_poly_div_divconquer(zmod_poly_t Q, zmod_poly_t A, zmod_poly_t B)
    217 
    218     #
    219     # Newton Inversion
    220     #
    221 
    222     cdef void zmod_poly_newton_invert_basecase(zmod_poly_t Q_inv, zmod_poly_t Q, unsigned long n)
    223     cdef void zmod_poly_newton_invert(zmod_poly_t Q_inv, zmod_poly_t Q, unsigned long n)
    224 
    225     #
    226     # Newton Division
    227     #
    228 
    229     cdef void zmod_poly_div_series(zmod_poly_t Q, zmod_poly_t A, zmod_poly_t B, unsigned long n)
    230     cdef void zmod_poly_div_newton(zmod_poly_t Q, zmod_poly_t A, zmod_poly_t B)
    231     cdef void zmod_poly_divrem_newton(zmod_poly_t Q, zmod_poly_t R, zmod_poly_t A, zmod_poly_t B)
    232 
    233     cdef void zmod_poly_divrem(zmod_poly_t Q, zmod_poly_t R, zmod_poly_t A, zmod_poly_t B)
    234 
    235     cdef void zmod_poly_div(zmod_poly_t Q, zmod_poly_t A, zmod_poly_t B)
    236 
    237     #
    238     # Resultant
    239     #
    240 
    241     cdef unsigned long zmod_poly_resultant_euclidean(zmod_poly_t a, zmod_poly_t b)
    242 
    243     cdef unsigned long zmod_poly_resultant(zmod_poly_t a, zmod_poly_t b)
    244 
    245     #
    246     # GCD
    247     #
    248 
    249     cdef void zmod_poly_gcd(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2)
    250     cdef int zmod_poly_gcd_invert(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2)
    251     cdef void zmod_poly_xgcd(zmod_poly_t res, zmod_poly_t s, zmod_poly_t t, zmod_poly_t poly1, zmod_poly_t poly2)
    252 
    253 
    254 
    255     # Composition / evaluation
    256 
    257     cdef unsigned long zmod_poly_evaluate(zmod_poly_t, unsigned long)
    258     cdef void zmod_poly_compose_horner(zmod_poly_t, zmod_poly_t, zmod_poly_t)
    259 
    260     # Factorization
    261 
    262     cdef bint zmod_poly_isirreducible(zmod_poly_t p)
    263 
    264     ctypedef struct zmod_poly_factors_struct:
    265         unsigned long num_factors
    266         unsigned long* exponents
    267         zmod_poly_t* factors
    268 
    269     ctypedef zmod_poly_factors_struct* zmod_poly_factor_t
    270 
    271     cdef void zmod_poly_factor_init(zmod_poly_factor_t)
    272     cdef void zmod_poly_factor_clear(zmod_poly_factor_t)
    273     cdef unsigned long zmod_poly_factor(zmod_poly_factor_t, zmod_poly_t)
    274     cdef void zmod_poly_factor_square_free(zmod_poly_factor_t, zmod_poly_t)
    275     cdef void zmod_poly_factor_berlekamp(zmod_poly_factor_t factors, zmod_poly_t f)
    276 
    277     cdef void zmod_poly_factor_add(zmod_poly_factor_t fac, zmod_poly_t poly)
    278     cdef void zmod_poly_factor_concat(zmod_poly_factor_t res, zmod_poly_factor_t fac)
    279     cdef void zmod_poly_factor_print(zmod_poly_factor_t fac)
    280     cdef void zmod_poly_factor_pow(zmod_poly_factor_t fac, unsigned long exp)
    281 
    282     #
    283     # Differentiation
    284     #
    285  
    286     cdef void zmod_poly_derivative(zmod_poly_t res, zmod_poly_t poly)
    287    
    288     #
    289     # Arithmetic modulo a polynomial
    290     #
    291    
    292     cdef void zmod_poly_mulmod(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, zmod_poly_t f)
    293     cdef void zmod_poly_powmod(zmod_poly_t res,zmod_poly_t pol, long exp, zmod_poly_t f)
     174    cdef void nmod_poly_factor_clear(nmod_poly_factor_t fac)
     175    cdef void nmod_poly_factor_init(nmod_poly_factor_t fac)
     176    cdef void nmod_poly_factor_insert(nmod_poly_factor_t fac, nmod_poly_t poly, unsigned long exp)
     177    cdef void nmod_poly_factor_print(nmod_poly_factor_t fac)
     178    cdef void nmod_poly_factor_concat(nmod_poly_factor_t res, nmod_poly_factor_t fac)
     179    cdef void nmod_poly_factor_pow(nmod_poly_factor_t fac, unsigned long exp)
     180    cdef unsigned long nmod_poly_remove(nmod_poly_t f, nmod_poly_t p)
     181    cdef int nmod_poly_is_irreducible(nmod_poly_t f)
     182    cdef int nmod_poly_is_squarefree(nmod_poly_t f)
     183    cdef void nmod_poly_factor_cantor_zassenhaus(nmod_poly_factor_t res, nmod_poly_t f)
     184    cdef void nmod_poly_factor_berlekamp(nmod_poly_factor_t factors, nmod_poly_t f)
     185    cdef void nmod_poly_factor_squarefree(nmod_poly_factor_t res, nmod_poly_t f)
     186    cdef mp_limb_t nmod_poly_factor_with_berlekamp(nmod_poly_factor_t result, nmod_poly_t input)
     187    cdef mp_limb_t nmod_poly_factor_with_cantor_zassenhaus(nmod_poly_factor_t result, nmod_poly_t input)
     188    cdef mp_limb_t nmod_poly_factor(nmod_poly_factor_t result, nmod_poly_t input)
  • sage/libs/flint/zmod_poly_linkage.pxi

    diff --git a/sage/libs/flint/zmod_poly_linkage.pxi b/sage/libs/flint/zmod_poly_linkage.pxi
    a b  
    1717#                  http://www.gnu.org/licenses/
    1818#*****************************************************************************
    1919
    20 from sage.libs.flint.zmod_poly cimport *, zmod_poly_t
    21 from sage.libs.flint.long_extras cimport *
     20from sage.libs.flint.zmod_poly cimport *, nmod_poly_t
     21from sage.libs.flint.ulong_extras cimport *
    2222
    2323include "../../ext/stdsage.pxi"
    2424
    2525cdef inline celement *celement_new(unsigned long n):
    26     cdef celement *g = <celement *>sage_malloc(sizeof(zmod_poly_t))
    27     zmod_poly_init(g, n)
     26    cdef celement *g = <celement *>sage_malloc(sizeof(nmod_poly_t))
     27    nmod_poly_init(g, n)
    2828    return g
    2929
    30 cdef inline int celement_delete(zmod_poly_t e, unsigned long n):
    31     zmod_poly_clear(e)
     30cdef inline int celement_delete(nmod_poly_t e, unsigned long n):
     31    nmod_poly_clear(e)
    3232    sage_free(e)
    3333
    34 cdef inline int celement_construct(zmod_poly_t e, unsigned long n):
     34cdef inline int celement_construct(nmod_poly_t e, unsigned long n):
    3535    """
    3636    EXAMPLE:
    3737        sage: P.<x> = GF(32003)[]
    3838
    3939        sage: Q.<x> = GF(7)[]
    4040    """
    41     zmod_poly_init(e, n)
     41    nmod_poly_init(e, n)
    4242
    43 cdef inline int celement_destruct(zmod_poly_t e, unsigned long n):
     43cdef inline int celement_destruct(nmod_poly_t e, unsigned long n):
    4444    """
    4545    EXAMPLE:
    4646        sage: P.<x> = GF(32003)[]
     
    4949        sage: Q.<x> = GF(7)[]
    5050        sage: del x
    5151    """
    52     zmod_poly_clear(e)
     52    nmod_poly_clear(e)
    5353
    54 cdef inline int celement_gen(zmod_poly_t e, long i, unsigned long n) except -2:
     54cdef inline int celement_gen(nmod_poly_t e, long i, unsigned long n) except -2:
    5555    """
    5656    EXAMPLE:
    5757        sage: P.<x> = GF(32003)[]
    5858
    5959        sage: Q.<x> = GF(7)[]
    6060    """
    61     zmod_poly_zero(e)
    62     zmod_poly_set_coeff_ui(e, 1, 1)
     61    nmod_poly_zero(e)
     62    nmod_poly_set_coeff_ui(e, 1, 1)
    6363
    64 cdef object celement_repr(zmod_poly_t e, unsigned long n):
     64cdef object celement_repr(nmod_poly_t e, unsigned long n):
    6565    raise NotImplementedError
    6666
    67 cdef inline int celement_set(zmod_poly_t res, zmod_poly_t a, unsigned long n) except -2:
     67cdef inline int celement_set(nmod_poly_t res, nmod_poly_t a, unsigned long n) except -2:
    6868    """
    6969    EXAMPLE:
    7070        sage: P.<x> = GF(32003)[]
     
    8989        6*x
    9090    """
    9191    cdef unsigned long i
    92     if a.p <= n:
    93         zmod_poly_set(res, a)
     92    if a.mod.n <= n:
     93        nmod_poly_set(res, a)
    9494    else:
    95         zmod_poly_zero(res)
     95        nmod_poly_zero(res)
    9696        for i from 0 <= i < a.length:
    97             zmod_poly_set_coeff_ui(res, i, zmod_poly_get_coeff_ui(a, i) % n)
     97            nmod_poly_set_coeff_ui(res, i, nmod_poly_get_coeff_ui(a, i) % n)
    9898
    99 cdef inline int celement_set_si(zmod_poly_t res, long i, unsigned long n) except -2:
     99cdef inline int celement_set_si(nmod_poly_t res, long i, unsigned long n) except -2:
    100100    """
    101101    EXAMPLE:
    102102        sage: P.<x> = GF(32003)[]
     
    117117    """
    118118    while i < 0:
    119119        i += n
    120     zmod_poly_zero(res)
     120    nmod_poly_zero(res)
    121121    if i:
    122         zmod_poly_set_coeff_ui(res, 0, <unsigned long>i)
     122        nmod_poly_set_coeff_ui(res, 0, <unsigned long>i)
    123123
    124 cdef inline long celement_get_si(zmod_poly_t res, unsigned long n) except -2:
     124cdef inline long celement_get_si(nmod_poly_t res, unsigned long n) except -2:
    125125    raise NotImplementedError
    126126   
    127 cdef inline bint celement_is_zero(zmod_poly_t a, unsigned long n) except -2:
     127cdef inline bint celement_is_zero(nmod_poly_t a, unsigned long n) except -2:
    128128    """
    129129    EXAMPLE:
    130130        sage: P.<x> = GF(32003)[]
     
    139139        sage: Q(0).is_zero()
    140140        True
    141141    """
    142     # is_zero doesn't exist
    143     return zmod_poly_degree(a) == -1
     142    return nmod_poly_is_zero(a)
    144143
    145 cdef inline bint celement_is_one(zmod_poly_t a, unsigned long n) except -2:
     144cdef inline bint celement_is_one(nmod_poly_t a, unsigned long n) except -2:
    146145    """
    147146    EXAMPLE:
    148147        sage: P.<x> = GF(32003)[]
     
    158157        False
    159158    """
    160159
    161     return zmod_poly_is_one(a)
     160    return nmod_poly_is_one(a)
    162161
    163 cdef inline bint celement_equal(zmod_poly_t a, zmod_poly_t b, unsigned long n) except -2:
     162cdef inline bint celement_equal(nmod_poly_t a, nmod_poly_t b, unsigned long n) except -2:
    164163    """
    165164    EXAMPLE:
    166165        sage: P.<x> = GF(32003)[]
     
    179178        sage: (3*2)*x + 7 == 3*(2*x) + 1 + 1
    180179        False
    181180    """
    182     return zmod_poly_equal(a, b)
     181    return nmod_poly_equal(a, b)
    183182
    184 cdef inline int celement_cmp(zmod_poly_t l, zmod_poly_t r, unsigned long n) except -2:
     183cdef inline int celement_cmp(nmod_poly_t l, nmod_poly_t r, unsigned long n) except -2:
    185184    """
    186185    EXAMPLE:
    187186        sage: P.<x> = GF(32003)[]
     
    215214        sage: f < g
    216215        True
    217216    """
    218     cdef int deg_right = zmod_poly_degree(r)
    219     cdef int degdiff = deg_right - zmod_poly_degree(l)
     217    cdef int deg_right = nmod_poly_degree(r)
     218    cdef int degdiff = deg_right - nmod_poly_degree(l)
    220219    cdef int i
    221220    cdef unsigned long rcoeff, lcoeff
    222221    if degdiff > 0:
     
    224223    elif degdiff < 0:
    225224        return 1
    226225    else:
    227         if zmod_poly_equal(l, r):
     226        if nmod_poly_equal(l, r):
    228227            return 0
    229228        i = deg_right
    230         rcoeff = zmod_poly_get_coeff_ui(r, i)
    231         lcoeff = zmod_poly_get_coeff_ui(l, i)
     229        rcoeff = nmod_poly_get_coeff_ui(r, i)
     230        lcoeff = nmod_poly_get_coeff_ui(l, i)
    232231        while rcoeff == lcoeff and i > 0:
    233232            i -= 1
    234             rcoeff = zmod_poly_get_coeff_ui(r, i)
    235             lcoeff = zmod_poly_get_coeff_ui(l, i)
     233            rcoeff = nmod_poly_get_coeff_ui(r, i)
     234            lcoeff = nmod_poly_get_coeff_ui(l, i)
    236235        return cmp(lcoeff, rcoeff)
    237236
    238 cdef long celement_len(zmod_poly_t a, unsigned long n) except -2:
     237cdef long celement_len(nmod_poly_t a, unsigned long n) except -2:
    239238    """
    240239    EXAMPLE:
    241240        sage: P.<x> = GF(32003)[]
     
    254253        sage: P(0).degree()
    255254        -1
    256255    """
    257     return <long>zmod_poly_length(a)
     256    return <long>nmod_poly_length(a)
    258257
    259 cdef inline int celement_add(zmod_poly_t res, zmod_poly_t a, zmod_poly_t b, unsigned long n) except -2:
     258cdef inline int celement_add(nmod_poly_t res, nmod_poly_t a, nmod_poly_t b, unsigned long n) except -2:
    260259    """
    261260    EXAMPLE:
    262261        sage: P.<x> = GF(32003)[]
     
    267266        sage: x + 1
    268267        x + 1
    269268    """
    270     zmod_poly_add(res, a, b)
     269    nmod_poly_add(res, a, b)
    271270
    272 cdef inline int celement_sub(zmod_poly_t res, zmod_poly_t a, zmod_poly_t b, unsigned long n) except -2:
     271cdef inline int celement_sub(nmod_poly_t res, nmod_poly_t a, nmod_poly_t b, unsigned long n) except -2:
    273272    """
    274273    EXAMPLE:
    275274        sage: P.<x> = GF(32003)[]
     
    280279        sage: x - 1
    281280        x + 6
    282281    """
    283     zmod_poly_sub(res, a, b)
     282    nmod_poly_sub(res, a, b)
    284283
    285 cdef inline int celement_neg(zmod_poly_t res, zmod_poly_t a, unsigned long n) except -2:
     284cdef inline int celement_neg(nmod_poly_t res, nmod_poly_t a, unsigned long n) except -2:
    286285    """
    287286    EXAMPLE:
    288287        sage: P.<x> = GF(32003)[]
     
    293292        sage: -(x + 2)
    294293        6*x + 5       
    295294    """
    296     zmod_poly_neg(res, a)
     295    nmod_poly_neg(res, a)
    297296
    298 cdef inline int celement_mul_scalar(zmod_poly_t res, zmod_poly_t p,
     297cdef inline int celement_mul_scalar(nmod_poly_t res, nmod_poly_t p,
    299298        object c, unsigned long n) except -2:
    300299    """
    301300    TESTS::
     
    307306        sage: p*983
    308307        29561*x^2 + 18665*x + 17051
    309308    """   
    310     zmod_poly_scalar_mul(res, p, (<unsigned long>c)%n)
     309    nmod_poly_scalar_mul_nmod(res, p, (<unsigned long>c)%n)
    311310
    312 cdef inline int celement_mul(zmod_poly_t res, zmod_poly_t a, zmod_poly_t b, unsigned long n) except -2:
     311cdef inline int celement_mul(nmod_poly_t res, nmod_poly_t a, nmod_poly_t b, unsigned long n) except -2:
    313312    """
    314313    EXAMPLE:
    315314        sage: P.<x> = GF(32003)[]
     
    320319        sage: (x + 1) * (x + 2)
    321320        x^2 + 3*x + 2
    322321    """
    323     zmod_poly_mul(res, a, b)
     322    nmod_poly_mul(res, a, b)
    324323
    325 cdef inline int celement_div(zmod_poly_t res, zmod_poly_t a, zmod_poly_t b, unsigned long n) except -2:
     324cdef inline int celement_div(nmod_poly_t res, nmod_poly_t a, nmod_poly_t b, unsigned long n) except -2:
    326325    raise NotImplementedError
    327326
    328 cdef inline int celement_floordiv(zmod_poly_t res, zmod_poly_t a, zmod_poly_t b, unsigned long n) except -2:
     327cdef inline int celement_floordiv(nmod_poly_t res, nmod_poly_t a, nmod_poly_t b, unsigned long n) except -2:
    329328    """
    330329    EXAMPLE:
    331330        sage: P.<x> = GF(32003)[]
     
    348347        sage: (x^2 + 3*x + 2)//(x + 2)
    349348        x + 1
    350349    """
    351     zmod_poly_div(res, a, b)
     350    nmod_poly_div(res, a, b)
    352351
    353 cdef inline int celement_mod(zmod_poly_t res, zmod_poly_t a, zmod_poly_t b, unsigned long n) except -2:
     352cdef inline int celement_mod(nmod_poly_t res, nmod_poly_t a, nmod_poly_t b, unsigned long n) except -2:
    354353    """
    355354    EXAMPLE:
    356355        sage: P.<x> = GF(32003)[]
     
    381380        ...
    382381        ZeroDivisionError
    383382    """
    384     cdef zmod_poly_t q
     383    cdef nmod_poly_t q
    385384    cdef unsigned long leadcoeff, modulus
    386385   
    387     zmod_poly_init(q, n)
    388     leadcoeff = zmod_poly_get_coeff_ui(b, zmod_poly_degree(b))
    389     modulus = zmod_poly_modulus(b)
    390     if (leadcoeff > 1 and z_gcd(modulus,leadcoeff) != 1):
     386    nmod_poly_init(q, n)
     387    leadcoeff = nmod_poly_get_coeff_ui(b, nmod_poly_degree(b))
     388    modulus = nmod_poly_modulus(b)
     389    if (leadcoeff > 1 and n_gcd(modulus,leadcoeff) != 1):
    391390        raise ValueError("Leading coefficient of a must be invertible.")
    392391
    393     zmod_poly_divrem(q, res, a, b)
    394     zmod_poly_clear(q)
     392    nmod_poly_divrem(q, res, a, b)
     393    nmod_poly_clear(q)
    395394
    396 cdef inline int celement_quorem(zmod_poly_t q, zmod_poly_t r, zmod_poly_t a, zmod_poly_t b, unsigned long n) except -2:
     395cdef inline int celement_quorem(nmod_poly_t q, nmod_poly_t r, nmod_poly_t a, nmod_poly_t b, unsigned long n) except -2:
    397396    """
    398397    EXAMPLES:
    399398        sage: R.<x> = Integers(125)[]
     
    418417    """
    419418    cdef unsigned long leadcoeff, modulus
    420419
    421     leadcoeff = zmod_poly_get_coeff_ui(b, zmod_poly_degree(b))
    422     modulus = zmod_poly_modulus(b)
    423     if (leadcoeff > 1 and z_gcd(modulus,leadcoeff) != 1):
     420    leadcoeff = nmod_poly_get_coeff_ui(b, nmod_poly_degree(b))
     421    modulus = nmod_poly_modulus(b)
     422    if (leadcoeff > 1 and n_gcd(modulus,leadcoeff) != 1):
    424423        raise ValueError("Leading coefficient of a must be invertible.")
    425424
    426     zmod_poly_divrem(q, r, a, b)
     425    nmod_poly_divrem(q, r, a, b)
    427426
    428 cdef inline int celement_inv(zmod_poly_t res, zmod_poly_t a, unsigned long n) except -2:
     427cdef inline int celement_inv(nmod_poly_t res, nmod_poly_t a, unsigned long n) except -2:
    429428    raise NotImplementedError
    430429
    431 cdef inline int celement_pow(zmod_poly_t res, zmod_poly_t x, long e, zmod_poly_t modulus, unsigned long n) except -2:
     430cdef inline int celement_pow(nmod_poly_t res, nmod_poly_t x, long e, nmod_poly_t modulus, unsigned long n) except -2:
    432431    """
    433432    EXAMPLE:
    434433        sage: P.<x> = GF(32003)[]
     
    480479        sage: f^5 % g
    481480        7231*x + 17274
    482481    """
    483     cdef zmod_poly_t pow2
    484     cdef zmod_poly_t q
    485     cdef zmod_poly_t tmp
     482    cdef nmod_poly_t pow2
     483    cdef nmod_poly_t q
     484    cdef nmod_poly_t tmp
    486485
    487     zmod_poly_init(q, n)
    488     zmod_poly_init(tmp, n)
     486    nmod_poly_init(q, n)
     487    nmod_poly_init(tmp, n)
    489488
    490     if zmod_poly_degree(x) == 1 and zmod_poly_get_coeff_ui(x,0) == 0 and zmod_poly_get_coeff_ui(x,1) == 1:
    491         zmod_poly_zero(res)
    492         zmod_poly_set_coeff_ui(res,e,1)
     489    if nmod_poly_degree(x) == 1 and nmod_poly_get_coeff_ui(x,0) == 0 and nmod_poly_get_coeff_ui(x,1) == 1:
     490        nmod_poly_zero(res)
     491        nmod_poly_set_coeff_ui(res,e,1)
    493492    elif e == 0:
    494         zmod_poly_zero(res)
    495         zmod_poly_set_coeff_ui(res,0,1)
     493        nmod_poly_zero(res)
     494        nmod_poly_set_coeff_ui(res,0,1)
    496495    elif e == 1:
    497         zmod_poly_set(res, x)
     496        nmod_poly_set(res, x)
    498497    elif e == 2:
    499         zmod_poly_sqr(res, x)
     498        nmod_poly_mul(res, x, x)
    500499    else:
    501500        if res == x:
    502             zmod_poly_set(tmp, x)
     501            nmod_poly_set(tmp, x)
    503502            x = tmp
    504         zmod_poly_init(pow2, n)
    505         zmod_poly_set(pow2, x)
     503        nmod_poly_init(pow2, n)
     504        nmod_poly_set(pow2, x)
    506505        if e % 2:
    507             zmod_poly_set(res, x)
     506            nmod_poly_set(res, x)
    508507        else:
    509             zmod_poly_zero(res)
    510             zmod_poly_set_coeff_ui(res, 0, 1)
     508            nmod_poly_zero(res)
     509            nmod_poly_set_coeff_ui(res, 0, 1)
    511510        e = e >> 1
    512511        while(e != 0):
    513             zmod_poly_sqr(pow2, pow2)
     512            nmod_poly_mul(pow2, pow2, pow2)
    514513            if e % 2:
    515                 zmod_poly_mul(res, res, pow2)
     514                nmod_poly_mul(res, res, pow2)
    516515            e = e >> 1
    517516            if modulus != NULL:
    518                 zmod_poly_divrem(q, res, res, modulus)
    519         zmod_poly_clear(pow2)
     517                nmod_poly_divrem(q, res, res, modulus)
     518        nmod_poly_clear(pow2)
    520519
    521520    if modulus != NULL:
    522         zmod_poly_divrem(q, res, res, modulus)
    523     zmod_poly_clear(q)
    524     zmod_poly_clear(tmp)
     521        nmod_poly_divrem(q, res, res, modulus)
     522    nmod_poly_clear(q)
     523    nmod_poly_clear(tmp)
    525524
    526 cdef inline int celement_gcd(zmod_poly_t res, zmod_poly_t a, zmod_poly_t b, unsigned long n) except -2:
     525cdef inline int celement_gcd(nmod_poly_t res, nmod_poly_t a, nmod_poly_t b, unsigned long n) except -2:
    527526    """
    528527    EXAMPLE:
    529528        sage: P.<x> = GF(32003)[]
     
    563562        True
    564563    """
    565564    if celement_is_zero(b, n):
    566         zmod_poly_set(res, a)
     565        nmod_poly_set(res, a)
    567566        return 0
    568567
    569     zmod_poly_gcd(res, a, b)
    570     cdef unsigned long leadcoeff = zmod_poly_get_coeff_ui(res, zmod_poly_degree(res))
    571     cdef unsigned long modulus = zmod_poly_modulus(res)
    572     if z_gcd(modulus,leadcoeff) == 1:
    573         zmod_poly_make_monic(res, res)
     568    nmod_poly_gcd(res, a, b)
     569    cdef unsigned long leadcoeff = nmod_poly_get_coeff_ui(res, nmod_poly_degree(res))
     570    cdef unsigned long modulus = nmod_poly_modulus(res)
     571    if n_gcd(modulus,leadcoeff) == 1:
     572        nmod_poly_make_monic(res, res)
    574573
    575 cdef inline int celement_xgcd(zmod_poly_t res, zmod_poly_t s, zmod_poly_t t, zmod_poly_t a, zmod_poly_t b, unsigned long n) except -2:
     574cdef inline int celement_xgcd(nmod_poly_t res, nmod_poly_t s, nmod_poly_t t, nmod_poly_t a, nmod_poly_t b, unsigned long n) except -2:
    576575    """
    577576    EXAMPLE:
    578577        sage: P.<x> = GF(32003)[]
     
    611610        sage: (G//d)*d == G
    612611        True
    613612    """
    614     zmod_poly_xgcd(res, s, t, a, b)
     613    nmod_poly_xgcd(res, s, t, a, b)
    615614
    616615
    617616cdef factor_helper(Polynomial_zmod_flint poly, bint squarefree=False):
     
    624623        sage: (prod(P.random_element()^i for i in range(5))).squarefree_decomposition()
    625624        (54) * (x^2 + 55*x + 839) * (x^2 + 48*x + 496)^2 * (x^2 + 435*x + 104)^3 * (x^2 + 176*x + 156)^4
    626625    """
    627     cdef zmod_poly_factor_t factors_c
    628     zmod_poly_factor_init(factors_c)
     626    cdef nmod_poly_factor_t factors_c
     627    nmod_poly_factor_init(factors_c)
    629628
    630629    if squarefree:
    631         zmod_poly_factor_square_free(factors_c, &poly.x)
     630        nmod_poly_factor_squarefree(factors_c, &poly.x)
    632631    else:
    633         zmod_poly_factor(factors_c, &poly.x)
     632        nmod_poly_factor(factors_c, &poly.x)
    634633
    635634    factor_list = []
    636635    cdef Polynomial_zmod_flint t
    637     for i in range(factors_c.num_factors):
     636    for i in range(factors_c.num):
    638637        t = poly._new()
    639         zmod_poly_swap(&t.x, factors_c.factors[i])
    640         factor_list.append((t, factors_c.exponents[i]))
     638        nmod_poly_swap(&t.x, &factors_c.p[i])
     639        factor_list.append((t, factors_c.exp[i]))
    641640
    642     zmod_poly_factor_clear(factors_c)
     641    nmod_poly_factor_clear(factors_c)
    643642
    644643    return Factorization(factor_list, unit=poly.leading_coefficient(),
    645644            sort=(not squarefree))
  • sage/modular/modsym/apply.pyx

    diff --git a/sage/modular/modsym/apply.pyx b/sage/modular/modsym/apply.pyx
    a b  
    5454        fmpz_poly_set_coeff_si(self.g, 1, c)
    5555
    5656        # h = (f**i)*(g**(j-i))
    57         fmpz_poly_power(self.ff, self.f, i)
    58         fmpz_poly_power(self.gg, self.g, j-i)
     57        fmpz_poly_pow(self.ff, self.f, i)
     58        fmpz_poly_pow(self.gg, self.g, j-i)
    5959        fmpz_poly_mul(ans, self.ff, self.gg)
    6060
    6161        return 0
  • sage/rings/fraction_field_FpT.pxd

    diff --git a/sage/rings/fraction_field_FpT.pxd b/sage/rings/fraction_field_FpT.pxd
    a b  
    77from sage.categories.map cimport Section
    88
    99cdef class FpTElement(RingElement):
    10     cdef zmod_poly_t _numer, _denom
     10    cdef nmod_poly_t _numer, _denom
    1111    cdef bint initalized
    1212    cdef long p
    1313
     
    2323    cdef parent
    2424    cdef long degree
    2525    cdef FpTElement cur
    26     cdef zmod_poly_t g
     26    cdef nmod_poly_t g
  • sage/rings/fraction_field_FpT.pyx

    diff --git a/sage/rings/fraction_field_FpT.pyx b/sage/rings/fraction_field_FpT.pyx
    a b  
    105105            numer = parent.poly_ring(numer)
    106106            denom = parent.poly_ring(denom)
    107107        self.p = parent.p
    108         zmod_poly_init(self._numer, self.p)
    109         zmod_poly_init(self._denom, self.p)
     108        nmod_poly_init(self._numer, self.p)
     109        nmod_poly_init(self._denom, self.p)
    110110        self.initalized = True
    111111        cdef long n
    112112        for n, a in enumerate(numer):
    113             zmod_poly_set_coeff_ui(self._numer, n, a)
     113            nmod_poly_set_coeff_ui(self._numer, n, a)
    114114        for n, a in enumerate(denom):
    115             zmod_poly_set_coeff_ui(self._denom, n, a)
     115            nmod_poly_set_coeff_ui(self._denom, n, a)
    116116        if reduce:
    117117            normalize(self._numer, self._denom, self.p)
    118118
     
    127127            sage: del t # indirect doctest
    128128        """
    129129        if self.initalized:
    130             zmod_poly_clear(self._numer)
    131             zmod_poly_clear(self._denom)
     130            nmod_poly_clear(self._numer)
     131            nmod_poly_clear(self._denom)
    132132   
    133133    def __reduce__(self):
    134134        """
     
    152152        cdef FpTElement x = <FpTElement>PY_NEW(FpTElement)
    153153        x._parent = self._parent
    154154        x.p = self.p
    155         zmod_poly_init_precomp(x._numer, x.p, self._numer.p_inv)
    156         zmod_poly_init_precomp(x._denom, x.p, self._numer.p_inv)
     155        nmod_poly_init_preinv(x._numer, x.p, self._numer.mod.ninv)
     156        nmod_poly_init_preinv(x._denom, x.p, self._numer.mod.ninv)
    157157        x.initalized = True
    158158        return x
    159159   
     
    164164        cdef FpTElement x = <FpTElement>PY_NEW(FpTElement)
    165165        x._parent = self._parent
    166166        x.p = self.p
    167         zmod_poly_init2_precomp(x._numer, x.p, self._numer.p_inv, self._numer.length)
    168         zmod_poly_init2_precomp(x._denom, x.p, self._denom.p_inv, self._denom.length)
    169         zmod_poly_set(x._numer, self._numer)
    170         zmod_poly_set(x._denom, self._denom)
     167        nmod_poly_init2_preinv(x._numer, x.p, self._numer.mod.ninv, self._numer.length)
     168        nmod_poly_init2_preinv(x._denom, x.p, self._denom.mod.ninv, self._denom.length)
     169        nmod_poly_set(x._numer, self._numer)
     170        nmod_poly_set(x._denom, self._denom)
    171171        x.initalized = True
    172172        return x
    173173   
     
    196196            t^6 + 3*t^4 + 10*t^3 + 3*t^2 + 1
    197197        """
    198198        cdef Polynomial_zmod_flint res = <Polynomial_zmod_flint>PY_NEW(Polynomial_zmod_flint)
    199         zmod_poly_init2_precomp(&res.x, self.p, self._numer.p_inv, self._numer.length)
    200         zmod_poly_set(&res.x, self._numer)
     199        nmod_poly_init2_preinv(&res.x, self.p, self._numer.mod.ninv, self._numer.length)
     200        nmod_poly_set(&res.x, self._numer)
    201201        res._parent = self._parent.poly_ring
    202202        res._cparent = get_cparent(self._parent.poly_ring)
    203203        return res
     
    227227            t^3
    228228        """
    229229        cdef Polynomial_zmod_flint res = <Polynomial_zmod_flint>PY_NEW(Polynomial_zmod_flint)
    230         zmod_poly_init2_precomp(&res.x, self.p, self._denom.p_inv, self._denom.length)
    231         zmod_poly_set(&res.x, self._denom)
     230        nmod_poly_init2_preinv(&res.x, self.p, self._denom.mod.ninv, self._denom.length)
     231        nmod_poly_set(&res.x, self._denom)
    232232        res._parent = self._parent.poly_ring
    233233        res._cparent = get_cparent(self._parent.poly_ring)
    234234        return res
     
    315315            sage: (1-t)/t
    316316            (16*t + 1)/t
    317317        """
    318         if zmod_poly_degree(self._denom) == 0 and zmod_poly_get_coeff_ui(self._denom, 0) == 1:
     318        if nmod_poly_degree(self._denom) == 0 and nmod_poly_get_coeff_ui(self._denom, 0) == 1:
    319319            return repr(self.numer())
    320320        else:
    321321            numer_s = repr(self.numer())
     
    338338            sage: latex((t + 1)/(t-1))
    339339            \frac{t + 1}{t + 6}
    340340        """
    341         if zmod_poly_degree(self._denom) == 0 and zmod_poly_get_coeff_ui(self._denom, 0) == 1:
     341        if nmod_poly_degree(self._denom) == 0 and nmod_poly_get_coeff_ui(self._denom, 0) == 1:
    342342            return self.numer()._latex_()
    343343        else:
    344344            return "\\frac{%s}{%s}" % (self.numer()._latex_(), self.denom()._latex_())
     
    391391            False
    392392        """
    393393        # They are normalized.
    394         cdef int j = sage_cmp_zmod_poly_t(self._numer, (<FpTElement>other)._numer)
     394        cdef int j = sage_cmp_nmod_poly_t(self._numer, (<FpTElement>other)._numer)
    395395        if j: return j
    396         return sage_cmp_zmod_poly_t(self._denom, (<FpTElement>other)._denom)
     396        return sage_cmp_nmod_poly_t(self._denom, (<FpTElement>other)._denom)
    397397   
    398398    def __hash__(self):
    399399        """
     
    428428            (4*t^2 + 3)/(t + 4)
    429429        """
    430430        cdef FpTElement x = self._copy_c()
    431         zmod_poly_neg(x._numer, x._numer)
     431        nmod_poly_neg(x._numer, x._numer)
    432432        return x
    433433       
    434434    def __invert__(self):
     
    442442            sage: ~a # indirect doctest
    443443            (t + 4)/(t^2 + 2)
    444444        """
    445         if zmod_poly_degree(self._numer) == -1:
     445        if nmod_poly_degree(self._numer) == -1:
    446446            raise ZeroDivisionError
    447447        cdef FpTElement x = self._copy_c()
    448         zmod_poly_swap(x._numer, x._denom)
     448        nmod_poly_swap(x._numer, x._denom)
    449449        return x
    450450
    451451    cpdef ModuleElement _add_(self, ModuleElement _other):
     
    469469        """
    470470        cdef FpTElement other = <FpTElement>_other
    471471        cdef FpTElement x = self._new_c()
    472         zmod_poly_mul(x._numer, self._numer, other._denom)
    473         zmod_poly_mul(x._denom, self._denom, other._numer) # use x._denom as a temp
    474         zmod_poly_add(x._numer, x._numer, x._denom)
    475         zmod_poly_mul(x._denom, self._denom, other._denom)
     472        nmod_poly_mul(x._numer, self._numer, other._denom)
     473        nmod_poly_mul(x._denom, self._denom, other._numer) # use x._denom as a temp
     474        nmod_poly_add(x._numer, x._numer, x._denom)
     475        nmod_poly_mul(x._denom, self._denom, other._denom)
    476476        normalize(x._numer, x._denom, self.p)
    477477        return x
    478478
     
    491491        """
    492492        cdef FpTElement other = <FpTElement>_other
    493493        cdef FpTElement x = self._new_c()
    494         zmod_poly_mul(x._numer, self._numer, other._denom)
    495         zmod_poly_mul(x._denom, self._denom, other._numer) # use x._denom as a temp
    496         zmod_poly_sub(x._numer, x._numer, x._denom)
    497         zmod_poly_mul(x._denom, self._denom, other._denom)
     494        nmod_poly_mul(x._numer, self._numer, other._denom)
     495        nmod_poly_mul(x._denom, self._denom, other._numer) # use x._denom as a temp
     496        nmod_poly_sub(x._numer, x._numer, x._denom)
     497        nmod_poly_mul(x._denom, self._denom, other._denom)
    498498        normalize(x._numer, x._denom, self.p)
    499499        return x
    500500
     
    513513        """
    514514        cdef FpTElement other = <FpTElement>_other
    515515        cdef FpTElement x = self._new_c()
    516         zmod_poly_mul(x._numer, self._numer, other._numer)
    517         zmod_poly_mul(x._denom, self._denom, other._denom)
     516        nmod_poly_mul(x._numer, self._numer, other._numer)
     517        nmod_poly_mul(x._denom, self._denom, other._denom)
    518518        normalize(x._numer, x._denom, self.p)
    519519        return x
    520520
     
    534534            t + 1
    535535        """
    536536        cdef FpTElement other = <FpTElement>_other
    537         if zmod_poly_degree(other._numer) == -1:
     537        if nmod_poly_degree(other._numer) == -1:
    538538            raise ZeroDivisionError
    539539        cdef FpTElement x = self._new_c()
    540         zmod_poly_mul(x._numer, self._numer, other._denom)
    541         zmod_poly_mul(x._denom, self._denom, other._numer)
     540        nmod_poly_mul(x._numer, self._numer, other._denom)
     541        nmod_poly_mul(x._denom, self._denom, other._numer)
    542542        normalize(x._numer, x._denom, self.p)
    543543        return x
    544544   
     
    603603        """
    604604        cdef FpTElement next = self._copy_c()
    605605        cdef long a, lead
    606         cdef zmod_poly_t g
    607         if zmod_poly_degree(self._numer) == -1:
     606        cdef nmod_poly_t g
     607        if nmod_poly_degree(self._numer) == -1:
    608608            # self should be normalized, so denom == 1
    609             zmod_poly_set_coeff_ui(next._numer, 0, 1)
     609            nmod_poly_set_coeff_ui(next._numer, 0, 1)
    610610            return next
    611         lead = zmod_poly_leading(next._numer)
     611        lead = nmod_poly_leading(next._numer)
    612612        if lead < self.p - 1:
    613613            a = mod_inverse_int(lead, self.p)
    614614            # no overflow since self.p < 2^16
    615615            a = a * (lead + 1) % self.p
    616             zmod_poly_scalar_mul(next._numer, next._numer, a)
     616            nmod_poly_scalar_mul_nmod(next._numer, next._numer, a)
    617617        else:
    618618            a = mod_inverse_int(lead, self.p)
    619             zmod_poly_scalar_mul(next._numer, next._numer, a)
     619            nmod_poly_scalar_mul_nmod(next._numer, next._numer, a)
    620620            # now both next._numer and next._denom are monic.  We figure out which is lexicographically bigger:
    621             a = zmod_poly_cmp(next._numer, next._denom)
     621            a = nmod_poly_cmp(next._numer, next._denom)
    622622            if a == 0:
    623623                # next._numer and next._denom are relatively prime, so they're both 1.
    624                 zmod_poly_inc(next._denom, True)
     624                nmod_poly_inc(next._denom, True)
    625625                return next
    626             zmod_poly_set(next._denom, next._numer)
    627             zmod_poly_set(next._numer, self._denom)
     626            nmod_poly_set(next._denom, next._numer)
     627            nmod_poly_set(next._numer, self._denom)
    628628            if a < 0:
    629629                # since next._numer is smaller, we flip and return the inverse.
    630630                return next
    631631            elif a > 0:
    632632                # since next._numer is bigger, we're in the flipped phase.  We flip back, and increment the numerator (until we reach the denominator).
    633                 zmod_poly_init(g, self.p)
     633                nmod_poly_init(g, self.p)
    634634                try:
    635635                    while True:
    636                         zmod_poly_inc(next._numer, True)
    637                         if zmod_poly_equal(next._numer, next._denom):
     636                        nmod_poly_inc(next._numer, True)
     637                        if nmod_poly_equal(next._numer, next._denom):
    638638                            # Since we've reached the denominator, we reset the numerator to 1 and increment the denominator.
    639                             zmod_poly_inc(next._denom, True)
    640                             zmod_poly_zero(next._numer)
    641                             zmod_poly_set_coeff_ui(next._numer, 0, 1)
     639                            nmod_poly_inc(next._denom, True)
     640                            nmod_poly_zero(next._numer)
     641                            nmod_poly_set_coeff_ui(next._numer, 0, 1)
    642642                            break
    643643                        else:
    644644                            # otherwise, we keep incrementing until we have a relatively prime numerator.
    645                             zmod_poly_gcd(g, next._numer, next._denom)
    646                             if zmod_poly_is_one(g):
     645                            nmod_poly_gcd(g, next._numer, next._denom)
     646                            if nmod_poly_is_one(g):
    647647                                break
    648648                finally:
    649                     zmod_poly_clear(g)
     649                    nmod_poly_clear(g)
    650650        return next
    651651       
    652652    cpdef _sqrt_or_None(self):
     
    686686            []
    687687
    688688        """
    689         if zmod_poly_degree(self._numer) == -1:
     689        if nmod_poly_degree(self._numer) == -1:
    690690            return self
    691         if not zmod_poly_sqrt_check(self._numer) or not zmod_poly_sqrt_check(self._denom):
    692             return None
    693         cdef zmod_poly_t numer
    694         cdef zmod_poly_t denom
     691        cdef nmod_poly_t numer
     692        cdef nmod_poly_t denom
    695693        cdef FpTElement res
    696         try:
    697             zmod_poly_init(denom, self.p)
    698             zmod_poly_init(numer, self.p)
    699             if not zmod_poly_sqrt0(numer, self._numer):
    700                 return None
    701             if not zmod_poly_sqrt0(denom, self._denom):
    702                 return None
     694
     695        nmod_poly_init(denom, self.p)
     696        nmod_poly_init(numer, self.p)
     697
     698        if nmod_poly_sqrt(numer, self._numer) and nmod_poly_sqrt(denom, self._denom):
    703699            res = self._new_c()
    704             zmod_poly_swap(numer, res._numer)
    705             zmod_poly_swap(denom, res._denom)
     700            nmod_poly_swap(numer, res._numer)
     701            nmod_poly_swap(denom, res._denom)
    706702            return res
    707         finally:
    708             zmod_poly_clear(numer)
    709             zmod_poly_clear(denom)
    710    
     703        else:
     704            nmod_poly_clear(numer)
     705            nmod_poly_clear(denom)
     706            return None   
     707
    711708    cpdef bint is_square(self):
    712709        """
    713710        Returns True if this element is the square of another element of the fraction field.
     
    792789        assert dummy is None
    793790        cdef FpTElement x = self._new_c()
    794791        if e >= 0:
    795             zmod_poly_pow(x._numer, self._numer, e)
    796             zmod_poly_pow(x._denom, self._denom, e)
     792            nmod_poly_pow(x._numer, self._numer, e)
     793            nmod_poly_pow(x._denom, self._denom, e)
    797794        else:
    798             zmod_poly_pow(x._denom, self._numer, -e)
    799             zmod_poly_pow(x._numer, self._denom, -e)
    800             if zmod_poly_leading(x._denom) != 1:
    801                 a = mod_inverse_int(zmod_poly_leading(x._denom), self.p)
    802                 zmod_poly_scalar_mul(x._numer, x._numer, a)
    803                 zmod_poly_scalar_mul(x._denom, x._denom, a)
     795            nmod_poly_pow(x._denom, self._numer, -e)
     796            nmod_poly_pow(x._numer, self._denom, -e)
     797            if nmod_poly_leading(x._denom) != 1:
     798                a = mod_inverse_int(nmod_poly_leading(x._denom), self.p)
     799                nmod_poly_scalar_mul_nmod(x._numer, x._numer, a)
     800                nmod_poly_scalar_mul_nmod(x._denom, x._denom, a)
    804801        return x
    805802
    806803       
     
    878875            sage: I
    879876            <sage.rings.fraction_field_FpT.FpT_iter object at ...>
    880877        """
    881         zmod_poly_init(self.g, parent.characteristic())
     878        nmod_poly_init(self.g, parent.characteristic())
    882879   
    883880    def __dealloc__(self):
    884881        """
     
    891888            sage: I = FpT_iter(K, 3)
    892889            sage: del I # indirect doctest
    893890        """
    894         zmod_poly_clear(self.g)
     891        nmod_poly_clear(self.g)
    895892   
    896893    def __iter__(self):
    897894        """
     
    971968            next = self.cur._copy_c()
    972969            sig_on()
    973970            while True:
    974                 zmod_poly_inc(next._numer, False)
    975                 if zmod_poly_degree(next._numer) > self.degree:
    976                     zmod_poly_inc(next._denom, True)
    977                     if zmod_poly_degree(next._denom) > self.degree:
     971                nmod_poly_inc(next._numer, False)
     972                if nmod_poly_degree(next._numer) > self.degree:
     973                    nmod_poly_inc(next._denom, True)
     974                    if nmod_poly_degree(next._denom) > self.degree:
    978975                        sig_off()
    979976                        raise StopIteration
    980                     zmod_poly_zero(next._numer)
    981                     zmod_poly_set_coeff_ui(next._numer, 0, 1)
    982                 zmod_poly_gcd(self.g, next._numer, next._denom)
    983                 if zmod_poly_is_one(self.g):
     977                    nmod_poly_zero(next._numer)
     978                    nmod_poly_set_coeff_ui(next._numer, 0, 1)
     979                nmod_poly_gcd(self.g, next._numer, next._denom)
     980                if nmod_poly_is_one(self.g):
    984981                    break
    985982            sig_off()
    986983            self.cur = next
    987984#            self.cur = self.cur.next()
    988 #            if zmod_poly_degree(self.cur._numer) > self.degree:
     985#            if nmod_poly_degree(self.cur._numer) > self.degree:
    989986#                raise StopIteration
    990987        return self.cur
    991988
     
    10351032        cdef FpTElement ans = <FpTElement>PY_NEW(FpTElement)
    10361033        ans._parent = self._codomain
    10371034        ans.p = self.p
    1038         zmod_poly_init(ans._numer, ans.p)
    1039         zmod_poly_init(ans._denom, ans.p)
    1040         zmod_poly_set(ans._numer, &x.x)
    1041         zmod_poly_set_coeff_ui(ans._denom, 0, 1)
     1035        nmod_poly_init(ans._numer, ans.p)
     1036        nmod_poly_init(ans._denom, ans.p)
     1037        nmod_poly_set(ans._numer, &x.x)
     1038        nmod_poly_set_coeff_ui(ans._denom, 0, 1)
    10421039        ans.initalized = True
    10431040        return ans
    10441041
     
    10641061        cdef FpTElement ans = <FpTElement>PY_NEW(FpTElement)
    10651062        ans._parent = self._codomain
    10661063        ans.p = self.p
    1067         zmod_poly_init(ans._numer, ans.p)
    1068         zmod_poly_init(ans._denom, ans.p)
     1064        nmod_poly_init(ans._numer, ans.p)
     1065        nmod_poly_init(ans._denom, ans.p)
    10691066        cdef long r
    1070         zmod_poly_set(ans._numer, &x.x)
     1067        nmod_poly_set(ans._numer, &x.x)
    10711068        if len(args) == 0:
    1072             zmod_poly_set_coeff_ui(ans._denom, 0, 1)
     1069            nmod_poly_set_coeff_ui(ans._denom, 0, 1)
    10731070        elif len(args) == 1:
    10741071            y = args[0]
    10751072            if PY_TYPE_CHECK(y, Integer):
    10761073                r = mpz_fdiv_ui((<Integer>y).value, self.p)
    10771074                if r == 0:
    10781075                    raise ZeroDivisionError
    1079                 zmod_poly_set_coeff_ui(ans._denom, 0, r)
     1076                nmod_poly_set_coeff_ui(ans._denom, 0, r)
    10801077            else:
    10811078                # could use the coerce keyword being set to False to not check this...
    10821079                if not (PY_TYPE_CHECK(y, Element) and y.parent() is self._domain):
    10831080                    # We could special case integers and GF(p) elements here.
    10841081                    y = self._domain(y)
    1085                 zmod_poly_set(ans._denom, &((<Polynomial_zmod_flint?>y).x))
     1082                nmod_poly_set(ans._denom, &((<Polynomial_zmod_flint?>y).x))
    10861083        else:
    10871084            raise ValueError, "FpT only supports two positional arguments"
    10881085        if not (kwds.has_key('reduce') and not kwds['reduce']):
     
    11691166        """
    11701167        cdef FpTElement x = <FpTElement?>_x
    11711168        cdef Polynomial_zmod_flint ans
    1172         if zmod_poly_degree(x._denom) != 0:
     1169        if nmod_poly_degree(x._denom) != 0:
    11731170            normalize(x._numer, x._denom, self.p)
    1174             if zmod_poly_degree(x._denom) != 0:
     1171            if nmod_poly_degree(x._denom) != 0:
    11751172                raise ValueError, "not integral"
    11761173        ans = PY_NEW(Polynomial_zmod_flint)
    1177         if zmod_poly_get_coeff_ui(x._denom, 0) != 1:
     1174        if nmod_poly_get_coeff_ui(x._denom, 0) != 1:
    11781175            normalize(x._numer, x._denom, self.p)
    1179         zmod_poly_init(&ans.x, self.p)
    1180         zmod_poly_set(&ans.x, x._numer)
     1176        nmod_poly_init(&ans.x, self.p)
     1177        nmod_poly_set(&ans.x, x._numer)
    11811178        ans._parent = self._codomain
    11821179        ans._cparent = get_cparent(self._codomain)
    11831180        return ans
     
    12281225        cdef FpTElement ans = <FpTElement>PY_NEW(FpTElement)
    12291226        ans._parent = self._codomain
    12301227        ans.p = self.p
    1231         zmod_poly_init(ans._numer, ans.p)
    1232         zmod_poly_init(ans._denom, ans.p)
    1233         zmod_poly_set_coeff_ui(ans._numer, 0, x.ivalue)
    1234         zmod_poly_set_coeff_ui(ans._denom, 0, 1)
     1228        nmod_poly_init(ans._numer, ans.p)
     1229        nmod_poly_init(ans._denom, ans.p)
     1230        nmod_poly_set_coeff_ui(ans._numer, 0, x.ivalue)
     1231        nmod_poly_set_coeff_ui(ans._denom, 0, 1)
    12351232        ans.initalized = True
    12361233        return ans
    12371234
     
    12571254        cdef FpTElement ans = <FpTElement>PY_NEW(FpTElement)
    12581255        ans._parent = self._codomain
    12591256        ans.p = self.p
    1260         zmod_poly_init(ans._numer, ans.p)
    1261         zmod_poly_init(ans._denom, ans.p)
     1257        nmod_poly_init(ans._numer, ans.p)
     1258        nmod_poly_init(ans._denom, ans.p)
    12621259        cdef long r
    1263         zmod_poly_set_coeff_ui(ans._numer, 0, x.ivalue)
     1260        nmod_poly_set_coeff_ui(ans._numer, 0, x.ivalue)
    12641261        if len(args) == 0:
    1265             zmod_poly_set_coeff_ui(ans._denom, 0, 1)
     1262            nmod_poly_set_coeff_ui(ans._denom, 0, 1)
    12661263        if len(args) == 1:
    12671264            y = args[0]
    12681265            if PY_TYPE_CHECK(y, Integer):
    12691266                r = mpz_fdiv_ui((<Integer>y).value, self.p)
    12701267                if r == 0:
    12711268                    raise ZeroDivisionError
    1272                 zmod_poly_set_coeff_ui(ans._denom, 0, r)
     1269                nmod_poly_set_coeff_ui(ans._denom, 0, r)
    12731270            else:
    12741271                R = self._codomain.ring_of_integers()
    12751272                # could use the coerce keyword being set to False to not check this...
    12761273                if not (PY_TYPE_CHECK(y, Element) and y.parent() is R):
    12771274                    # We could special case integers and GF(p) elements here.
    12781275                    y = R(y)
    1279                 zmod_poly_set(ans._denom, &((<Polynomial_zmod_flint?>y).x))
     1276                nmod_poly_set(ans._denom, &((<Polynomial_zmod_flint?>y).x))
    12801277        else:
    12811278            raise ValueError, "FpT only supports two positional arguments"
    12821279        if not (kwds.has_key('reduce') and not kwds['reduce']):
     
    13731370        """
    13741371        cdef FpTElement x = <FpTElement?>_x
    13751372        cdef IntegerMod_int ans
    1376         if zmod_poly_degree(x._denom) != 0 or zmod_poly_degree(x._numer) > 0:
     1373        if nmod_poly_degree(x._denom) != 0 or nmod_poly_degree(x._numer) > 0:
    13771374            normalize(x._numer, x._denom, self.p)
    1378             if zmod_poly_degree(x._denom) != 0:
     1375            if nmod_poly_degree(x._denom) != 0:
    13791376                raise ValueError, "not integral"
    1380             if zmod_poly_degree(x._numer) > 0:
     1377            if nmod_poly_degree(x._numer) > 0:
    13811378                raise ValueError, "not constant"
    13821379        ans = PY_NEW(IntegerMod_int)
    13831380        ans.__modulus = self._codomain._pyx_order
    1384         if zmod_poly_get_coeff_ui(x._denom, 0) != 1:
     1381        if nmod_poly_get_coeff_ui(x._denom, 0) != 1:
    13851382            normalize(x._numer, x._denom, self.p)
    1386         ans.ivalue = zmod_poly_get_coeff_ui(x._numer, 0)
     1383        ans.ivalue = nmod_poly_get_coeff_ui(x._numer, 0)
    13871384        ans._parent = self._codomain
    13881385        return ans
    13891386
     
    14331430        cdef FpTElement ans = <FpTElement>PY_NEW(FpTElement)
    14341431        ans._parent = self._codomain
    14351432        ans.p = self.p
    1436         zmod_poly_init(ans._numer, ans.p)
    1437         zmod_poly_init(ans._denom, ans.p)
    1438         zmod_poly_set_coeff_ui(ans._numer, 0, mpz_fdiv_ui(x.value, self.p))
    1439         zmod_poly_set_coeff_ui(ans._denom, 0, 1)
     1433        nmod_poly_init(ans._numer, ans.p)
     1434        nmod_poly_init(ans._denom, ans.p)
     1435        nmod_poly_set_coeff_ui(ans._numer, 0, mpz_fdiv_ui(x.value, self.p))
     1436        nmod_poly_set_coeff_ui(ans._denom, 0, 1)
    14401437        ans.initalized = True
    14411438        return ans
    14421439
     
    14641461        cdef FpTElement ans = <FpTElement>PY_NEW(FpTElement)
    14651462        ans._parent = self._codomain
    14661463        ans.p = self.p
    1467         zmod_poly_init(ans._numer, ans.p)
    1468         zmod_poly_init(ans._denom, ans.p)
     1464        nmod_poly_init(ans._numer, ans.p)
     1465        nmod_poly_init(ans._denom, ans.p)
    14691466        cdef long r