1 /*
2 Copyright (C) 2015, Elias Tsigaridas
3 Copyright (C) 2016, Vincent Delecroix
4
5 This file is part of FLINT.
6
7 The implementation follows the function slv_poly_root_upper_bound_2exp in
8 the library SLV version 0.5 by Elias Tsigaridas (in the file poly_ops.c
9 lines 67-125).
10
11 FLINT is free software: you can redistribute it and/or modify it under
12 the terms of the GNU Lesser General Public License (LGPL) as published
13 by the Free Software Foundation; either version 2.1 of the License, or
14 (at your option) any later version. See <http://www.gnu.org/licenses/>.
15 */
16
17 #include <assert.h>
18
19 #include <arb.h>
20 #include <flint/fmpz_poly.h>
21 #include "../e-antic/fmpz_poly_extra.h"
22
_fmpz_poly_positive_root_upper_bound_2exp_local_max(const fmpz * pol,slong len)23 slong _fmpz_poly_positive_root_upper_bound_2exp_local_max(const fmpz * pol, slong len)
24 {
25 slong b, b0, bmin;
26 slong i, j, jmin = -1;
27 fmpz_t tmp;
28
29 fmpz_init(tmp);
30
31 assert(len >= 0 && "len must be non-negative");
32 slong * coeffs = flint_malloc((ulong)len * sizeof(slong));
33 for (i = 0; i < len; i++)
34 coeffs[i] = 1;
35
36 b0 = WORD_MIN;
37 int sgn = fmpz_sgn(pol + len - 1);
38 for (i = len - 2; i >= 0; i--)
39 {
40 if (fmpz_sgn(pol + i) == 0 || fmpz_sgn(pol + i) == sgn) continue;
41
42 bmin = WORD_MAX;
43 for (j = i + 1; j < len; j++)
44 {
45 /* compare the current bound with the log (in base 2) of */
46 /* (- 2^coeffs[j] * p_i / p_j) ^ (1 / (j - i)) */
47 /* which equals */
48 /* (coeffs[j] + log(|p_i|) + log(|p_j|)) / (j - i) */
49 b = coeffs[j];
50
51 fmpz_set(tmp, pol + i);
52 fmpz_abs(tmp, tmp);
53 b += fmpz_clog_ui(tmp, 2);
54
55 fmpz_set(tmp, pol + j);
56 fmpz_abs(tmp, tmp);
57 b -= fmpz_flog_ui(tmp, 2);
58
59 b = (b + j - i - 1) / (j - i);
60
61 if (b < bmin)
62 {
63 jmin = j;
64 bmin = b;
65 if (bmin < b0)
66 break;
67 }
68 }
69
70 b0 = FLINT_MAX(b0, bmin);
71
72 assert(jmin >= 0);
73 coeffs[jmin] ++;
74 }
75
76 fmpz_clear(tmp);
77 flint_free(coeffs);
78
79 return b0;
80 }
81
_fmpz_poly_positive_root_upper_bound_2exp(const fmpz * pol,slong len)82 slong _fmpz_poly_positive_root_upper_bound_2exp(const fmpz * pol, slong len)
83 {
84 return _fmpz_poly_positive_root_upper_bound_2exp_local_max(pol, len);
85 }
86
fmpz_poly_positive_root_upper_bound_2exp(const fmpz_poly_t pol)87 slong fmpz_poly_positive_root_upper_bound_2exp(const fmpz_poly_t pol)
88 {
89 slong i0;
90
91 if (fmpz_poly_is_zero(pol)) return 0;
92
93 i0 = 0;
94 while (fmpz_is_zero(pol->coeffs + i0)) i0++;
95
96 return _fmpz_poly_positive_root_upper_bound_2exp_local_max(
97 pol->coeffs + i0,
98 fmpz_poly_length(pol) - i0);
99 }
100