1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
|
/*
Copyright (C) 2017 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "arb_poly.h"
void
_arb_poly_product_roots_complex(arb_ptr poly, arb_srcptr r, slong rn,
acb_srcptr c, slong cn, slong prec)
{
if (rn == 0 && cn == 0)
{
arb_one(poly);
}
else if (rn == 1 && cn == 0)
{
arb_neg(poly, r);
arb_one(poly + 1);
}
else if (rn == 2 && cn == 0)
{
arb_mul(poly, r + 0, r + 1, prec);
arb_add(poly + 1, r + 0, r + 1, prec);
arb_neg(poly + 1, poly + 1);
arb_one(poly + 2);
}
else if (rn == 3 && cn == 0)
{
arb_mul(poly + 1, r, r + 1, prec);
arb_mul(poly, poly + 1, r + 2, prec);
arb_neg(poly, poly);
arb_add(poly + 2, r, r + 1, prec);
arb_addmul(poly + 1, poly + 2, r + 2, prec);
arb_add(poly + 2, poly + 2, r + 2, prec);
arb_neg(poly + 2, poly + 2);
arb_one(poly + 3);
}
else if (rn == 0 && cn == 1)
{
arb_mul(poly, acb_realref(c), acb_realref(c), prec);
arb_addmul(poly, acb_imagref(c), acb_imagref(c), prec);
arb_mul_2exp_si(poly + 1, acb_realref(c), 1);
arb_neg(poly + 1, poly + 1);
arb_one(poly + 2);
}
else if (rn == 1 && cn == 1)
{
arb_mul(poly + 1, acb_realref(c), acb_realref(c), prec);
arb_addmul(poly + 1, acb_imagref(c), acb_imagref(c), prec);
arb_mul(poly, poly + 1, r, prec);
arb_neg(poly, poly);
arb_mul_2exp_si(poly + 2, acb_realref(c), 1);
arb_addmul(poly + 1, poly + 2, r, prec);
arb_add(poly + 2, poly + 2, r, prec);
arb_neg(poly + 2, poly + 2);
arb_one(poly + 3);
}
else
{
slong rm, cm, rm2, cm2;
arb_ptr tmp, tmp2;
rm = (rn + 1) / 2;
cm = cn / 2;
rm2 = rn - rm;
cm2 = cn - cm;
tmp = _arb_vec_init(rn + 2 * cn + 2);
tmp2 = tmp + rm + (2 * cm) + 1;
_arb_poly_product_roots_complex(tmp, r, rm, c, cm, prec);
_arb_poly_product_roots_complex(tmp2, r + rm, rm2, c + cm, cm2, prec);
_arb_poly_mul_monic(poly, tmp, rm + 2 * cm + 1, tmp2, rm2 + 2 * cm2 + 1, prec);
_arb_vec_clear(tmp, rn + 2 * cn + 2);
}
}
void
arb_poly_product_roots_complex(arb_poly_t poly,
arb_srcptr r, slong rn, acb_srcptr c, slong cn, slong prec)
{
arb_poly_fit_length(poly, rn + 2 * cn + 1);
_arb_poly_product_roots_complex(poly->coeffs, r, rn, c, cn, prec);
_arb_poly_set_length(poly, rn + 2 * cn + 1);
}
|
