.. _nmod-poly:

**nmod_poly.h** -- univariate polynomials over integers mod n (word-size n)
===============================================================================

The :type:`nmod_poly_t` data type represents elements of
`\mathbb{Z}/n\mathbb{Z}[x]` for a fixed modulus `n`. The ``nmod_poly``
module provides routines for memory management, basic arithmetic and
some higher level functions such as GCD, etc.

Each coefficient of an :type:`nmod_poly_t` is of type ``ulong``
and represents an integer reduced modulo the fixed modulus `n`.

Unless otherwise specified, all functions in this section permit
aliasing between their input arguments and between their input and
output arguments.

The :type:`nmod_poly_t` type is a typedef for an array of length 1 of
:type:`nmod_poly_struct`'s. This permits passing parameters of type
:type:`nmod_poly_t` by reference.

In reality one never deals directly with the ``struct`` and simply
deals with objects of type :type:`nmod_poly_t`. For simplicity we will
think of an :type:`nmod_poly_t` as a ``struct``, though in practice to
access fields of this ``struct``, one needs to dereference first,
e.g.\ to access the ``length`` field of an :type:`nmod_poly_t` called
``poly1`` one writes ``poly1->length``.

An :type:`nmod_poly_t` is said to be *normalised* if either ``length``
is zero, or if the leading coefficient of the polynomial is non-zero.
All ``nmod_poly`` functions expect their inputs to be normalised and
for all coefficients to be reduced modulo `n` and unless otherwise
specified they produce output that is normalised with coefficients
reduced modulo `n`.

It is recommended that users do not access the fields of an
:type:`nmod_poly_t` or its coefficient data directly, but make use of
the functions designed for this purpose, detailed below.

Functions in ``nmod_poly`` do all the memory management for the user.
One does not need to specify the maximum length in advance before
using a polynomial object. FLINT reallocates space automatically as
the computation proceeds, if more space is required.

Simple example
--------------

The following example computes the square of the polynomial `5x^3 + 6`
in `\mathbb{Z}/7\mathbb{Z}[x]`.

.. code:: c

   #include "nmod_poly.h"
   int main()
   {
       nmod_poly_t x, y;
       nmod_poly_init(x, 7);
       nmod_poly_init(y, 7);
       nmod_poly_set_coeff_ui(x, 3, 5);
       nmod_poly_set_coeff_ui(x, 0, 6);
       nmod_poly_mul(y, x, x);
       nmod_poly_print(x); flint_printf("\n");
       nmod_poly_print(y); flint_printf("\n");
       nmod_poly_clear(x);
       nmod_poly_clear(y);
   }

The output is:

::

   4 7  6 0 0 5
   7 7  1 0 0 4 0 0 4

Types, macros and constants
-------------------------------------------------------------------------------

.. type:: nmod_poly_struct

.. type:: nmod_poly_t


Memory management
--------------------------------------------------------------------------------

.. function:: void nmod_poly_init(nmod_poly_t poly, ulong n)

    Initialises ``poly``. It will have coefficients modulo `n`.

.. function:: void nmod_poly_init_preinv(nmod_poly_t poly, ulong n, ulong ninv)

    Initialises ``poly``. It will have coefficients modulo `n`.
    The caller supplies a precomputed inverse limb generated by
    :func:`n_preinvert_limb`.

.. function:: void nmod_poly_init_mod(nmod_poly_t poly, const nmod_t mod)

    Initialises ``poly`` using an already initialised modulus ``mod``.

.. function:: void nmod_poly_init2(nmod_poly_t poly, ulong n, slong alloc)

    Initialises ``poly``. It will have coefficients modulo `n`.
    Up to ``alloc`` coefficients may be stored in ``poly``.

.. function:: void nmod_poly_init2_preinv(nmod_poly_t poly, ulong n, ulong ninv, slong alloc)

    Initialises ``poly``. It will have coefficients modulo `n`.
    The caller supplies a precomputed inverse limb generated by
    :func:`n_preinvert_limb`. Up to ``alloc`` coefficients may
    be stored in ``poly``.

.. function:: void nmod_poly_realloc(nmod_poly_t poly, slong alloc)

    Reallocates ``poly`` to the given length. If the current
    length is less than ``alloc``, the polynomial is truncated
    and normalised.  If ``alloc`` is zero, the polynomial is
    cleared.

.. function:: void nmod_poly_clear(nmod_poly_t poly)

    Clears the polynomial and releases any memory it used. The polynomial
    cannot be used again until it is initialised.

.. function:: void nmod_poly_fit_length(nmod_poly_t poly, slong alloc)

    Ensures ``poly`` has space for at least ``alloc`` coefficients.
    This function only ever grows the allocated space, so no data loss can
    occur.

.. function:: void _nmod_poly_normalise(nmod_poly_t poly)

    Internal function for normalising a polynomial so that the top
    coefficient, if there is one at all, is not zero.


Polynomial properties
--------------------------------------------------------------------------------


.. function:: slong nmod_poly_length(const nmod_poly_t poly)

    Returns the length of the polynomial ``poly``. The zero polynomial
    has length zero.

.. function:: slong nmod_poly_degree(const nmod_poly_t poly)

    Returns the degree of the polynomial ``poly``. The zero polynomial
    is deemed to have degree `-1`.

.. function:: ulong nmod_poly_modulus(const nmod_poly_t poly)

    Returns the modulus of the polynomial ``poly``. This will be a
    positive integer.

.. function:: flint_bitcnt_t nmod_poly_max_bits(const nmod_poly_t poly)

    Returns the maximum number of bits of any coefficient of ``poly``.

.. function:: int nmod_poly_is_unit(const nmod_poly_t poly)

   Returns `1` if the polynomial is a nonzero constant (in the case of prime
   modulus, this is equivalent to being a unit), otherwise `0`.

.. function:: int nmod_poly_is_monic(const nmod_poly_t poly)

   Returns `1` if the polynomial is monic, i.e. nonzero with leading
   coefficient `1`, otherwise `0`.


Assignment and basic manipulation
--------------------------------------------------------------------------------


.. function:: void nmod_poly_set(nmod_poly_t a, const nmod_poly_t b)

    Sets ``a`` to a copy of ``b``.

.. function:: void nmod_poly_swap(nmod_poly_t poly1, nmod_poly_t poly2)

    Efficiently swaps ``poly1`` and ``poly2`` by swapping pointers
    internally.

.. function:: void nmod_poly_zero(nmod_poly_t res)

    Sets ``res`` to the zero polynomial.

.. function:: void nmod_poly_truncate(nmod_poly_t poly, slong len)

    Truncates ``poly`` to the given length and normalises it.
    If ``len`` is greater than the current length of ``poly``,
    then nothing happens.

.. function:: void nmod_poly_set_trunc(nmod_poly_t res, const nmod_poly_t poly, slong len)

    Notionally truncate ``poly`` to length ``len`` and set ``res`` to the
    result. The result is normalised.

.. function:: void _nmod_poly_reverse(nn_ptr output, nn_srcptr input, slong len, slong m)

    Sets ``output`` to the reverse of ``input``, which is of length
    ``len``, but thinking of it as a polynomial of length ``m``,
    notionally zero-padded if necessary. The length ``m`` must be
    non-negative, but there are no other restrictions. The polynomial
    ``output`` must have space for ``m`` coefficients. Supports
    aliasing of ``output`` and ``input``, but the behaviour is
    undefined in case of partial overlap.

.. function:: void nmod_poly_reverse(nmod_poly_t output, const nmod_poly_t input, slong m)

    Sets ``output`` to the reverse of ``input``, thinking of it as
    a polynomial of length ``m``, notionally zero-padded if necessary).
    The length ``m`` must be non-negative, but there are no other
    restrictions. The output polynomial will be set to length ``m``
    and then normalised.


Randomization
--------------------------------------------------------------------------------


.. function:: void nmod_poly_randtest(nmod_poly_t poly, flint_rand_t state, slong len)

    Generates a random polynomial with length up to ``len``.

.. function:: void nmod_poly_randtest_monic(nmod_poly_t poly, flint_rand_t state, slong len)

    Generates a random monic polynomial with length ``len``.

.. function:: void nmod_poly_randtest_trinomial(nmod_poly_t poly, flint_rand_t state, slong len)

    Generates a random monic trinomial of length ``len``.

.. function:: void nmod_poly_randtest_pentomial(nmod_poly_t poly, flint_rand_t state, slong len)

    Generates a random monic pentomial of length ``len``.

Construction of irreducible polynomials
--------------------------------------------------------------------------------

The following functions assume a prime modulus.

.. function:: void nmod_poly_minimal_irreducible(nmod_poly_t res, ulong n)

    Generates a monic irreducible polynomial of degree ``n`` with minimal
    weight (minimal number of nonzero terms). We generate a binomial
    if possible, otherwise a trinomial, etc.
    It is conjectured that one never needs more than a pentanomial
    modulo `p = 2` and a tetranomial modulo `p > 2`.

    More specifically, this function returns the first among all minimal-weight
    polynomials in the following ordering.
    Firstly, for trinomials, `x^n + a x^k + b` comes before
    its monic reversal `x^n + a' x^{n-k} + b'` if `k < n - k`.
    Secondly, writing
    `f = x^n + a_1 x^{k_1} + a_2 x^{k_2} + \ldots + a_t x^{k_t}`
    with `n > k_1 > k_2 \ldots > k_t`, we order tuples
    `(a_1, \ldots, a_{t}, k_1, \ldots, k_t)` lexicographically.
    We thus favor polynomials with smaller coefficients
    (all 1 if possible), and secondly with smaller degrees for the
    middle terms.

.. function:: void nmod_poly_randtest_irreducible(nmod_poly_t poly, flint_rand_t state, slong len)

    Generates a random irreducible polynomial with length up to ``len``.

.. function:: void nmod_poly_randtest_monic_irreducible(nmod_poly_t poly, flint_rand_t state, slong len)

    Generates a random monic irreducible polynomial with length ``len``.

.. function:: void nmod_poly_randtest_monic_primitive(nmod_poly_t poly, flint_rand_t state, slong len)

    Generates a random monic irreducible primitive polynomial with
    length ``len``.

.. function:: int nmod_poly_randtest_trinomial_irreducible(nmod_poly_t poly, flint_rand_t state, slong len, slong max_attempts)

    Attempts to set ``poly`` to a monic irreducible trinomial of
    length ``len``.  It will generate up to ``max_attempts``
    trinomials in attempt to find an irreducible one.  If
    ``max_attempts`` is ``0``, then it will keep generating
    trinomials until an irreducible one is found.  Returns `1` if one
    is found and `0` otherwise.

.. function:: int nmod_poly_randtest_pentomial_irreducible(nmod_poly_t poly, flint_rand_t state, slong len, slong max_attempts)

    Attempts to set ``poly`` to a monic irreducible pentomial of
    length ``len``.  It will generate up to ``max_attempts``
    pentomials in attempt to find an irreducible one.  If
    ``max_attempts`` is ``0``, then it will keep generating
    pentomials until an irreducible one is found.  Returns `1` if one
    is found and `0` otherwise.

.. function:: void nmod_poly_randtest_sparse_irreducible(nmod_poly_t poly, flint_rand_t state, slong len)

    Attempts to set ``poly`` to a sparse, monic irreducible polynomial
    with length ``len``.  It attempts to find an irreducible
    trinomial.  If that does not succeed, it attempts to find a
    irreducible pentomial.  If that fails, then ``poly`` is just
    set to a random monic irreducible polynomial.


Getting and setting coefficients
--------------------------------------------------------------------------------


.. function:: ulong nmod_poly_get_coeff_ui(const nmod_poly_t poly, slong j)

    Returns the coefficient of ``poly`` at index ``j``, where
    coefficients are numbered with zero being the constant coefficient,
    and returns it as an ``ulong``. If ``j`` refers to a
    coefficient beyond the end of ``poly``, zero is returned.

.. function:: void nmod_poly_set_coeff_ui(nmod_poly_t poly, slong j, ulong c)

    Sets the coefficient of ``poly`` at index ``j``, where
    coefficients are numbered with zero being the constant coefficient,
    to the value ``c`` reduced modulo the modulus of ``poly``.
    If ``j`` refers to a coefficient beyond the current end of ``poly``,
    the polynomial is first resized, with intervening coefficients being
    set to zero.


Input and output
--------------------------------------------------------------------------------


.. function:: char * nmod_poly_get_str(const nmod_poly_t poly)

    Writes ``poly`` to a string representation. The format is as
    described for :func:`nmod_poly_print`. The string must be freed by the
    user when finished. For this it is sufficient to call :func:`flint_free`.

.. function:: char * nmod_poly_get_str_pretty(const nmod_poly_t poly, const char * x)

    Writes ``poly`` to a pretty string representation. The format is as
    described for :func:`nmod_poly_print_pretty`. The string must be freed
    by the user when finished. For this it is sufficient to call
    :func:`flint_free`.

    It is assumed that the top coefficient is non-zero.

.. function:: int nmod_poly_set_str(nmod_poly_t poly, const char * s)

    Reads ``poly`` from a string ``s``. The format is as described
    for :func:`nmod_poly_print`. If a polynomial in the correct format
    is read, a positive value is returned, otherwise a non-positive value
    is returned.

.. function:: int nmod_poly_print(const nmod_poly_t a)

    Prints the polynomial to ``stdout``. The length is printed,
    followed by a space, then the modulus. If the length is zero this is
    all that is printed, otherwise two spaces followed by a space
    separated list of coefficients is printed, beginning with the constant
    coefficient.

    In case of success, returns a positive value.  In case of failure,
    returns a non-positive value.

.. function:: int nmod_poly_print_pretty(const nmod_poly_t a, const char * x)

    Prints the polynomial to ``stdout`` using the string ``x`` to
    represent the indeterminate.

    It is assumed that the top coefficient is non-zero.

    In case of success, returns a positive value.  In case of failure,
    returns a non-positive value.

.. function:: int nmod_poly_fread(FILE * f, nmod_poly_t poly)

    Reads ``poly`` from the file stream ``f``. If this is a file
    that has just been written, the file should be closed then opened
    again. The format is as described for :func:`nmod_poly_print`. If a
    polynomial in the correct format is read, a positive value is returned,
    otherwise a non-positive value is returned.

.. function:: int nmod_poly_fprint(FILE * f, const nmod_poly_t poly)

    Writes a polynomial to the file stream ``f``. If this is a file
    then the file should be closed and reopened before being read.
    The format is as described for :func:`nmod_poly_print`. If the
    polynomial is written correctly, a positive value is returned,
    otherwise a non-positive value is returned.

    In case of success, returns a positive value.  In case of failure,
    returns a non-positive value.

.. function:: int nmod_poly_fprint_pretty(FILE * f, const nmod_poly_t poly, const char * x)

    Writes a polynomial to the file stream ``f``. If this is a file
    then the file should be closed and reopened before being read.
    The format is as described for :func:`nmod_poly_print_pretty`. If the
    polynomial is written correctly, a positive value is returned,
    otherwise a non-positive value is returned.

    It is assumed that the top coefficient is non-zero.

    In case of success, returns a positive value.  In case of failure,
    returns a non-positive value.

.. function:: int nmod_poly_read(nmod_poly_t poly)

    Read ``poly`` from ``stdin``. The format is as described for
    :func:`nmod_poly_print`. If a polynomial in the correct format is read, a
    positive value is returned, otherwise a non-positive value is returned.


Comparison
--------------------------------------------------------------------------------


.. function:: int nmod_poly_equal(const nmod_poly_t a, const nmod_poly_t b)

    Returns `1` if the polynomials are equal, otherwise `0`.

.. function:: int nmod_poly_equal_nmod(const nmod_poly_t poly, ulong cst)

    Returns `1` if the polynomial ``poly`` is constant, equal to ``cst``,
    otherwise `0`.
    ``cst`` is assumed to be already reduced, i.e. less than the modulus of
    ``poly``.

.. function:: int nmod_poly_equal_ui(const nmod_poly_t poly, ulong cst)

    Returns `1` if the polynomial ``poly`` is constant and equal to ``cst`` up to
    reduction modulo the modulus of ``poly``, otherwise returns `0`.

.. function:: int nmod_poly_equal_trunc(const nmod_poly_t poly1, const nmod_poly_t poly2, slong n)

    Notionally truncate ``poly1`` and ``poly2`` to length `n` and return
    `1` if the truncations are equal, otherwise return `0`.

.. function:: int nmod_poly_is_zero(const nmod_poly_t poly)

    Returns `1` if the polynomial ``poly`` is the zero polynomial,
    otherwise returns `0`.

.. function:: int nmod_poly_is_one(const nmod_poly_t poly)

    Returns `1` if the polynomial ``poly`` is the constant polynomial 1,
    otherwise returns `0`.

