1 /*
2 Copyright (C) 2016 Pascal Molin
3 Copyright (C) 2016 Fredrik Johansson
4
5 This file is part of Arb.
6
7 Arb is free software: you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License (LGPL) as published
9 by the Free Software Foundation; either version 2.1 of the License, or
10 (at your option) any later version. See <http://www.gnu.org/licenses/>.
11 */
12
13 #include "acb.h"
14
15 void
_acb_vec_unit_roots(acb_ptr res,slong n,slong len,slong prec)16 _acb_vec_unit_roots(acb_ptr res, slong n, slong len, slong prec)
17 {
18 int conj = 0;
19 slong k, len1, wp;
20 acb_t t;
21
22 if (len <= 0)
23 return;
24 if (n == 0)
25 {
26 flint_printf("\n_acb_vec_unit_roots: need order != 0\n");
27 abort();
28 }
29
30 if (n < 0)
31 {
32 n = -n;
33 conj = 1;
34 }
35
36 if (n % 4 == 0)
37 len1 = FLINT_MIN(len, n / 8 + 1);
38 else if (n % 2 == 0)
39 len1 = FLINT_MIN(len, n / 4 + 1);
40 else
41 len1 = FLINT_MIN(len, n / 2 + 1);
42
43 wp = prec + 6 + 2 * FLINT_BIT_COUNT(len1);
44
45 acb_init(t);
46 acb_unit_root(t, n, prec);
47 _acb_vec_set_powers(res, t, len1, wp);
48 acb_clear(t);
49 _acb_vec_set_round(res, res, len1, prec);
50
51 if (n % 4 == 0)
52 {
53 for (k = n / 8 + 1; k <= n / 4 && k < len; k++)
54 {
55 arb_set(acb_realref(res + k), acb_imagref(res + n / 4 - k));
56 arb_set(acb_imagref(res + k), acb_realref(res + n / 4 - k));
57 }
58
59 for (k = n / 4 + 1; k <= n / 2 && k < len; k++)
60 acb_mul_onei(res + k, res + k - n / 4);
61 }
62 else if (n % 2 == 0)
63 {
64 for (k = n / 4 + 1; k <= n / 2 && k < len; k++)
65 {
66 acb_set(res + k, res + n / 2 - k);
67 arb_neg(acb_realref(res + k), acb_realref(res + k));
68 }
69 }
70
71 for (k = n / 2 + 1; k < len && k < n; k++)
72 acb_conj(res + k, res + n - k);
73
74 for (k = n; k < len; k++)
75 acb_set(res + k, res - n + k);
76
77 if (conj)
78 {
79 for (k = 1; k < len; k++)
80 acb_conj(res + k, res + k);
81 }
82
83 }
84