companion matrices

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);
 

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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);
+}

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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;
+}

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);
 

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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);
+}

arb_mat/indeterminate.c

@@ -9,7 +9,7 @@
     (at your option) any later version.  See <http://www.gnu.org/licenses/>.
 */
 
-#include "acb_mat.h"
+#include "arb_mat.h"
 
 void
 arb_mat_indeterminate(arb_mat_t A)

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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;
+}

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
 -------------------------------------------------------------------------------
 

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
 -------------------------------------------------------------------------------