hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 09:21:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

add arb_fmpz_poly_complex_roots

Browse source
This commit is contained in:
Fredrik Johansson 2017-06-23 15:14:30 +02:00
parent 84c79a3a10
commit 7d0f75f38e
7 changed files with 541 additions and 208 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

7
arb_fmpz_poly.h
View file

@ -26,6 +26,8 @@
#include "arb_poly.h"
#include "acb_poly.h"
#define ARB_FMPZ_POLY_ROOTS_VERBOSE 1
void _arb_fmpz_poly_evaluate_acb_horner(acb_t res, const fmpz * f, slong len, const acb_t x, slong prec);
void arb_fmpz_poly_evaluate_acb_horner(acb_t res, const fmpz_poly_t f, const acb_t a, slong prec);
void _arb_fmpz_poly_evaluate_acb_rectangular(acb_t res, const fmpz * f, slong len, const acb_t x, slong prec);
@ -40,6 +42,11 @@ void arb_fmpz_poly_evaluate_arb_rectangular(arb_t res, const fmpz_poly_t f, cons
void _arb_fmpz_poly_evaluate_arb(arb_t res, const fmpz * f, slong len, const arb_t x, slong prec);
void arb_fmpz_poly_evaluate_arb(arb_t res, const fmpz_poly_t f, const arb_t a, slong prec);
void arb_fmpz_poly_deflate(fmpz_poly_t result, const fmpz_poly_t input, ulong deflation);
ulong arb_fmpz_poly_deflation(const fmpz_poly_t input);
void arb_fmpz_poly_complex_roots(acb_ptr roots, const fmpz_poly_t poly, int flags, slong target_prec);
void arb_fmpz_poly_gauss_period_minpoly(fmpz_poly_t res, ulong q, ulong n);
#endif

190
arb_fmpz_poly/complex_roots.c Normal file
View file

@ -0,0 +1,190 @@
/*
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_fmpz_poly.h"
#include "flint/profiler.h"
int check_accuracy(acb_ptr vec, slong len, slong prec)
{
slong i;
for (i = 0; i < len; i++)
{
if (acb_rel_accuracy_bits(vec + i) < prec)
return 0;
}
return 1;
}
void
arb_fmpz_poly_complex_roots(acb_ptr roots, const fmpz_poly_t poly, int flags, slong target_prec)
{
slong i, j, prec, deg, deg_deflated, isolated, maxiter, deflation;
slong initial_prec, num_real;
acb_poly_t cpoly, cpoly_deflated;
fmpz_poly_t poly_deflated;
acb_ptr roots_deflated;
int removed_zero;
if (fmpz_poly_degree(poly) < 1)
return;
initial_prec = 32;
fmpz_poly_init(poly_deflated);
acb_poly_init(cpoly);
acb_poly_init(cpoly_deflated);
/* try to write poly as poly_deflated(x^deflation), possibly multiplied by x */
removed_zero = fmpz_is_zero(poly->coeffs);
if (removed_zero)
fmpz_poly_shift_right(poly_deflated, poly, 1);
else
fmpz_poly_set(poly_deflated, poly);
deflation = arb_fmpz_poly_deflation(poly_deflated);
arb_fmpz_poly_deflate(poly_deflated, poly_deflated, deflation);
deg = fmpz_poly_degree(poly);
deg_deflated = fmpz_poly_degree(poly_deflated);
if (flags & ARB_FMPZ_POLY_ROOTS_VERBOSE)
{
flint_printf("searching for %wd roots, %wd deflated\n", deg, deg_deflated);
}
/* we only need deg_deflated entries, but the remainder will be useful
as scratch space */
roots_deflated = _acb_vec_init(deg);
for (prec = initial_prec; ; prec *= 2)
{
acb_poly_set_fmpz_poly(cpoly_deflated, poly_deflated, prec);
maxiter = FLINT_MIN(FLINT_MAX(deg_deflated, 32), prec);
if (flags & ARB_FMPZ_POLY_ROOTS_VERBOSE)
{
TIMEIT_ONCE_START
flint_printf("prec=%wd: ", prec);
isolated = acb_poly_find_roots(roots_deflated, cpoly_deflated,
prec == initial_prec ? NULL : roots_deflated, maxiter, prec);
flint_printf("%wd isolated roots | ", isolated);
TIMEIT_ONCE_STOP
}
else
{
isolated = acb_poly_find_roots(roots_deflated, cpoly_deflated,
prec == initial_prec ? NULL : roots_deflated, maxiter, prec);
}
if (isolated == deg_deflated)
{
if (!check_accuracy(roots_deflated, deg_deflated, target_prec))
continue;
if (deflation == 1)
{
_acb_vec_set(roots, roots_deflated, deg_deflated);
}
else /* compute all nth roots */
{
acb_t w, w2;
acb_init(w);
acb_init(w2);
acb_unit_root(w, deflation, prec);
acb_unit_root(w2, 2 * deflation, prec);
for (i = 0; i < deg_deflated; i++)
{
if (arf_sgn(arb_midref(acb_realref(roots_deflated + i))) > 0)
{
acb_root_ui(roots + i * deflation,
roots_deflated + i, deflation, prec);
}
else
{
acb_neg(roots + i * deflation, roots_deflated + i);
acb_root_ui(roots + i * deflation,
roots + i * deflation, deflation, prec);
acb_mul(roots + i * deflation,
roots + i * deflation, w2, prec);
}
for (j = 1; j < deflation; j++)
{
acb_mul(roots + i * deflation + j,
roots + i * deflation + j - 1, w, prec);
}
}
acb_clear(w);
acb_clear(w2);
}
/* by assumption that poly is squarefree, must be just one */
if (removed_zero)
acb_zero(roots + deg_deflated * deflation);
if (!check_accuracy(roots, deg, target_prec))
continue;
acb_poly_set_fmpz_poly(cpoly, poly, prec);
if (!acb_poly_validate_real_roots(roots, cpoly, prec))
continue;
for (i = 0; i < deg; i++)
{
if (arb_contains_zero(acb_imagref(roots + i)))
arb_zero(acb_imagref(roots + i));
}
if (flags & ARB_FMPZ_POLY_ROOTS_VERBOSE)
flint_printf("done!\n");
break;
}
}
_acb_vec_sort_pretty(roots, deg);
/* pair conjugates */
num_real = 0;
for (i = 0; i < deg; i++)
if (acb_is_real(roots + i))
num_real++;
if (deg != num_real)
{
for (i = num_real, j = 0; i < deg; i++)
{
if (arb_is_positive(acb_imagref(roots + i)))
{
acb_swap(roots_deflated + j, roots + i);
j++;
}
}
for (i = 0; i < (deg - num_real) / 2; i++)
{
acb_swap(roots + num_real + 2 * i, roots_deflated + i);
acb_conj(roots + num_real + 2 * i + 1, roots + num_real + 2 * i);
}
}
fmpz_poly_clear(poly_deflated);
acb_poly_clear(cpoly);
acb_poly_clear(cpoly_deflated);
_acb_vec_clear(roots_deflated, deg);
}

