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(acb_ptr res,acb_srcptr poly1,slong len1,acb_srcptr poly2,slong len2,slong n,slong prec)15 _acb_poly_mullow_transpose(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;
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 /* (e+fi) = (a+bi)(c+di) = (ac - bd) + (ad + bc)i */
35 t = _arb_vec_init(n);
36
37 for (i = 0; i < len1; i++)
38 {
39 a[i] = *acb_realref(poly1 + i);
40 b[i] = *acb_imagref(poly1 + i);
41 }
42
43 for (i = 0; i < len2; i++)
44 {
45 c[i] = *acb_realref(poly2 + i);
46 d[i] = *acb_imagref(poly2 + i);
47 }
48
49 for (i = 0; i < n; i++)
50 {
51 e[i] = *acb_realref(res + i);
52 f[i] = *acb_imagref(res + i);
53 }
54
55 _arb_poly_mullow(e, a, len1, c, len2, n, prec);
56 _arb_poly_mullow(t, b, len1, d, len2, n, prec);
57 _arb_vec_sub(e, e, t, n, prec);
58
59 _arb_poly_mullow(f, a, len1, d, len2, n, prec);
60 /* squaring */
61 if (poly1 == poly2 && len1 == len2)
62 {
63 _arb_vec_scalar_mul_2exp_si(f, f, n, 1);
64 }
65 else
66 {
67 _arb_poly_mullow(t, b, len1, c, len2, n, prec);
68 _arb_vec_add(f, f, t, n, prec);
69 }
70
71 for (i = 0; i < n; i++)
72 {
73 *acb_realref(res + i) = e[i];
74 *acb_imagref(res + i) = f[i];
75 }
76
77 _arb_vec_clear(t, n);
78 flint_free(w);
79 }
80
81 void
acb_poly_mullow_transpose(acb_poly_t res,const acb_poly_t poly1,const acb_poly_t poly2,slong n,slong prec)82 acb_poly_mullow_transpose(acb_poly_t res, const acb_poly_t poly1,
83 const acb_poly_t poly2,
84 slong n, slong prec)
85 {
86 slong len1, len2;
87
88 len1 = poly1->length;
89 len2 = poly2->length;
90
91 if (len1 == 0 || len2 == 0 || n == 0)
92 {
93 acb_poly_zero(res);
94 return;
95 }
96
97 n = FLINT_MIN((len1 + len2 - 1), n);
98 len1 = FLINT_MIN(len1, n);
99 len2 = FLINT_MIN(len2, n);
100
101 if (res == poly1 || res == poly2)
102 {
103 acb_poly_t t;
104 acb_poly_init2(t, n);
105 _acb_poly_mullow_transpose(t->coeffs, poly1->coeffs, len1,
106 poly2->coeffs, len2, n, prec);
107 acb_poly_swap(res, t);
108 acb_poly_clear(t);
109 }
110 else
111 {
112 acb_poly_fit_length(res, n);
113 _acb_poly_mullow_transpose(res->coeffs, poly1->coeffs, len1,
114 poly2->coeffs, len2, n, prec);
115 }
116
117 _acb_poly_set_length(res, n);
118 _acb_poly_normalise(res);
119 }
120
121