add arb_fmpz_poly module with evaluation, plus a function to construct minimal polynomials of Gaussian periods

Makefile.in

@@ -15,6 +15,7 @@ AT=@
 BUILD_DIRS = fmpr arf mag arb arb_mat arb_poly arb_calc acb acb_mat acb_poly \
    acb_calc acb_hypgeom acb_elliptic acb_modular dirichlet acb_dirichlet \
    arb_hypgeom bernoulli hypgeom fmpz_extras bool_mat partitions dlog \
+   arb_fmpz_poly \
    $(EXTRA_BUILD_DIRS)
 
 TEMPLATE_DIRS = 

arb_fmpz_poly.h

@@ -0,0 +1,45 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef ARB_FMPZ_POLY_H
+#define ARB_FMPZ_POLY_H
+
+#ifdef ARB_FMPZ_POLY_INLINES_C
+#define ARB_FMPZ_POLY_INLINE
+#else
+#define ARB_FMPZ_POLY_INLINE static __inline__
+#endif
+
+#include "flint/ulong_extras.h"
+#include "flint/fmpz.h"
+#include "flint/fmpz_poly.h"
+#include "arb.h"
+#include "acb.h"
+#include "arb_poly.h"
+#include "acb_poly.h"
+
+void _arb_fmpz_poly_evaluate_acb_horner(acb_t res, const fmpz * f, slong len, const acb_t x, slong prec);
+void arb_fmpz_poly_evaluate_acb_horner(acb_t res, const fmpz_poly_t f, const acb_t a, slong prec);
+void _arb_fmpz_poly_evaluate_acb_rectangular(acb_t res, const fmpz * f, slong len, const acb_t x, slong prec);
+void arb_fmpz_poly_evaluate_acb_rectangular(acb_t res, const fmpz_poly_t f, const acb_t a, slong prec);
+void _arb_fmpz_poly_evaluate_acb(acb_t res, const fmpz * f, slong len, const acb_t x, slong prec);
+void arb_fmpz_poly_evaluate_acb(acb_t res, const fmpz_poly_t f, const acb_t a, slong prec);
+
+void _arb_fmpz_poly_evaluate_arb_horner(arb_t res, const fmpz * f, slong len, const arb_t x, slong prec);
+void arb_fmpz_poly_evaluate_arb_horner(arb_t res, const fmpz_poly_t f, const arb_t a, slong prec);
+void _arb_fmpz_poly_evaluate_arb_rectangular(arb_t res, const fmpz * f, slong len, const arb_t x, slong prec);
+void arb_fmpz_poly_evaluate_arb_rectangular(arb_t res, const fmpz_poly_t f, const arb_t a, slong prec);
+void _arb_fmpz_poly_evaluate_arb(arb_t res, const fmpz * f, slong len, const arb_t x, slong prec);
+void arb_fmpz_poly_evaluate_arb(arb_t res, const fmpz_poly_t f, const arb_t a, slong prec);
+
+void arb_fmpz_poly_gauss_period_minpoly(fmpz_poly_t res, ulong q, ulong n);
+
+#endif

arb_fmpz_poly/evaluate_acb.c

