Merge pull request #67 from argriffing/add-trace

ENH: add trace

acb_mat.h

@@ -322,6 +322,8 @@ 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_trace(acb_t trace, const acb_mat_t mat, slong prec);
+
 #ifdef __cplusplus
 }
 #endif

acb_mat/test/t-trace.c

@@ -0,0 +1,142 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2015 Arb authors
+
+******************************************************************************/
+
+#include "acb_mat.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("trace....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* check that the acb trace contains the fmpq trace */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        fmpq_mat_t Q;
+        fmpq_t Qtrace;
+        acb_mat_t A;
+        acb_t Atrace;
+        slong n, qbits, prec;
+
+        n = n_randint(state, 8);
+        qbits = 1 + n_randint(state, 100);
+        prec = 2 + n_randint(state, 200);
+
+        fmpq_mat_init(Q, n, n);
+        fmpq_init(Qtrace);
+
+        acb_mat_init(A, n, n);
+        acb_init(Atrace);
+
+        fmpq_mat_randtest(Q, state, qbits);
+        fmpq_mat_trace(Qtrace, Q);
+
+        acb_mat_set_fmpq_mat(A, Q, prec);
+        acb_mat_trace(Atrace, A, prec);
+
+        if (!acb_contains_fmpq(Atrace, Qtrace))
+        {
+            flint_printf("FAIL (containment, iter = %wd)\n", iter);
+            flint_printf("n = %wd, prec = %wd\n", n, prec);
+            flint_printf("\n");
+
+            flint_printf("Q = \n"); fmpq_mat_print(Q); flint_printf("\n\n");
+            flint_printf("Qtrace = \n"); fmpq_print(Qtrace); flint_printf("\n\n");
+
+            flint_printf("A = \n"); acb_mat_printd(A, 15); flint_printf("\n\n");
+            flint_printf("Atrace = \n"); acb_printd(Atrace, 15); flint_printf("\n\n");
+            flint_printf("Atrace = \n"); acb_print(Atrace); flint_printf("\n\n");
+
+            abort();
+        }
+
+        fmpq_mat_clear(Q);
+        fmpq_clear(Qtrace);
+        acb_mat_clear(A);
+        acb_clear(Atrace);
+    }
+
+    /* check trace(A*B) = trace(B*A) */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        slong m, n, prec;
+        acb_mat_t a, b, ab, ba;
+        acb_t trab, trba;
+
+        prec = 2 + n_randint(state, 200);
+
+        m = n_randint(state, 10);
+        n = n_randint(state, 10);
+
+        acb_mat_init(a, m, n);
+        acb_mat_init(b, n, m);
+        acb_mat_init(ab, m, m);
+        acb_mat_init(ba, n, n);
+
+        acb_init(trab);
+        acb_init(trba);
+
+        acb_mat_randtest(a, state, 2 + n_randint(state, 100), 10);
+        acb_mat_randtest(b, state, 2 + n_randint(state, 100), 10);
+
+        acb_mat_mul(ab, a, b, prec);
+        acb_mat_mul(ba, b, a, prec);
+
+        acb_mat_trace(trab, ab, prec);
+        acb_mat_trace(trba, ba, prec);
+
+        if (!acb_overlaps(trab, trba))
+        {
+            flint_printf("FAIL (overlap, iter = %wd)\n", iter);
+            flint_printf("m = %wd, n = %wd, prec = %wd\n", m, n, prec);
+            flint_printf("\n");
+
+            flint_printf("a = \n"); acb_mat_printd(a, 15); flint_printf("\n\n");
+            flint_printf("b = \n"); acb_mat_printd(b, 15); flint_printf("\n\n");
+            flint_printf("ab = \n"); acb_mat_printd(ab, 15); flint_printf("\n\n");
+            flint_printf("ba = \n"); acb_mat_printd(ba, 15); flint_printf("\n\n");
+
+            flint_printf("trace(ab) = \n"); acb_printd(trab, 15); flint_printf("\n\n");
+            flint_printf("trace(ba) = \n"); acb_printd(trba, 15); flint_printf("\n\n");
+        }
+
+        acb_clear(trab);
+        acb_clear(trba);
+
+        acb_mat_clear(a);
+        acb_mat_clear(b);
+        acb_mat_clear(ab);
+        acb_mat_clear(ba);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

acb_mat/trace.c

@@ -0,0 +1,45 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2015 Arb authors
+
+******************************************************************************/
+
+#include "acb_mat.h"
+
+void
+acb_mat_trace(acb_t trace, const acb_mat_t mat, slong prec)
+{
+    slong i, n = acb_mat_nrows(mat);
+
+    if (n == 0)
+    {
+        acb_zero(trace);
+    }
+    else
+    {
+        acb_set(trace, acb_mat_entry(mat, 0, 0));
+        for (i = 1; i < n; i++)
+        {
+            acb_add(trace, trace, acb_mat_entry(mat, i, i), prec);
+        }
+    }
+}

arb_mat.h

@@ -290,6 +290,8 @@ 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_trace(arb_t trace, const arb_mat_t mat, slong prec);
+
 #ifdef __cplusplus
 }
 #endif

arb_mat/test/t-trace.c

