From 84c79a3a1049a9e596541c489dd93e9954bbe426 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Fri, 23 Jun 2017 13:08:45 +0200 Subject: [PATCH] add arb_fmpz_poly module with evaluation, plus a function to construct minimal polynomials of Gaussian periods --- Makefile.in | 1 + arb_fmpz_poly.h | 45 ++++++ arb_fmpz_poly/evaluate_acb.c | 34 ++++ arb_fmpz_poly/evaluate_acb_horner.c | 59 +++++++ arb_fmpz_poly/evaluate_acb_rectangular.c | 63 ++++++++ arb_fmpz_poly/evaluate_arb.c | 39 +++++ arb_fmpz_poly/evaluate_arb_horner.c | 59 +++++++ arb_fmpz_poly/evaluate_arb_rectangular.c | 63 ++++++++ arb_fmpz_poly/gauss_period_minpoly.c | 165 ++++++++++++++++++++ arb_fmpz_poly/test/t-evaluate_acb.c | 111 +++++++++++++ arb_fmpz_poly/test/t-evaluate_arb.c | 111 +++++++++++++ arb_fmpz_poly/test/t-gauss_period_minpoly.c | 95 +++++++++++ doc/source/arb_fmpz_poly.rst | 76 +++++++++ doc/source/index.rst | 1 + 14 files changed, 922 insertions(+) create mode 100644 arb_fmpz_poly.h create mode 100644 arb_fmpz_poly/evaluate_acb.c create mode 100644 arb_fmpz_poly/evaluate_acb_horner.c create mode 100644 arb_fmpz_poly/evaluate_acb_rectangular.c create mode 100644 arb_fmpz_poly/evaluate_arb.c create mode 100644 arb_fmpz_poly/evaluate_arb_horner.c create mode 100644 arb_fmpz_poly/evaluate_arb_rectangular.c create mode 100644 arb_fmpz_poly/gauss_period_minpoly.c create mode 100644 arb_fmpz_poly/test/t-evaluate_acb.c create mode 100644 arb_fmpz_poly/test/t-evaluate_arb.c create mode 100644 arb_fmpz_poly/test/t-gauss_period_minpoly.c create mode 100644 doc/source/arb_fmpz_poly.rst diff --git a/Makefile.in b/Makefile.in index 1f8d5a17..7ff5b63d 100644 --- a/Makefile.in +++ b/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 = diff --git a/arb_fmpz_poly.h b/arb_fmpz_poly.h new file mode 100644 index 00000000..4bcd1c8f --- /dev/null +++ b/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 . +*/ + +#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 diff --git a/arb_fmpz_poly/evaluate_acb.c b/arb_fmpz_poly/evaluate_acb.c new file mode 100644 index 00000000..442ff324 --- /dev/null +++ b/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 . +*/ + +#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); +} + diff --git a/arb_fmpz_poly/evaluate_acb_horner.c b/arb_fmpz_poly/evaluate_acb_horner.c new file mode 100644 index 00000000..fc43ba27 --- /dev/null +++ b/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 . +*/ + +#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); +} + diff --git a/arb_fmpz_poly/evaluate_acb_rectangular.c b/arb_fmpz_poly/evaluate_acb_rectangular.c new file mode 100644 index 00000000..b5147f78 --- /dev/null +++ b/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 . +*/ + +#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); +} + diff --git a/arb_fmpz_poly/evaluate_arb.c b/arb_fmpz_poly/evaluate_arb.c new file mode 100644 index 00000000..80fa2d33 --- /dev/null +++ b/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 . +*/ + +#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); +} + diff --git a/arb_fmpz_poly/evaluate_arb_horner.c b/arb_fmpz_poly/evaluate_arb_horner.c new file mode 100644 index 00000000..8e3e8d02 --- /dev/null +++ b/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 . +*/ + +#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); +} + diff --git a/arb_fmpz_poly/evaluate_arb_rectangular.c b/arb_fmpz_poly/evaluate_arb_rectangular.c new file mode 100644 index 00000000..dfd9f4d3 --- /dev/null +++ b/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 . +*/ + +#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); +} + diff --git a/arb_fmpz_poly/gauss_period_minpoly.c b/arb_fmpz_poly/gauss_period_minpoly.c new file mode 100644 index 00000000..67d0611b --- /dev/null +++ b/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 . +*/ + +#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); +} + diff --git a/arb_fmpz_poly/test/t-evaluate_acb.c b/arb_fmpz_poly/test/t-evaluate_acb.c new file mode 100644 index 00000000..d18abc11 --- /dev/null +++ b/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 . +*/ + +#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; +} + diff --git a/arb_fmpz_poly/test/t-evaluate_arb.c b/arb_fmpz_poly/test/t-evaluate_arb.c new file mode 100644 index 00000000..64e622f6 --- /dev/null +++ b/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 . +*/ + +#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; +} + diff --git a/arb_fmpz_poly/test/t-gauss_period_minpoly.c b/arb_fmpz_poly/test/t-gauss_period_minpoly.c new file mode 100644 index 00000000..a387850d --- /dev/null +++ b/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 . +*/ + +#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; +} + diff --git a/doc/source/arb_fmpz_poly.rst b/doc/source/arb_fmpz_poly.rst new file mode 100644 index 00000000..af8034db --- /dev/null +++ b/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]_. + diff --git a/doc/source/index.rst b/doc/source/index.rst index bab1c7e3..f16bcd5e 100644 --- a/doc/source/index.rst +++ b/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 ::::::::::::::::::::::::::::::::::::