38
arb_fmpz_poly/deflate.c Normal file
View file

@ -0,0 +1,38 @@
/*
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_fmpz_poly.h"
void
arb_fmpz_poly_deflate(fmpz_poly_t result, const fmpz_poly_t input, ulong deflation)
{
slong res_length, i;
if (deflation == 0)
{
flint_printf("Exception (fmpz_poly_deflate). Division by zero.\n");
flint_abort();
}
if (input->length <= 1 || deflation == 1)
{
fmpz_poly_set(result, input);
return;
}
res_length = (input->length - 1) / deflation + 1;
fmpz_poly_fit_length(result, res_length);
for (i = 0; i < res_length; i++)
fmpz_set(result->coeffs + i, input->coeffs + i*deflation);
result->length = res_length;
}

42
arb_fmpz_poly/deflation.c Normal file
View file

@ -0,0 +1,42 @@
/*
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_fmpz_poly.h"
ulong arb_fmpz_poly_deflation(const fmpz_poly_t input)
{
slong i, coeff, deflation;
if (input->length <= 1)
return input->length;
coeff = 1;
while (fmpz_is_zero(input->coeffs + coeff))
coeff++;
deflation = n_gcd(input->length - 1, coeff);
while ((deflation > 1) && (coeff + deflation < input->length))
{
for (i = 0; i < deflation - 1; i++)
{
coeff++;
if (!fmpz_is_zero(input->coeffs + coeff))
deflation = n_gcd(coeff, deflation);
}
if (i == deflation - 1)
coeff++;
}
return deflation;
}

192
arb_fmpz_poly/test/t-complex_roots.c Normal file
View file

@ -0,0 +1,192 @@
/*
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 "flint/arith.h"
#include "arb_fmpz_poly.h"
void
check_roots(const fmpz_poly_t poly, acb_srcptr roots, slong prec)
{
arb_ptr real;
acb_ptr upper;
arb_poly_t rpoly;
arb_t lead;
slong i, num_real, num_upper, deg;
deg = fmpz_poly_degree(poly);
num_real = 0;
for (i = 0; i < deg; i++)
if (acb_is_real(roots + i))
num_real++;
num_upper = (deg - num_real) / 2;
real = _arb_vec_init(num_real);
upper = _acb_vec_init(num_upper);
arb_poly_init(rpoly);
arb_init(lead);
for (i = 0; i < num_real; i++)
arb_set(real + i, acb_realref(roots + i));
for (i = 0; i < num_upper; i++)
acb_set(upper + i, roots + num_real + 2 * i);
arb_poly_product_roots_complex(rpoly, real, num_real, upper, num_upper, prec);
arb_set_fmpz(lead, poly->coeffs + deg);
arb_poly_scalar_mul(rpoly, rpoly, lead, prec);
if (!arb_poly_contains_fmpz_poly(rpoly, poly))
{
flint_printf("FAIL!\n");
flint_printf("deg = %wd, num_real = %wd, num_upper = %wd\n\n", deg, num_real, num_upper);
for (i = 0; i < deg; i++)
{
acb_printn(roots + i, 30, 0);
flint_printf("\n");
}
flint_printf("\npoly:\n");
fmpz_poly_print(poly); flint_printf("\n\n");
flint_printf("rpoly:\n");
arb_poly_printd(rpoly, 30); flint_printf("\n\n");
flint_abort();
}
_arb_vec_clear(real, num_real);
_acb_vec_clear(upper, num_upper);
arb_poly_clear(rpoly);
arb_clear(lead);
}
int main()
{
slong iter;
flint_rand_t state;
flint_printf("complex_roots....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 500 * arb_test_multiplier(); iter++)
{
fmpz_poly_t f, g;
fmpq_poly_t h;
fmpz_poly_factor_t fac;
fmpz_t t;
acb_ptr roots;
slong i, j, n, deg, prec, num_factors;
int flags;
prec = 20 + n_randint(state, 1000);
flags = 0; /* ARB_FMPZ_POLY_ROOTS_VERBOSE; */
fmpz_poly_init(f);
fmpz_poly_init(g);
fmpq_poly_init(h);
fmpz_init(t);
fmpz_poly_one(f);
num_factors = 1 + n_randint(state, 3);
for (i = 0; i < num_factors; i++)
{
n = n_randint(state, 18);
switch (n_randint(state, 12))
{
case 0:
fmpz_poly_zero(g);
for (j = 0; j <= n; j++)
fmpz_poly_set_coeff_ui(g, j, j+1);
break;
case 1:
arith_chebyshev_t_polynomial(g, n);
break;
case 2:
arith_chebyshev_u_polynomial(g, n);
break;
case 3:
arith_legendre_polynomial(h, n);
fmpq_poly_get_numerator(g, h);
break;
case 4:
arith_cyclotomic_polynomial(g, n);
break;
case 5:
arith_swinnerton_dyer_polynomial(g, n % 4);
break;
case 6:
arith_bernoulli_polynomial(h, n);
fmpq_poly_get_numerator(g, h);
break;
case 7:
fmpz_poly_zero(g);
fmpz_poly_fit_length(g, n+2);
arith_stirling_number_1_vec(g->coeffs, n+1, n+2);
_fmpz_poly_set_length(g, n+2);
fmpz_poly_shift_right(g, g, 1);
break;
case 8:
fmpq_poly_zero(h);
fmpq_poly_set_coeff_si(h, 0, 0);
fmpq_poly_set_coeff_si(h, 1, 1);
fmpq_poly_exp_series(h, h, n + 1);
fmpq_poly_get_numerator(g, h);
break;
case 9:
fmpz_poly_zero(g);
fmpz_poly_set_coeff_ui(g, 0, 1);
fmpz_poly_set_coeff_ui(g, 1, 100);
fmpz_poly_pow(g, g, n_randint(state, 5));
fmpz_poly_set_coeff_ui(g, n, 1);
break;
default:
fmpz_poly_randtest(g, state, 1 + n, 1 + n_randint(state, 300));
break;
}
fmpz_poly_mul(f, f, g);
}
if (!fmpz_poly_is_zero(f))
{
fmpz_poly_factor_init(fac);
fmpz_poly_factor_squarefree(fac, f);
for (i = 0; i < fac->num; i++)
{
deg = fmpz_poly_degree(fac->p + i);
roots = _acb_vec_init(deg);
arb_fmpz_poly_complex_roots(roots, fac->p + i, flags, prec);
check_roots(fac->p + i, roots, prec);
_acb_vec_clear(roots, deg);
}
fmpz_poly_factor_clear(fac);
}
fmpz_poly_clear(f);
fmpz_poly_clear(g);
fmpq_poly_clear(h);
fmpz_clear(t);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

