diff --git a/acb_mat.h b/acb_mat.h index ef55e9a2..65cb4066 100644 --- a/acb_mat.h +++ b/acb_mat.h @@ -432,9 +432,10 @@ void acb_mat_exp_taylor_sum(acb_mat_t S, const acb_mat_t A, slong N, slong prec) void acb_mat_exp(acb_mat_t B, const acb_mat_t A, slong prec); -void _acb_mat_charpoly(acb_ptr cp, const acb_mat_t mat, slong prec); - -void acb_mat_charpoly(acb_poly_t cp, const acb_mat_t mat, slong prec); +void _acb_mat_charpoly(acb_ptr poly, const acb_mat_t mat, slong prec); +void acb_mat_charpoly(acb_poly_t poly, const acb_mat_t mat, slong prec); +void _acb_mat_companion(acb_mat_t mat, acb_srcptr poly, slong prec); +void acb_mat_companion(acb_mat_t mat, const acb_poly_t poly, slong prec); void acb_mat_trace(acb_t trace, const acb_mat_t mat, slong prec); diff --git a/acb_mat/companion.c b/acb_mat/companion.c new file mode 100644 index 00000000..143e95f6 --- /dev/null +++ b/acb_mat/companion.c @@ -0,0 +1,49 @@ +/* + Copyright (C) 2018 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 . +*/ + +#include "acb_mat.h" + +void +_acb_mat_companion(acb_mat_t A, acb_srcptr poly, slong prec) +{ + slong i, j, n; + acb_t c; + + n = acb_mat_nrows(A); + + if (n == 0) + return; + + for (i = 0; i < n - 1; i++) + for (j = 0; j < n; j++) + acb_set_ui(acb_mat_entry(A, i, j), (i + 1) == j); + + acb_init(c); + acb_inv(c, poly + n, prec); + acb_neg(c, c); + for (j = 0; j < n; j++) + acb_mul(acb_mat_entry(A, n - 1, j), poly + j, c, prec); + acb_clear(c); +} + +void +acb_mat_companion(acb_mat_t A, const acb_poly_t poly, slong prec) +{ + slong n = acb_mat_nrows(A); + + if (n != acb_poly_degree(poly) || n != acb_mat_ncols(A)) + { + flint_printf("acb_mat_companion: incompatible dimensions!\n"); + flint_abort(); + } + + _acb_mat_companion(A, poly->coeffs, prec); +} diff --git a/acb_mat/test/t-companion.c b/acb_mat/test/t-companion.c new file mode 100644 index 00000000..a62098f7 --- /dev/null +++ b/acb_mat/test/t-companion.c @@ -0,0 +1,65 @@ +/* + Copyright (C) 2018 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 . +*/ + +#include "acb_mat.h" + +int +main(void) +{ + slong iter; + flint_rand_t state; + + flint_printf("companion...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100 * arb_test_multiplier(); iter++) + { + acb_mat_t A; + acb_poly_t f, g; + slong n, prec; + + acb_poly_init(f); + acb_poly_init(g); + + do { + acb_poly_randtest(f, state, 1 + n_randint(state, 8), 1 + n_randint(state, 1000), 10); + } while (acb_poly_degree(f) < 0); + + n = acb_poly_degree(f); + prec = 2 + n_randint(state, 200); + + acb_mat_init(A, n, n); + acb_mat_randtest(A, state, 1 + n_randint(state, 1000), 10); + acb_mat_companion(A, f, prec); + acb_mat_charpoly(g, A, prec); + acb_poly_scalar_mul(g, g, acb_poly_get_coeff_ptr(f, n), prec); + + if (!acb_poly_contains(g, f)) + { + flint_printf("FAIL\n"); + flint_printf("A:\n"), acb_mat_printd(A, 15), flint_printf("\n"); + flint_printf("f:\n"), acb_poly_printd(f, 15), flint_printf("\n"); + flint_printf("g:\n"), acb_poly_printd(g, 15), flint_printf("\n"); + flint_abort(); + } + + acb_mat_clear(A); + acb_poly_clear(f); + acb_poly_clear(g); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return 0; +} diff --git a/arb_mat.h b/arb_mat.h index a31437b6..5c06ab4c 100644 --- a/arb_mat.h +++ b/arb_mat.h @@ -397,9 +397,10 @@ void arb_mat_exp_taylor_sum(arb_mat_t S, const arb_mat_t A, slong N, slong prec) void arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec); -void _arb_mat_charpoly(arb_ptr cp, const arb_mat_t mat, slong prec); - -void arb_mat_charpoly(arb_poly_t cp, const arb_mat_t mat, slong prec); +void _arb_mat_charpoly(arb_ptr poly, const arb_mat_t mat, slong prec); +void arb_mat_charpoly(arb_poly_t poly, const arb_mat_t mat, slong prec); +void _arb_mat_companion(arb_mat_t mat, arb_srcptr poly, slong prec); +void arb_mat_companion(arb_mat_t mat, const arb_poly_t poly, slong prec); void arb_mat_trace(arb_t trace, const arb_mat_t mat, slong prec); diff --git a/arb_mat/companion.c b/arb_mat/companion.c new file mode 100644 index 00000000..e3bd5068 --- /dev/null +++ b/arb_mat/companion.c @@ -0,0 +1,49 @@ +/* + Copyright (C) 2018 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 . +*/ + +#include "arb_mat.h" + +void +_arb_mat_companion(arb_mat_t A, arb_srcptr poly, slong prec) +{ + slong i, j, n; + arb_t c; + + n = arb_mat_nrows(A); + + if (n == 0) + return; + + for (i = 0; i < n - 1; i++) + for (j = 0; j < n; j++) + arb_set_ui(arb_mat_entry(A, i, j), (i + 1) == j); + + arb_init(c); + arb_inv(c, poly + n, prec); + arb_neg(c, c); + for (j = 0; j < n; j++) + arb_mul(arb_mat_entry(A, n - 1, j), poly + j, c, prec); + arb_clear(c); +} + +void +arb_mat_companion(arb_mat_t A, const arb_poly_t poly, slong prec) +{ + slong n = arb_mat_nrows(A); + + if (n != arb_poly_degree(poly) || n != arb_mat_ncols(A)) + { + flint_printf("arb_mat_companion: incompatible dimensions!\n"); + flint_abort(); + } + + _arb_mat_companion(A, poly->coeffs, prec); +} diff --git a/arb_mat/indeterminate.c b/arb_mat/indeterminate.c index 24d02982..b7207f41 100644 --- a/arb_mat/indeterminate.c +++ b/arb_mat/indeterminate.c @@ -9,7 +9,7 @@ (at your option) any later version. See . */ -#include "acb_mat.h" +#include "arb_mat.h" void arb_mat_indeterminate(arb_mat_t A) diff --git a/arb_mat/test/t-companion.c b/arb_mat/test/t-companion.c new file mode 100644 index 00000000..3b9e5115 --- /dev/null +++ b/arb_mat/test/t-companion.c @@ -0,0 +1,65 @@ +/* + Copyright (C) 2018 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 . +*/ + +#include "arb_mat.h" + +int +main(void) +{ + slong iter; + flint_rand_t state; + + flint_printf("companion...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100 * arb_test_multiplier(); iter++) + { + arb_mat_t A; + arb_poly_t f, g; + slong n, prec; + + arb_poly_init(f); + arb_poly_init(g); + + do { + arb_poly_randtest(f, state, 1 + n_randint(state, 8), 1 + n_randint(state, 1000), 10); + } while (arb_poly_degree(f) < 0); + + n = arb_poly_degree(f); + prec = 2 + n_randint(state, 200); + + arb_mat_init(A, n, n); + arb_mat_randtest(A, state, 1 + n_randint(state, 1000), 10); + arb_mat_companion(A, f, prec); + arb_mat_charpoly(g, A, prec); + arb_poly_scalar_mul(g, g, arb_poly_get_coeff_ptr(f, n), prec); + + if (!arb_poly_contains(g, f)) + { + flint_printf("FAIL\n"); + flint_printf("A:\n"), arb_mat_printd(A, 15), flint_printf("\n"); + flint_printf("f:\n"), arb_poly_printd(f, 15), flint_printf("\n"); + flint_printf("g:\n"), arb_poly_printd(g, 15), flint_printf("\n"); + flint_abort(); + } + + arb_mat_clear(A); + arb_poly_clear(f); + arb_poly_clear(g); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return 0; +} diff --git a/doc/source/acb_mat.rst b/doc/source/acb_mat.rst index b4ea2236..06383407 100644 --- a/doc/source/acb_mat.rst +++ b/doc/source/acb_mat.rst @@ -484,18 +484,26 @@ Gaussian elimination and solving output matrices are set to the approximate floating-point results with zeroed error bounds. -Characteristic polynomial +Characteristic polynomial and companion matrix ------------------------------------------------------------------------------- -.. function:: void _acb_mat_charpoly(acb_ptr cp, const acb_mat_t mat, slong prec) +.. function:: void _acb_mat_charpoly(acb_ptr poly, const acb_mat_t mat, slong prec) -.. function:: void acb_mat_charpoly(acb_poly_t cp, const acb_mat_t mat, slong prec) +.. function:: void acb_mat_charpoly(acb_poly_t poly, const acb_mat_t mat, slong prec) - Sets *cp* to the characteristic polynomial of *mat* which must be + Sets *poly* to the characteristic polynomial of *mat* which must be a square matrix. If the matrix has *n* rows, the underscore method requires space for `n + 1` output coefficients. Employs a division-free algorithm using `O(n^4)` operations. +.. function:: void _acb_mat_companion(acb_mat_t mat, acb_srcptr poly, slong prec) + +.. function:: void acb_mat_companion(acb_mat_t mat, const acb_poly_t poly, slong prec) + + Sets the *n* by *n* matrix *mat* to the companion matrix of the polynomial + *poly* which must have degree *n*. + The underscore method reads `n + 1` input coefficients. + Special functions ------------------------------------------------------------------------------- diff --git a/doc/source/arb_mat.rst b/doc/source/arb_mat.rst index fbd4626e..a6042016 100644 --- a/doc/source/arb_mat.rst +++ b/doc/source/arb_mat.rst @@ -614,18 +614,26 @@ Cholesky decomposition and solving The inverse is calculated using the method of [Kri2013]_ which is more efficient than solving `AX = I` with :func:`arb_mat_solve_ldl_precomp`. -Characteristic polynomial +Characteristic polynomial and companion matrix ------------------------------------------------------------------------------- -.. function:: void _arb_mat_charpoly(arb_ptr cp, const arb_mat_t mat, slong prec) +.. function:: void _arb_mat_charpoly(arb_ptr poly, const arb_mat_t mat, slong prec) -.. function:: void arb_mat_charpoly(arb_poly_t cp, const arb_mat_t mat, slong prec) +.. function:: void arb_mat_charpoly(arb_poly_t poly, const arb_mat_t mat, slong prec) - Sets *cp* to the characteristic polynomial of *mat* which must be + Sets *poly* to the characteristic polynomial of *mat* which must be a square matrix. If the matrix has *n* rows, the underscore method requires space for `n + 1` output coefficients. Employs a division-free algorithm using `O(n^4)` operations. +.. function:: void _arb_mat_companion(arb_mat_t mat, arb_srcptr poly, slong prec) + +.. function:: void arb_mat_companion(arb_mat_t mat, const arb_poly_t poly, slong prec) + + Sets the *n* by *n* matrix *mat* to the companion matrix of the polynomial + *poly* which must have degree *n*. + The underscore method reads `n + 1` input coefficients. + Special functions -------------------------------------------------------------------------------