.. function:: int nmod_poly_is_gen(const nmod_poly_t poly)

   Returns `1` if the polynomial is the generating indeterminate (i.e. has
   degree `1`, constant coefficient `0`, and leading coefficient `1`), otherwise
   returns `0`.


Shifting
--------------------------------------------------------------------------------


.. function:: void _nmod_poly_shift_left(nn_ptr res, nn_srcptr poly, slong len, slong k)

    Sets ``(res, len + k)`` to ``(poly, len)`` shifted left by
    ``k`` coefficients. Assumes that ``res`` has space for
    ``len + k`` coefficients.

.. function:: void nmod_poly_shift_left(nmod_poly_t res, const nmod_poly_t poly, slong k)

    Sets ``res`` to ``poly`` shifted left by ``k`` coefficients,
    i.e. multiplied by `x^k`.

.. function:: void _nmod_poly_shift_right(nn_ptr res, nn_srcptr poly, slong len, slong k)

    Sets ``(res, len - k)`` to ``(poly, len)`` shifted left by
    ``k`` coefficients. It is assumed that ``k <= len`` and that
    ``res`` has space for at least ``len - k`` coefficients.

.. function:: void nmod_poly_shift_right(nmod_poly_t res, const nmod_poly_t poly, slong k)

    Sets ``res`` to ``poly`` shifted right by ``k`` coefficients,
    i.e. divide by `x^k` and throw away the remainder. If ``k`` is
    greater than or equal to the length of ``poly``, the result is the
    zero polynomial.


Addition and subtraction
--------------------------------------------------------------------------------


.. function:: void _nmod_poly_add(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, nmod_t mod)

    Sets ``res`` to the sum of ``(poly1, len1)`` and
    ``(poly2, len2)``. There are no restrictions on the lengths.

.. function:: void nmod_poly_add(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2)

    Sets ``res`` to the sum of ``poly1`` and ``poly2``.

.. function:: void nmod_poly_add_series(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2, slong n)

    Notionally truncate ``poly1`` and ``poly2`` to length `n` and set
    ``res`` to the sum.

.. function:: void _nmod_poly_sub(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, nmod_t mod)

    Sets ``res`` to the difference of ``(poly1, len1)`` and
    ``(poly2, len2)``. There are no restrictions on the lengths.

.. function:: void nmod_poly_sub(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2)

    Sets ``res`` to the difference of ``poly1`` and ``poly2``.

.. function:: void nmod_poly_sub_series(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2, slong n)

    Notionally truncate ``poly1`` and ``poly2`` to length `n` and set
    ``res`` to the difference.

.. function:: void nmod_poly_neg(nmod_poly_t res, const nmod_poly_t poly)

    Sets ``res`` to the negation of ``poly``.


Scalar multiplication and division
--------------------------------------------------------------------------------


.. function:: void nmod_poly_scalar_mul_nmod(nmod_poly_t res, const nmod_poly_t poly, ulong c)

    Sets ``res`` to ``poly`` multiplied by `c`. The element `c` is assumed
    to be less than the modulus of ``poly``.

.. function:: void nmod_poly_scalar_addmul_nmod(nmod_poly_t res, const nmod_poly_t poly, ulong c)

    Adds ``poly`` multiplied by `c` to ``res``. The element `c` is assumed
    to be less than the modulus of ``poly``.

.. function:: void _nmod_poly_make_monic(nn_ptr res, nn_srcptr poly, slong len, nmod_t mod)

    Requires that ``res`` and ``poly`` have length at least ``len``, with ``len
    > 0``, and that ``poly[len-1]`` is invertible modulo ``mod.n``. Sets
    ``res[i]`` to the modular product of `c` and ``poly[i]`` for `i` from `0`
    to ``len-1``, where `c` is the inverse of ``poly[len-1]``.

.. function:: void nmod_poly_make_monic(nmod_poly_t res, const nmod_poly_t poly)

    Sets ``res`` to be the scalar multiple of ``poly`` with leading coefficient
    one. If ``poly`` is zero, an exception is raised.


Bit packing and unpacking
--------------------------------------------------------------------------------


.. function:: void _nmod_poly_bit_pack(nn_ptr res, nn_srcptr poly, slong len, flint_bitcnt_t bits)

    Packs ``len`` coefficients of ``poly`` into fields of the given
    number of bits in the large integer ``res``, i.e. evaluates
    ``poly`` at ``2^bits`` and store the result in ``res``.
    Assumes ``len > 0`` and ``bits > 0``. Also assumes that no
    coefficient of ``poly`` is bigger than ``bits/2`` bits. We
    also assume ``bits < 3 * FLINT_BITS``.

.. function:: void _nmod_poly_bit_unpack(nn_ptr res, slong len, nn_srcptr mpn, ulong bits, nmod_t mod)

    Unpacks ``len`` coefficients stored in the big integer ``mpn``
    in bit fields of the given number of bits, reduces them modulo the
    given modulus, then stores them in the polynomial ``res``.
    We assume ``len > 0`` and ``3 * FLINT_BITS > bits > 0``.
    There are no restrictions on the size of the actual coefficients as
    stored within the bitfields.

.. function:: void nmod_poly_bit_pack(fmpz_t f, const nmod_poly_t poly, flint_bitcnt_t bit_size)

    Packs ``poly`` into bitfields of size ``bit_size``, writing the
    result to ``f``.

.. function:: void nmod_poly_bit_unpack(nmod_poly_t poly, const fmpz_t f, flint_bitcnt_t bit_size)

    Unpacks the polynomial from fields of size ``bit_size`` as
    represented by the integer ``f``.


.. function:: void _nmod_poly_KS2_pack1(nn_ptr res, nn_srcptr op, slong n, slong s, ulong b, ulong k, slong r)

    Same as ``_nmod_poly_KS2_pack``, but requires ``b <= FLINT_BITS``.

.. function:: void _nmod_poly_KS2_pack(nn_ptr res, nn_srcptr op, slong n, slong s, ulong b, ulong k, slong r)

    Bit packing routine used by KS2 and KS4 multiplication.

.. function:: void _nmod_poly_KS2_unpack1(nn_ptr res, nn_srcptr op, slong n, ulong b, ulong k)

    Same as ``_nmod_poly_KS2_unpack``, but requires ``b <= FLINT_BITS``
    (i.e. writes one word per coefficient).

.. function:: void _nmod_poly_KS2_unpack2(nn_ptr res, nn_srcptr op, slong n, ulong b, ulong k)

    Same as ``_nmod_poly_KS2_unpack``, but requires
    ``FLINT_BITS < b <= 2 * FLINT_BITS`` (i.e. writes two words per
    coefficient).

.. function:: void _nmod_poly_KS2_unpack3(nn_ptr res, nn_srcptr op, slong n, ulong b, ulong k)

    Same as ``_nmod_poly_KS2_unpack``, but requires
    ``2 * FLINT_BITS < b < 3 * FLINT_BITS`` (i.e. writes three words per
    coefficient).

.. function:: void _nmod_poly_KS2_unpack(nn_ptr res, nn_srcptr op, slong n, ulong b, ulong k)

    Bit unpacking code used by KS2 and KS4 multiplication.



KS2/KS4 Reduction
--------------------------------------------------------------------------------


.. function:: void _nmod_poly_KS2_reduce(nn_ptr res, slong s, nn_srcptr op, slong n, ulong w, nmod_t mod)

    Reduction code used by KS2 and KS4 multiplication.

.. function:: void _nmod_poly_KS2_recover_reduce1(nn_ptr res, slong s, nn_srcptr op1, nn_srcptr op2, slong n, ulong b, nmod_t mod)

    Same as ``_nmod_poly_KS2_recover_reduce``, but requires
    ``0 < 2 * b <= FLINT_BITS``.

.. function:: void _nmod_poly_KS2_recover_reduce2(nn_ptr res, slong s, nn_srcptr op1, nn_srcptr op2, slong n, ulong b, nmod_t mod)

    Same as ``_nmod_poly_KS2_recover_reduce``, but requires
    ``FLINT_BITS < 2 * b < 2*FLINT_BITS``.

.. function:: void _nmod_poly_KS2_recover_reduce2b(nn_ptr res, slong s, nn_srcptr op1, nn_srcptr op2, slong n, ulong b, nmod_t mod)

    Same as ``_nmod_poly_KS2_recover_reduce``, but requires
    ``b == FLINT_BITS``.

.. function:: void _nmod_poly_KS2_recover_reduce3(nn_ptr res, slong s, nn_srcptr op1, nn_srcptr op2, slong n, ulong b, nmod_t mod)

    Same as ``_nmod_poly_KS2_recover_reduce``, but requires
    ``2 * FLINT_BITS < 2 * b <= 3 * FLINT_BITS``.

.. function:: void _nmod_poly_KS2_recover_reduce(nn_ptr res, slong s, nn_srcptr op1, nn_srcptr op2, slong n, ulong b, nmod_t mod)

    Reduction code used by KS4 multiplication.



Multiplication
--------------------------------------------------------------------------------


.. function:: void _nmod_poly_mul_classical(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, nmod_t mod)

    Sets ``(res, len1 + len2 - 1)`` to the product of ``(poly1, len1)``
    and ``(poly2, len2)``. Assumes ``len1 >= len2 > 0``. Aliasing of
    inputs and output is not permitted.

.. function:: void nmod_poly_mul_classical(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2)

    Sets ``res`` to the product of ``poly1`` and ``poly2``.

.. function:: void _nmod_poly_mullow_classical(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong trunc, nmod_t mod)

    Sets ``res`` to the lower ``trunc`` coefficients of the product of
    ``(poly1, len1)`` and ``(poly2, len2)``. Assumes that
    ``len1 >= len2 > 0`` and ``trunc > 0``. Aliasing of inputs and
    output is not permitted.

.. function:: void nmod_poly_mullow_classical(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2, slong trunc)

    Sets ``res`` to the lower ``trunc`` coefficients of the product
    of ``poly1`` and ``poly2``.

.. function:: void _nmod_poly_mulhigh_classical(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong start, nmod_t mod)

    Computes the product of ``(poly1, len1)`` and ``(poly2, len2)``
    and writes the coefficients from ``start`` onwards into the high
    coefficients of ``res``, the remaining coefficients being arbitrary
    but reduced.  Assumes that ``len1 >= len2 > 0``. Aliasing of inputs
    and output is not permitted.

.. function:: void nmod_poly_mulhigh_classical(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2, slong start)

    Computes the product of ``poly1`` and ``poly2`` and writes the
    coefficients from ``start`` onwards into the high coefficients of
    ``res``, the remaining coefficients being arbitrary but reduced.

.. function:: void _nmod_poly_mul_KS(nn_ptr out, nn_srcptr in1, slong len1, nn_srcptr in2, slong len2, flint_bitcnt_t bits, nmod_t mod)

    Sets ``res`` to the product of ``in1`` and ``in2``
    assuming the output coefficients are at most the given number of
    bits wide. If ``bits`` is set to `0` an appropriate value is
    computed automatically.  Assumes that ``len1 >= len2 > 0``.

.. function:: void nmod_poly_mul_KS(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2, flint_bitcnt_t bits)

    Sets ``res`` to the product of ``poly1`` and ``poly2``
    assuming the output coefficients are at most the given number of
    bits wide. If ``bits`` is set to `0` an appropriate value
    is computed automatically.

.. function:: void _nmod_poly_mul_KS2(nn_ptr res, nn_srcptr op1, slong n1, nn_srcptr op2, slong n2, nmod_t mod)

    Sets ``res`` to the product of ``op1`` and ``op2``.
    Assumes that ``len1 >= len2 > 0``.

.. function:: void nmod_poly_mul_KS2(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2)

    Sets ``res`` to the product of ``poly1`` and ``poly2``.

.. function:: void _nmod_poly_mul_KS4(nn_ptr res, nn_srcptr op1, slong n1, nn_srcptr op2, slong n2, nmod_t mod)

    Sets ``res`` to the product of ``op1`` and ``op2``.
    Assumes that ``len1 >= len2 > 0``.

.. function:: void nmod_poly_mul_KS4(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2)

    Sets ``res`` to the product of ``poly1`` and ``poly2``.

.. function:: void _nmod_poly_mullow_KS(nn_ptr out, nn_srcptr in1, slong len1, nn_srcptr in2, slong len2, flint_bitcnt_t bits, slong n, nmod_t mod)

    Sets ``out`` to the low `n` coefficients of ``in1`` of length
    ``len1`` times ``in2`` of length ``len2``. The output must have
    space for ``n`` coefficients. We assume that ``len1 >= len2 > 0``
    and that ``0 < n <= len1 + len2 - 1``.

.. function:: void nmod_poly_mullow_KS(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2, flint_bitcnt_t bits, slong n)

    Set ``res`` to the low `n` coefficients of ``in1`` of length
    ``len1`` times ``in2`` of length ``len2``.

.. function:: void _nmod_poly_mul(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, nmod_t mod)

    Sets ``res`` to the product of ``poly1`` of length ``len1``
    and ``poly2`` of length ``len2``. Assumes ``len1 >= len2 > 0``.
    No aliasing is permitted between the inputs and the output.

.. function:: void nmod_poly_mul(nmod_poly_t res, const nmod_poly_t poly, const nmod_poly_t poly2)

    Sets ``res`` to the product of ``poly1`` and ``poly2``.

.. function:: void _nmod_poly_mullow(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong n, nmod_t mod)

    Sets ``res`` to the first ``n`` coefficients of the
    product of ``poly1`` of length ``len1`` and ``poly2`` of
    length ``len2``. It is assumed that ``0 < n <= len1 + len2 - 1``
    and that ``len1 >= len2 > 0``. No aliasing of inputs and output
    is permitted.

.. function:: void nmod_poly_mullow(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2, slong trunc)

    Sets ``res`` to the first ``trunc`` coefficients of the
    product of ``poly1`` and ``poly2``.

.. function:: void _nmod_poly_mulhigh(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong n, nmod_t mod)

    Sets all but the low `n` coefficients of ``res`` to the
    corresponding coefficients of the product of ``poly1`` of length
    ``len1`` and ``poly2`` of length ``len2``, the other
    coefficients being arbitrary. It is assumed that
    ``len1 >= len2 > 0`` and that ``0 < n <= len1 + len2 - 1``.
    Aliasing of inputs and output is not permitted.

.. function:: void nmod_poly_mulhigh(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2, slong n)

    Sets all but the low `n` coefficients of ``res`` to the
    corresponding coefficients of the product of ``poly1`` and
    ``poly2``, the remaining coefficients being arbitrary.

.. function:: int _nmod_poly_mullow_want_fft_small(slong len1, slong len2, slong n, int squaring, nmod_t mod)

    Estimate whether *fft_small* multiplication should be used instead of
    other multiplication algorithms, given inputs of length *len1* and *len2*
    and output truncation to length *n*.

.. function:: int _nmod_poly_mullow_fft_small_repack(nn_ptr z, nn_srcptr a, slong an, nn_srcptr b, slong bn, slong zn, nmod_t mod)

    Internal helper function for :func:`_nmod_poly_mullow_fft_small`: if the
    inputs are small enough to perform a repacked convolution of half the
    length, multiply and return 1, otherwise do nothing and return 0.
    The conditions on the arguments are the same as for :func:`_nmod_poly_mullow`.

.. function:: void _nmod_poly_mullow_fft_small(nn_ptr z, nn_srcptr a, slong an, nn_srcptr b, slong bn, slong zn, nmod_t mod)

    Low multiplication via the *fft_small* module. Throws an error
    if *fft_small* is not available. The conditions on the arguments
    are the same as for :func:`_nmod_poly_mullow`.

.. function:: void _nmod_poly_mulmod(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, nn_srcptr f, slong lenf, nmod_t mod)

    Sets ``res`` to the remainder of the product of ``poly1`` and
    ``poly2`` upon polynomial division by ``f``.

    It is required that ``len1 + len2 - lenf > 0``, which is equivalent
    to requiring that the result will actually be reduced. Otherwise, simply
    use ``_nmod_poly_mul`` instead.

    Aliasing of ``f`` and ``res`` is not permitted.

.. function:: void nmod_poly_mulmod(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2, const nmod_poly_t f)

    Sets ``res`` to the remainder of the product of ``poly1`` and
    ``poly2`` upon polynomial division by ``f``.

