1 /*
2     Copyright (C) 2016 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 void
_acb_poly_taylor_shift_convolution(acb_ptr p,const acb_t c,slong len,slong prec)15 _acb_poly_taylor_shift_convolution(acb_ptr p, const acb_t c, slong len, slong prec)
16 {
17     slong i, n = len - 1;
18     acb_t d;
19     arb_t f;
20     acb_ptr t, u;
21 
22     if (acb_is_zero(c) || len <= 1)
23         return;
24 
25     t = _acb_vec_init(len);
26     u = _acb_vec_init(len);
27 
28     arb_init(f);
29     acb_init(d);
30 
31     arb_one(f);
32     for (i = 2; i <= n; i++)
33     {
34         arb_mul_ui(f, f, i, prec);
35         acb_mul_arb(p + i, p + i, f, prec);
36     }
37 
38     _acb_poly_reverse(p, p, len, len);
39 
40     acb_one(t + n);
41     for (i = n; i > 0; i--)
42         acb_mul_ui(t + i - 1, t + i, i, prec);
43 
44     if (acb_equal_si(c, -1))
45     {
46         for (i = 1; i <= n; i += 2)
47             acb_neg(t + i, t + i);
48     }
49     else if (!acb_is_one(c))
50     {
51         acb_set(d, c);
52 
53         for (i = 1; i <= n; i++)
54         {
55             acb_mul(t + i, t + i, d, prec);
56             acb_mul(d, d, c, prec);
57         }
58     }
59 
60     _acb_poly_mullow(u, p, len, t, len, len, prec);
61 
62     arb_mul(f, f, f, prec);
63 
64     if (arb_bits(f) > 0.25 * prec)
65     {
66         arb_inv(f, f, prec);
67     }
68     else
69     {
70         for (i = 0; i <= n; i++)
71             acb_div_arb(u + i, u + i, f, prec);
72 
73         arb_one(f);
74     }
75 
76     for (i = n; i >= 0; i--)
77     {
78         acb_mul_arb(p + i, u + n - i, f, prec);
79         arb_mul_ui(f, f, (i == 0) ? 1 : i, prec);
80     }
81 
82     _acb_vec_clear(t, len);
83     _acb_vec_clear(u, len);
84 
85     arb_clear(f);
86     acb_clear(d);
87 }
88 
89 void
acb_poly_taylor_shift_convolution(acb_poly_t g,const acb_poly_t f,const acb_t c,slong prec)90 acb_poly_taylor_shift_convolution(acb_poly_t g, const acb_poly_t f,
91         const acb_t c, slong prec)
92 {
93     if (f != g)
94         acb_poly_set_round(g, f, prec);
95 
96     _acb_poly_taylor_shift_convolution(g->coeffs, c, g->length, prec);
97 }
98 
99