@@ -0,0 +1,34 @@
+/*
+    Copyright (C) 2013 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_fmpz_poly.h"
+
+void
+_arb_fmpz_poly_evaluate_acb(acb_t res, const fmpz * f, slong len,
+                           const acb_t x, slong prec)
+{
+    if (acb_is_real(x))
+    {
+        _arb_fmpz_poly_evaluate_arb(acb_realref(res), f, len, acb_realref(x), prec);
+        arb_zero(acb_imagref(res));
+    }
+    else
+    {
+        _arb_fmpz_poly_evaluate_acb_rectangular(res, f, len, x, prec);
+    }
+}
+
+void
+arb_fmpz_poly_evaluate_acb(acb_t res, const fmpz_poly_t f, const acb_t a, slong prec)
+{
+    _arb_fmpz_poly_evaluate_acb(res, f->coeffs, f->length, a, prec);
+}
+

arb_fmpz_poly/evaluate_acb_horner.c

@@ -0,0 +1,59 @@
+/*
+    Copyright (C) 2010 Sebastian Pancratz
+    Copyright (C) 2012 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_fmpz_poly.h"
+
+void
+_arb_fmpz_poly_evaluate_acb_horner(acb_t y, const fmpz * f, slong len,
+                           const acb_t x, slong prec)
+{
+    if (len == 0)
+    {
+        acb_zero(y);
+    }
+    else if (len == 1 || acb_is_zero(x))
+    {
+        acb_set_round_fmpz(y, f, prec);
+    }
+    else if (len == 2)
+    {
+        acb_mul_fmpz(y, x, f + 1, prec);
+        acb_add_fmpz(y, y, f + 0, prec);
+    }
+    else
+    {
+        slong i = len - 1;
+        acb_t t, u;
+
+        acb_init(t);
+        acb_init(u);
+        acb_set_fmpz(u, f + i);
+
+        for (i = len - 2; i >= 0; i--)
+        {
+            acb_mul(t, u, x, prec);
+            acb_add_fmpz(u, t, f + i, prec);
+        }
+
+        acb_swap(y, u);
+
+        acb_clear(t);
+        acb_clear(u);
+    }
+}
+
+void
+arb_fmpz_poly_evaluate_acb_horner(acb_t res, const fmpz_poly_t f, const acb_t a, slong prec)
+{
+    _arb_fmpz_poly_evaluate_acb_horner(res, f->coeffs, f->length, a, prec);
+}
+

arb_fmpz_poly/evaluate_acb_rectangular.c

@@ -0,0 +1,63 @@
+/*
+    Copyright (C) 2013 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_fmpz_poly.h"
+
+void
+_arb_fmpz_poly_evaluate_acb_rectangular(acb_t y, const fmpz * poly,
+    slong len, const acb_t x, slong prec)
+{
+    slong i, j, m, r;
+    acb_ptr xs;
+    acb_t s, t, c;
+
+    if (len < 3)
+    {
+        _arb_fmpz_poly_evaluate_acb_horner(y, poly, len, x, prec);
+        return;
+    }
+
+    m = n_sqrt(len) + 1;
+    r = (len + m - 1) / m;
+
+    xs = _acb_vec_init(m + 1);
+    acb_init(s);
+    acb_init(t);
+    acb_init(c);
+
+    _acb_vec_set_powers(xs, x, m + 1, prec);
+
+    acb_set_fmpz(y, poly + (r - 1) * m);
+    for (j = 1; (r - 1) * m + j < len; j++)
+        acb_addmul_fmpz(y, xs + j, poly + (r - 1) * m + j, prec);
+
+    for (i = r - 2; i >= 0; i--)
+    {
+        acb_set_fmpz(s, poly + i * m);
+        for (j = 1; j < m; j++)
+            acb_addmul_fmpz(s, xs + j, poly + i * m + j, prec);
+
+        acb_mul(y, y, xs + m, prec);
+        acb_add(y, y, s, prec);
+    }
+
+    _acb_vec_clear(xs, m + 1);
+    acb_clear(s);
+    acb_clear(t);
+    acb_clear(c);
+}
+
+void
+arb_fmpz_poly_evaluate_acb_rectangular(acb_t res, const fmpz_poly_t f, const acb_t a, slong prec)
+{
+    _arb_fmpz_poly_evaluate_acb_rectangular(res, f->coeffs, f->length, a, prec);
+}
+

arb_fmpz_poly/evaluate_arb.c

@@ -0,0 +1,39 @@
+/*
+    Copyright (C) 2013 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_fmpz_poly.h"
+
+void
+_arb_fmpz_poly_evaluate_arb(arb_t res, const fmpz * f, slong len,
+                           const arb_t x, slong prec)
+{
+    if ((prec >= 1024) && (len >= 5 + 20000 / prec))
+    {
+        slong fbits;
+
+        fbits = _fmpz_vec_max_bits(f, len);
+
+        if (fbits <= prec / 2)
+        {
+            _arb_fmpz_poly_evaluate_arb_rectangular(res, f, len, x, prec);
+            return;
+        }
+    }
+
+    _arb_fmpz_poly_evaluate_arb_horner(res, f, len, x, prec);
+}
+
+void
+arb_fmpz_poly_evaluate_arb(arb_t res, const fmpz_poly_t f, const arb_t a, slong prec)
+{
+    _arb_fmpz_poly_evaluate_arb(res, f->coeffs, f->length, a, prec);
+}
+

arb_fmpz_poly/evaluate_arb_horner.c

@@ -0,0 +1,59 @@
+/*
+    Copyright (C) 2010 Sebastian Pancratz
+    Copyright (C) 2012 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_fmpz_poly.h"
+
+void
+_arb_fmpz_poly_evaluate_arb_horner(arb_t y, const fmpz * f, slong len,
+                           const arb_t x, slong prec)
+{
+    if (len == 0)
+    {
+        arb_zero(y);
+    }
+    else if (len == 1 || arb_is_zero(x))
+    {
+        arb_set_round_fmpz(y, f, prec);
+    }
+    else if (len == 2)
+    {
+        arb_mul_fmpz(y, x, f + 1, prec);
+        arb_add_fmpz(y, y, f + 0, prec);
+    }
+    else
+    {
+        slong i = len - 1;
+        arb_t t, u;
+
+        arb_init(t);
+        arb_init(u);
+        arb_set_fmpz(u, f + i);
+
+        for (i = len - 2; i >= 0; i--)
+        {
+            arb_mul(t, u, x, prec);
+            arb_add_fmpz(u, t, f + i, prec);
+        }
+
+        arb_swap(y, u);
+
+        arb_clear(t);
+        arb_clear(u);
+    }
+}
+
+void
+arb_fmpz_poly_evaluate_arb_horner(arb_t res, const fmpz_poly_t f, const arb_t a, slong prec)
+{
+    _arb_fmpz_poly_evaluate_arb_horner(res, f->coeffs, f->length, a, prec);
+}
+

arb_fmpz_poly/evaluate_arb_rectangular.c

@@ -0,0 +1,63 @@
+/*
+    Copyright (C) 2013 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_fmpz_poly.h"
+
+void
+_arb_fmpz_poly_evaluate_arb_rectangular(arb_t y, const fmpz * poly,
+    slong len, const arb_t x, slong prec)
+{
+    slong i, j, m, r;
+    arb_ptr xs;
+    arb_t s, t, c;
+
+    if (len < 3)
+    {
+        _arb_fmpz_poly_evaluate_arb_horner(y, poly, len, x, prec);
+        return;
+    }
+
+    m = n_sqrt(len) + 1;
+    r = (len + m - 1) / m;
+
+    xs = _arb_vec_init(m + 1);
+    arb_init(s);
+    arb_init(t);
+    arb_init(c);
+
+    _arb_vec_set_powers(xs, x, m + 1, prec);
+
+    arb_set_fmpz(y, poly + (r - 1) * m);
+    for (j = 1; (r - 1) * m + j < len; j++)
+        arb_addmul_fmpz(y, xs + j, poly + (r - 1) * m + j, prec);
+
+    for (i = r - 2; i >= 0; i--)
+    {
+        arb_set_fmpz(s, poly + i * m);
+        for (j = 1; j < m; j++)
+            arb_addmul_fmpz(s, xs + j, poly + i * m + j, prec);
+
+        arb_mul(y, y, xs + m, prec);
+        arb_add(y, y, s, prec);
+    }
+
+    _arb_vec_clear(xs, m + 1);
+    arb_clear(s);
+    arb_clear(t);
+    arb_clear(c);
+}
+
+void
+arb_fmpz_poly_evaluate_arb_rectangular(arb_t res, const fmpz_poly_t f, const arb_t a, slong prec)
+{
+    _arb_fmpz_poly_evaluate_arb_rectangular(res, f->coeffs, f->length, a, prec);
+}
+

arb_fmpz_poly/gauss_period_minpoly.c

@@ -0,0 +1,165 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_fmpz_poly.h"
+#include "acb_dirichlet.h"
+
+void
+arb_fmpz_poly_gauss_period_minpoly(fmpz_poly_t res, ulong q, ulong n)
+{
+    ulong k, d, e, g, gk;
+    ulong * es;
+    slong prec, initial_prec;
+    int done, real;
+    int * lower_plane;
+
+    if (n == 0 || !n_is_prime(q) || ((q - 1) % n) != 0 ||
+            n_gcd_full(n, (q - 1) / n) != 1 || q >= ULONG_MAX / 2)
+    {
+        fmpz_poly_zero(res);
+        return;
+    }
+
+    d = (q - 1) / n;
+
+    /* this is much faster */
+    if (d == 1)
+    {
+        fmpz_poly_cyclotomic(res, q);
+        return;
+    }
+
+    g = n_primitive_root_prime(q);
+
+    es = flint_malloc(sizeof(ulong) * d);
+    lower_plane = flint_calloc(n, sizeof(int));
+
+    for (e = 0; e < d; e++)
+        es[e] = n_powmod2(g, n * e, 2 * q);
+
+    /* either all roots are real, or all roots are complex */
+    real = (n % 2) == 1;
+
+    /* first estimate precision crudely based on d and n */
+    initial_prec = n * log(2 * d) * 1.4426950408889 * 1.01 + 20;
+    initial_prec = FLINT_MAX(initial_prec, 48);
+
+    /* if high, start lower to get a good estimate */
+    if (initial_prec > 200)
+        initial_prec = 48;
+
+    for (prec = initial_prec, done = 0; !done; )
+    {
+        acb_dirichlet_roots_t zeta;
+        arb_poly_t pz;
+        arb_ptr roots;
+        acb_ptr croots;
+        acb_t t, u;
+        slong root_index;
+
+        acb_dirichlet_roots_init(zeta, q, n * d, prec);
+        roots = _arb_vec_init(n);
+        croots = (acb_ptr) roots;
+
+        acb_init(t);
+        acb_init(u);
+        arb_poly_init(pz);
+
+        root_index = 0;
+
+        for (k = 0; k < n; k++)
+        {
+            if (lower_plane[k])
+                continue;
+
+            gk = n_powmod2(g, k, 2 * q);
+            acb_zero(u);
+
+            if (real)
+            {
+                for (e = 0; e < d / 2; e++)
+                {
+                    acb_dirichlet_root(t, zeta, n_mulmod2(gk, es[e], 2 * q), prec);
+                    acb_mul_2exp_si(t, t, 1);  /* compute conjugates */
+                    acb_add(u, u, t, prec);
+                }
+
+                arb_set(roots + k, acb_realref(u));
+            }
+            else
+            {
+                for (e = 0; e < d; e++)
+                {
+                    acb_dirichlet_root(t, zeta, n_mulmod2(gk, es[e], 2 * q), prec);
+                    acb_add(u, u, t, prec);
+                }
+
+                if (arb_is_negative(acb_imagref(u)))
+                {
+                    lower_plane[k] = 1;
+                }
+                else if (arb_contains_zero(acb_imagref(u)))
+                {
+                    /* todo: could increase precision */
+                    flint_printf("fail! imaginary part should be nonzero\n");
+                    flint_abort();
+                }
+                else
+                {
+                    acb_set(croots + root_index, u);
+                    root_index++;
+                }
+            }
+        }
+
+        if (real)
+            arb_poly_product_roots(pz, roots, n, prec);
+        else
+            arb_poly_product_roots_complex(pz, NULL, 0, croots, root_index, prec);
+
+        done = arb_poly_get_unique_fmpz_poly(res, pz);
+
+        if (!done && prec == initial_prec)
+        {
+            mag_t m, mmax;
+            mag_init(m);
+            mag_init(mmax);
+
+            for (k = 0; k < n; k++)
+            {
+                arb_get_mag(m, pz->coeffs + k);
+                mag_max(mmax, mmax, m);
+            }
+
+            prec = mag_get_d_log2_approx(mmax) * 1.01 + 20;
+
+            if (prec < 2 * initial_prec)
+                prec = 2 * initial_prec;
+
+            mag_clear(m);
+            mag_clear(mmax);
+        }
+        else if (!done)
+        {
+            prec *= 2;
+        }
+
+        acb_dirichlet_roots_clear(zeta);
+        _arb_vec_clear(roots, n);
+        acb_clear(t);
+        acb_clear(u);
+        arb_poly_clear(pz);
+    }
+
+    flint_free(es);
+    flint_free(lower_plane);
+}
+

