diff --git a/acb/test/t-vec_nth_roots.c b/acb/test/t-vec_nth_roots.c new file mode 100644 index 00000000..f278d04f --- /dev/null +++ b/acb/test/t-vec_nth_roots.c @@ -0,0 +1,63 @@ +/* + Copyright (C) 2016 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 "acb.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("vec_nth_roots...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100; iter++) + { + acb_ptr vec; + fmpq_t q; + acb_t t; + slong k; + slong prec = 53; + + vec = _acb_vec_init(iter); + acb_init(t); + fmpq_init(q); + + _acb_vec_nth_roots(vec, iter, prec); + + for (k = 0; k < iter; k++) + { + fmpq_set_si(q, 2 * k, iter); + arb_sin_cos_pi_fmpq(acb_imagref(t), acb_realref(t), q, prec); + + if (!acb_overlaps(vec + k, t)) + { + flint_printf("FAIL: overlap\n\n"); + flint_printf("n = %wu k = %wd\n\n", iter, k); + flint_printf("vec = "); acb_printn(vec + k, 30, 0); flint_printf("\n\n"); + flint_printf("t = "); acb_printn(t, 30, 0); flint_printf("\n\n"); + abort(); + } + } + + _acb_vec_clear(vec, iter); + acb_clear(t); + fmpq_clear(q); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/acb/vec_nth_root.c b/acb/vec_nth_root.c deleted file mode 100644 index 10b22621..00000000 --- a/acb/vec_nth_root.c +++ /dev/null @@ -1,91 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - - 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 "acb_dirichlet.h" - -/* assume xs has size >= len * step */ -void -_acb_vec_set_powers_step(acb_ptr xs, slong n, slong len, slong step, slong prec) -{ - slong i, j; - prec += n_clog(n, 2); - - for (i = 0, j = 0; i < len; i++, j += step) - { - if (i == 0) - acb_one(xs + j); - else if (i == 1) - acb_nth_root(xs + j, n, prec); - else if (i % 2 == 0) - acb_mul(xs + j, xs + j / 2, xs + j / 2, prec); - else - acb_mul(xs + j, xs + j - step, xs + step, prec); - } -} - -/* assume len = p^e and z has size >= len * step */ -void -_acb_vec_nth_roots_pe(acb_ptr z, slong p, slong e, slong len, slong step, slong prec) -{ - if (e <= 1) - { - _acb_vec_set_powers_step(z, p, len, step, prec); - } - else - { - slong q, r; - _acb_vec_nth_roots_pe(z, p, e - 1, len / p, step * p, prec); - _acb_vec_set_powers_step(z, n_pow(p, e), p, step, prec); - - for (q = p; q < len; q += p) - for (r = 1; r < p; r++) - acb_mul(z + (q + r) * step, z + q * step, z + r * step, prec); - } -} - -void -_acb_vec_nth_roots(acb_ptr z, slong len, slong prec) -{ - slong i, q; - n_factor_t fac; - - acb_one(z + 0); - - n_factor_init(&fac); - n_factor(&fac, len, 0); - q = 1; - - for (i = 0; i < fac.num; i++) - { - slong p, e, pe, mp, mq; - p = fac.p[i]; - e = fac.exp[i]; - pe = n_pow(p, e); - mp = len / pe; - mq = len / q; - - _acb_vec_nth_roots_pe(z, p, e, pe, mp, prec); - - if (i > 0) - { - slong k, l; - - for (k = mp; k < len; k += mp) - { - for (l = mq; l < len - k; l += mq) - acb_mul(z + k + l, z + k, z + l, prec); - for (; l < len; l += mq) - acb_mul(z + k + l - len, z + k, z + l, prec); - } - } - q *= pe; - } -} diff --git a/acb/vec_nth_roots.c b/acb/vec_nth_roots.c new file mode 100644 index 00000000..79d28f87 --- /dev/null +++ b/acb/vec_nth_roots.c @@ -0,0 +1,62 @@ +/* + Copyright (C) 2016 Pascal Molin + Copyright (C) 2016 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 "acb_dirichlet.h" + +void +_acb_vec_nth_roots(acb_ptr res, slong n, slong prec) +{ + slong k, len, wp; + acb_t t; + + if (n <= 0) + return; + + if (n % 4 == 0) + len = n / 8 + 1; + else if (n % 2 == 0) + len = n / 4 + 1; + else + len = n / 2 + 1; + + wp = prec + 6 + 2 * FLINT_BIT_COUNT(len); + + acb_init(t); + acb_nth_root(t, n, wp); + _acb_vec_set_powers(res, t, len, wp); + acb_clear(t); + _acb_vec_set_round(res, res, len, prec); + + if (n % 4 == 0) + { + for (k = n / 8 + 1; k <= n / 4; k++) + { + arb_set(acb_realref(res + k), acb_imagref(res + n / 4 - k)); + arb_set(acb_imagref(res + k), acb_realref(res + n / 4 - k)); + } + + for (k = n / 4 + 1; k <= n / 2; k++) + acb_mul_onei(res + k, res + k - n / 4); + } + else if (n % 2 == 0) + { + for (k = n / 4 + 1; k <= n / 2; k++) + { + acb_set(res + k, res + n / 2 - k); + arb_neg(acb_realref(res + k), acb_realref(res + k)); + } + } + + for (k = n / 2 + 1; k < n; k++) + acb_conj(res + k, res + n - k); +} +