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 void
_acb_poly_compose_series_horner(acb_ptr res,acb_srcptr poly1,slong len1,acb_srcptr poly2,slong len2,slong n,slong prec)15 _acb_poly_compose_series_horner(acb_ptr res, acb_srcptr poly1, slong len1,
16 acb_srcptr poly2, slong len2, slong n, slong prec)
17 {
18 if (n == 1)
19 {
20 acb_set(res, poly1);
21 }
22 else
23 {
24 slong i = len1 - 1;
25 slong lenr;
26
27 acb_ptr t = _acb_vec_init(n);
28
29 lenr = len2;
30 _acb_vec_scalar_mul(res, poly2, len2, poly1 + i, prec);
31 i--;
32 acb_add(res, res, poly1 + i, prec);
33
34 while (i > 0)
35 {
36 i--;
37 if (lenr + len2 - 1 < n)
38 {
39 _acb_poly_mul(t, res, lenr, poly2, len2, prec);
40 lenr = lenr + len2 - 1;
41 }
42 else
43 {
44 _acb_poly_mullow(t, res, lenr, poly2, len2, n, prec);
45 lenr = n;
46 }
47 _acb_poly_add(res, t, lenr, poly1 + i, 1, prec);
48 }
49
50 _acb_vec_zero(res + lenr, n - lenr);
51 _acb_vec_clear(t, n);
52 }
53 }
54
55 void
acb_poly_compose_series_horner(acb_poly_t res,const acb_poly_t poly1,const acb_poly_t poly2,slong n,slong prec)56 acb_poly_compose_series_horner(acb_poly_t res,
57 const acb_poly_t poly1,
58 const acb_poly_t poly2, slong n, slong prec)
59 {
60 slong len1 = poly1->length;
61 slong len2 = poly2->length;
62 slong lenr;
63
64 if (len2 != 0 && !acb_is_zero(poly2->coeffs))
65 {
66 flint_printf("exception: compose_series: inner "
67 "polynomial must have zero constant term\n");
68 flint_abort();
69 }
70
71 if (len1 == 0 || n == 0)
72 {
73 acb_poly_zero(res);
74 return;
75 }
76
77 if (len2 == 0 || len1 == 1)
78 {
79 acb_poly_set_acb(res, poly1->coeffs);
80 return;
81 }
82
83 lenr = FLINT_MIN((len1 - 1) * (len2 - 1) + 1, n);
84 len1 = FLINT_MIN(len1, lenr);
85 len2 = FLINT_MIN(len2, lenr);
86
87 if ((res != poly1) && (res != poly2))
88 {
89 acb_poly_fit_length(res, lenr);
90 _acb_poly_compose_series_horner(res->coeffs, poly1->coeffs, len1,
91 poly2->coeffs, len2, lenr, prec);
92 _acb_poly_set_length(res, lenr);
93 _acb_poly_normalise(res);
94 }
95 else
96 {
97 acb_poly_t t;
98 acb_poly_init2(t, lenr);
99 _acb_poly_compose_series_horner(t->coeffs, poly1->coeffs, len1,
100 poly2->coeffs, len2, lenr, prec);
101 _acb_poly_set_length(t, lenr);
102 _acb_poly_normalise(t);
103 acb_poly_swap(res, t);
104 acb_poly_clear(t);
105 }
106 }
107