.. function:: void _nmod_poly_mulmod_preinv(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, nn_srcptr f, slong lenf, nn_srcptr finv, slong lenfinv, nmod_t mod)

    Sets ``res`` to the remainder of the product of ``poly1`` and
    ``poly2`` upon polynomial division by ``f``.

    It is required that ``finv`` is the inverse of the reverse of ``f``
    mod ``x^lenf``. It is required that ``len1 + len2 - lenf > 0``,
    which is equivalent to requiring that the result will actually be reduced.
    It is required that ``len1 < lenf`` and ``len2 < lenf``.
    Otherwise, simply use ``_nmod_poly_mul`` instead.

    Aliasing of ```res`` with any of the inputs is not permitted.

.. function:: void nmod_poly_mulmod_preinv(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2, const nmod_poly_t f, const nmod_poly_t finv)

    Sets ``res`` to the remainder of the product of ``poly1`` and
    ``poly2`` upon polynomial division by ``f``. ``finv`` is the
    inverse of the reverse of ``f``. It is required that ``poly1`` and
    ``poly2`` are reduced modulo ``f``.

Preconditioned modular multiplication
--------------------------------------------------------------------------------

.. type:: nmod_poly_mulmod_precond_struct
          nmod_poly_mulmod_precond_t

    Stores precomputed data for evaluating `ab \bmod d` where both `a`
    and `d` are fixed.

.. function:: void _nmod_poly_mulmod_precond_init_method(nmod_poly_mulmod_precond_t precond, nn_srcptr a, slong alen, nn_srcptr d, slong dlen, nn_srcptr dinv, slong lendinv, int method, nmod_t mod)
              void nmod_poly_mulmod_precond_init_method(nmod_poly_mulmod_precond_t precond, const nmod_poly_t a, const nmod_poly_t d, const nmod_poly_t dinv, int method)
              void _nmod_poly_mulmod_precond_init_num(nmod_poly_mulmod_precond_t precond, nn_srcptr a, slong alen, nn_srcptr d, slong dlen, nn_srcptr dinv, slong lendinv, slong num, nmod_t mod)
              void nmod_poly_mulmod_precond_init_num(nmod_poly_mulmod_precond_t precond, const nmod_poly_t a, const nmod_poly_t d, const nmod_poly_t dinv, slong num)

    Initialize ``precond`` for computing  `ab \bmod d`.
    It is assumed that `a` is already reduced modulo `d`.
    The *method* parameter must be one of the following:

    * ``NMOD_POLY_MULMOD_PRECOND_NONE`` (no precomputation; multiplication will simply delegate to :func:`_nmod_poly_mulmod_preinv`)

    * ``NMOD_POLY_MULMOD_PRECOND_SHOUP`` (use Shoup multiplication)

    * ``NMOD_POLY_MULMOD_PRECOND_MATRIX`` (use the matrix algorithm)

    The *num* versions of these functions attempt to choose the optimal
    method automatically assuming that one intends to perform *num*
    multiplications.

    Shallow references to ``a``, ``d`` and ``dinv`` may be stored
    in ``precond``; the original objects must therefore be kept alive
    without modification as long as ``precond`` is used.
    The user must supply the precomputed inverse of ``d``, with the same
    meaning as in :func:`_nmod_poly_mulmod_preinv` and :func:`nmod_poly_mulmod_preinv`.

.. function:: void nmod_poly_mulmod_precond_clear(nmod_poly_mulmod_precond_t precond)

    Clears ``precond``, freeing any allocated memory.

.. function:: void _nmod_poly_mulmod_precond(nn_ptr res, const nmod_poly_mulmod_precond_t precond, nn_srcptr b, slong blen, nmod_t mod)
              void nmod_poly_mulmod_precond(nmod_poly_t res, const nmod_poly_mulmod_precond_t precond, const nmod_poly_t b)

    Compute `ab \bmod d` where both `a` and `d` are fixed and represented by
    the ``precond`` object. We require that `b` is already reduced modulo `d`.
    The underscore method requires nonzero lengths and does not allow aliasing
    between the output and any inputs (including ``a`` and ``d``).
    The non-underscore method allows aliasing between ``b`` and ``res``.


Powering
--------------------------------------------------------------------------------


.. function:: void _nmod_poly_pow_binexp(nn_ptr res, nn_srcptr poly, slong len, ulong e, nmod_t mod)

    Raises ``poly`` of length ``len`` to the power ``e`` and sets
    ``res`` to the result. We require that ``res`` has enough space
    for ``(len - 1)*e + 1`` coefficients. Assumes that ``len > 0``,
    ``e > 1``. Aliasing is not permitted. Uses the binary exponentiation
    method.

.. function:: void nmod_poly_pow_binexp(nmod_poly_t res, const nmod_poly_t poly, ulong e)

    Raises ``poly`` to the power ``e`` and sets ``res`` to the
    result. Uses the binary exponentiation method.

.. function:: void _nmod_poly_pow(nn_ptr res, nn_srcptr poly, slong len, ulong e, nmod_t mod)

    Raises ``poly`` of length ``len`` to the power ``e`` and sets
    ``res`` to the result. We require that ``res`` has enough space
    for ``(len - 1)*e + 1`` coefficients. Assumes that ``len > 0``,
    ``e > 1``. Aliasing is not permitted.

.. function:: void nmod_poly_pow(nmod_poly_t res, const nmod_poly_t poly, ulong e)

    Raises ``poly`` to the power ``e`` and sets ``res`` to the
    result.

.. function:: void _nmod_poly_pow_trunc_binexp(nn_ptr res, nn_srcptr poly, ulong e, slong trunc, nmod_t mod)

    Sets ``res`` to the low ``trunc`` coefficients of ``poly``
    (assumed to be zero padded if necessary to length ``trunc``) to
    the power ``e``. This is equivalent to doing a powering followed
    by a truncation. We require that ``res`` has enough space for
    ``trunc`` coefficients, that ``trunc > 0`` and that
    ``e > 1``. Aliasing is not permitted. Uses the binary
    exponentiation method.

.. function:: void nmod_poly_pow_trunc_binexp(nmod_poly_t res, const nmod_poly_t poly, ulong e, slong trunc)

    Sets ``res`` to the low ``trunc`` coefficients of ``poly``
    to the power ``e``. This is equivalent to doing a powering
    followed by a truncation. Uses the binary exponentiation method.

.. function:: void _nmod_poly_pow_trunc(nn_ptr res, nn_srcptr poly, ulong e, slong trunc, nmod_t mod)

    Sets ``res`` to the low ``trunc`` coefficients of ``poly``
    (assumed to be zero padded if necessary to length ``trunc``) to
    the power ``e``. This is equivalent to doing a powering followed
    by a truncation. We require that ``res`` has enough space for
    ``trunc`` coefficients, that ``trunc > 0`` and that
    ``e > 1``. Aliasing is not permitted.

.. function:: void nmod_poly_pow_trunc(nmod_poly_t res, const nmod_poly_t poly, ulong e, slong trunc)

    Sets ``res`` to the low ``trunc`` coefficients of ``poly``
    to the power ``e``. This is equivalent to doing a powering
    followed by a truncation.

.. function:: void _nmod_poly_powmod_ui_binexp(nn_ptr res, nn_srcptr poly, ulong e, nn_srcptr f, slong lenf, nmod_t mod)

    Sets ``res`` to ``poly`` raised to the power ``e``
    modulo ``f``, using binary exponentiation. We require ``e > 0``.

    We require ``lenf > 1``. It is assumed that ``poly`` is already
    reduced modulo ``f`` and zero-padded as necessary to have length
    exactly ``lenf - 1``. The output ``res`` must have room for
    ``lenf - 1`` coefficients.

.. function:: void nmod_poly_powmod_ui_binexp(nmod_poly_t res, const nmod_poly_t poly, ulong e, const nmod_poly_t f)

    Sets ``res`` to ``poly`` raised to the power ``e``
    modulo ``f``, using binary exponentiation. We require ``e >= 0``.

.. function:: void _nmod_poly_powmod_fmpz_binexp(nn_ptr res, nn_srcptr poly, fmpz_t e, nn_srcptr f, slong lenf, nmod_t mod)

    Sets ``res`` to ``poly`` raised to the power ``e``
    modulo ``f``, using binary exponentiation. We require ``e > 0``.

    We require ``lenf > 1``. It is assumed that ``poly`` is already
    reduced modulo ``f`` and zero-padded as necessary to have length
    exactly ``lenf - 1``. The output ``res`` must have room for ``lenf - 1`` coefficients.

.. function:: void nmod_poly_powmod_fmpz_binexp(nmod_poly_t res, const nmod_poly_t poly, fmpz_t e, const nmod_poly_t f)

    Sets ``res`` to ``poly`` raised to the power ``e``
    modulo ``f``, using binary exponentiation. We require ``e >= 0``.

.. function:: void _nmod_poly_powmod_ui_binexp_preinv (nn_ptr res, nn_srcptr poly, ulong e, nn_srcptr f, slong lenf, nn_srcptr finv, slong lenfinv, nmod_t mod)

    Sets ``res`` to ``poly`` raised to the power ``e``
    modulo ``f``, using binary exponentiation. We require ``e > 0``.
    We require ``finv`` to be the inverse of the reverse of ``f``.

    We require ``lenf > 1``. It is assumed that ``poly`` is already
    reduced modulo ``f`` and zero-padded as necessary to have length
    exactly ``lenf - 1``. The output ``res`` must have room for
    ``lenf - 1`` coefficients.

.. function:: void nmod_poly_powmod_ui_binexp_preinv(nmod_poly_t res, const nmod_poly_t poly, ulong e, const nmod_poly_t f, const nmod_poly_t finv)

    Sets ``res`` to ``poly`` raised to the power ``e``
    modulo ``f``, using binary exponentiation. We require ``e >= 0``.
    We require ``finv`` to be the inverse of the reverse of ``f``.

.. function:: void _nmod_poly_powmod_fmpz_binexp_preinv (nn_ptr res, nn_srcptr poly, fmpz_t e, nn_srcptr f, slong lenf, nn_srcptr finv, slong lenfinv, nmod_t mod)

    Sets ``res`` to ``poly`` raised to the power ``e``
    modulo ``f``, using binary exponentiation. We require ``e > 0``.
    We require ``finv`` to be the inverse of the reverse of ``f``.

    We require ``lenf > 1``. It is assumed that ``poly`` is already
    reduced modulo ``f`` and zero-padded as necessary to have length
    exactly ``lenf - 1``. The output ``res`` must have room for
    ``lenf - 1`` coefficients.

.. function:: void nmod_poly_powmod_fmpz_binexp_preinv(nmod_poly_t res, const nmod_poly_t poly, fmpz_t e, const nmod_poly_t f, const nmod_poly_t finv)

    Sets ``res`` to ``poly`` raised to the power ``e``
    modulo ``f``, using binary exponentiation. We require ``e >= 0``.
    We require ``finv`` to be the inverse of the reverse of ``f``.

.. function:: void _nmod_poly_powmod_x_ui_preinv (nn_ptr res, ulong e, nn_srcptr f, slong lenf, nn_srcptr finv, slong lenfinv, nmod_t mod)

    Sets ``res`` to ``x`` raised to the power ``e`` modulo ``f``,
    using sliding window exponentiation. We require ``e > 0``.
    We require ``finv`` to be the inverse of the reverse of ``f``.

    We require ``lenf > 2``. The output ``res`` must have room for
    ``lenf - 1`` coefficients.

.. function:: void nmod_poly_powmod_x_ui_preinv(nmod_poly_t res, ulong e, const nmod_poly_t f, const nmod_poly_t finv)

    Sets ``res`` to ``x`` raised to the power ``e``
    modulo ``f``, using sliding window exponentiation. We require
    ``e >= 0``. We require ``finv`` to be the inverse of the reverse of
    ``f``.

.. function:: void _nmod_poly_powmod_x_fmpz_preinv (nn_ptr res, fmpz_t e, nn_srcptr f, slong lenf, nn_srcptr finv, slong lenfinv, nmod_t mod)

    Sets ``res`` to ``x`` raised to the power ``e`` modulo ``f``,
    using sliding window exponentiation. We require ``e > 0``.
    We require ``finv`` to be the inverse of the reverse of ``f``.

    We require ``lenf > 2``. The output ``res`` must have room for
    ``lenf - 1`` coefficients.

.. function:: void nmod_poly_powmod_x_fmpz_preinv(nmod_poly_t res, fmpz_t e, const nmod_poly_t f, const nmod_poly_t finv)

    Sets ``res`` to ``x`` raised to the power ``e``
    modulo ``f``, using sliding window exponentiation. We require
    ``e >= 0``. We require ``finv`` to be the inverse of the reverse of
    ``f``.

.. function:: void _nmod_poly_powers_mod_preinv_naive(nn_ptr * res, nn_srcptr f, slong flen, slong n, nn_srcptr g, slong glen, nn_srcptr ginv, slong ginvlen, const nmod_t mod)

    Compute ``f^0, f^1, ..., f^(n-1) mod g``, where ``g`` has length ``glen``
    and ``f`` is reduced mod ``g`` and has length ``flen`` (possibly zero
    spaced). Assumes ``res`` is an array of ``n`` arrays each with space for
    at least ``glen - 1`` coefficients and that ``flen > 0``. We require that
    ``ginv`` of length ``ginvlen`` is set to the power series inverse of the
    reverse of ``g``.

.. function:: void nmod_poly_powers_mod_naive(nmod_poly_struct * res, const nmod_poly_t f, slong n, const nmod_poly_t g)

    Set the entries of the array ``res`` to ``f^0, f^1, ..., f^(n-1) mod g``.
    No aliasing is permitted between the entries of ``res`` and either of the
    inputs.

.. function:: void _nmod_poly_powers_mod_preinv_threaded_pool(nn_ptr * res, nn_srcptr f, slong flen, slong n, nn_srcptr g, slong glen, nn_srcptr ginv, slong ginvlen, const nmod_t mod, thread_pool_handle * threads, slong num_threads)

    Compute ``f^0, f^1, ..., f^(n-1) mod g``, where ``g`` has length ``glen``
    and ``f`` is reduced mod ``g`` and has length ``flen`` (possibly zero
    spaced). Assumes ``res`` is an array of ``n`` arrays each with space for
    at least ``glen - 1`` coefficients and that ``flen > 0``. We require that
    ``ginv`` of length ``ginvlen`` is set to the power series inverse of the
    reverse of ``g``.

.. function:: void _nmod_poly_powers_mod_preinv_threaded(nn_ptr * res, nn_srcptr f, slong flen, slong n, nn_srcptr g, slong glen, nn_srcptr ginv, slong ginvlen, const nmod_t mod)

    Compute ``f^0, f^1, ..., f^(n-1) mod g``, where ``g`` has length ``glen``
    and ``f`` is reduced mod ``g`` and has length ``flen`` (possibly zero
    spaced). Assumes ``res`` is an array of ``n`` arrays each with space for
    at least ``glen - 1`` coefficients and that ``flen > 0``. We require that
    ``ginv`` of length ``ginvlen`` is set to the power series inverse of the
    reverse of ``g``.

.. function:: void nmod_poly_powers_mod_bsgs(nmod_poly_struct * res, const nmod_poly_t f, slong n, const nmod_poly_t g)

    Set the entries of the array ``res`` to ``f^0, f^1, ..., f^(n-1) mod g``.
    No aliasing is permitted between the entries of ``res`` and either of the
    inputs.

Division
--------------------------------------------------------------------------------


.. function:: void _nmod_poly_divrem_basecase(nn_ptr Q, nn_ptr R, nn_srcptr A, slong A_len, nn_srcptr B, slong B_len, nmod_t mod)

    Finds `Q` and `R` such that `A = B Q + R` with `\operatorname{len}(R) < \operatorname{len}(B)`.
    If `\operatorname{len}(B) = 0` an exception is raised. We require that ``W``
    is temporary space of ``NMOD_DIVREM_BC_ITCH(A_len, B_len, mod)``
    coefficients.

.. function:: void nmod_poly_divrem_basecase(nmod_poly_t Q, nmod_poly_t R, const nmod_poly_t A, const nmod_poly_t B)

    Finds `Q` and `R` such that `A = B Q + R` with `\operatorname{len}(R) < \operatorname{len}(B)`.
    If `\operatorname{len}(B) = 0` an exception is raised.

.. function:: void _nmod_poly_divrem(nn_ptr Q, nn_ptr R, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nmod_t mod)

    Computes `Q` and `R` such that `A = BQ + R` with `\operatorname{len}(R)` less than
    ``lenB``, where ``A`` is of length ``lenA`` and ``B`` is of
    length ``lenB``. We require that ``Q`` have space for
    ``lenA - lenB + 1`` coefficients.

.. function:: void nmod_poly_divrem(nmod_poly_t Q, nmod_poly_t R, const nmod_poly_t A, const nmod_poly_t B)

    Computes `Q` and `R` such that `A = BQ + R` with `\operatorname{len}(R) < \operatorname{len}(B)`.

.. function:: void _nmod_poly_div(nn_ptr Q, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nmod_t mod)

    Notionally computes polynomials `Q` and `R` such that `A = BQ + R` with
    `\operatorname{len}(R)` less than ``lenB``, where ``A`` is of length ``lenA``
    and ``B`` is of length ``lenB``, but returns only ``Q``. We
    require that ``Q`` have space for ``lenA - lenB + 1`` coefficients.

.. function:: void nmod_poly_div(nmod_poly_t Q, const nmod_poly_t A, const nmod_poly_t B)

    Computes the quotient `Q` on polynomial division of `A` and `B`.

.. function:: void _nmod_poly_rem_q1(nn_ptr R, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nmod_t mod)

.. function:: void _nmod_poly_rem(nn_ptr R, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nmod_t mod)

    Computes the remainder `R` on polynomial division of `A` by `B`.

.. function:: void nmod_poly_rem(nmod_poly_t R, const nmod_poly_t A, const nmod_poly_t B)

    Computes the remainder `R` on polynomial division of `A` by `B`.

.. function:: void _nmod_poly_divexact(nn_ptr Q, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nmod_t mod)
              void nmod_poly_divexact(nmod_poly_t Q, const nmod_poly_t A, const nmod_poly_t B)

    Computes the quotient `Q` of `A` and `B` assuming that the division
    is exact.

.. function:: void _nmod_poly_inv_series_basecase(nn_ptr Qinv, nn_srcptr Q, slong Qlen, slong n, nmod_t mod)

    Given ``Q`` of length ``Qlen`` whose leading coefficient is invertible
    modulo the given modulus, finds a polynomial ``Qinv`` of length ``n``
    such that the top ``n`` coefficients of the product ``Q * Qinv`` is
    `x^{n - 1}`. Requires that ``n > 0``. This function can be viewed as
    inverting a power series.

.. function:: void nmod_poly_inv_series_basecase(nmod_poly_t Qinv, const nmod_poly_t Q, slong n)

    Given ``Q`` of length at least ``n`` find ``Qinv`` of length
    ``n`` such that the top ``n`` coefficients of the product
    ``Q * Qinv`` is `x^{n - 1}`. An exception is raised if ``n = 0``
    or if the length of ``Q`` is less than ``n``. The leading
    coefficient of ``Q`` must be invertible modulo the modulus of
    ``Q``. This function can be viewed as inverting a power series.

.. function:: void _nmod_poly_inv_series_newton(nn_ptr Qinv, nn_srcptr Q, slong Qlen, slong n, nmod_t mod)

    Given ``Q`` of length ``Qlen`` whose constant coefficient is invertible
    modulo the given modulus, find a polynomial ``Qinv`` of length ``n``
    such that ``Q * Qinv`` is ``1`` modulo `x^n`. Requires ``n > 0``.
    This function can be viewed as inverting a power series via Newton
    iteration.

.. function:: void nmod_poly_inv_series_newton(nmod_poly_t Qinv, const nmod_poly_t Q, slong n)

    Given ``Q`` find ``Qinv`` such that ``Q * Qinv`` is ``1``
    modulo `x^n`. The constant coefficient of ``Q`` must be invertible
    modulo the modulus of ``Q``. An exception is raised if this is not
    the case or if ``n = 0``. This function can be viewed as inverting
    a power series via Newton iteration.

.. function:: void _nmod_poly_inv_series(nn_ptr Qinv, nn_srcptr Q, slong Qlen, slong n, nmod_t mod)

    Given ``Q`` of length ``Qlenn`` whose constant coefficient is invertible
    modulo the given modulus, find a polynomial ``Qinv`` of length ``n``
    such that ``Q * Qinv`` is ``1`` modulo `x^n`. Requires ``n > 0``.
    This function can be viewed as inverting a power series.

.. function:: void nmod_poly_inv_series(nmod_poly_t Qinv, const nmod_poly_t Q, slong n)

    Given ``Q`` find ``Qinv`` such that ``Q * Qinv`` is ``1``
    modulo `x^n`. The constant coefficient of ``Q`` must be invertible
    modulo the modulus of ``Q``. An exception is raised if this is not
    the case or if ``n = 0``. This function can be viewed as inverting
    a power series.

.. function:: void _nmod_poly_div_series_basecase(nn_ptr Q, nn_srcptr A, slong Alen, nn_srcptr B, slong Blen, slong n, nmod_t mod)

    Given polynomials ``A`` and ``B`` of length ``Alen`` and
    ``Blen``, finds the
    polynomial ``Q`` of length ``n`` such that ``Q * B = A``
    modulo `x^n`. We assume ``n > 0`` and that the constant coefficient
    of ``B`` is invertible modulo the given modulus. The polynomial
    ``Q`` must have space for ``n`` coefficients.

.. function:: void nmod_poly_div_series_basecase(nmod_poly_t Q, const nmod_poly_t A, const nmod_poly_t B, slong n)

    Given polynomials ``A`` and ``B`` considered modulo ``n``,
    finds the polynomial ``Q`` of length at most ``n`` such that
    ``Q * B = A`` modulo `x^n`. We assume ``n > 0`` and that the
    constant coefficient of ``B`` is invertible modulo the modulus.
    An exception is raised if ``n == 0`` or the constant coefficient
    of ``B`` is zero.

.. function:: void _nmod_poly_div_series(nn_ptr Q, nn_srcptr A, slong Alen, nn_srcptr B, slong Blen, slong n, nmod_t mod)

    Given polynomials ``A`` and ``B`` of length ``Alen`` and
    ``Blen``, finds the
    polynomial ``Q`` of length ``n`` such that ``Q * B = A``
    modulo `x^n`. We assume ``n > 0`` and that the constant coefficient
    of ``B`` is invertible modulo the given modulus. The polynomial
    ``Q`` must have space for ``n`` coefficients.

.. function:: void nmod_poly_div_series(nmod_poly_t Q, const nmod_poly_t A, const nmod_poly_t B, slong n)

    Given polynomials ``A`` and ``B`` considered modulo ``n``,
    finds the polynomial ``Q`` of length at most ``n`` such that
    ``Q * B = A`` modulo `x^n`. We assume ``n > 0`` and that the
    constant coefficient of ``B`` is invertible modulo the modulus.
    An exception is raised if ``n == 0`` or the constant coefficient
    of ``B`` is zero.

.. function:: void _nmod_poly_div_newton_n_preinv (nn_ptr Q, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nn_srcptr Binv, slong lenBinv, nmod_t mod)

    Notionally computes polynomials `Q` and `R` such that `A = BQ + R` with
    `\operatorname{len}(R)` less than ``lenB``, where ``A`` is of length ``lenA``
    and ``B`` is of length ``lenB``, but return only `Q`.

    We require that `Q` have space for ``lenA - lenB + 1`` coefficients
    and assume that the leading coefficient of `B` is a unit. Furthermore, we
    assume that `Binv` is the inverse of the reverse of `B` mod `x^{\operatorname{len}(B)}`.

    The algorithm used is to reverse the polynomials and divide the
    resulting power series, then reverse the result.

.. function:: void nmod_poly_div_newton_n_preinv (nmod_poly_t Q, const nmod_poly_t A, const nmod_poly_t B, const nmod_poly_t Binv)

    Notionally computes `Q` and `R` such that `A = BQ + R` with
    `\operatorname{len}(R) < \operatorname{len}(B)`, but returns only `Q`.

    We assume that the leading coefficient of `B` is a unit and that `Binv` is
    the inverse of the reverse of `B` mod `x^{\operatorname{len}(B)}`.

    It is required that the length of `A` is less than or equal to
    2*the length of `B` - 2.

    The algorithm used is to reverse the polynomials and divide the
    resulting power series, then reverse the result.

.. function:: void _nmod_poly_divrem_newton_n_preinv (nn_ptr Q, nn_ptr R, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nn_srcptr Binv, slong lenBinv, nmod_t mod)

    Computes `Q` and `R` such that `A = BQ + R` with `\operatorname{len}(R)` less than
    ``lenB``, where `A` is of length ``lenA`` and `B` is of length
    ``lenB``. We require that `Q` have space for ``lenA - lenB + 1``
    coefficients. Furthermore, we assume that `Binv` is the inverse of the
    reverse of `B` mod `x^{\operatorname{len}(B)}`. The algorithm used is to call
    :func:`div_newton_n_preinv` and then multiply out and compute
    the remainder.

.. function:: void nmod_poly_divrem_newton_n_preinv(nmod_poly_t Q, nmod_poly_t R, const nmod_poly_t A, const nmod_poly_t B, const nmod_poly_t Binv)

    Computes `Q` and `R` such that `A = BQ + R` with `\operatorname{len}(R) < \operatorname{len}(B)`.
    We assume `Binv` is the inverse of the reverse of `B` mod `x^{\operatorname{len}(B)}`.

    It is required that the length of `A` is less than or equal to
    2*the length of `B` - 2.

    The algorithm used is to call :func:`div_newton_n` and then multiply out
    and compute the remainder.

.. function:: ulong _nmod_poly_div_root(nn_ptr Q, nn_srcptr A, slong len, ulong c, nmod_t mod)

    Sets ``(Q, len-1)`` to the quotient of ``(A, len)`` on division
    by `(x - c)`, and returns the remainder, equal to the value of `A`
    evaluated at `c`. `A` and `Q` are allowed to be the same, but may
    not overlap partially in any other way.

.. function:: ulong nmod_poly_div_root(nmod_poly_t Q, const nmod_poly_t A, ulong c)

    Sets `Q` to the quotient of `A` on division by `(x - c)`, and returns
    the remainder, equal to the value of `A` evaluated at `c`.


Divisibility testing
--------------------------------------------------------------------------------


.. function:: int _nmod_poly_divides_classical(nn_ptr Q, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nmod_t mod)

    Returns `1` if `(B, lenB)` divides `(A, lenA)` and sets
    `(Q, lenA - lenB + 1)` to the quotient. Otherwise, returns `0` and sets
    `(Q, lenA - lenB + 1)` to zero. We require that `lenA >= lenB > 0`.

.. function:: int nmod_poly_divides_classical(nmod_poly_t Q, const nmod_poly_t A, const nmod_poly_t B)

    Returns `1` if `B` divides `A` and sets `Q` to the quotient. Otherwise
    returns `0` and sets `Q` to zero.

.. function:: int _nmod_poly_divides(nn_ptr Q, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nmod_t mod)

    Returns `1` if `(B, lenB)` divides `(A, lenA)` and sets
    `(Q, lenA - lenB + 1)` to the quotient. Otherwise, returns `0` and sets
    `(Q, lenA - lenB + 1)` to zero. We require that `lenA >= lenB > 0`.

.. function:: int nmod_poly_divides(nmod_poly_t Q, const nmod_poly_t A, const nmod_poly_t B)

    Returns `1` if `B` divides `A` and sets `Q` to the quotient. Otherwise
    returns `0` and sets `Q` to zero.

.. function:: ulong nmod_poly_remove(nmod_poly_t f, const nmod_poly_t p)

    Removes the highest possible power of ``p`` from ``f`` and
    returns the exponent.


Derivative and integral
--------------------------------------------------------------------------------


.. function:: void _nmod_poly_derivative(nn_ptr x_prime, nn_srcptr x, slong len, nmod_t mod)

    Sets the first ``len - 1`` coefficients of ``x_prime`` to the
    derivative of ``x`` which is assumed to be of length ``len``.
    It is assumed that ``len > 0``.

.. function:: void nmod_poly_derivative(nmod_poly_t x_prime, const nmod_poly_t x)

    Sets ``x_prime`` to the derivative of ``x``.

.. function:: void _nmod_poly_integral(nn_ptr x_int, nn_srcptr x, slong len, nmod_t mod)

    Set the first ``len`` coefficients of ``x_int`` to the
    integral of ``x`` which is assumed to be of length ``len - 1``.
    The constant term of ``x_int`` is set to zero.
    It is assumed that ``len > 0``. The result is only well-defined
    if the modulus is a prime number strictly larger than the degree of
    ``x``. Supports aliasing between the two polynomials.

.. function:: void nmod_poly_integral(nmod_poly_t x_int, const nmod_poly_t x)

    Set ``x_int`` to the indefinite integral of ``x`` with constant
    term zero. The result is only well-defined if the modulus
    is a prime number strictly larger than the degree of ``x``.



Evaluation
--------------------------------------------------------------------------------


.. function:: ulong _nmod_poly_evaluate_nmod_precomp(nn_srcptr poly, slong len, ulong c, ulong c_precomp, ulong n)

    Evaluates ``poly`` at the value ``c`` and reduces modulo the given modulus
    ``modn``. The value ``c`` should be reduced modulo the modulus, and the
    modulus must be less than `2^{\mathtt{FLINT\_BITS} - 1}`. The algorithm
    used is Horner's method, with multiplications done via
    :func:`n_mulmod_shoup` using the precomputed ``c_precomp`` obtained via
    :func:`n_mulmod_precomp_shoup`.

.. function:: ulong _nmod_poly_evaluate_nmod_precomp_lazy(nn_srcptr poly, slong len, ulong c, ulong c_precomp, ulong n)

    Evaluates ``poly`` at the value ``c`` modulo ``n``, with lazy reductions
    modulo `n`. Precisely, if all coefficients of ``poly`` are less than `m`,
    the input requirement is `m \le 2^{\mathtt{FLINT\_BITS}} - 2n + 1`, and the
    output value is in `[0, m+2n-1)` and equal to the sought evaluation modulo
    `n`. In particular the coefficients of ``poly`` need not be reduced modulo
    ``n``, and the output may not be either. However, the value ``c`` should be
    reduced modulo `n`.

    In the case where `m = n` (coefficients of ``poly`` are reduced modulo
    `n`), then the above leads to the requirement `3n-1 \le
    2^{\mathtt{FLINT\_BITS}}` (this is `n \le 6148914691236517205` for 64 bits,
    and `n \le 1431655765` for 32 bits), and reducing the output just amounts
    to subtracting `n` or `2n`. The algorithm used is Horner's method, with
    multiplications done as in :func:`n_mulmod_shoup` using the precomputed
    ``c_precomp`` obtained via :func:`n_mulmod_precomp_shoup`.

.. function:: ulong _nmod_poly_evaluate_nmod(nn_srcptr poly, slong len, ulong c, nmod_t mod)

    Evaluates ``poly`` at the value ``c`` and reduces modulo the given modulus
    of ``poly``. The value ``c`` should be reduced modulo the modulus. The
    algorithm used is Horner's method, with multiplications done via
    :func:`nmod_mul`.

.. function:: ulong nmod_poly_evaluate_nmod(const nmod_poly_t poly, ulong c)

    Evaluates ``poly`` at the value ``c`` and reduces modulo the modulus of
    ``poly``. The value ``c`` should be reduced modulo the modulus. The
    algorithm used is Horner's method, with multiplications and additions done
    differently depending on the modulus ``poly->mod`` and on the degree (calls
    one of :func:`_nmod_poly_evaluate_nmod`,
    :func:`_nmod_poly_evaluate_nmod_precomp`,
    :func:`_nmod_poly_evaluate_nmod_precomp_lazy`).

.. function:: void nmod_poly_evaluate_mat_horner(nmod_mat_t dest, const nmod_poly_t poly, const nmod_mat_t c)

    Evaluates ``poly`` with matrix as an argument at the value ``c``
    and stores the result in ``dest``. The dimension and modulus of
    ``dest`` is assumed to be same as that of ``c``. ``dest`` and
    ``c`` may be aliased. Horner's Method is used to compute the result.

.. function:: void nmod_poly_evaluate_mat_paterson_stockmeyer(nmod_mat_t dest, const nmod_poly_t poly, const nmod_mat_t c)

    Evaluates ``poly`` with matrix as an argument at the value ``c``
    and stores the result in ``dest``. The dimension and modulus of
    ``dest`` is assumed to be same as that of ``c``. ``dest`` and
    ``c`` may be aliased. Paterson-Stockmeyer algorithm is used to compute
    the result. The algorithm is described in [Paterson1973]_.

.. function:: void nmod_poly_evaluate_mat(nmod_mat_t dest, const nmod_poly_t poly, const nmod_mat_t c)

    Evaluates ``poly`` with matrix as an argument at the value ``c``
    and stores the result in ``dest``. The dimension and modulus of
    ``dest`` is assumed to be same as that of ``c``. ``dest`` and
    ``c`` may be aliased. This function automatically switches between
    Horner's method and the Paterson-Stockmeyer algorithm.


Multipoint evaluation
--------------------------------------------------------------------------------


.. function:: void _nmod_poly_evaluate_nmod_vec_iter(nn_ptr ys, nn_srcptr poly, slong len, nn_srcptr xs, slong n, nmod_t mod)

    Evaluates (``coeffs``, ``len``) at the ``n`` values
    given in the vector ``xs``, writing the output values
    to ``ys``. The values in ``xs`` should be reduced
    modulo the modulus.

    Uses Horner's method iteratively.

.. function:: void nmod_poly_evaluate_nmod_vec_iter(nn_ptr ys, const nmod_poly_t poly, nn_srcptr xs, slong n)

    Evaluates ``poly`` at the ``n`` values given in the vector
    ``xs``, writing the output values to ``ys``. The values in
    ``xs`` should be reduced modulo the modulus.

    Uses Horner's method iteratively.

.. function:: void _nmod_poly_evaluate_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr poly, slong plen, const nn_ptr * tree, slong len, nmod_t mod)

    Evaluates (``poly``, ``plen``) at the ``len`` values given
    by the precomputed subproduct tree ``tree``.

.. function:: void _nmod_poly_evaluate_nmod_vec_fast(nn_ptr ys, nn_srcptr poly, slong len, nn_srcptr xs, slong n, nmod_t mod)

    Evaluates (``coeffs``, ``len``) at the ``n`` values
    given in the vector ``xs``, writing the output values
    to ``ys``. The values in ``xs`` should be reduced
    modulo the modulus.

    Uses fast multipoint evaluation, building a temporary subproduct tree.

.. function:: void nmod_poly_evaluate_nmod_vec_fast(nn_ptr ys, const nmod_poly_t poly, nn_srcptr xs, slong n)

    Evaluates ``poly`` at the ``n`` values given in the vector
    ``xs``, writing the output values to ``ys``. The values in
    ``xs`` should be reduced modulo the modulus.

    Uses fast multipoint evaluation, building a temporary subproduct tree.

.. function:: void _nmod_poly_evaluate_nmod_vec(nn_ptr ys, nn_srcptr poly, slong len, nn_srcptr xs, slong n, nmod_t mod)

    Evaluates (``poly``, ``len``) at the ``n`` values
    given in the vector ``xs``, writing the output values
    to ``ys``. The values in ``xs`` should be reduced
    modulo the modulus.

.. function:: void nmod_poly_evaluate_nmod_vec(nn_ptr ys, const nmod_poly_t poly, nn_srcptr xs, slong n)

    Evaluates ``poly`` at the ``n`` values given in the vector
    ``xs``, writing the output values to ``ys``. The values in
    ``xs`` should be reduced modulo the modulus.

.. function:: void _nmod_poly_evaluate_geometric_nmod_vec_iter(nn_ptr ys, nn_srcptr coeffs, slong len, ulong r, slong n, nmod_t mod)

    Evaluates (``coeffs``, ``len``) at the first ``n`` powers
    of the square of ``r``, writing the output values
    to ``ys``. The value of ``r`` should be reduced
    modulo the modulus.

    Uses Horner's method iteratively.

.. function:: void nmod_poly_evaluate_geometric_nmod_vec_iter(nn_ptr ys, const nmod_poly_t poly, ulong r, slong n)

    Evaluates ``poly`` at the first ``n`` powers
    of the square of ``r``, writing the output values
    to ``ys``. The value of ``r`` should be reduced
    modulo the modulus.

    Uses Horner's method iteratively.

.. function:: void _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr poly, slong plen, const nmod_geometric_progression_t G, slong len)

    Evaluates (``poly``, ``plen``) at the ``len`` values given
    by the precomputed geometric progression ``G``. The value of
    ``len`` should be less than or equal to the precomputation size parameter ``G->len``.

.. function:: void _nmod_poly_evaluate_geometric_nmod_vec_fast(nn_ptr ys, nn_srcptr coeffs, slong len, ulong r, slong n, nmod_t mod)

    Evaluates (``coeffs``, ``len``) at the first ``n`` powers
    of the square of ``r``, writing the output values to ``ys``. 
    The value of ``r`` should be reduced modulo the modulus ``mod``
    and of sufficient multiplicative order such that none of 
    the first `n` powers of `r^2` is one.

    Uses fast geometric multipoint evaluation, building a temporary geometric progression precomputation.

.. function:: void nmod_poly_evaluate_geometric_nmod_vec_fast(nn_ptr ys, const nmod_poly_t poly, ulong r, slong n)

    Evaluates ``poly``  at the first ``n`` powers
    of the square of ``r``, writing the output values to ``ys``. 
    The value of ``r`` should be reduced modulo the modulus of the polynomial
    and of sufficient multiplicative order such that none of 
    the first `n` powers of `r^2` is one.

    Uses fast geometric multipoint evaluation, building a temporary geometric progression precomputation.

Interpolation
--------------------------------------------------------------------------------


.. function:: void _nmod_poly_interpolate_nmod_vec(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong n, nmod_t mod)

    Sets ``poly`` to the unique polynomial of length at most ``n``
    that interpolates the ``n`` given evaluation points ``xs`` and
    values ``ys``. If the interpolating polynomial is shorter than
    length ``n``, the leading coefficients are set to zero.

    The values in ``xs`` and ``ys`` should be reduced modulo the
    modulus, and all ``xs`` must be distinct. Aliasing between
    ``poly`` and ``xs`` or ``ys`` is not allowed.

.. function:: void nmod_poly_interpolate_nmod_vec(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong n)

    Sets ``poly`` to the unique polynomial of length ``n`` that
    interpolates the ``n`` given evaluation points ``xs`` and
    values ``ys``. The values in ``xs`` and ``ys`` should be
    reduced modulo the modulus, and all ``xs`` must be distinct.

.. function:: void _nmod_poly_interpolation_weights(nn_ptr w, const nn_ptr * tree, slong len, nmod_t mod)

    Sets ``w`` to the barycentric interpolation weights for fast
    Lagrange interpolation with respect to a given subproduct tree.

.. function:: void _nmod_poly_interpolate_nmod_vec_fast_precomp(nn_ptr poly, nn_srcptr ys, const nn_ptr * tree, nn_srcptr weights, slong len, nmod_t mod)

    Performs interpolation using the fast Lagrange interpolation
    algorithm, generating a temporary subproduct tree.

    The function values are given as ``ys``. The function takes
    a precomputed subproduct tree ``tree`` and barycentric
    interpolation weights ``weights`` corresponding to the
    roots.

.. function:: void _nmod_poly_interpolate_nmod_vec_fast(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong n, nmod_t mod)

    Performs interpolation using the fast Lagrange interpolation
    algorithm, generating a temporary subproduct tree.

.. function:: void nmod_poly_interpolate_nmod_vec_fast(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong n)

    Performs interpolation using the fast Lagrange interpolation algorithm,
    generating a temporary subproduct tree.

.. function:: void _nmod_poly_interpolate_nmod_vec_newton(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong n, nmod_t mod)

    Forms the interpolating polynomial in the Newton basis using
    the method of divided differences and then converts it to
    monomial form.

.. function:: void nmod_poly_interpolate_nmod_vec_newton(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong n)

    Forms the interpolating polynomial in the Newton basis using
    the method of divided differences and then converts it to
    monomial form.

.. function:: void _nmod_poly_interpolate_nmod_vec_barycentric(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong n, nmod_t mod)

    Forms the interpolating polynomial using a naive implementation
    of the barycentric form of Lagrange interpolation.

.. function:: void nmod_poly_interpolate_nmod_vec_barycentric(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong n)

    Forms the interpolating polynomial using a naive implementation
    of the barycentric form of Lagrange interpolation.

.. function:: void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, nn_srcptr v, const nmod_geometric_progression_t G, slong len, nmod_t mod)
              void nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nmod_poly_t poly, nn_srcptr v, const nmod_geometric_progression_t G, slong len)

    Performs interpolation using the geometric progression precomputation ``G``.

    Sets ``poly`` to the unique polynomial of length at most ``len``
    that interpolates according to the parameter set of ``G``.
    The value of ``len`` should be equal to the precomputation size parameter ``G->len``.

    Uses fast geometric multipoint interpolation using a supplied geometric progression precomputation.

.. function:: void nmod_poly_interpolate_geometric_nmod_vec_fast(nmod_poly_t poly, ulong r, nn_srcptr ys, slong n)

    Sets ``poly`` to the unique polynomial of length at most ``n``
    that interpolates the first ``n`` powers of ``r`` and
    values ``ys``.

    The values ``ys`` and ``r`` should be reduced modulo the
    modulus, and all ``r`` should be of sufficient order such that
    none of the first `n` powers of `r^2` is one. Aliasing between
    ``poly`` and ``ys`` is not allowed.

    Uses fast geometric multipoint interpolation, building a temporary geometric progression precomputation.

    
Composition
--------------------------------------------------------------------------------


.. function:: void _nmod_poly_compose_horner(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, nmod_t mod)

    Composes ``poly1`` of length ``len1`` with ``poly2`` of length
    ``len2`` and sets ``res`` to the result, i.e. evaluates
    ``poly1`` at ``poly2``. The algorithm used is Horner's algorithm.
    We require that ``res`` have space for ``(len1 - 1)*(len2 - 1) + 1``
    coefficients. It is assumed that ``len1 > 0`` and ``len2 > 0``.

.. function:: void nmod_poly_compose_horner(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2)

    Composes ``poly1`` with ``poly2`` and sets ``res`` to the result,
    i.e. evaluates ``poly1`` at ``poly2``. The algorithm used is
    Horner's algorithm.

.. function:: void _nmod_poly_compose_divconquer(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, nmod_t mod)

    Composes ``poly1`` of length ``len1`` with ``poly2`` of length
    ``len2`` and sets ``res`` to the result, i.e. evaluates
    ``poly1`` at ``poly2``. The algorithm used is the divide and
    conquer algorithm. We require that ``res`` have space for
    ``(len1 - 1)*(len2 - 1) + 1`` coefficients. It is assumed that
    ``len1 > 0`` and ``len2 > 0``.

.. function:: void nmod_poly_compose_divconquer(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2)

    Composes ``poly1`` with ``poly2`` and sets ``res`` to the result,
    i.e. evaluates ``poly1`` at ``poly2``. The algorithm used is
    the divide and conquer algorithm.

.. function:: void _nmod_poly_compose(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, nmod_t mod)

    Composes ``poly1`` of length ``len1`` with ``poly2`` of length
    ``len2`` and sets ``res`` to the result, i.e. evaluates ``poly1``
    at ``poly2``. We require that ``res`` have space for
    ``(len1 - 1)*(len2 - 1) + 1`` coefficients. It is assumed that
    ``len1 > 0`` and ``len2 > 0``.

.. function:: void nmod_poly_compose(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2)

    Composes ``poly1`` with ``poly2`` and sets ``res`` to the result,
    that is, evaluates ``poly1`` at ``poly2``.


Taylor shift
--------------------------------------------------------------------------------


.. function:: void _nmod_poly_taylor_shift_horner(nn_ptr poly, ulong c, slong len, nmod_t mod)

    Performs the Taylor shift composing ``poly`` by `x+c` in-place.
    Uses an efficient version Horner's rule.

.. function:: void nmod_poly_taylor_shift_horner(nmod_poly_t g, const nmod_poly_t f, ulong c)

    Performs the Taylor shift composing ``f`` by `x+c`.

.. function:: void _nmod_poly_taylor_shift_convolution(nn_ptr poly, ulong c, slong len, nmod_t mod)

    Performs the Taylor shift composing ``poly`` by `x+c` in-place.
    Writes the composition as a single convolution with cost `O(M(n))`.
    We require that the modulus is a prime at least as large as the length.

.. function:: void nmod_poly_taylor_shift_convolution(nmod_poly_t g, const nmod_poly_t f, ulong c)

    Performs the Taylor shift composing ``f`` by `x+c`.
    Writes the composition as a single convolution with cost `O(M(n))`.
    We require that the modulus is a prime at least as large as the length.

.. function:: void _nmod_poly_taylor_shift(nn_ptr poly, ulong c, slong len, nmod_t mod)

    Performs the Taylor shift composing ``poly`` by `x+c` in-place.
    We require that the modulus is a prime.

.. function:: void nmod_poly_taylor_shift(nmod_poly_t g, const nmod_poly_t f, ulong c)

    Performs the Taylor shift composing ``f`` by `x+c`.
    We require that the modulus is a prime.


Modular composition
--------------------------------------------------------------------------------


.. function:: void _nmod_poly_compose_mod_horner(nn_ptr res, nn_srcptr f, slong lenf, nn_srcptr g, nn_srcptr h, slong lenh, nmod_t mod)

    Sets ``res`` to the composition `f(g)` modulo `h`. We require that
    `h` is nonzero and that the length of `g` is one less than the
    length of `h` (possibly with zero padding). The output is not allowed
    to be aliased with any of the inputs.

    The algorithm used is Horner's rule.

.. function:: void nmod_poly_compose_mod_horner(nmod_poly_t res, const nmod_poly_t f, const nmod_poly_t g, const nmod_poly_t h)

    Sets ``res`` to the composition `f(g)` modulo `h`. We require that
    `h` is nonzero. The algorithm used is Horner's rule.

.. function:: void _nmod_poly_compose_mod_brent_kung(nn_ptr res, nn_srcptr f, slong lenf, nn_srcptr g, nn_srcptr h, slong lenh, nmod_t mod)

    Sets ``res`` to the composition `f(g)` modulo `h`. We require that
    `h` is nonzero and that the length of `g` is one less than the
    length of `h` (possibly with zero padding). We also require that
    the length of `f` is less than the length of `h`. The output is not allowed
    to be aliased with any of the inputs.

    The algorithm used is the Brent-Kung matrix algorithm.

.. function:: void nmod_poly_compose_mod_brent_kung(nmod_poly_t res, const nmod_poly_t f, const nmod_poly_t g, const nmod_poly_t h)

    Sets ``res`` to the composition `f(g)` modulo `h`. We require that
    `h` is nonzero and that `f` has smaller degree than `h`.
    The algorithm used is the Brent-Kung matrix algorithm.

.. function:: void _nmod_poly_compose_mod_brent_kung_preinv(nn_ptr res, nn_srcptr f, slong lenf, nn_srcptr g, nn_srcptr h, slong lenh, nn_srcptr hinv, slong lenhinv, nmod_t mod)

    Sets ``res`` to the composition `f(g)` modulo `h`. We require that
    `h` is nonzero and that the length of `g` is one less than the
    length of `h` (possibly with zero padding). We also require that
    the length of `f` is less than the length of `h`. Furthermore, we require
    ``hinv`` to be the inverse of the reverse of ``h``.
    The output is not allowed to be aliased with any of the inputs.

    The algorithm used is the Brent-Kung matrix algorithm.

.. function:: void nmod_poly_compose_mod_brent_kung_preinv(nmod_poly_t res, const nmod_poly_t f, const nmod_poly_t g, const nmod_poly_t h, const nmod_poly_t hinv)

    Sets ``res`` to the composition `f(g)` modulo `h`. We require that
    `h` is nonzero and that `f` has smaller degree than `h`. Furthermore,
    we require ``hinv`` to be the inverse of the reverse of ``h``.
    The algorithm used is the Brent-Kung matrix algorithm.

.. function:: void _nmod_poly_mod_matrix_rows_evaluate(nn_ptr res, const nmod_mat_t A, nn_srcptr h, slong n, nn_srcptr poly3, slong len3, nn_srcptr poly3inv, slong len3inv, nmod_t mod)

    Set ``res`` to the polynomial `\sum_{i=0}^{d-1} A_i h^i` modulo ``poly3`` where `A_i`
    denotes the `i`-th row of the `d \times n` matrix ``A``
    and ``res`` and ``h`` have length ``n``. The length of ``poly3`` must
    be equal to ``n + 1``.

.. function:: void _nmod_poly_reduce_matrix_mod_poly (nmod_mat_t A, const nmod_mat_t B, const nmod_poly_t f)

    Sets the ith row of ``A`` to the reduction of the ith row of `B` modulo
    `f` for `i=1,\ldots,\sqrt{\deg(f)}`. We require `B` to be at least
    a `\sqrt{\deg(f)}\times \deg(f)` matrix and `f` to be nonzero.

.. function:: void _nmod_poly_precompute_matrix_worker (void * arg_ptr)

    Worker function version of ``_nmod_poly_precompute_matrix``.
    Input/output is stored in ``nmod_poly_matrix_precompute_arg_t``.

.. function:: void _nmod_poly_precompute_matrix (nmod_mat_t A, nn_srcptr f, nn_srcptr g, slong leng, nn_srcptr ginv, slong lenginv, nmod_t mod)

    Sets the ith row of ``A`` to `f^i` modulo `g` for
    `i=1,\ldots,\sqrt{\deg(g)}`. We require `A` to be
    a `\sqrt{\deg(g)}\times \deg(g)` matrix. We require
    ``ginv`` to be the inverse of the reverse of ``g`` and `g` to be
    nonzero. ``f`` has to be reduced modulo ``g`` and of length one less
    than ``leng`` (possibly with zero padding).

.. function:: void nmod_poly_precompute_matrix (nmod_mat_t A, const nmod_poly_t f, const nmod_poly_t g, const nmod_poly_t ginv)

    Sets the ith row of ``A`` to `f^i` modulo `g` for
    `i=1,\ldots,\sqrt{\deg(g)}`. We require `A` to be
    a `\sqrt{\deg(g)}\times \deg(g)` matrix. We require
    ``ginv`` to be the inverse of the reverse of ``g``.

.. function:: void _nmod_poly_compose_mod_brent_kung_precomp_preinv_worker(void * arg_ptr)

    Worker function version of
    ``_nmod_poly_compose_mod_brent_kung_precomp_preinv``.
    Input/output is stored in
    ``nmod_poly_compose_mod_precomp_preinv_arg_t``.

.. function:: void _nmod_poly_compose_mod_brent_kung_precomp_preinv(nn_ptr res, nn_srcptr f, slong lenf, const nmod_mat_t A, nn_srcptr h, slong lenh, nn_srcptr hinv, slong lenhinv, nmod_t mod)

    Sets ``res`` to the composition `f(g)` modulo `h`. We require that
    `h` is nonzero. We require that the ith row of `A` contains `g^i` for
    `i=1,\ldots,\sqrt{\deg(h)}`, i.e. `A` is a
    `\sqrt{\deg(h)}\times \deg(h)` matrix. We also require that
    the length of `f` is less than the length of `h`. Furthermore, we require
    ``hinv`` to be the inverse of the reverse of ``h``.
    The output is not allowed to be aliased with any of the inputs.

    The algorithm used is the Brent-Kung matrix algorithm.

.. function:: void nmod_poly_compose_mod_brent_kung_precomp_preinv(nmod_poly_t res, const nmod_poly_t f, const nmod_mat_t A, const nmod_poly_t h, const nmod_poly_t hinv)

    Sets ``res`` to the composition `f(g)` modulo `h`. We require that the
    ith row of `A` contains `g^i` for `i=1,\ldots,\sqrt{\deg(h)}`, i.e. `A` is a
    `\sqrt{\deg(h)}\times \deg(h)` matrix. We require that `h` is nonzero and
    that `f` has smaller degree than `h`. Furthermore, we require ``hinv`` to
    be the inverse of the reverse of ``h``. This version of Brent-Kung
    modular composition is particularly useful if one has to perform several
    modular composition of the form `f(g)` modulo `h` for fixed `g` and `h`.

.. function:: void _nmod_poly_compose_mod_brent_kung_vec_preinv(nmod_poly_struct * res, const nmod_poly_struct * polys, slong len1, slong l, nn_srcptr g, slong leng, nn_srcptr h, slong lenh, nn_srcptr hinv, slong lenhinv, nmod_t mod)

    Sets ``res`` to the composition `f_i(g)` modulo `h` for `1\leq i \leq l`,
    where `f_i` are the first ``l`` elements of ``polys``. We require that `h`
    is nonzero and that the length of `g` is less than the length of `h`. We
    also require that the length of `f_i` is less than the length of `h`. We
    require ``res`` to have enough memory allocated to hold ``l``
    ``nmod_poly_struct``'s. The entries of ``res`` need to be initialised and
    ``l`` needs to be less than ``len1`` Furthermore, we require ``hinv`` to
    be the inverse of the reverse of ``h``. The output is not allowed to be
    aliased with any of the inputs.

    The algorithm used is the Brent-Kung matrix algorithm.

.. function:: void nmod_poly_compose_mod_brent_kung_vec_preinv(nmod_poly_struct * res, const nmod_poly_struct * polys, slong len1, slong n, const nmod_poly_t g, const nmod_poly_t h, const nmod_poly_t hinv)

    Sets ``res`` to the composition `f_i(g)` modulo `h` for `1\leq i \leq n`
    where `f_i` are the first ``n`` elements of ``polys``. We require ``res``
    to have enough memory allocated to hold ``n`` ``nmod_poly_struct``. The
    entries of ``res`` need to be initialised and ``n`` needs to be less than
    ``len1``. We require that `h` is nonzero and that `f_i` and `g` have
    smaller degree than `h`. Furthermore, we require ``hinv`` to be the inverse
    of the reverse of ``h``. No aliasing of ``res`` and ``polys`` is allowed.
    The algorithm used is the Brent-Kung matrix algorithm.

.. function:: void _nmod_poly_compose_mod_brent_kung_vec_preinv_threaded_pool(nmod_poly_struct * res, const nmod_poly_struct * polys, slong lenpolys, slong l, nn_srcptr g, slong glen, nn_srcptr poly, slong len, nn_srcptr polyinv, slong leninv, nmod_t mod, thread_pool_handle * threads, slong num_threads)

    Multithreaded version of
    :func:`_nmod_poly_compose_mod_brent_kung_vec_preinv`. Distributing the
    Horner evaluations across :func:`flint_get_num_threads` threads.

.. function:: void nmod_poly_compose_mod_brent_kung_vec_preinv_threaded_pool(nmod_poly_struct * res, const nmod_poly_struct * polys, slong len1, slong n, const nmod_poly_t g, const nmod_poly_t poly, const nmod_poly_t polyinv, thread_pool_handle * threads, slong num_threads)

    Multithreaded version of
    :func:`nmod_poly_compose_mod_brent_kung_vec_preinv`. Distributing the
    Horner evaluations across :func:`flint_get_num_threads` threads.

.. function:: void nmod_poly_compose_mod_brent_kung_vec_preinv_threaded(nmod_poly_struct * res, const nmod_poly_struct * polys, slong len1, slong n, const nmod_poly_t g, const nmod_poly_t poly, const nmod_poly_t polyinv)

    Multithreaded version of
    :func:`nmod_poly_compose_mod_brent_kung_vec_preinv`. Distributing the
    Horner evaluations across :func:`flint_get_num_threads` threads.

.. function:: void _nmod_poly_compose_mod(nn_ptr res, nn_srcptr f, slong lenf, nn_srcptr g, nn_srcptr h, slong lenh, nmod_t mod)

    Sets ``res`` to the composition `f(g)` modulo `h`. We require that
    `h` is nonzero and that the length of `g` is one less than the
    length of `h` (possibly with zero padding). The output is not allowed
    to be aliased with any of the inputs.

.. function:: void nmod_poly_compose_mod(nmod_poly_t res, const nmod_poly_t f, const nmod_poly_t g, const nmod_poly_t h)

    Sets ``res`` to the composition `f(g)` modulo `h`. We require that
    `h` is nonzero.



Greatest common divisor
--------------------------------------------------------------------------------


.. function:: slong _nmod_poly_gcd_euclidean(nn_ptr G, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nmod_t mod)

    Computes the GCD of `A` of length ``lenA`` and `B` of length
    ``lenB``, where ``lenA >= lenB > 0``. The length of the GCD `G`
    is returned by the function. No attempt is made to make the GCD monic. It
    is required that `G` have space for ``lenB`` coefficients.

.. function:: void nmod_poly_gcd_euclidean(nmod_poly_t G, const nmod_poly_t A, const nmod_poly_t B)

    Computes the GCD of `A` and `B`. The GCD of zero polynomials is
    defined to be zero, whereas the GCD of the zero polynomial and some other
    polynomial `P` is defined to be `P`. Except in the case where
    the GCD is zero, the GCD `G` is made monic.

.. function:: slong _nmod_poly_hgcd(nn_ptr * M, slong * lenM, nn_ptr A, slong * lenA, nn_ptr B, slong * lenB, nn_srcptr a, slong lena, nn_srcptr b, slong lenb, nmod_t mod)

    Computes the HGCD of `a` and `b`, that is, a matrix `M`, a sign `\sigma`
    and two polynomials `A` and `B` such that

    .. math::


        (A,B)^t = M^{-1} (a,b)^t, \sigma = \det(M),


    and `A` and `B` are consecutive remainders in the Euclidean remainder
    sequence for the division of `a` by `b` satisfying \deg(A) \ge \frac{\deg(a)}{2} > \deg(B).
    Furthermore, `M` will be the product of ``[[q 1][1 0]]`` for the quotients ``q`` generated by such a remainder sequence.
    Assumes that `\operatorname{len}(a) > \operatorname{len}(b) > 0`, i.e. `\deg(a) > `deg(b) > 1`.

    Assumes that `A` and `B` have space of size at least `\operatorname{len}(a)`
    and `\operatorname{len}(b)`, respectively.  On exit, ``*lenA`` and ``*lenB``
    will contain the correct lengths of `A` and `B`.

    Assumes that ``M[0]``, ``M[1]``, ``M[2]``, and ``M[3]``
    each point to a vector of size at least `\operatorname{len}(a)`.

