Merge pull request #2189 from lgoettgens/lg/fmpz_mod_mat_pow

Add `gr_mat_pow_ui` and `gr_mat_pow_fmpz`, use it in added `fmpz_mod_mat_pow_ui`

doc/source/fmpz_mod_mat.rst

@@ -234,6 +234,11 @@ Matrix multiplication
 
     Set ``B`` to ``A^2``. The matrix ``A`` must be square.
 
+.. function:: void fmpz_mod_mat_pow_ui(fmpz_mod_mat_t B, const fmpz_mod_mat_t A, ulong e, const fmpz_mod_ctx_t ctx)
+
+    Sets ``B`` to the matrix ``A`` raised to the power ``e``,
+    where ``A`` must be a square matrix. Aliasing is allowed.
+
 .. function:: void fmpz_mod_mat_mul_fmpz_vec(fmpz * c, const fmpz_mod_mat_t A, const fmpz * b, slong blen, const fmpz_mod_ctx_t ctx)
               void fmpz_mod_mat_mul_fmpz_vec_ptr(fmpz * const * c, const fmpz_mod_mat_t A, const fmpz * const * b, slong blen, const fmpz_mod_ctx_t ctx)
 

doc/source/gr_mat.rst

@@ -290,6 +290,8 @@ Arithmetic
     exact division by two.
 
 .. function:: int gr_mat_sqr(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
+.. function:: int gr_mat_pow_ui(gr_mat_t res, const gr_mat_t mat, ulong e, gr_ctx_t ctx)
+.. function:: int gr_mat_pow_fmpz(gr_mat_t res, const gr_mat_t mat, fmpz_t e, gr_ctx_t ctx)
 
 .. function:: int gr_mat_add_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t ctx)
               int gr_mat_scalar_add(gr_mat_t res, gr_srcptr x, const gr_mat_t mat, gr_ctx_t ctx)

src/fmpz_mat/pow.c

@@ -46,7 +46,7 @@ fmpz_mat_pow(fmpz_mat_t B, const fmpz_mat_t A, ulong exp)
         fmpz_mat_init_set(T, A);
         fmpz_mat_init(U, d, d);
 
