ENH: arb_mat_sqr

arb_mat.h

@@ -137,6 +137,10 @@ void arb_mat_mul_classical(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, sl
 
 void arb_mat_mul_threaded(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec);
 
+void arb_mat_sqr(arb_mat_t B, const arb_mat_t A, slong prec);
+
+void arb_mat_sqr_classical(arb_mat_t B, const arb_mat_t A, slong prec);
+
 void arb_mat_pow_ui(arb_mat_t B, const arb_mat_t A, ulong exp, slong prec);
 
 /* Scalar arithmetic */

arb_mat/exp.c

@@ -69,7 +69,7 @@ _arb_mat_exp_taylor(arb_mat_t S, const arb_mat_t A, slong N, slong prec)
     {
         arb_mat_t T;
         arb_mat_init(T, arb_mat_nrows(A), arb_mat_nrows(A));
-        arb_mat_mul(T, A, A, prec);
+        arb_mat_sqr(T, A, prec);
         arb_mat_scalar_mul_2exp_si(T, T, -1);
         arb_mat_add(S, A, T, prec);
         arb_mat_one(T);
@@ -203,7 +203,7 @@ arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec)
 
         for (i = 0; i < r; i++)
         {
-            arb_mat_mul(T, B, B, wp);
+            arb_mat_sqr(T, B, wp);
             arb_mat_swap(T, B);
         }
 

arb_mat/pow_ui.c

@@ -47,7 +47,7 @@ arb_mat_pow_ui(arb_mat_t B, const arb_mat_t A, ulong exp, slong prec)
         }
         else if (exp == 2)
         {
-            arb_mat_mul(B, A, A, prec);   /* todo: sqr */
+            arb_mat_sqr(B, A, prec);
         }
     }
     else
@@ -61,7 +61,7 @@ arb_mat_pow_ui(arb_mat_t B, const arb_mat_t A, ulong exp, slong prec)
 
         for (i = ((slong) FLINT_BIT_COUNT(exp)) - 2; i >= 0; i--)
         {
-            arb_mat_mul(U, T, T, prec);   /* todo: sqr */
+            arb_mat_sqr(U, T, prec);
 
             if (exp & (WORD(1) << i))
                 arb_mat_mul(T, U, A, prec);

arb_mat/sqr.c

@@ -0,0 +1,42 @@
+/*=============================================================================
+
+    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_sqr(arb_mat_t B, const arb_mat_t A, slong prec)
+{
+    double d = (double) arb_mat_nrows(A);
+
+    if (flint_get_num_threads() > 1 && d*d*d * (double) prec > 100000)
+    {
+        arb_mat_mul_threaded(B, A, A, prec);   /* todo: sqr threaded */
+    }
+    else
+    {
+        arb_mat_sqr_classical(B, A, prec);
+    }
+}
+

arb_mat/sqr_classical.c

