add fmpz_mod_mat_det

doc/source/fmpz_mod_mat.rst

@@ -262,6 +262,10 @@ Gaussian elimination
 --------------------------------------------------------------------------------
 
 
+.. function:: void fmpz_mod_mat_det(fmpz_t res, const fmpz_mod_mat_t mat, const fmpz_mod_ctx_t ctx)
+
+    Set ``res`` to the determinant of the matrix ``mat``.
+
 .. function:: slong fmpz_mod_mat_rref(fmpz_mod_mat_t res, const fmpz_mod_mat_t mat, const fmpz_mod_ctx_t ctx)
 
     Sets ``res`` to the reduced row echelon form of ``mat``

src/fmpz_mod_mat.h

@@ -223,6 +223,8 @@ void fmpz_mod_mat_trace(fmpz_t trace, const fmpz_mod_mat_t mat, const fmpz_mod_c
 
 /* Gaussian elimination *********************************************/
 
+void fmpz_mod_mat_det(fmpz_t res, const fmpz_mod_mat_t mat, const fmpz_mod_ctx_t ctx);
+
 slong fmpz_mod_mat_rref(fmpz_mod_mat_t res, const fmpz_mod_mat_t mat, const fmpz_mod_ctx_t ctx);
 
 slong fmpz_mod_mat_reduce_row(fmpz_mod_mat_t A, slong * P, slong * L, slong m, const fmpz_mod_ctx_t ctx);

src/fmpz_mod_mat/det.c

@@ -0,0 +1,42 @@
+/*
+    Copyright (C) 2024 Fredrik Johansson
+
+    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 2.1 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include "fmpz_mod.h"
+#include "fmpz_mod_mat.h"
+#include "gr.h"
+#include "gr_mat.h"
+
+void fmpz_mod_mat_det(fmpz_t res, const fmpz_mod_mat_t mat, const fmpz_mod_ctx_t ctx)
+{
+    gr_ctx_t gr_ctx;
+    slong n = mat->r;
+
+    if (!fmpz_mod_mat_is_square(mat, ctx))
+    {
+        flint_throw(FLINT_ERROR, "Exception (fmpz_mod_mat_charpoly_berkowitz). Non-square matrix.\n");
+    }
+
+    _gr_ctx_init_fmpz_mod_from_ref(gr_ctx, ctx);
+
+    if (n <= 4)
+    {
+        GR_MUST_SUCCEED(gr_mat_det_cofactor(res, (const gr_mat_struct *) mat, gr_ctx));
+    }
+    else
+    {
+        if (gr_mat_det_lu(res, (const gr_mat_struct *) mat, gr_ctx) != GR_SUCCESS)
+        {
+            /* Fall back on division-free algorithm if we encountered an impossible inverse */
+            /* Could try something else here: faddeev_bsgs (O(n^3.5)) or Howell form. */
+            GR_MUST_SUCCEED(gr_mat_det_berkowitz(res, (const gr_mat_struct *) mat, gr_ctx));
+        }
+    }
+}

src/fmpz_mod_mat/test/main.c

@@ -17,6 +17,7 @@
 #include "t-add_sub_neg.c"
 #include "t-can_solve.c"
 #include "t-charpoly.c"
+#include "t-det.c"
 #include "t-fmpz_vec_mul.c"
 #include "t-get_set_fmpz_mat.c"
 #include "t-howell_form.c"
@@ -47,6 +48,7 @@ test_struct tests[] =
     TEST_FUNCTION(fmpz_mod_mat_add_sub_neg),
     TEST_FUNCTION(fmpz_mod_mat_can_solve),
     TEST_FUNCTION(fmpz_mod_mat_charpoly),
+    TEST_FUNCTION(fmpz_mod_mat_det),
     TEST_FUNCTION(fmpz_mod_mat_fmpz_vec_mul),
     TEST_FUNCTION(fmpz_mod_mat_get_set_fmpz_mat),
     TEST_FUNCTION(fmpz_mod_mat_howell_form),

src/fmpz_mod_mat/test/t-det.c

