1 /*
2 Copyright (C) 2012 Fredrik Johansson
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "acb_poly.h"
13
14 /* This gives some speedup for small lengths. */
15 static __inline__ void
_acb_poly_rem_2(acb_ptr r,acb_srcptr a,slong al,acb_srcptr b,slong bl,slong prec)16 _acb_poly_rem_2(acb_ptr r, acb_srcptr a, slong al,
17 acb_srcptr b, slong bl, slong prec)
18 {
19 if (al == 2)
20 {
21 acb_mul(r + 0, a + 1, b + 0, prec);
22 acb_sub(r + 0, a + 0, r + 0, prec);
23 }
24 else
25 {
26 _acb_poly_rem(r, a, al, b, bl, prec);
27 }
28 }
29
30 void
_acb_poly_evaluate_vec_fast_precomp(acb_ptr vs,acb_srcptr poly,slong plen,acb_ptr * tree,slong len,slong prec)31 _acb_poly_evaluate_vec_fast_precomp(acb_ptr vs, acb_srcptr poly,
32 slong plen, acb_ptr * tree, slong len, slong prec)
33 {
34 slong height, i, j, pow, left;
35 slong tree_height;
36 slong tlen;
37 acb_ptr t, u, swap, pa, pb, pc;
38
39 /* avoid worrying about some degenerate cases */
40 if (len < 2 || plen < 2)
41 {
42 if (len == 1)
43 {
44 acb_t tmp;
45 acb_init(tmp);
46 acb_neg(tmp, tree[0] + 0);
47 _acb_poly_evaluate(vs + 0, poly, plen, tmp, prec);
48 acb_clear(tmp);
49 }
50 else if (len != 0 && plen == 0)
51 {
52 _acb_vec_zero(vs, len);
53 }
54 else if (len != 0 && plen == 1)
55 {
56 for (i = 0; i < len; i++)
57 acb_set(vs + i, poly + 0);
58 }
59 return;
60 }
61
62 t = _acb_vec_init(len);
63 u = _acb_vec_init(len);
64
65 left = len;
66
67 /* Initial reduction. We allow the polynomial to be larger
68 or smaller than the number of points. */
69 height = FLINT_BIT_COUNT(plen - 1) - 1;
70 tree_height = FLINT_CLOG2(len);
71 while (height >= tree_height)
72 height--;
73 pow = WORD(1) << height;
74
75 for (i = j = 0; i < len; i += pow, j += (pow + 1))
76 {
77 tlen = ((i + pow) <= len) ? pow : len % pow;
78 _acb_poly_rem(t + i, poly, plen, tree[height] + j, tlen + 1, prec);
79 }
80
81 for (i = height - 1; i >= 0; i--)
82 {
83 pow = WORD(1) << i;
84 left = len;
85 pa = tree[i];
86 pb = t;
87 pc = u;
88
89 while (left >= 2 * pow)
90 {
91 _acb_poly_rem_2(pc, pb, 2 * pow, pa, pow + 1, prec);
92 _acb_poly_rem_2(pc + pow, pb, 2 * pow, pa + pow + 1, pow + 1, prec);
93
94 pa += 2 * pow + 2;
95 pb += 2 * pow;
96 pc += 2 * pow;
97 left -= 2 * pow;
98 }
99
100 if (left > pow)
101 {
102 _acb_poly_rem(pc, pb, left, pa, pow + 1, prec);
103 _acb_poly_rem(pc + pow, pb, left, pa + pow + 1, left - pow + 1, prec);
104 }
105 else if (left > 0)
106 _acb_vec_set(pc, pb, left);
107
108 swap = t;
109 t = u;
110 u = swap;
111 }
112
113 _acb_vec_set(vs, t, len);
114 _acb_vec_clear(t, len);
115 _acb_vec_clear(u, len);
116 }
117
_acb_poly_evaluate_vec_fast(acb_ptr ys,acb_srcptr poly,slong plen,acb_srcptr xs,slong n,slong prec)118 void _acb_poly_evaluate_vec_fast(acb_ptr ys, acb_srcptr poly, slong plen,
119 acb_srcptr xs, slong n, slong prec)
120 {
121 acb_ptr * tree;
122
123 tree = _acb_poly_tree_alloc(n);
124 _acb_poly_tree_build(tree, xs, n, prec);
125 _acb_poly_evaluate_vec_fast_precomp(ys, poly, plen, tree, n, prec);
126 _acb_poly_tree_free(tree, n);
127 }
128
129 void
acb_poly_evaluate_vec_fast(acb_ptr ys,const acb_poly_t poly,acb_srcptr xs,slong n,slong prec)130 acb_poly_evaluate_vec_fast(acb_ptr ys,
131 const acb_poly_t poly, acb_srcptr xs, slong n, slong prec)
132 {
133 _acb_poly_evaluate_vec_fast(ys, poly->coeffs,
134 poly->length, xs, n, prec);
135 }
136