@@ -0,0 +1,116 @@
+/*=============================================================================
+
+    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_sqr_classical(arb_mat_t B, const arb_mat_t A, slong prec)
+{
+    slong n, i, j, k;
+    arb_t p, s;
+
+    n = arb_mat_nrows(A);
+
+    if (arb_mat_ncols(A) != n || arb_mat_nrows(B) != n || arb_mat_ncols(B) != n)
+    {
+        flint_printf("arb_mat_sqr: incompatible dimensions\n");
+        abort();
+    }
+
+    if (n == 0)
+    {
+        arb_mat_zero(B);
+        return;
+    }
+
+    if (n == 1)
+    {
+        arb_mul(arb_mat_entry(B, 0, 0),
+                  arb_mat_entry(A, 0, 0),
+                  arb_mat_entry(A, 0, 0), prec);
+        return;
+    }
+
+    if (A == B)
+    {
+        arb_mat_t T;
+        arb_mat_init(T, n, n);
+        arb_mat_sqr(T, A, prec);
+        arb_mat_swap(T, B);
+        arb_mat_clear(T);
+        return;
+    }
+
+    arb_init(p);
+    arb_init(s);
+
+    /* contribution of diagonal of A to diagonal of B */
+    for (i = 0; i < n; i++)
+    {
+        arb_mul(arb_mat_entry(B, i, i),
+                  arb_mat_entry(A, i, i),
+                  arb_mat_entry(A, i, i), prec);
+    }
+
+    for (i = 0; i < n; i++)
+    {
+        for (j = 0; j < i; j++)
+        {
+            /* contribution of off-diagonal of A to diagonal of B */
+            arb_mul(p, arb_mat_entry(A, i, j), arb_mat_entry(A, j, i), prec);
+            arb_add(arb_mat_entry(B, i, i), arb_mat_entry(B, i, i), p, prec);
+            arb_add(arb_mat_entry(B, j, j), arb_mat_entry(B, j, j), p, prec);
+
+            /* contribution of diagonal of A to off-diagonal of B */
+            arb_add(s, arb_mat_entry(A, i, i), arb_mat_entry(A, j, j), prec);
+            arb_mul(arb_mat_entry(B, i, j), arb_mat_entry(A, i, j), s, prec);
+            arb_mul(arb_mat_entry(B, j, i), arb_mat_entry(A, j, i), s, prec);
+        }
+    }
+
+    /* contribution of off-diagonal of A to off-diagonal of B */
+    if (n > 2)
+    {
+        for (i = 0; i < n; i++)
+        {
+            for (j = 0; j < n; j++)
+            {
+                for (k = 0; k < n; k++)
+                {
+                    if (i != j && j != k && k != i)
+                    {
+                        arb_addmul(arb_mat_entry(B, i, j),
+                                     arb_mat_entry(A, i, k),
+                                     arb_mat_entry(A, k, j), prec);
+                    }
+                }
+            }
+        }
+    }
+
+    arb_clear(p);
+    arb_clear(s);
+}
+

arb_mat/test/t-sqr.c

@@ -0,0 +1,99 @@
+/*=============================================================================
+
+    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("sqr....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000; iter++)
+    {
+        slong n, qbits1, rbits1, rbits2;
+        fmpq_mat_t A, B;
+        arb_mat_t a, b, c;
+
+        qbits1 = 2 + n_randint(state, 200);
+        rbits1 = 2 + n_randint(state, 200);
+        rbits2 = 2 + n_randint(state, 200);
+
+        n = n_randint(state, 10);
+
+        fmpq_mat_init(A, n, n);
+        fmpq_mat_init(B, n, n);
+
+        arb_mat_init(a, n, n);
+        arb_mat_init(b, n, n);
+        arb_mat_init(c, n, n);
+
+        fmpq_mat_randtest(A, state, qbits1);
+        fmpq_mat_mul(B, A, A);
+
+        arb_mat_set_fmpq_mat(a, A, rbits1);
+        arb_mat_sqr(b, a, rbits2);
+
+        if (!arb_mat_contains_fmpq_mat(b, B))
+        {
+            flint_printf("FAIL\n\n");
+            flint_printf("n = %wd, bits2 = %wd\n", n, rbits2);
+
+            flint_printf("A = "); fmpq_mat_print(A); flint_printf("\n\n");
+            flint_printf("B = "); fmpq_mat_print(B); flint_printf("\n\n");
+
+            flint_printf("a = "); arb_mat_printd(a, 15); flint_printf("\n\n");
+            flint_printf("b = "); arb_mat_printd(b, 15); flint_printf("\n\n");
+
+            abort();
+        }
+
+        /* test aliasing */
+        arb_mat_set(c, a);
+        arb_mat_sqr(c, c, rbits2);
+        if (!arb_mat_equal(c, b))
+        {
+            flint_printf("FAIL (aliasing 1)\n\n");
+            abort();
+        }
+
+        fmpq_mat_clear(A);
+        fmpq_mat_clear(B);
+
+        arb_mat_clear(a);
+        arb_mat_clear(b);
+        arb_mat_clear(c);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

doc/source/arb_mat.rst

@@ -171,6 +171,13 @@ Arithmetic
     if the matrices are sufficiently large and more than one thread
     can be used.
 
+.. function:: void arb_mat_sqr_classical(arb_mat_t B, const arb_mat_t A, slong prec)
+
+.. function:: void arb_mat_sqr(arb_mat_t res, const arb_mat_t mat, slong prec)
+
+   Sets *res* to the matrix square of *mat*. The operands must both be square
+   with the same dimensions.
+
 .. function:: void arb_mat_pow_ui(arb_mat_t res, const arb_mat_t mat, ulong exp, slong prec)
 
     Sets *res* to *mat* raised to the power *exp*. Requires that *mat*