50
doc/source/arb_fmpz_poly.rst
View file

@ -49,6 +49,56 @@ Evaluation
at the given real or complex number, respectively using Horner's rule, rectangular
splitting, or a default algorithm choice.
Utility methods
-------------------------------------------------------------------------------
.. function:: ulong arb_fmpz_poly_deflation(const fmpz_poly_t poly)
Finds the maximal exponent by which *poly* can be deflated.
.. function:: void arb_fmpz_poly_deflate(fmpz_poly_t res, const fmpz_poly_t poly, ulong deflation)
Sets *res* to a copy of *poly* deflated by the exponent *deflation*.
Polynomial roots
-------------------------------------------------------------------------------
.. function:: void arb_fmpz_poly_complex_roots(acb_ptr roots, const fmpz_poly_t poly, int flags, slong prec)
Writes to *roots* all the real and complex roots of the polynomial *poly*,
computed to *prec* accurate bits.
The real roots are written first in ascending order (with
the imaginary parts set exactly to zero). The following
nonreal roots are written in arbitrary order, but with conjugate pairs
grouped together (the root in the upper plane leading
the root in the lower plane).
The input polynomial *must* be squarefree. For a general polynomial,
do a squarefree factorization (which also gives the correct multiplicities
of the roots)::
fmpz_poly_factor_init(fac);
fmpz_poly_factor_squarefree(fac, poly);
for (i = 0; i < fac->num; i++)
{
deg = fmpz_poly_degree(fac->p + i);
flint_printf("%wd roots of multiplicity %wd\n", deg, fac->exp[i]);
roots = _acb_vec_init(deg);
arb_fmpz_poly_complex_roots(roots, fac->p + i, 0, prec);
_acb_vec_clear(roots, deg);
}
fmpz_poly_factor_clear(fac);
All roots are refined to a relative accuracy of at least *prec* bits.
The output values will generally have higher actual precision,
depending on the precision used internally the algorithm.
This implementation should be adequate for general use, but it is not
currently competitive with state-of-the-art isolation
methods for finding real roots alone.
Special polynomials
-------------------------------------------------------------------------------