.. function:: slong _nmod_poly_gcd_hgcd(nn_ptr G, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nmod_t mod)

    Computes the monic GCD of `A` and `B`, assuming that
    `\operatorname{len}(A) \geq \operatorname{len}(B) > 0`.

    Assumes that `G` has space for `\operatorname{len}(B)` coefficients and
    returns the length of `G` on output.

.. function:: void nmod_poly_gcd_hgcd(nmod_poly_t G, const nmod_poly_t A, const nmod_poly_t B)

    Computes the monic GCD of `A` and `B` using the HGCD algorithm.

    As a special case, the GCD of two zero polynomials is defined to be
    the zero polynomial.

    The time complexity of the algorithm is `\mathcal{O}(n \log^2 n)`.
    For further details, see [ThullYap1990]_.

.. function:: slong _nmod_poly_gcd(nn_ptr G, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nmod_t mod)

    Computes the GCD of `A` of length ``lenA`` and `B` of length
    ``lenB``, where ``lenA >= lenB > 0``. The length of the GCD `G`
    is returned by the function. No attempt is made to make the GCD monic. It
    is required that `G` have space for ``lenB`` coefficients.

.. function:: void nmod_poly_gcd(nmod_poly_t G, const nmod_poly_t A, const nmod_poly_t B)

    Computes the GCD of `A` and `B`. The GCD of zero polynomials is
    defined to be zero, whereas the GCD of the zero polynomial and some other
    polynomial `P` is defined to be `P`. Except in the case where
    the GCD is zero, the GCD `G` is made monic.