@@ -0,0 +1,56 @@
+/*
+    Copyright (C) 2024 Fredrik Johansson
+
+    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 2.1 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include "test_helpers.h"
+#include "fmpz.h"
+#include "fmpz_mod.h"
+#include "fmpz_mod_mat.h"
+#include "fmpz_mod_poly.h"
+#include "fmpz_mat.h"
+
+TEST_FUNCTION_START(fmpz_mod_mat_det, state)
+{
+    fmpz_mod_ctx_t ctx;
+    fmpz_mod_mat_t A;
+    fmpz_mat_t B;
+    fmpz_t d1, d2;
+    slong i, n;
+
+    for (i = 0; i < 1000 * flint_test_multiplier(); i++)
+    {
+        n = n_randint(state, 7);
+        fmpz_mod_ctx_init_rand_bits(ctx, state, 200);
+
+        fmpz_init(d1);
+        fmpz_init(d2);
+
+        fmpz_mod_mat_init(A, n, n, ctx);
+        fmpz_mat_init(B, n, n);
+
+        fmpz_mod_mat_randtest(A, state, ctx);
+        fmpz_mod_mat_get_fmpz_mat(B, A, ctx);
+
+        fmpz_mod_mat_det(d1, A, ctx);
+        fmpz_mat_det(d2, B);
+        fmpz_mod_set_fmpz(d2, d2, ctx);
+
+        FLINT_TEST(fmpz_equal(d1, d2));
+
+        fmpz_mod_mat_clear(A, ctx);
+        fmpz_mat_clear(B);
+        fmpz_clear(d1);
+        fmpz_clear(d2);
+
+        fmpz_mod_ctx_clear(ctx);
+    }
+
+    TEST_FUNCTION_END(state);
+}

src/gr/fmpz_mod.c

@@ -613,7 +613,7 @@ _gr_fmpz_mod_roots_gr_poly(gr_vec_t roots, gr_vec_t mult, const fmpz_mod_poly_t
 }
 
 int
-_gr_fmpz_mod_mat_mul(fmpz_mat_t res, const fmpz_mat_t x, const fmpz_mat_t y, gr_ctx_t ctx)
+_gr_fmpz_mod_mat_mul(fmpz_mod_mat_t res, const fmpz_mod_mat_t x, const fmpz_mod_mat_t y, gr_ctx_t ctx)
 {
     fmpz_mat_mul(res, x, y);
     _fmpz_mod_mat_reduce(res, FMPZ_MOD_CTX(ctx));
@@ -623,11 +623,18 @@ _gr_fmpz_mod_mat_mul(fmpz_mat_t res, const fmpz_mat_t x, const fmpz_mat_t y, gr_
 /* todo: tune cutoff for different bit sizes */
 /* also tune cutoff for triangular solving */
 int
-_gr_fmpz_mod_mat_lu(slong * rank, slong * P, fmpz_mat_t LU, const fmpz_mat_t A, int rank_check, gr_ctx_t ctx)
+_gr_fmpz_mod_mat_lu(slong * rank, slong * P, fmpz_mod_mat_t LU, const fmpz_mod_mat_t A, int rank_check, gr_ctx_t ctx)
 {
     return gr_mat_lu_recursive(rank, P, (gr_mat_struct *) LU, (const gr_mat_struct *) A, rank_check, 8, ctx);
 }
 
+int
+_gr_fmpz_mod_mat_det(fmpz_t res, const fmpz_mod_mat_t mat, gr_ctx_t ctx)
+{
+    fmpz_mod_mat_det(res, mat, FMPZ_MOD_CTX(ctx));
+    return GR_SUCCESS;
+}
+
 int _fmpz_mod_methods_initialized = 0;
 
 gr_static_method_table _fmpz_mod_methods;
@@ -694,6 +701,7 @@ gr_method_tab_input _fmpz_mod_methods_input[] =
     {GR_METHOD_POLY_ROOTS,      (gr_funcptr) _gr_fmpz_mod_roots_gr_poly},
     {GR_METHOD_MAT_MUL,         (gr_funcptr) _gr_fmpz_mod_mat_mul},
     {GR_METHOD_MAT_LU,          (gr_funcptr) _gr_fmpz_mod_mat_lu},
+    {GR_METHOD_MAT_DET,         (gr_funcptr) _gr_fmpz_mod_mat_det},
     {0,                         (gr_funcptr) NULL},
 };
 

src/qadic/norm_resultant.c

@@ -10,88 +10,8 @@
 */
 
 #include "qadic.h"
