add fmprb_mat_pow_ui, contains_fmpz_mat

doc/source/fmprb_mat.rst

@@ -78,6 +78,8 @@ Comparisons
     Returns nonzero iff the matrices have the same dimensions
     and identical entries.
 
+.. function:: int fmprb_mat_contains_fmpz_mat(const fmprb_mat_t mat1, const fmpz_mat_t mat2)
+
 .. function:: int fmprb_mat_contains_fmpq_mat(const fmprb_mat_t mat1, const fmpq_mat_t mat2)
 
     Returns nonzero iff the matrices have the same dimensions and each entry
@@ -119,6 +121,10 @@ Arithmetic
     Sets *res* to the matrix product of *mat1* and *mat2*. The operands must have
     compatible dimensions for matrix multiplication.
 
+.. function:: void fmprb_mat_pow_ui(fmprb_mat_t res, const fmprb_mat_t mat, ulong exp, long prec)
+
+    Sets *res* to *mat* raised to the power *exp*. Requires that *mat*
+    is a square matrix.
 
 Gaussian elimination and solving
 -------------------------------------------------------------------------------

fmprb_mat.h

@@ -79,6 +79,8 @@ int fmprb_mat_equal(const fmprb_mat_t mat1, const fmprb_mat_t mat2);
 
 int fmprb_mat_contains_fmpq_mat(const fmprb_mat_t mat1, const fmpq_mat_t mat2);
 
+int fmprb_mat_contains_fmpz_mat(const fmprb_mat_t mat1, const fmpz_mat_t mat2);
+
 /* Special matrices */
 
 void fmprb_mat_zero(fmprb_mat_t mat);
@@ -95,6 +97,8 @@ void fmprb_mat_sub(fmprb_mat_t res, const fmprb_mat_t mat1, const fmprb_mat_t ma
 
 void fmprb_mat_mul(fmprb_mat_t res, const fmprb_mat_t mat1, const fmprb_mat_t mat2, long prec);
 
+void fmprb_mat_pow_ui(fmprb_mat_t B, const fmprb_mat_t A, ulong exp, long prec);
+
 /* Solving */
 
 static __inline__ void

fmprb_mat/contains_fmpz_mat.c

@@ -0,0 +1,44 @@
+/*=============================================================================
+
+    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) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmprb_mat.h"
+
+int
+fmprb_mat_contains_fmpz_mat(const fmprb_mat_t mat1, const fmpz_mat_t mat2)
+{
+    long i, j;
+
+    if ((fmprb_mat_nrows(mat1) != fmprb_mat_nrows(mat2)) ||
+        (fmprb_mat_ncols(mat1) != fmprb_mat_ncols(mat2)))
+        return 0;
+
+    for (i = 0; i < fmprb_mat_nrows(mat1); i++)
+        for (j = 0; j < fmprb_mat_ncols(mat1); j++)
+            if (!fmprb_contains_fmpz(fmprb_mat_entry(mat1, i, j),
+                fmpz_mat_entry(mat2, i, j)))
+                return 0;
+
+    return 1;
+}

fmprb_mat/pow_ui.c

@@ -0,0 +1,77 @@
+/*=============================================================================
+
+    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) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmprb_mat.h"
+
+void
+fmprb_mat_pow_ui(fmprb_mat_t B, const fmprb_mat_t A, ulong exp, long prec)
+{
+    long d = fmprb_mat_nrows(A);
+
+    if (exp <= 2 || d <= 1)
+    {
+        if (exp == 0 || d == 0)
+        {
+            fmprb_mat_one(B);
+        }
+        else if (d == 1)
+        {
+            fmprb_pow_ui(fmprb_mat_entry(B, 0, 0),
+                 fmprb_mat_entry(A, 0, 0), exp, prec);
+        }
+        else if (exp == 1)
+        {
+            fmprb_mat_set(B, A);
+        }
+        else if (exp == 2)
+        {
+            fmprb_mat_mul(B, A, A, prec);   /* todo: sqr */
+        }
+    }
+    else
+    {
+        fmprb_mat_t T, U;
+        long i;
+
+        fmprb_mat_init(T, d, d);
+        fmprb_mat_set(T, A);
+        fmprb_mat_init(U, d, d);
+
+        for (i = ((long) FLINT_BIT_COUNT(exp)) - 2; i >= 0; i--)
+        {
+            fmprb_mat_mul(U, T, T, prec);   /* todo: sqr */
+
+            if (exp & (1L << i))
+                fmprb_mat_mul(T, U, A, prec);
+            else
+                fmprb_mat_swap(T, U);
+        }
+
+        fmprb_mat_swap(B, T);
+        fmprb_mat_clear(T);
+        fmprb_mat_clear(U);
+    }
+}
+