.. function:: slong _nmod_poly_xgcd_euclidean(nn_ptr G, nn_ptr S, nn_ptr T, nn_srcptr A, slong A_len, nn_srcptr B, slong B_len, nmod_t mod)

    Computes the GCD of `A` and `B` together with cofactors `S` and `T`
    such that `S A + T B = G`.  Returns the length of `G`.

    Assumes that `\operatorname{len}(A) \geq \operatorname{len}(B) \geq 1` and
    `(\operatorname{len}(A),\operatorname{len}(B)) \neq (1,1)`.

    No attempt is made to make the GCD monic.

    Requires that `G` have space for `\operatorname{len}(B)` coefficients.  Writes
    `\operatorname{len}(B)-1` and `\operatorname{len}(A)-1` coefficients to `S` and `T`, respectively.
    Note that, in fact, `\operatorname{len}(S) \leq \max(\operatorname{len}(B) - \operatorname{len}(G), 1)` and
    `\operatorname{len}(T) \leq \max(\operatorname{len}(A) - \operatorname{len}(G), 1)`.

    No aliasing of input and output operands is permitted.

.. function:: void nmod_poly_xgcd_euclidean(nmod_poly_t G, nmod_poly_t S, nmod_poly_t T, const nmod_poly_t A, const nmod_poly_t B)

    Computes the GCD of `A` and `B`. The GCD of zero polynomials is
    defined to be zero, whereas the GCD of the zero polynomial and some other
    polynomial `P` is defined to be `P`. Except in the case where
    the GCD is zero, the GCD `G` is made monic.

    Polynomials ``S`` and ``T`` are computed such that
    ``S*A + T*B = G``. The length of ``S`` will be at most
    ``lenB`` and the length of ``T`` will be at most ``lenA``.