230
examples/poly_roots.c
View file

@ -3,219 +3,18 @@
#include <string.h>
#include <ctype.h>
#include "acb.h"
#include "acb_poly.h"
#include "arb_fmpz_poly.h"
#include "flint/arith.h"
#include "flint/profiler.h"
int check_accuracy(acb_ptr vec, slong len, slong prec)
{
slong i;
for (i = 0; i < len; i++)
{
if (mag_cmp_2exp_si(arb_radref(acb_realref(vec + i)), -prec) >= 0
|| mag_cmp_2exp_si(arb_radref(acb_imagref(vec + i)), -prec) >= 0)
return 0;
}
return 1;
}
slong fmpz_poly_deflation(const fmpz_poly_t input)
{
slong i, coeff, deflation;
if (input->length <= 1)
return input->length;
coeff = 1;
while (fmpz_is_zero(input->coeffs + coeff))
coeff++;
deflation = n_gcd(input->length - 1, coeff);
while ((deflation > 1) && (coeff + deflation < input->length))
{
for (i = 0; i < deflation - 1; i++)
{
coeff++;
if (!fmpz_is_zero(input->coeffs + coeff))
deflation = n_gcd(coeff, deflation);
}
if (i == deflation - 1)
coeff++;
}
return deflation;
}
void
fmpz_poly_deflate(fmpz_poly_t result, const fmpz_poly_t input, ulong deflation)
{
slong res_length, i;
if (deflation == 0)
{
flint_printf("Exception (fmpz_poly_deflate). Division by zero.\n");
flint_abort();
}
if (input->length <= 1 || deflation == 1)
{
fmpz_poly_set(result, input);
return;
}
res_length = (input->length - 1) / deflation + 1;
fmpz_poly_fit_length(result, res_length);
for (i = 0; i < res_length; i++)
fmpz_set(result->coeffs + i, input->coeffs + i*deflation);
result->length = res_length;
}
void
fmpz_poly_complex_roots_squarefree(const fmpz_poly_t poly,
slong initial_prec,
slong target_prec,
slong print_digits)
{
slong i, j, prec, deg, deg_deflated, isolated, maxiter, deflation;
acb_poly_t cpoly, cpoly_deflated;
fmpz_poly_t poly_deflated;
acb_ptr roots, roots_deflated;
int removed_zero;
if (fmpz_poly_degree(poly) < 1)
return;
fmpz_poly_init(poly_deflated);
acb_poly_init(cpoly);
acb_poly_init(cpoly_deflated);
/* try to write poly as poly_deflated(x^deflation), possibly multiplied by x */
removed_zero = fmpz_is_zero(poly->coeffs);
if (removed_zero)
fmpz_poly_shift_right(poly_deflated, poly, 1);
else
fmpz_poly_set(poly_deflated, poly);
deflation = fmpz_poly_deflation(poly_deflated);
fmpz_poly_deflate(poly_deflated, poly_deflated, deflation);
deg = fmpz_poly_degree(poly);
deg_deflated = fmpz_poly_degree(poly_deflated);
flint_printf("searching for %wd roots, %wd deflated\n", deg, deg_deflated);
roots = _acb_vec_init(deg);
roots_deflated = _acb_vec_init(deg_deflated);
for (prec = initial_prec; ; prec *= 2)
{
acb_poly_set_fmpz_poly(cpoly_deflated, poly_deflated, prec);
maxiter = FLINT_MIN(FLINT_MAX(deg_deflated, 32), prec);
TIMEIT_ONCE_START
flint_printf("prec=%wd: ", prec);
isolated = acb_poly_find_roots(roots_deflated, cpoly_deflated,
prec == initial_prec ? NULL : roots_deflated, maxiter, prec);
flint_printf("%wd isolated roots | ", isolated);
TIMEIT_ONCE_STOP
if (isolated == deg_deflated)
{
if (!check_accuracy(roots_deflated, deg_deflated, target_prec))
continue;
if (deflation == 1)
{
_acb_vec_set(roots, roots_deflated, deg_deflated);
}
else /* compute all nth roots */
{
acb_t w, w2;
acb_init(w);
acb_init(w2);
acb_unit_root(w, deflation, prec);
acb_unit_root(w2, 2 * deflation, prec);
for (i = 0; i < deg_deflated; i++)
{
if (arf_sgn(arb_midref(acb_realref(roots_deflated + i))) > 0)
{
acb_root_ui(roots + i * deflation,
roots_deflated + i, deflation, prec);
}
else
{
acb_neg(roots + i * deflation, roots_deflated + i);
acb_root_ui(roots + i * deflation,
roots + i * deflation, deflation, prec);
acb_mul(roots + i * deflation,
roots + i * deflation, w2, prec);
}
for (j = 1; j < deflation; j++)
{
acb_mul(roots + i * deflation + j,
roots + i * deflation + j - 1, w, prec);
}
}
acb_clear(w);
acb_clear(w2);
}
/* by assumption that poly is squarefree, must be just one */
if (removed_zero)
acb_zero(roots + deg_deflated * deflation);
if (!check_accuracy(roots, deg, target_prec))
continue;
acb_poly_set_fmpz_poly(cpoly, poly, prec);
if (!acb_poly_validate_real_roots(roots, cpoly, prec))
continue;
for (i = 0; i < deg; i++)
{
if (arb_contains_zero(acb_imagref(roots + i)))
arb_zero(acb_imagref(roots + i));
}
flint_printf("done!\n");
break;
}
}
if (print_digits != 0)
{
_acb_vec_sort_pretty(roots, deg);
for (i = 0; i < deg; i++)
{
acb_printn(roots + i, print_digits, 0);
flint_printf("\n");
}
}
fmpz_poly_clear(poly_deflated);
acb_poly_clear(cpoly);
acb_poly_clear(cpoly_deflated);
_acb_vec_clear(roots, deg);
_acb_vec_clear(roots_deflated, deg_deflated);
}
int main(int argc, char *argv[])
{
fmpz_poly_t f, g;
fmpz_poly_factor_t fac;
fmpz_t t;
slong compd, printd, i, j;
acb_ptr roots;
slong compd, printd, i, j, deg;
int flags;
if (argc < 2)
{
@ -252,6 +51,7 @@ int main(int argc, char *argv[])
compd = 0;
printd = 0;
flags = ARB_FMPZ_POLY_ROOTS_VERBOSE;
fmpz_poly_init(f);
fmpz_poly_init(g);
@ -387,9 +187,23 @@ int main(int argc, char *argv[])
TIMEIT_ONCE_START
for (i = 0; i < fac->num; i++)
{
flint_printf("roots with multiplicity %wd\n", fac->exp[i]);
fmpz_poly_complex_roots_squarefree(fac->p + i,
32, compd * 3.32193 + 2, printd);
deg = fmpz_poly_degree(fac->p + i);
flint_printf("%wd roots with multiplicity %wd\n", deg, fac->exp[i]);
roots = _acb_vec_init(deg);
arb_fmpz_poly_complex_roots(roots, fac->p + i, flags, compd * 3.32193 + 2);
if (printd)
{
for (i = 0; i < deg; i++)
{
acb_printn(roots + i, printd, 0);
flint_printf("\n");
}
}
_acb_vec_clear(roots, deg);
}
TIMEIT_ONCE_STOP