-        for (i = ((slong) FLINT_BIT_COUNT(exp)) - 2; i >= 0; i--)
+        for (i = FLINT_BIT_COUNT(exp) - 2; i >= 0; i--)
         {
             fmpz_mat_sqr(U, T);
 

src/fmpz_mod_mat.h

@@ -206,6 +206,8 @@ void fmpz_mod_mat_mul_classical_threaded(fmpz_mod_mat_t C,
 
 void fmpz_mod_mat_sqr(fmpz_mod_mat_t B, const fmpz_mod_mat_t A, const fmpz_mod_ctx_t ctx);
 
+void fmpz_mod_mat_pow_ui(fmpz_mod_mat_t B, const fmpz_mod_mat_t A, ulong exp, const fmpz_mod_ctx_t ctx);
+
 void fmpz_mod_mat_submul(fmpz_mod_mat_t D, const fmpz_mod_mat_t C,
                                const fmpz_mod_mat_t A, const fmpz_mod_mat_t B, const fmpz_mod_ctx_t ctx);
 

src/fmpz_mod_mat/pow.c

@@ -0,0 +1,28 @@
+/*
+    Copyright (C) 2025 Lars Göttgens
+
+    This file is part of FLINT.
+
+    FLINT 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 3 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include "longlong.h"
+#include "fmpz_mod.h"
+#include "fmpz_mod_mat.h"
+#include "gr.h"
+#include "gr_mat.h"
+
+void
+fmpz_mod_mat_pow_ui(fmpz_mod_mat_t B, const fmpz_mod_mat_t A, ulong exp,
+                    const fmpz_mod_ctx_t ctx)
+{
+    gr_ctx_t gr_ctx;
+
+    _gr_ctx_init_fmpz_mod_from_ref(gr_ctx, ctx);
+    GR_MUST_SUCCEED(gr_mat_pow_ui
+                    ((gr_mat_struct *) B, (const gr_mat_struct *) A, exp,
+                     gr_ctx));
+}

src/fmpz_mod_mat/test/main.c

@@ -26,6 +26,7 @@
 #include "t-mul_classical_threaded.c"
 #include "t-mul_fmpz_vec.c"
 #include "t-nullspace.c"
+#include "t-pow.c"
 #include "t-rank.c"
 #include "t-rref.c"
 #include "t-scalar_mul_fmpz.c"
@@ -57,6 +58,7 @@ test_struct tests[] =
     TEST_FUNCTION(fmpz_mod_mat_mul_classical_threaded),
     TEST_FUNCTION(fmpz_mod_mat_mul_fmpz_vec),
     TEST_FUNCTION(fmpz_mod_mat_nullspace),
+    TEST_FUNCTION(fmpz_mod_mat_pow),
     TEST_FUNCTION(fmpz_mod_mat_rank),
     TEST_FUNCTION(fmpz_mod_mat_rref),
     TEST_FUNCTION(fmpz_mod_mat_scalar_mul_fmpz),

src/fmpz_mod_mat/test/t-pow.c

@@ -0,0 +1,70 @@
+/*
+    Copyright (C) 2025 Lars Göttgens
+
+    This file is part of FLINT.
+
+    FLINT 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 3 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include "test_helpers.h"
+#include "fmpz_mod_mat.h"
+
+TEST_FUNCTION_START(fmpz_mod_mat_pow, state)
+{
+    slong i;
+
+    for (i = 0; i < 100 * flint_test_multiplier(); i++)
+    {
+        fmpz_mod_mat_t A, B, C;
+        slong n, i;
+        ulong e;
+        fmpz_mod_ctx_t ctx;
+
+        n = n_randint(state, 10);
+        e = n_randint(state, 20);
+
+        fmpz_mod_ctx_init_rand_bits(ctx, state, 200);
+
+        fmpz_mod_mat_init(A, n, n, ctx);
+        fmpz_mod_mat_init(B, n, n, ctx);
+        fmpz_mod_mat_init(C, n, n, ctx);
+
+        fmpz_mod_mat_randtest(A, state, ctx);
+        fmpz_mod_mat_randtest(B, state, ctx);
+
+        /* Make sure noise in the output is ok */
+        fmpz_mod_mat_randtest(B, state, ctx);
+
+        fmpz_mod_mat_pow_ui(B, A, e, ctx);
+
+        fmpz_mod_mat_one(C, ctx);
+        for (i = 0; i < e; i++)
+            fmpz_mod_mat_mul(C, C, A, ctx);
+
+        if (!fmpz_mod_mat_equal(C, B, ctx))
+        {
+            flint_printf("FAIL: results not equal\n");
+            fflush(stdout);
+            flint_abort();
+        }
+
+        fmpz_mod_mat_pow_ui(A, A, e, ctx);
+
+        if (!fmpz_mod_mat_equal(A, B, ctx))
+        {
+            flint_printf("FAIL: aliasing failed\n");
+            fflush(stdout);
+            flint_abort();
+        }
+
+        fmpz_mod_mat_clear(A, ctx);
+        fmpz_mod_mat_clear(B, ctx);
+        fmpz_mod_mat_clear(C, ctx);
+        fmpz_mod_ctx_clear(ctx);
+    }
+
+    TEST_FUNCTION_END(state);
+}

src/gr_mat.h

@@ -199,6 +199,9 @@ gr_mat_sqr(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
     return gr_mat_mul(res, mat, mat, ctx);
 }
 