.. function:: slong _nmod_poly_xgcd_hgcd(nn_ptr G, nn_ptr S, nn_ptr T, nn_srcptr A, slong A_len, nn_srcptr B, slong B_len, nmod_t mod)

    Computes the GCD of `A` and `B`, where `\operatorname{len}(A) \geq \operatorname{len}(B) > 0`,
    together with cofactors `S` and `T` such that `S A + T B = G`. Returns
    the length of `G`.

    No attempt is made to make the GCD monic.

    Requires that `G` have space for `\operatorname{len}(B)` coefficients.  Writes
    `\operatorname{len}(B) - 1` and `\operatorname{len}(A) - 1` coefficients to `S` and `T`,
    respectively.  Note that, in fact, `\operatorname{len}(S) \leq \operatorname{len}(B) - \operatorname{len}(G)`
    and `\operatorname{len}(T) \leq \operatorname{len}(A) - \operatorname{len}(G)`.

    Both `S` and `T` must have space for at least `2` coefficients.

    No aliasing of input and output operands is permitted.

.. function:: void nmod_poly_xgcd_hgcd(nmod_poly_t G, nmod_poly_t S, nmod_poly_t T, const nmod_poly_t A, const nmod_poly_t B)

    Computes the GCD of `A` and `B`. The GCD of zero polynomials is
    defined to be zero, whereas the GCD of the zero polynomial and some other
    polynomial `P` is defined to be `P`. Except in the case where
    the GCD is zero, the GCD `G` is made monic.

    Polynomials ``S`` and ``T`` are computed such that
    ``S*A + T*B = G``. The length of ``S`` will be at most
    ``lenB`` and the length of ``T`` will be at most ``lenA``.

.. function:: slong _nmod_poly_xgcd(nn_ptr G, nn_ptr S, nn_ptr T, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nmod_t mod)

    Computes the GCD of `A` and `B`, where `\operatorname{len}(A) \geq \operatorname{len}(B) > 0`,
    together with cofactors `S` and `T` such that `S A + T B = G`. Returns
    the length of `G`.

    No attempt is made to make the GCD monic.

    Requires that `G` have space for `\operatorname{len}(B)` coefficients.  Writes
    `\operatorname{len}(B) - 1` and `\operatorname{len}(A) - 1` coefficients to `S` and `T`,
    respectively.  Note that, in fact, `\operatorname{len}(S) \leq \operatorname{len}(B) - \operatorname{len}(G)`
    and `\operatorname{len}(T) \leq \operatorname{len}(A) - \operatorname{len}(G)`.

    No aliasing of input and output operands is permitted.

.. function:: void nmod_poly_xgcd(nmod_poly_t G, nmod_poly_t S, nmod_poly_t T, const nmod_poly_t A, const nmod_poly_t B)

    Computes the GCD of `A` and `B`. The GCD of zero polynomials is
    defined to be zero, whereas the GCD of the zero polynomial and some other
    polynomial `P` is defined to be `P`. Except in the case where
    the GCD is zero, the GCD `G` is made monic.

    The polynomials ``S`` and ``T`` are set such that
    ``S*A + T*B = G``. The length of ``S`` will be at most
    ``lenB`` and the length of ``T`` will be at most ``lenA``.

.. function:: ulong _nmod_poly_resultant_euclidean(nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, nmod_t mod)

    Returns the resultant of ``(poly1, len1)`` and
    ``(poly2, len2)`` using the Euclidean algorithm.

    Assumes that ``len1 >= len2 > 0``.

    Assumes that the modulus is prime.

.. function:: ulong nmod_poly_resultant_euclidean(const nmod_poly_t f, const nmod_poly_t g)

    Computes the resultant of `f` and `g` using the Euclidean algorithm.

    For two non-zero polynomials `f(x) = a_m x^m + \dotsb + a_0` and
    `g(x) = b_n x^n + \dotsb + b_0` of degrees `m` and `n`, the resultant
    is defined to be

    .. math::
            a_m^n b_n^m \prod_{(x, y) : f(x) = g(y) = 0} (x - y).


    For convenience, we define the resultant to be equal to zero if either
    of the two polynomials is zero.

