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_mullow_transpose_gauss(acb_ptr res,acb_srcptr poly1,slong len1,acb_srcptr poly2,slong len2,slong n,slong prec)15 _acb_poly_mullow_transpose_gauss(acb_ptr res,
16 acb_srcptr poly1, slong len1,
17 acb_srcptr poly2, slong len2, slong n, slong prec)
18 {
19 arb_ptr a, b, c, d, e, f, w;
20 arb_ptr t, u, v;
21 slong i;
22
23 len1 = FLINT_MIN(len1, n);
24 len2 = FLINT_MIN(len2, n);
25
26 w = flint_malloc(sizeof(arb_struct) * (2 * (len1 + len2 + n)));
27 a = w;
28 b = a + len1;
29 c = b + len1;
30 d = c + len2;
31 e = d + len2;
32 f = e + n;
33
34 t = _arb_vec_init(n);
35 u = _arb_vec_init(n);
36 v = _arb_vec_init(n);
37
38 for (i = 0; i < len1; i++)
39 {
40 a[i] = *acb_realref(poly1 + i);
41 b[i] = *acb_imagref(poly1 + i);
42 }
43
44 for (i = 0; i < len2; i++)
45 {
46 c[i] = *acb_realref(poly2 + i);
47 d[i] = *acb_imagref(poly2 + i);
48 }
49
50 for (i = 0; i < n; i++)
51 {
52 e[i] = *acb_realref(res + i);
53 f[i] = *acb_imagref(res + i);
54 }
55
56 _arb_vec_add(t, a, b, len1, prec);
57 _arb_vec_add(u, c, d, len2, prec);
58
59 _arb_poly_mullow(v, t, len1, u, len2, n, prec);
60 _arb_poly_mullow(t, a, len1, c, len2, n, prec);
61 _arb_poly_mullow(u, b, len1, d, len2, n, prec);
62
63 _arb_vec_sub(e, t, u, n, prec);
64 _arb_vec_sub(f, v, t, n, prec);
65 _arb_vec_sub(f, f, u, n, prec);
66
67 for (i = 0; i < n; i++)
68 {
69 *acb_realref(res + i) = e[i];
70 *acb_imagref(res + i) = f[i];
71 }
72
73 _arb_vec_clear(t, n);
74 _arb_vec_clear(u, n);
75 _arb_vec_clear(v, n);
76
77 flint_free(w);
78 }
79
80 void
acb_poly_mullow_transpose_gauss(acb_poly_t res,const acb_poly_t poly1,const acb_poly_t poly2,slong n,slong prec)81 acb_poly_mullow_transpose_gauss(acb_poly_t res, const acb_poly_t poly1,
82 const acb_poly_t poly2,
83 slong n, slong prec)
84 {
85 slong len1, len2;
86
87 len1 = poly1->length;
88 len2 = poly2->length;
89
90 if (len1 == 0 || len2 == 0 || n == 0)
91 {
92 acb_poly_zero(res);
93 return;
94 }
95
96 n = FLINT_MIN((len1 + len2 - 1), n);
97 len1 = FLINT_MIN(len1, n);
98 len2 = FLINT_MIN(len2, n);
99
100 if (res == poly1 || res == poly2)
101 {
102 acb_poly_t t;
103 acb_poly_init2(t, n);
104 _acb_poly_mullow_transpose_gauss(t->coeffs, poly1->coeffs, len1,
105 poly2->coeffs, len2, n, prec);
106 acb_poly_swap(res, t);
107 acb_poly_clear(t);
108 }
109 else
110 {
111 acb_poly_fit_length(res, n);
112 _acb_poly_mullow_transpose_gauss(res->coeffs, poly1->coeffs, len1,
113 poly2->coeffs, len2, n, prec);
114 }
115
116 _acb_poly_set_length(res, n);
117 _acb_poly_normalise(res);
118 }
119
120