+WARN_UNUSED_RESULT int gr_mat_pow_ui(gr_mat_t res, const gr_mat_t mat, ulong exp, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_pow_fmpz(gr_mat_t res, const gr_mat_t mat, fmpz_t exp, gr_ctx_t ctx);
+
 WARN_UNUSED_RESULT int _gr_mat_gr_poly_evaluate(gr_mat_t y, gr_srcptr poly, slong len, const gr_mat_t x, gr_ctx_t ctx);
 WARN_UNUSED_RESULT int gr_mat_gr_poly_evaluate(gr_mat_t res, const gr_poly_t f, const gr_mat_t a, gr_ctx_t ctx);
 

src/gr_mat/pow.c

@@ -0,0 +1,107 @@
+/*
+    Copyright (C) 2025 Lars Göttgens
+
+    This file is part of FLINT.
+
+    FLINT 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 3 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include "longlong.h"
+#include "fmpz.h"
+#include "gr.h"
+#include "gr_mat.h"
+
+int
+gr_mat_pow_ui(gr_mat_t res, const gr_mat_t mat, ulong exp, gr_ctx_t ctx)
+{
+    int status;
+    slong sz = ctx->sizeof_elem;
+    slong d;
+
+    d = gr_mat_nrows(res, ctx);
+
+    if (d != gr_mat_ncols(res, ctx) || d != gr_mat_nrows(mat, ctx)
+        || d != gr_mat_ncols(mat, ctx))
+    {
+        return GR_DOMAIN;
+    }
+
+    status = GR_SUCCESS;
+
+    if (exp <= 2 || d <= 1)
+    {
+        if (exp == 0 || d == 0)
+        {
+            status |= gr_mat_one(res, ctx);
+        }
+        else if (d == 1)
+        {
+            status |= gr_pow_ui(GR_MAT_ENTRY(res, 0, 0, sz),
+                                GR_MAT_ENTRY(mat, 0, 0, sz), exp, ctx);
+        }
+        else if (exp == 1)
+        {
+            status |= gr_mat_set(res, mat, ctx);
+        }
+        else if (exp == 2)
+        {
+            status |= gr_mat_sqr(res, mat, ctx);
+        }
+    }
+    else
+    {
+        gr_ctx_t mctx;
+
+        gr_ctx_init_matrix_ring(mctx, ctx, d);
+        status |= gr_pow_ui(res, mat, exp, mctx);
+        gr_ctx_clear(mctx);
+    }
+
+    return status;
+}
+
+int
+gr_mat_pow_fmpz(gr_mat_t res, const gr_mat_t mat, fmpz_t exp, gr_ctx_t ctx)
+{
+    int status;
+    slong sz = ctx->sizeof_elem;
+    slong d;
+
+    d = gr_mat_nrows(res, ctx);
+
+    if (d != gr_mat_ncols(res, ctx) || d != gr_mat_nrows(mat, ctx)
+        || d != gr_mat_ncols(mat, ctx))
+    {
+        return GR_DOMAIN;
+    }
+
+    status = GR_SUCCESS;
+
+
+    if (fmpz_is_zero(exp) || d == 0)
+    {
+        status |= gr_mat_one(res, ctx);
+    }
+    else if (d == 1)
+    {
+        status |= gr_pow_fmpz(GR_MAT_ENTRY(res, 0, 0, sz),
+                              GR_MAT_ENTRY(mat, 0, 0, sz), exp, ctx);
+    }
+    else if (fmpz_is_one(exp))
+    {
+        status |= gr_mat_set(res, mat, ctx);
+    }
+    else
+    {
+        gr_ctx_t mctx;
+
+        gr_ctx_init_matrix_ring(mctx, ctx, d);
+        status |= gr_pow_fmpz(res, mat, exp, mctx);
+        gr_ctx_clear(mctx);
+    }
+
+    return status;
+}

src/gr_mat/test/main.c

@@ -38,6 +38,7 @@
 #include "t-mul_strassen.c"
 #include "t-mul_waksman.c"
 #include "t-nullspace.c"
+#include "t-pow.c"
 #include "t-properties.c"
 #include "t-randrank.c"
 #include "t-rank.c"
@@ -88,6 +89,7 @@ test_struct tests[] =
     TEST_FUNCTION(gr_mat_mul_strassen),
     TEST_FUNCTION(gr_mat_mul_waksman),
     TEST_FUNCTION(gr_mat_nullspace),
+    TEST_FUNCTION(gr_mat_pow),
     TEST_FUNCTION(gr_mat_properties),
     TEST_FUNCTION(gr_mat_randrank),
     TEST_FUNCTION(gr_mat_rank),