.. function:: ulong _nmod_poly_resultant_hgcd(nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, nmod_t mod)

    Returns the resultant of ``(poly1, len1)`` and
    ``(poly2, len2)`` using the half-gcd algorithm.

    This algorithm computes the half-gcd as per :func:`_nmod_poly_gcd_hgcd`
    but additionally updates the resultant every time a division occurs. The
    half-gcd algorithm computes the GCD recursively. Given inputs `a` and `b`
    it lets ``m = len(a)/2`` and (recursively) performs all quotients in
    the Euclidean algorithm which do not require the low `m` coefficients of
    `a` and `b`.

    This performs quotients in exactly the same order as the ordinary
    Euclidean algorithm except that the low `m` coefficients of the polynomials
    in the remainder sequence are not computed. A correction step after hgcd
    has been called computes these low `m` coefficients (by matrix
    multiplication by a transformation matrix also computed by hgcd).

    This means that from the point of view of the resultant, all but the last
    quotient performed by a recursive call to hgcd is an ordinary quotient as
    per the usual Euclidean algorithm. However, the final quotient may give
    a remainder of less than `m + 1` coefficients, which won't be corrected
    until the hgcd correction step is performed afterwards.

    To compute the adjustments to the resultant coming from this corrected
    quotient, we save the relevant information in an :type:`nmod_poly_res_t`
    struct at the time the quotient is performed so that when the correction
    step is performed later, the adjustments to the resultant can be computed
    at that time also.

    The only time an adjustment to the resultant is not required after a
    call to hgcd is if hgcd does nothing (the remainder may already have had
    less than `m + 1` coefficients when hgcd was called).

    Assumes that ``len1 >= len2 > 0``.

    Assumes that the modulus is prime.

.. function:: ulong nmod_poly_resultant_hgcd(const nmod_poly_t f, const nmod_poly_t g)

    Computes the resultant of `f` and `g` using the half-gcd algorithm.

    For two non-zero polynomials `f(x) = a_m x^m + \dotsb + a_0` and
    `g(x) = b_n x^n + \dotsb + b_0` of degrees `m` and `n`, the resultant
    is defined to be

    .. math::


            a_m^n b_n^m \prod_{(x, y) : f(x) = g(y) = 0} (x - y).


    For convenience, we define the resultant to be equal to zero if either
    of the two polynomials is zero.

.. function:: ulong _nmod_poly_resultant(nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, nmod_t mod)

    Returns the resultant of ``(poly1, len1)`` and
    ``(poly2, len2)``.

    Assumes that ``len1 >= len2 > 0``.

    Assumes that the modulus is prime.

.. function:: ulong nmod_poly_resultant(const nmod_poly_t f, const nmod_poly_t g)

    Computes the resultant of `f` and `g`.

    For two non-zero polynomials `f(x) = a_m x^m + \dotsb + a_0` and
    `g(x) = b_n x^n + \dotsb + b_0` of degrees `m` and `n`, the resultant
    is defined to be

    .. math::


            a_m^n b_n^m \prod_{(x, y) : f(x) = g(y) = 0} (x - y).


    For convenience, we define the resultant to be equal to zero if either
    of the two polynomials is zero.

.. function:: slong _nmod_poly_gcdinv(ulong * G, ulong * S, const ulong * A, slong lenA, const ulong * B, slong lenB, const nmod_t mod)

    Computes ``(G, lenA)``, ``(S, lenB-1)`` such that
    `G \cong S A \pmod{B}`, returning the actual length of `G`.

    Assumes that `0 < \operatorname{len}(A) < \operatorname{len}(B)`.

.. function:: void nmod_poly_gcdinv(nmod_poly_t G, nmod_poly_t S, const nmod_poly_t A, const nmod_poly_t B)

    Computes polynomials `G` and `S`, both reduced modulo `B`,
    such that `G \cong S A \pmod{B}`, where `B` is assumed to
    have `\operatorname{len}(B) \geq 2`.

    In the case that `A = 0 \pmod{B}`, returns `G = S = 0`.

.. function:: int _nmod_poly_invmod(ulong * A, const ulong * B, slong lenB, const ulong * P, slong lenP, const nmod_t mod)

    Attempts to set ``(A, lenP-1)`` to the inverse of ``(B, lenB)``
    modulo the polynomial ``(P, lenP)``.  Returns `1` if ``(B, lenB)``
    is invertible and `0` otherwise.

    Assumes that `0 < \operatorname{len}(B) < \operatorname{len}(P)`, and hence also `\operatorname{len}(P) \geq 2`,
    but supports zero-padding in ``(B, lenB)``.

    Does not support aliasing.

    Assumes that `mod` is a prime number.

.. function:: int nmod_poly_invmod(nmod_poly_t A, const nmod_poly_t B, const nmod_poly_t P)

    Attempts to set `A` to the inverse of `B` modulo `P` in the polynomial
    ring `(\mathbf{Z}/p\mathbf{Z})[X]`, where we assume that `p` is a prime
    number.

    If `\operatorname{len}(P) < 2`, raises an exception.

    If the greatest common divisor of `B` and `P` is `1`, returns `1` and
    sets `A` to the inverse of `B`.  Otherwise, returns `0` and the value
    of `A` on exit is undefined.



Discriminant
--------------------------------------------------------------------------------


.. function:: ulong _nmod_poly_discriminant(nn_srcptr poly, slong len, nmod_t mod)

    Return the discriminant of ``(poly, len)``. Assumes ``len > 1``.