arb_fmpz_poly/test/t-evaluate_acb.c

@@ -0,0 +1,111 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_fmpz_poly.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("evaluate_acb....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 3000 * arb_test_multiplier(); iter++)
+    {
+        fmpz_poly_t f, g, h;
+        acb_t x, fx, gx, hx, fxgx;
+        slong prec1, prec2, prec3;
+
+        fmpz_poly_init(f);
+        fmpz_poly_init(g);
+        fmpz_poly_init(h);
+        acb_init(x);
+        acb_init(fx);
+        acb_init(gx);
+        acb_init(hx);
+        acb_init(fxgx);
+
+        fmpz_poly_randtest(f, state, 1 + n_randint(state, 50), 1 + n_randint(state, 1000));
+        fmpz_poly_randtest(g, state, 1 + n_randint(state, 50), 1 + n_randint(state, 1000));
+        fmpz_poly_add(h, f, g);
+
+        acb_randtest(x, state, 1 + n_randint(state, 2000), 1 + n_randint(state, 100));
+        acb_randtest(fx, state, 1 + n_randint(state, 2000), 1 + n_randint(state, 100));
+        acb_randtest(gx, state, 1 + n_randint(state, 2000), 1 + n_randint(state, 100));
+        acb_randtest(hx, state, 1 + n_randint(state, 2000), 1 + n_randint(state, 100));
+        acb_randtest(fxgx, state, 1 + n_randint(state, 2000), 1 + n_randint(state, 100));
+
+        prec1 = 2 + n_randint(state, 2000);
+        prec2 = 2 + n_randint(state, 2000);
+        prec3 = 2 + n_randint(state, 2000);
+
+        switch (n_randint(state, 6))
+        {
+            case 0:
+                arb_fmpz_poly_evaluate_acb_horner(fx, f, x, prec1);
+                break;
+            case 1:
+                arb_fmpz_poly_evaluate_acb_rectangular(fx, f, x, prec1);
+                break;
+            case 2:
+                arb_fmpz_poly_evaluate_acb(fx, f, x, prec1);
+                break;
+            case 3:
+                acb_set(fx, x);
+                arb_fmpz_poly_evaluate_acb_horner(fx, f, fx, prec1);
+                break;
+            case 4:
+                acb_set(fx, x);
+                arb_fmpz_poly_evaluate_acb_rectangular(fx, f, fx, prec1);
+                break;
+            default:
+                acb_set(fx, x);
+                arb_fmpz_poly_evaluate_acb(fx, f, fx, prec1);
+                break;
+        }
+
+        arb_fmpz_poly_evaluate_acb(gx, g, x, prec2);
+        arb_fmpz_poly_evaluate_acb(hx, h, x, prec3);
+        acb_add(fxgx, fx, gx, prec3);
+
+        if (!acb_overlaps(fxgx, hx))
+        {
+            flint_printf("FAIL\n");
+            fmpz_poly_print(f); flint_printf("\n\n");
+            fmpz_poly_print(g); flint_printf("\n\n");
+            fmpz_poly_print(h); flint_printf("\n\n");
+            acb_printd(x, 30); flint_printf("\n\n");
+            acb_printd(fx, 30); flint_printf("\n\n");
+            acb_printd(gx, 30); flint_printf("\n\n");
+            acb_printd(hx, 30); flint_printf("\n\n");
+            acb_printd(fxgx, 30); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        fmpz_poly_clear(f);
+        fmpz_poly_clear(g);
+        fmpz_poly_clear(h);
+        acb_clear(x);
+        acb_clear(fx);
+        acb_clear(gx);
+        acb_clear(hx);
+        acb_clear(fxgx);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

arb_fmpz_poly/test/t-evaluate_arb.c

@@ -0,0 +1,111 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_fmpz_poly.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("evaluate_arb....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 3000 * arb_test_multiplier(); iter++)
+    {
+        fmpz_poly_t f, g, h;
+        arb_t x, fx, gx, hx, fxgx;
+        slong prec1, prec2, prec3;
+
+        fmpz_poly_init(f);
+        fmpz_poly_init(g);
+        fmpz_poly_init(h);
+        arb_init(x);
+        arb_init(fx);
+        arb_init(gx);
+        arb_init(hx);
+        arb_init(fxgx);
+
+        fmpz_poly_randtest(f, state, 1 + n_randint(state, 50), 1 + n_randint(state, 1000));
+        fmpz_poly_randtest(g, state, 1 + n_randint(state, 50), 1 + n_randint(state, 1000));
+        fmpz_poly_add(h, f, g);
+
+        arb_randtest(x, state, 1 + n_randint(state, 2000), 1 + n_randint(state, 100));
+        arb_randtest(fx, state, 1 + n_randint(state, 2000), 1 + n_randint(state, 100));
+        arb_randtest(gx, state, 1 + n_randint(state, 2000), 1 + n_randint(state, 100));
+        arb_randtest(hx, state, 1 + n_randint(state, 2000), 1 + n_randint(state, 100));
+        arb_randtest(fxgx, state, 1 + n_randint(state, 2000), 1 + n_randint(state, 100));
+
+        prec1 = 2 + n_randint(state, 2000);
+        prec2 = 2 + n_randint(state, 2000);
+        prec3 = 2 + n_randint(state, 2000);
+
+        switch (n_randint(state, 6))
+        {
+            case 0:
+                arb_fmpz_poly_evaluate_arb_horner(fx, f, x, prec1);
+                break;
+            case 1:
+                arb_fmpz_poly_evaluate_arb_rectangular(fx, f, x, prec1);
+                break;
+            case 2:
+                arb_fmpz_poly_evaluate_arb(fx, f, x, prec1);
+                break;
+            case 3:
+                arb_set(fx, x);
+                arb_fmpz_poly_evaluate_arb_horner(fx, f, fx, prec1);
+                break;
+            case 4:
+                arb_set(fx, x);
+                arb_fmpz_poly_evaluate_arb_rectangular(fx, f, fx, prec1);
+                break;
+            default:
+                arb_set(fx, x);
+                arb_fmpz_poly_evaluate_arb(fx, f, fx, prec1);
+                break;
+        }
+
+        arb_fmpz_poly_evaluate_arb(gx, g, x, prec2);
+        arb_fmpz_poly_evaluate_arb(hx, h, x, prec3);
+        arb_add(fxgx, fx, gx, prec3);
+
+        if (!arb_overlaps(fxgx, hx))
+        {
+            flint_printf("FAIL\n");
+            fmpz_poly_print(f); flint_printf("\n\n");
+            fmpz_poly_print(g); flint_printf("\n\n");
+            fmpz_poly_print(h); flint_printf("\n\n");
+            arb_printd(x, 30); flint_printf("\n\n");
+            arb_printd(fx, 30); flint_printf("\n\n");
+            arb_printd(gx, 30); flint_printf("\n\n");
+            arb_printd(hx, 30); flint_printf("\n\n");
+            arb_printd(fxgx, 30); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        fmpz_poly_clear(f);
+        fmpz_poly_clear(g);
+        fmpz_poly_clear(h);
+        arb_clear(x);
+        arb_clear(fx);
+        arb_clear(gx);
+        arb_clear(hx);
+        arb_clear(fxgx);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

arb_fmpz_poly/test/t-gauss_period_minpoly.c

@@ -0,0 +1,95 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_fmpz_poly.h"
+#include "acb_dirichlet.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("gauss_period_minpoly....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    {
+        slong prec;
+        ulong n, q;
+        fmpz_poly_t pol;
+        fmpz_poly_init(pol);
+
+        for (q = 0; q < 1000; q++)
+        {
+            acb_dirichlet_roots_t zeta;
+            prec = 100 + n_randint(state, 500);
+
+            if (n_is_prime(q))
+                acb_dirichlet_roots_init(zeta, q, 30, prec);
+
+            for (n = 0; n < 1000; n++)
+            {
+                arb_fmpz_poly_gauss_period_minpoly(pol, q, n);
+
+                if (!fmpz_poly_is_zero(pol))
+                {
+                    ulong k, g, gk, e, d;
+                    acb_t t, u;
+
+                    acb_init(t);
+                    acb_init(u);
+                    d = (q - 1) / n;
+                    g = n_primitive_root_prime(q);
+
+                    for (iter = 0; iter < 3; iter++)
+                    {
+                        k = n_randint(state, n);
+                        gk = n_powmod2(g, k, 2 * q);
+                        acb_zero(u);
+
+                        for (e = 0; e < d; e++)
+                        {
+                            acb_dirichlet_root(t, zeta, n_mulmod2(gk, n_powmod2(g, n * e, 2 * q), 2 * q), prec);
+                            acb_add(u, u, t, prec);
+                        }
+
+                        arb_fmpz_poly_evaluate_acb(t, pol, u, prec);
+
+                        if (!acb_contains_zero(t) || fmpz_poly_degree(pol) != n)
+                        {
+                            flint_printf("FAIL\n");
+                            flint_printf("q = %wu, n = %wu, k = %wu\n\n", q, n, k);
+                            fmpz_poly_print(pol); flint_printf("\n\n");
+                            acb_printn(u, 30, 0); flint_printf("\n\n");
+                            acb_printn(t, 30, 0); flint_printf("\n\n");
+                            flint_abort();
+                        }
+                    }
+
+                    acb_clear(t);
+                    acb_clear(u);
+                }
+            }
+
+            if (n_is_prime(q))
+                acb_dirichlet_roots_clear(zeta);
+        }
+
+        fmpz_poly_clear(pol);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

doc/source/arb_fmpz_poly.rst

@@ -0,0 +1,76 @@
+.. _arb_fmpz_poly:
+
+**arb_fmpz_poly.h** -- extra methods for integer polynomials
+=========================================================================================
+
+This module provides methods for FLINT polynomials
+with integer and rational coefficients
+(:type:`fmpz_poly_t`) and (:type:`fmpq_poly_t`)
+requiring use of Arb real or complex numbers.
+
+Some methods output real or complex numbers while others
+use real and complex numbers internally to produce an exact result.
+This module also contains some useful helper functions not specifically
+related to real and complex numbers that are useful elsewhere in Arb.
+
+Note that methods that combine Arb *polynomials* and FLINT polynomials
+are found in the respective Arb polynomial modules, such as
+:func:`arb_poly_set_fmpz_poly`
+and :func:`arb_poly_get_unique_fmpz_poly`.
+
+Evaluation
+-------------------------------------------------------------------------------
+
+.. function:: void _arb_fmpz_poly_evaluate_arb_horner(arb_t res, const fmpz * poly, slong len, const arb_t x, slong prec)
+
+.. function:: void arb_fmpz_poly_evaluate_arb_horner(arb_t res, const fmpz_poly_t poly, const arb_t x, slong prec)
+
+.. function:: void _arb_fmpz_poly_evaluate_arb_rectangular(arb_t res, const fmpz * poly, slong len, const arb_t x, slong prec)
+
+.. function:: void arb_fmpz_poly_evaluate_arb_rectangular(arb_t res, const fmpz_poly_t poly, const arb_t x, slong prec)
+
+.. function:: void _arb_fmpz_poly_evaluate_arb(arb_t res, const fmpz * poly, slong len, const arb_t x, slong prec)
+
+.. function:: void arb_fmpz_poly_evaluate_arb(arb_t res, const fmpz_poly_t poly, const arb_t x, slong prec)
+
+.. function:: void _arb_fmpz_poly_evaluate_acb_horner(acb_t res, const fmpz * poly, slong len, const acb_t x, slong prec)
+
+.. function:: void arb_fmpz_poly_evaluate_acb_horner(acb_t res, const fmpz_poly_t poly, const acb_t x, slong prec)
+
+.. function:: void _arb_fmpz_poly_evaluate_acb_rectangular(acb_t res, const fmpz * poly, slong len, const acb_t x, slong prec)
+
+.. function:: void arb_fmpz_poly_evaluate_acb_rectangular(acb_t res, const fmpz_poly_t poly, const acb_t x, slong prec)
+
+.. function:: void _arb_fmpz_poly_evaluate_acb(acb_t res, const fmpz * poly, slong len, const acb_t x, slong prec)
+
+.. function:: void arb_fmpz_poly_evaluate_acb(acb_t res, const fmpz_poly_t poly, const acb_t x, slong prec)
+
+    Evaluates *poly* (given by a polynomial object or an array with *len* coefficients)
+    at the given real or complex number, respectively using Horner's rule, rectangular
+    splitting, or a default algorithm choice.
+
+Special polynomials
+-------------------------------------------------------------------------------
+
+Note: see also the methods available in FLINT (e.g. for cyclotomic polynomials).
+
+.. function:: void arb_fmpz_poly_gauss_period_minpoly(fmpz_poly_t res, ulong q, ulong n)
+
+    Sets *res* to the minimal polynomial of the Gaussian periods
+    `\sum_{a \in H} \zeta^a` where `\zeta = \exp(2 \pi i / q)`
+    and *H* are the cosets of the subgroups of order `d = (q - 1) / n` of
+    `(\mathbb{Z}/q\mathbb{Z})^{\times}`.
+    The resulting polynomial has degree *n*.
+    When `d = 1`, the result is the cyclotomic polynomial `\Phi_q`.
+
+    The implementation assumes that *q* is prime, and that *n* is a divisor of
+    `q - 1` such that *n* is coprime with *d*. If any condition is not met,
+    *res* is set to the zero polynomial.
+
+    This method provides a fast (in practice) way to
+    construct finite field extensions of prescribed degree.
+    If *q* satisfies the conditions stated above and `(q-1)/f` additionally
+    is coprime with *n*, where *f* is the multiplicative order of *p* mod *q*, then
+    the Gaussian period minimal polynomial is irreducible over
+    `\operatorname{GF}(p)` [CP2005]_.
+

doc/source/index.rst

@@ -121,6 +121,7 @@ on polynomials, without introducing a separate power series type.
 
    arb_poly.rst
    acb_poly.rst
+   arb_fmpz_poly.rst
 
 Matrices
 ::::::::::::::::::::::::::::::::::::