-
-/*
-    Computes the characteristic polynomial of the $n \times n$ matrix $M$
-    modulo \code{pN} using a division-free algorithm in $O(n^4)$ ring
-    operations.
-
-    Only returns the determinant.
-
-    Assumes that $n$ is at least $2$.
- */
-
-static
-void _fmpz_mod_mat_det(fmpz_t rop, const fmpz *M, slong n, const fmpz_t pN)
-{
-    fmpz *F;
-    fmpz *a;
-    fmpz *A;
-    fmpz_t s;
-    slong t, i, j, p, k;
-
-    F = _fmpz_vec_init(n);
-    a = _fmpz_vec_init((n-1) * n);
-    A = _fmpz_vec_init(n);
-
-    fmpz_init(s);
-
-    fmpz_neg(F + 0, M + 0*n + 0);
-
-    for (t = 1; t < n; t++)
-    {
-        for (i = 0; i <= t; i++)
-            fmpz_set(a + 0*n + i, M + i*n + t);
-
-        fmpz_set(A + 0, M + t*n + t);
-
-        for (p = 1; p < t; p++)
-        {
-            for (i = 0; i <= t; i++)
-            {
-                fmpz_zero(s);
-                for (j = 0; j <= t; j++)
-                    fmpz_addmul(s, M + i*n + j, a + (p-1)*n + j);
-                fmpz_mod(a + p*n + i, s, pN);
-            }
-
-            fmpz_set(A + p, a + p*n + t);
-        }
-
-        fmpz_zero(s);
-        for (j = 0; j <= t; j++)
-            fmpz_addmul(s, M + t*n + j, a + (t-1)*n + j);
-        fmpz_mod(A + t, s, pN);
-
-        for (p = 0; p <= t; p++)
-        {
-            fmpz_sub(F + p, F + p, A + p);
-            for (k = 0; k < p; k++)
-                fmpz_submul(F + p, A + k, F + (p-k-1));
-            fmpz_mod(F + p, F + p, pN);
-        }
-    }
-
-    /*
-        Now [F{n-1}, F{n-2}, ..., F{0}, 1] is the
-        characteristic polynomial of the matrix M.
-     */
-
-    if (n % WORD(2) == 0)
-    {
-        fmpz_set(rop, F + (n-1));
-    }
-    else
-    {
-        fmpz_neg(rop, F + (n-1));
-        fmpz_mod(rop, rop, pN);
-    }
-
-    _fmpz_vec_clear(F, n);
-    _fmpz_vec_clear(a, (n-1)*n);
-    _fmpz_vec_clear(A, n);
-    fmpz_clear(s);
-}
+#include "fmpz_mod.h"
+#include "fmpz_mod_mat.h"
 
 void _qadic_norm_resultant(fmpz_t rop, const fmpz *op, slong len,
                            const fmpz *a, const slong *j, slong lena,
@@ -110,32 +30,25 @@ void _qadic_norm_resultant(fmpz_t rop, const fmpz *op, slong len,
     }
     else  /* len >= 2 */
     {
-        {
-            const slong n = d + len - 1;
-            slong i, k;
-            fmpz *M;
-
-            M = flint_calloc(n * n, sizeof(fmpz));
-
-            for (k = 0; k < len-1; k++)
-            {
-                for (i = 0; i < lena; i++)
-                {
-                    M[k * n + k + (d - j[i])] = a[i];
-                }
-            }
-            for (k = 0; k < d; k++)
-            {
-                for (i = 0; i < len; i++)
-                {
-                    M[(len-1 + k) * n + k + (len-1 - i)] = op[i];
-                }
-            }
-
-            _fmpz_mod_mat_det(rop, M, n, pN);
-
-            flint_free(M);
-        }
+        fmpz_mod_ctx_t ctx;
+        fmpz_mod_mat_t M;
+        slong n = d + len - 1, i, k;
+
+        /* todo: should use fmpz_mod_poly_resultant when that exists */
+        fmpz_mod_ctx_init(ctx, pN);
+        fmpz_mod_mat_init(M, n, n, ctx);
+
+        for (k = 0; k < len-1; k++)
+            for (i = 0; i < lena; i++)
+                fmpz_mod_set_fmpz(fmpz_mod_mat_entry(M, k, k + d - j[i]), a + i, ctx);
+        for (k = 0; k < d; k++)
+            for (i = 0; i < len; i++)
+                fmpz_mod_set_fmpz(fmpz_mod_mat_entry(M, len - 1 + k, k + len -1 - i), op + i, ctx);
+
+        fmpz_mod_mat_det(rop, M, ctx);
+
+        fmpz_mod_mat_clear(M, ctx);
+        fmpz_mod_ctx_clear(ctx);
 
         /*
             XXX:  This part of the code is currently untested as the Conway