.. function:: ulong nmod_poly_discriminant(const nmod_poly_t f)

    Return the discriminant of `f`.
    We normalise the discriminant so that
    `\operatorname{disc}(f) = (-1)^{n(n-1)/2} \operatorname{res}(f, f') /
    \operatorname{lc}(f)^{n - m - 2}`, where ``n = len(f)`` and
    ``m = len(f')``. Thus `\operatorname{disc}(f) =
    \operatorname{lc}(f)^{2n - 2} \prod_{i < j} (r_i - r_j)^2`, where
    `\operatorname{lc}(f)` is the leading coefficient of `f` and `r_i` are the
    roots of `f`.



Power series composition
--------------------------------------------------------------------------------

.. function:: void _nmod_poly_compose_series(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong n, nmod_t mod)

    Sets ``res`` to the composition of ``poly1`` and ``poly2``
    modulo `x^n`, where the constant term of ``poly2`` is required
    to be zero.

    Assumes that ``len1, len2, n > 0``, that ``len1, len2 <= n``,
    and that ``(len1-1) * (len2-1) + 1 <= n``, and that ``res`` has
    space for ``n`` coefficients. Does not support aliasing between any
    of the inputs and the output.

    Wraps :func:`_gr_poly_compose_series` which chooses automatically
    between various algorithms.

.. function:: void nmod_poly_compose_series(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2, slong n)

    Sets ``res`` to the composition of ``poly1`` and ``poly2``
    modulo `x^n`, where the constant term of ``poly2`` is required
    to be zero.


Power series reversion
--------------------------------------------------------------------------------

.. function:: void _nmod_poly_revert_series(nn_ptr Qinv, nn_srcptr Q, slong Qlen, slong n, nmod_t mod)
              void nmod_poly_revert_series(nmod_poly_t Qinv, const nmod_poly_t Q, slong n)

    Sets ``Qinv`` to the compositional inverse or reversion of ``Q``
    as a power series, i.e. computes `Q^{-1}` such that
    `Q(Q^{-1}(x)) = Q^{-1}(Q(x)) = x \bmod x^n`.

    It is required that `Q_0 = 0` and that `Q_1` as well as the integers
    `1, 2, \ldots, n-1` are invertible modulo the modulus.

    Wraps :func:`_gr_poly_revert_series` which chooses automatically
    between various algorithms.

Square roots
--------------------------------------------------------------------------------

The series expansions for `\sqrt{h}` and `1/\sqrt{h}` are defined
by means of the generalised binomial theorem
``h^r = (1+y)^r =
\sum_{k=0}^{\infty} {r \choose k} y^k.``
It is assumed that `h` has constant term `1` and that the coefficients
`2^{-k}` exist in the coefficient ring (i.e. `2` must be invertible).

.. function:: void _nmod_poly_invsqrt_series(nn_ptr g, nn_srcptr h, slong hlen, slong n, nmod_t mod)

    Set the first `n` terms of `g` to the series expansion of `1/\sqrt{h}`.
    It is assumed that `n > 0`, that `h` has constant term 1. Aliasing is not permitted.

.. function:: void nmod_poly_invsqrt_series(nmod_poly_t g, const nmod_poly_t h, slong n)

    Set `g` to the series expansion of `1/\sqrt{h}` to order `O(x^n)`.
    It is assumed that `h` has constant term 1.

.. function:: void _nmod_poly_sqrt_series(nn_ptr g, nn_srcptr h, slong hlen, slong n, nmod_t mod)

    Set the first `n` terms of `g` to the series expansion of `\sqrt{h}`.
    It is assumed that `n > 0`, that `h` has constant term 1. Aliasing is not permitted.

.. function:: void nmod_poly_sqrt_series(nmod_poly_t g, const nmod_poly_t h, slong n)

    Set `g` to the series expansion of `\sqrt{h}` to order `O(x^n)`.
    It is assumed that `h` has constant term 1.

.. function:: int _nmod_poly_sqrt(nn_ptr s, nn_srcptr p, slong n, nmod_t mod)

    If ``(p, n)`` is a perfect square, sets ``(s, n / 2 + 1)``
    to a square root of `p` and returns 1. Otherwise returns 0.

.. function:: int nmod_poly_sqrt(nmod_poly_t s, const nmod_poly_t p)

    If `p` is a perfect square, sets `s` to a square root of `p`
    and returns 1. Otherwise returns 0.


Power sums
--------------------------------------------------------------------------------


.. function:: void _nmod_poly_power_sums_naive(nn_ptr res, nn_srcptr poly, slong len, slong n, nmod_t mod)

    Compute the (truncated) power sums series of the polynomial
    ``(poly,len)`` up to length `n` using Newton identities.

.. function:: void nmod_poly_power_sums_naive(nmod_poly_t res, const nmod_poly_t poly, slong n)

    Compute the (truncated) power sum series of the polynomial
    ``poly`` up to length `n` using Newton identities.

.. function:: void _nmod_poly_power_sums_schoenhage(nn_ptr res, nn_srcptr poly, slong len, slong n, nmod_t mod)

    Compute the (truncated) power sums series of the polynomial
    ``(poly,len)`` up to length `n` using a series expansion
    (a formula due to Schoenhage).

.. function:: void nmod_poly_power_sums_schoenhage(nmod_poly_t res, const nmod_poly_t poly, slong n)

    Compute the (truncated) power sums series of the polynomial
    ``poly`` up to length `n` using a series expansion
    (a formula due to Schoenhage).

.. function:: void _nmod_poly_power_sums(nn_ptr res, nn_srcptr poly, slong len, slong n, nmod_t mod)

    Compute the (truncated) power sums series of the polynomial
    ``(poly,len)`` up to length `n`.

.. function:: void nmod_poly_power_sums(nmod_poly_t res, const nmod_poly_t poly, slong n)

    Compute the (truncated) power sums series of the polynomial
    ``poly`` up to length `n`.

.. function:: void _nmod_poly_power_sums_to_poly_naive(nn_ptr res, nn_srcptr poly, slong len, nmod_t mod)

    Compute the (monic) polynomial given by its power sums series
    ``(poly,len)`` using Newton identities.

.. function:: void nmod_poly_power_sums_to_poly_naive(nmod_poly_t res, const nmod_poly_t Q)

    Compute the (monic) polynomial given by its power sums series
    ``Q`` using Newton identities.

.. function:: void _nmod_poly_power_sums_to_poly_schoenhage(nn_ptr res, nn_srcptr poly, slong len, nmod_t mod)

    Compute the (monic) polynomial given by its power sums series
    ``(poly,len)`` using series expansion (a formula due to Schoenhage).

.. function:: void nmod_poly_power_sums_to_poly_schoenhage(nmod_poly_t res, const nmod_poly_t Q)

    Compute the (monic) polynomial given by its power sums series
    ``Q`` using series expansion (a formula due to Schoenhage).

.. function:: void _nmod_poly_power_sums_to_poly(nn_ptr res, nn_srcptr poly, slong len, nmod_t mod)

    Compute the (monic) polynomial given by its power sums series
    ``(poly,len)``.

.. function:: void nmod_poly_power_sums_to_poly(nmod_poly_t res, const nmod_poly_t Q)

     Compute the (monic) polynomial given by its power sums series ``Q``.


Transcendental functions
--------------------------------------------------------------------------------

The elementary transcendental functions of a formal power series `h`
are defined as

`\exp(h(x)) = \sum_{k=0}^{\infty} \frac{(h(x))^k}{k!}`

`\log(h(x)) = \int_0^x \frac{h'(t)}{h(t)} dt`

`\operatorname{atan}(h(x)) = \int_0^x\frac{h'(t)}{1+(h(t))^2} dt`

`\operatorname{atanh}(h(x)) = \int_0^x\frac{h'(t)}{1-(h(t))^2} dt`

`\operatorname{asin}(h(x)) = \int_0^x\frac{h'(t)}{\sqrt{1-(h(t))^2}} dt`

`\operatorname{asinh}(h(x)) = \int_0^x\frac{h'(t)}{\sqrt{1+(h(t))^2}} dt`

The functions sin, cos, tan, etc. are defined using standard inverse
or functional relations.
The logarithm function assumes that `h` has constant term `1`. All
other functions assume that `h` has constant term `0`.
All functions assume that the coefficient `1/k` or `1/k!` exists
for all indices `k`. When computing to order `O(x^n)`, the modulus `p`
must therefore be a prime satisfying `p \ge n`. Further, we always
require that `p > 2` in order to be able to multiply by `1/2` for
internal purposes.
If the input does not satisfy all these conditions, results are undefined.
Except where otherwise noted, functions are implemented with optimal
(up to constants) complexity `O(M(n))`, where `M(n)` is the cost
of polynomial multiplication.

.. function:: void _nmod_poly_log_series(nn_ptr g, nn_srcptr h, slong hlen, slong n, nmod_t mod)

    Set `g = \log(h) + O(x^n)`. Assumes `n > 0` and ``hlen > 0``.
    Aliasing of `g` and `h` is allowed.

.. function:: void nmod_poly_log_series(nmod_poly_t g, const nmod_poly_t h, slong n)

    Set `g = \log(h) + O(x^n)`. The case `h = 1+cx^r` is automatically
    detected and handled efficiently.

.. function:: void _nmod_poly_exp_series(nn_ptr f, nn_srcptr h, slong hlen, slong n, nmod_t mod)

    Set `f = \exp(h) + O(x^n)` where ``h`` is a polynomial. Assume
    `n > 0`. Aliasing of `g` and `h` is not allowed.

    Uses Newton iteration (an improved version of the
    algorithm in [HanZim2004]_).
    For small `n`, falls back to the basecase algorithm.

.. function:: void  _nmod_poly_exp_expinv_series(nn_ptr f, nn_ptr g, nn_srcptr h, slong hlen, slong n, nmod_t mod)

    Set `f = \exp(h) + O(x^n)` and `g = \exp(-h) + O(x^n)`, more efficiently
    for large `n` than performing a separate inversion to obtain `g`.
    Assumes `n > 0` and that `h` is zero-padded
    as necessary to length `n`. Aliasing is not allowed.

    Uses Newton iteration (the version given in [HanZim2004]_).
    For small `n`, falls back to the basecase algorithm.

.. function:: void nmod_poly_exp_series(nmod_poly_t g, const nmod_poly_t h, slong n)

    Set `g = \exp(h) + O(x^n)`. The case `h = cx^r` is automatically
    detected and handled efficiently. Otherwise this function automatically
    uses the basecase algorithm for small `n` and Newton iteration otherwise.

.. function:: void _nmod_poly_atan_series(nn_ptr g, nn_srcptr h, slong hlen, slong n, nmod_t mod)

    Set `g = \operatorname{atan}(h) + O(x^n)`. Assumes `n > 0`.
    Aliasing of `g` and `h` is allowed.

.. function:: void nmod_poly_atan_series(nmod_poly_t g, const nmod_poly_t h, slong n)

    Set `g = \operatorname{atan}(h) + O(x^n)`.

.. function:: void _nmod_poly_atanh_series(nn_ptr g, nn_srcptr h, slong hlen, slong n, nmod_t mod)

    Set `g = \operatorname{atanh}(h) + O(x^n)`. Assumes `n > 0`.
    Aliasing of `g` and `h` is allowed.

.. function:: void nmod_poly_atanh_series(nmod_poly_t g, const nmod_poly_t h, slong n)

    Set `g = \operatorname{atanh}(h) + O(x^n)`.

.. function:: void _nmod_poly_asin_series(nn_ptr g, nn_srcptr h, slong hlen, slong n, nmod_t mod)

    Set `g = \operatorname{asin}(h) + O(x^n)`. Assumes `n > 0`.
    Aliasing of `g` and `h` is allowed.

.. function:: void nmod_poly_asin_series(nmod_poly_t g, const nmod_poly_t h, slong n)

    Set `g = \operatorname{asin}(h) + O(x^n)`.

.. function:: void _nmod_poly_asinh_series(nn_ptr g, nn_srcptr h, slong hlen, slong n, nmod_t mod)

    Set `g = \operatorname{asinh}(h) + O(x^n)`. Assumes `n > 0`.
    Aliasing of `g` and `h` is allowed.

.. function:: void nmod_poly_asinh_series(nmod_poly_t g, const nmod_poly_t h, slong n)

    Set `g = \operatorname{asinh}(h) + O(x^n)`.

.. function:: void _nmod_poly_sin_series(nn_ptr g, nn_srcptr h, slong n, nmod_t mod)

    Set `g = \operatorname{sin}(h) + O(x^n)`. Assumes `n > 0` and that `h`
    is zero-padded as necessary to length `n`. Aliasing of `g` and `h` is
    allowed. The value is computed using the identity
    `\sin(x) = 2 \tan(x/2)) / (1 + \tan^2(x/2)).`

.. function:: void nmod_poly_sin_series(nmod_poly_t g, const nmod_poly_t h, slong n)

    Set `g = \operatorname{sin}(h) + O(x^n)`.

.. function:: void _nmod_poly_cos_series(nn_ptr g, nn_srcptr h, slong n, nmod_t mod)

    Set `g = \operatorname{cos}(h) + O(x^n)`. Assumes `n > 0` and that `h`
    is zero-padded as necessary to length `n`. Aliasing of `g` and `h` is
    allowed. The value is computed using the identity
    `\cos(x) = (1-\tan^2(x/2)) / (1 + \tan^2(x/2)).`

.. function:: void nmod_poly_cos_series(nmod_poly_t g, const nmod_poly_t h, slong n)

    Set `g = \operatorname{cos}(h) + O(x^n)`.

.. function:: void _nmod_poly_tan_series(nn_ptr g, nn_srcptr h, slong hlen, slong n, nmod_t mod)

    Set `g = \operatorname{tan}(h) + O(x^n)`. Assumes `n > 0` and that `h`
    is zero-padded as necessary to length `n`. Aliasing of `g` and `h` is
    not allowed. Uses Newton iteration to invert the atan function.

.. function:: void nmod_poly_tan_series(nmod_poly_t g, const nmod_poly_t h, slong n)

    Set `g = \operatorname{tan}(h) + O(x^n)`.

.. function:: void _nmod_poly_sinh_series(nn_ptr g, nn_srcptr h, slong n, nmod_t mod)

    Set `g = \operatorname{sinh}(h) + O(x^n)`. Assumes `n > 0` and that `h`
    is zero-padded as necessary to length `n`. Aliasing of `g` and `h` is
    not allowed. Uses the identity `\sinh(x) = (e^x - e^{-x})/2`.

.. function:: void nmod_poly_sinh_series(nmod_poly_t g, const nmod_poly_t h, slong n)

    Set `g = \operatorname{sinh}(h) + O(x^n)`.

.. function:: void _nmod_poly_cosh_series(nn_ptr g, nn_srcptr h, slong n, nmod_t mod)

    Set `g = \operatorname{cos}(h) + O(x^n)`. Assumes `n > 0` and that `h`
    is zero-padded as necessary to length `n`. Aliasing of `g` and `h` is
    not allowed.
    Uses the identity `\cosh(x) = (e^x + e^{-x})/2`.

.. function:: void nmod_poly_cosh_series(nmod_poly_t g, const nmod_poly_t h, slong n)

    Set `g = \operatorname{cosh}(h) + O(x^n)`.

.. function:: void _nmod_poly_tanh_series(nn_ptr g, nn_srcptr h, slong n, nmod_t mod)

    Set `g = \operatorname{tanh}(h) + O(x^n)`. Assumes `n > 0` and that `h`
    is zero-padded as necessary to length `n`. Uses the identity
    `\tanh(x) = (e^{2x}-1)/(e^{2x}+1)`.

.. function:: void nmod_poly_tanh_series(nmod_poly_t g, const nmod_poly_t h, slong n)

    Set `g = \operatorname{tanh}(h) + O(x^n)`.


Special polynomials
--------------------------------------------------------------------------------

.. function:: int _nmod_poly_conway(nn_ptr op, ulong prime, slong deg)

    Sets ``op`` to the coefficients to the Conway polynomial `C_{p, d}`, where
    `p` is ``prime`` and `d` is ``deg``. This is done by checking against Frank
    Lübeck's database [Lüb2004]_, which has been compressed in FLINT. Returns
    `1` in case of success and returns `0` in case of failure.

.. function:: ulong _nmod_poly_conway_rand(slong * degree, flint_rand_t state, int type)

    Returns a pseudorandom prime and sets ``degree`` that when put into
    :func:`_nmod_poly_conway` will always succeed.

    Here, ``type`` can be the following values:

    * ``0`` for which there is a bijection between the image of this function
      and the database of Conway polynomials,
    * ``1`` returns a random prime found in the database and sets ``degree`` to
      some degree less than `15` along with some prime found in the database,
    * ``2`` returns a random prime less than `2^{10}` and sets ``degree`` to
      some random degree found in the database,
    * ``3`` returns a random prime less than `2^{10}` and sets ``degree`` to
      some random degree less than `15`.


Products
--------------------------------------------------------------------------------


.. function:: void _nmod_poly_product_roots_nmod_vec(nn_ptr poly, nn_srcptr xs, slong n, nmod_t mod)

    Sets ``(poly, n + 1)`` to the monic polynomial which is the product
    of `(x - x_0)(x - x_1) \cdots (x - x_{n-1})`, the roots `x_i` being
    given by ``xs``.

    Aliasing of the input and output is not allowed.

.. function:: void nmod_poly_product_roots_nmod_vec(nmod_poly_t poly, nn_srcptr xs, slong n)

    Sets ``poly`` to the monic polynomial which is the product
    of `(x - x_0)(x - x_1) \cdots (x - x_{n-1})`, the roots `x_i` being
    given by ``xs``.

.. function:: int nmod_poly_find_distinct_nonzero_roots(ulong * roots, const nmod_poly_t A)

    If ``A`` has `\deg(A)` distinct nonzero roots in `\mathbb{F}_p`, write these roots out to ``roots[0]`` to ``roots[deg(A) - 1]`` and return ``1``.
    Otherwise, return ``0``. It is assumed that ``A`` is nonzero and that the modulus of ``A`` is prime.
    This function uses Rabin's probabilistic method via gcd's with `(x + \delta)^{\frac{p-1}{2}} - 1`.


Subproduct trees
--------------------------------------------------------------------------------


.. function:: nn_ptr * _nmod_poly_tree_alloc(slong len)

    Allocates space for a subproduct tree of the given length, having
    linear factors at the lowest level.

    Entry `i` in the tree is a pointer to a single array of limbs,
    capable of storing `\lfloor n / 2^i \rfloor` subproducts of
    degree `2^i` adjacently, plus a trailing entry if `n / 2^i` is
    not an integer.

    For example, a tree of length 7 built from monic linear factors has
    the following structure, where spaces have been inserted
    for illustrative purposes::

        X1 X1 X1 X1 X1 X1 X1
        XX1   XX1   XX1   X1
        XXXX1       XX1   X1
        XXXXXXX1

.. function:: void _nmod_poly_tree_free(nn_ptr * tree, slong len)

    Free the allocated space for the subproduct.

.. function:: void _nmod_poly_tree_build(nn_ptr * tree, nn_srcptr roots, slong len, nmod_t mod)

    Builds a subproduct tree in the preallocated space from
    the ``len`` monic linear factors `(x-r_i)`. The top level
    product is not computed.



Geometric progression
--------------------------------------------------------------------------------


.. function:: void nmod_geometric_progression_init(nmod_geometric_progression_t G, ulong r, slong len, nmod_t mod)

    Builds a geometric progression multipoint evaluation / interpolation structure.

    The set of points used will be `1, r^2, r^4, \ldots, r^{2(len-1)}`.

    The value of ``r`` should be reduced modulo the modulus ``mod``
    and of sufficient multiplicative order such that none of
    the powers `r^2, r^4, \ldots, r^{2(len-1)}` is one.

    The value of ``len`` should be both greater than or equal to the number of evaluation points to be
    considered, and greater than or equal to the length of the polynomials to be evaluated / interpolated.
    This allocates vectors and polynomials for a total space of `8 len - 1` coefficients.

    If the modulus is not prime, this function will work under the additional
    assumption that all the used points `r^{2k}` as well as the axuiliary
    values `r^{2k} - 1` are invertible.

.. function:: void nmod_geometric_progression_clear(nmod_geometric_progression_t G)

    Clears the allocated polynomials and vectors used in the geometric progression precomputation ``G``.


Inflation and deflation
--------------------------------------------------------------------------------


.. function:: void nmod_poly_inflate(nmod_poly_t result, const nmod_poly_t input, slong inflation)

    Sets ``result`` to the inflated polynomial `p(x^n)` where
    `p` is given by ``input`` and `n` is given by ``deflation``.

.. function:: void nmod_poly_deflate(nmod_poly_t result, const nmod_poly_t input, slong deflation)

    Sets ``result`` to the deflated polynomial `p(x^{1/n})` where
    `p` is given by ``input`` and `n` is given by ``deflation``.
    Requires `n > 0`.

.. function:: slong nmod_poly_deflation(const nmod_poly_t input)

    Returns the largest integer by which ``input`` can be deflated.
    As special cases, returns 0 if ``input`` is the zero polynomial
    and 1 of ``input`` is a constant polynomial.


Chinese Remaindering
--------------------------------------------------------------------------------

    In all of these functions the moduli (mod.n) of all of the ``nmod_poly``'s involved is assumed to match and be prime.

.. function:: void nmod_poly_multi_crt_init(nmod_poly_multi_crt_t CRT)

    Initialize ``CRT`` for Chinese remaindering.

.. function:: int nmod_poly_multi_crt_precompute(nmod_poly_multi_crt_t CRT, const nmod_poly_struct * moduli, slong len)
              int nmod_poly_multi_crt_precompute_p(nmod_poly_multi_crt_t CRT, const nmod_poly_struct * const * moduli, slong len)

    Configure ``CRT`` for repeated Chinese remaindering of ``moduli``. The number of moduli, ``len``, should be positive.
    A return of ``0`` indicates that the compilation failed and future calls to :func:`nmod_poly_multi_crt_precomp` will leave the output undefined.
    A return of ``1`` indicates that the compilation was successful, which occurs if and only if either (1) ``len == 1`` and ``modulus + 0`` is nonzero, or (2) all of the moduli have positive degree and are pairwise relatively prime.

.. function:: void nmod_poly_multi_crt_precomp(nmod_poly_t output, const nmod_poly_multi_crt_t CRT, const nmod_poly_struct * values)
              void nmod_poly_multi_crt_precomp_p(nmod_poly_t output, const nmod_poly_multi_crt_t CRT, const nmod_poly_struct * const * values)

    Set ``output`` to the polynomial of lowest possible degree that is congruent to ``values + i`` modulo the ``moduli + i`` in :func:`nmod_poly_multi_crt_precompute`.
    The inputs ``values + 0, ..., values + len - 1`` where ``len`` was used in :func:`nmod_poly_multi_crt_precompute` are expected to be valid and have modulus matching the modulus of the moduli used in :func:`nmod_poly_multi_crt_precompute`.

.. function:: int nmod_poly_multi_crt(nmod_poly_t output, const nmod_poly_struct * moduli, const nmod_poly_struct * values, slong len)

    Perform the same operation as :func:`nmod_poly_multi_crt_precomp` while internally constructing and destroying the precomputed data.
    All of the remarks in :func:`nmod_poly_multi_crt_precompute` apply.

.. function:: void nmod_poly_multi_crt_clear(nmod_poly_multi_crt_t CRT)

    Free all space used by ``CRT``.

.. function:: slong _nmod_poly_multi_crt_local_size(const nmod_poly_multi_crt_t CRT)

    Return the required length of the output for :func:`_nmod_poly_multi_crt_run`.

.. function:: void _nmod_poly_multi_crt_run(nmod_poly_struct * outputs, const nmod_poly_multi_crt_t CRT, const nmod_poly_struct * inputs)
              void _nmod_poly_multi_crt_run_p(nmod_poly_struct * outputs, const nmod_poly_multi_crt_t CRT, const nmod_poly_struct * const * inputs)

    Perform the same operation as :func:`nmod_poly_multi_crt_precomp` using supplied temporary space.
    The actual output is placed in ``outputs + 0``, and ``outputs`` should contain space for all temporaries and should be at least as long as ``_nmod_poly_multi_crt_local_size(CRT)``.
    Of course the moduli of these temporaries should match the modulus of the inputs.


Berlekamp-Massey Algorithm
--------------------------------------------------------------------------------

    The nmod_berlekamp_massey_t manages an unlimited stream of points `a_1, a_2, \dots.`
    At any point in time, after, say, `n` points have been added, a call to :func:`nmod_berlekamp_massey_reduce` will
    calculate the polynomials `U`, `V` and `R` in the extended euclidean remainder sequence with

    .. math::

        U x^n + V (a_1 x^{n-1} + a_{n-1} x + \cdots + a_n) = R, \quad \deg(U) < \deg(V) \le n/2, \quad \deg(R) < n/2.

    The polynomials `V` and `R` may be obtained with :func:`nmod_berlekamp_massey_V_poly` and :func:`nmod_berlekamp_massey_R_poly`.
    This class differs from :func:`fmpz_mod_poly_minpoly` in the following respect. Let `v_i` denote the coefficient of `x^i` in `V`.
    :func:`fmpz_mod_poly_minpoly` will return a polynomial `V` of lowest degree that annihilates the whole sequence `a_1, \dots, a_n` as

    .. math::

        \sum_{i} v_i a_{j + i} = 0, \quad 1 \le j \le n - \deg(V).

    The cost is that a polynomial of degree `n-1` might be returned and the return is not generally uniquely determined by the input sequence.
    For the nmod_berlekamp_massey_t we have

    .. math::

        \sum_{i,j} v_i a_{j+i} x^{-j} = -U + \frac{R}{x^n}\text{,}

    and it can be seen that `\sum_{i} v_i a_{j + i}` is zero for `1 \le j < n - \deg(R)`. Thus whether or not `V` has annihilated the whole sequence may be checked by comparing the degrees of `V` and `R`.

.. function:: void nmod_berlekamp_massey_init(nmod_berlekamp_massey_t B, ulong p)

    Initialize ``B`` in characteristic ``p`` with an empty stream.

.. function:: void nmod_berlekamp_massey_clear(nmod_berlekamp_massey_t B)

    Free any space used by ``B``.

.. function:: void nmod_berlekamp_massey_start_over(nmod_berlekamp_massey_t B)

    Empty the stream of points in ``B``.

.. function:: void nmod_berlekamp_massey_set_prime(nmod_berlekamp_massey_t B, ulong p)

    Set the characteristic of the field and empty the stream of points in ``B``.

.. function:: void nmod_berlekamp_massey_add_points(nmod_berlekamp_massey_t B, const ulong * a, slong count)
              void nmod_berlekamp_massey_add_zeros(nmod_berlekamp_massey_t B, slong count)
              void nmod_berlekamp_massey_add_point(nmod_berlekamp_massey_t B, ulong a)

    Add point(s) to the stream processed by ``B``. The addition of any number of points will not update the `V` and `R` polynomial.

.. function:: int nmod_berlekamp_massey_reduce(nmod_berlekamp_massey_t B)

    Ensure that the polynomials `V` and `R` are up to date. The return value is ``1`` if this function changed `V` and ``0`` otherwise.
    For example, if this function is called twice in a row without adding any points in between, the return of the second call should be ``0``.
    As another example, suppose the object is emptied, the points `1, 1, 2, 3` are added, then reduce is called. This reduce should return ``1`` with `\deg(R) < \deg(V) = 2` because the Fibonacci sequence has been recognized. The further addition of the two points `5, 8` and a reduce will result in a return value of ``0``.

.. function:: slong nmod_berlekamp_massey_point_count(const nmod_berlekamp_massey_t B)

    Return the number of points stored in ``B``.

.. function:: const ulong * nmod_berlekamp_massey_points(const nmod_berlekamp_massey_t B)

    Return a pointer to the array of points stored in ``B``. This may be ``NULL`` if :func:`nmod_berlekamp_massey_point_count` returns ``0``.

.. function:: const nmod_poly_struct * nmod_berlekamp_massey_V_poly(const nmod_berlekamp_massey_t B)

    Return the polynomial `V` in ``B``.

.. function:: const nmod_poly_struct * nmod_berlekamp_massey_R_poly(const nmod_berlekamp_massey_t B)

    Return the polynomial `R` in ``B``.