_acb_vec_nth_roots: simpler, better algorithm

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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;
+}
+

acb/vec_nth_root.c

@@ -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 <http://www.gnu.org/licenses/>.
-*/
-
-#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;
-   }
-}

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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);
+}
+