@@ -0,0 +1,142 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2015 Arb authors
+
+******************************************************************************/
+
+#include "arb_mat.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("trace....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* check that the arb trace contains the fmpq trace */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        fmpq_mat_t Q;
+        fmpq_t Qtrace;
+        arb_mat_t A;
+        arb_t Atrace;
+        slong n, qbits, prec;
+
+        n = n_randint(state, 8);
+        qbits = 1 + n_randint(state, 100);
+        prec = 2 + n_randint(state, 200);
+
+        fmpq_mat_init(Q, n, n);
+        fmpq_init(Qtrace);
+
+        arb_mat_init(A, n, n);
+        arb_init(Atrace);
+
+        fmpq_mat_randtest(Q, state, qbits);
+        fmpq_mat_trace(Qtrace, Q);
+
+        arb_mat_set_fmpq_mat(A, Q, prec);
+        arb_mat_trace(Atrace, A, prec);
+
+        if (!arb_contains_fmpq(Atrace, Qtrace))
+        {
+            flint_printf("FAIL (containment, iter = %wd)\n", iter);
+            flint_printf("n = %wd, prec = %wd\n", n, prec);
+            flint_printf("\n");
+
+            flint_printf("Q = \n"); fmpq_mat_print(Q); flint_printf("\n\n");
+            flint_printf("Qtrace = \n"); fmpq_print(Qtrace); flint_printf("\n\n");
+
+            flint_printf("A = \n"); arb_mat_printd(A, 15); flint_printf("\n\n");
+            flint_printf("Atrace = \n"); arb_printd(Atrace, 15); flint_printf("\n\n");
+            flint_printf("Atrace = \n"); arb_print(Atrace); flint_printf("\n\n");
+
+            abort();
+        }
+
+        fmpq_mat_clear(Q);
+        fmpq_clear(Qtrace);
+        arb_mat_clear(A);
+        arb_clear(Atrace);
+    }
+
+    /* check trace(A*B) = trace(B*A) */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        slong m, n, prec;
+        arb_mat_t a, b, ab, ba;
+        arb_t trab, trba;
+
+        prec = 2 + n_randint(state, 200);
+
+        m = n_randint(state, 10);
+        n = n_randint(state, 10);
+
+        arb_mat_init(a, m, n);
+        arb_mat_init(b, n, m);
+        arb_mat_init(ab, m, m);
+        arb_mat_init(ba, n, n);
+
+        arb_init(trab);
+        arb_init(trba);
+
+        arb_mat_randtest(a, state, 2 + n_randint(state, 100), 10);
+        arb_mat_randtest(b, state, 2 + n_randint(state, 100), 10);
+
+        arb_mat_mul(ab, a, b, prec);
+        arb_mat_mul(ba, b, a, prec);
+
+        arb_mat_trace(trab, ab, prec);
+        arb_mat_trace(trba, ba, prec);
+
+        if (!arb_overlaps(trab, trba))
+        {
+            flint_printf("FAIL (overlap, iter = %wd)\n", iter);
+            flint_printf("m = %wd, n = %wd, prec = %wd\n", m, n, prec);
+            flint_printf("\n");
+
+            flint_printf("a = \n"); arb_mat_printd(a, 15); flint_printf("\n\n");
+            flint_printf("b = \n"); arb_mat_printd(b, 15); flint_printf("\n\n");
+            flint_printf("ab = \n"); arb_mat_printd(ab, 15); flint_printf("\n\n");
+            flint_printf("ba = \n"); arb_mat_printd(ba, 15); flint_printf("\n\n");
+
+            flint_printf("trace(ab) = \n"); arb_printd(trab, 15); flint_printf("\n\n");
+            flint_printf("trace(ba) = \n"); arb_printd(trba, 15); flint_printf("\n\n");
+        }
+
+        arb_clear(trab);
+        arb_clear(trba);
+
+        arb_mat_clear(a);
+        arb_mat_clear(b);
+        arb_mat_clear(ab);
+        arb_mat_clear(ba);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

arb_mat/trace.c

@@ -0,0 +1,45 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2015 Arb authors
+
+******************************************************************************/
+
+#include "arb_mat.h"
+
+void
+arb_mat_trace(arb_t trace, const arb_mat_t mat, slong prec)
+{
+    slong i, n = arb_mat_nrows(mat);
+
+    if (n == 0)
+    {
+        arb_zero(trace);
+    }
+    else
+    {
+        arb_set(trace, arb_mat_entry(mat, 0, 0));
+        for (i = 1; i < n; i++)
+        {
+            arb_add(trace, trace, arb_mat_entry(mat, i, i), prec);
+        }
+    }
+}

doc/source/acb_mat.rst

@@ -306,3 +306,7 @@ Special functions
     evaluated using rectangular splitting.
     Error bounds are computed as for :func:`arb_mat_exp`.
 
+.. function:: void acb_mat_trace(acb_t trace, const acb_mat_t mat, slong prec)
+
+    Sets *trace* to the trace of the matrix, i.e. the sum of entries on the
+    main diagonal of *mat*. The matrix is required to be square.

doc/source/arb_mat.rst

@@ -313,3 +313,7 @@ Special functions
 
     We bound the sum on the right using :func:`mag_exp_tail`.
 
+.. function:: void arb_mat_trace(arb_t trace, const arb_mat_t mat, slong prec)
+
+    Sets *trace* to the trace of the matrix, i.e. the sum of entries on the
+    main diagonal of *mat*. The matrix is required to be square.