Merge pull request #2179 from fredrik-johansson/matrix3

Use delayed canonicalization in fmpz_mod_mat_reduce_row

src/fmpz_mod_mat.h

@@ -231,7 +231,7 @@ void fmpz_mod_mat_det(fmpz_t res, const fmpz_mod_mat_t mat, const fmpz_mod_ctx_t
 
 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);
+int fmpz_mod_mat_reduce_row(slong * column, fmpz_mod_mat_t A, slong * P, slong * L, slong m, const fmpz_mod_ctx_t ctx);
 
 slong fmpz_mod_mat_lu(slong * P, fmpz_mod_mat_t A, int rank_check, const fmpz_mod_ctx_t ctx);
 

src/fmpz_mod_mat/reduce_row.c

@@ -1,6 +1,7 @@
 /*
     Copyright (C) 2015 William Hart
     Copyright (C) 2021 Daniel Schultz
+    Copyright (C) 2025 Fredrik Johansson
 
     This file is part of FLINT.
 
@@ -10,52 +11,73 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
+#include "gr.h"
 #include "fmpz.h"
 #include "fmpz_mod.h"
 #include "fmpz_mod_mat.h"
 
-slong fmpz_mod_mat_reduce_row(fmpz_mod_mat_t A, slong * P, slong * L,
-                                             slong m, const fmpz_mod_ctx_t ctx)
+int fmpz_mod_mat_reduce_row(slong * column, fmpz_mod_mat_t A, slong * P, slong * L,
+                                         slong m, const fmpz_mod_ctx_t ctx)
 {
-    slong n = fmpz_mod_mat_ncols(A, ctx), i, j, r;
-    fmpz_t h;
+    slong n = fmpz_mod_mat_ncols(A, ctx), i, j, r, res = -WORD(1);
+    fmpz_t d, h;
+    int status = GR_SUCCESS;
 
+    fmpz_init(d);
     fmpz_init(h);
 
     for (i = 0; i < n; i++)
     {
+        if (i != 0)
+            fmpz_mod_set_fmpz(fmpz_mod_mat_entry(A, m, i),
+                              fmpz_mod_mat_entry(A, m, i), ctx);
+
         if (!fmpz_is_zero(fmpz_mod_mat_entry(A, m, i)))
         {
             r = P[i];
+
             if (r != -WORD(1))
             {
                 for (j = i + 1; j < L[r]; j++)
                 {
-                    fmpz_mod_mul(h, fmpz_mod_mat_entry(A, r, j),
-                                    fmpz_mod_mat_entry(A, m, i), ctx);
-                    fmpz_mod_sub(fmpz_mod_mat_entry(A, m, j),
-                                 fmpz_mod_mat_entry(A, m, j), h, ctx);
+                    fmpz_submul(fmpz_mod_mat_entry(A, m, j),
+                        fmpz_mod_mat_entry(A, r, j),
+                        fmpz_mod_mat_entry(A, m, i));
                 }
 
                 fmpz_zero(fmpz_mod_mat_entry(A, m, i));
             }
             else
             {
-                fmpz_mod_inv(h, fmpz_mod_mat_entry(A, m, i), ctx);
-                fmpz_one(fmpz_mod_mat_entry(A, m, i));
+                /* fmpz_mod_inv(h, fmpz_mod_mat_entry(A, m, i), ctx); */
+                fmpz_gcdinv(d, h, fmpz_mod_mat_entry(A, m, i), ctx->n);
+                if (!fmpz_is_one(d))
+                {
+                    status = GR_DOMAIN;
+                    goto cleanup;
+                }
+
+                fmpz_mod_set_ui(fmpz_mod_mat_entry(A, m, i), 1, ctx);
 
                 for (j = i + 1; j < L[m]; j++)
+                {
+                    fmpz_mod_set_fmpz(fmpz_mod_mat_entry(A, m, j),
+                        fmpz_mod_mat_entry(A, m, j), ctx);
                     fmpz_mod_mul(fmpz_mod_mat_entry(A, m, j),
-                                 fmpz_mod_mat_entry(A, m, j), h, ctx);
+                        fmpz_mod_mat_entry(A, m, j), h, ctx);
+                }
 
                 P[i] = m;
+                res = i;
+                break;
 
-                fmpz_clear(h);
-                return i;
             }
         }
     }
 
+cleanup:
     fmpz_clear(h);
-    return -WORD(1);
+    fmpz_clear(d);
+    *column = res;
+    return status;
 }

src/gr/fmpz_mod.c

@@ -763,6 +763,11 @@ _gr_fmpz_mod_mat_charpoly(fmpz * res, const fmpz_mod_mat_t mat, gr_ctx_t ctx)
     return _gr_mat_charpoly_berkowitz(res, (const gr_mat_struct *) mat, ctx);
 }
 
+int
+_gr_fmpz_mod_mat_reduce_row(slong * column, fmpz_mod_mat_t mat, slong * P, slong * L, slong n, gr_ctx_t ctx)
+{
+    return fmpz_mod_mat_reduce_row(column, mat, P, L, n, FMPZ_MOD_CTX(ctx));
+}
 
 
 int _fmpz_mod_methods_initialized = 0;
@@ -856,6 +861,7 @@ gr_method_tab_input _fmpz_mod_methods_input[] =
     {GR_METHOD_MAT_LU,          (gr_funcptr) _gr_fmpz_mod_mat_lu},
     {GR_METHOD_MAT_DET,         (gr_funcptr) _gr_fmpz_mod_mat_det},
     {GR_METHOD_MAT_CHARPOLY,    (gr_funcptr) _gr_fmpz_mod_mat_charpoly},
+    {GR_METHOD_MAT_REDUCE_ROW,  (gr_funcptr) _gr_fmpz_mod_mat_reduce_row},
     {0,                         (gr_funcptr) NULL},
 };
 

src/gr/init_random.c

@@ -16,6 +16,7 @@
 #include "long_extras.h"
 #include "mpoly.h"
 #include "gr.h"
+#include "mpn_mod.h"
 
 /* For random composite rings, some base rings that don't require
    memory allocation. */
@@ -80,7 +81,7 @@ gr_ctx_init_random_ring_integers_mod(gr_ctx_t ctx, flint_rand_t state)
 {
     fmpz_t t;
 
-    switch (n_randint(state, 4))
+    switch (n_randint(state, 5))
     {
         case 0:
             gr_ctx_init_nmod8(ctx, n_randtest(state) % 255 + 1);
@@ -98,6 +99,16 @@ gr_ctx_init_random_ring_integers_mod(gr_ctx_t ctx, flint_rand_t state)
             gr_ctx_init_fmpz_mod(ctx, t);
             fmpz_clear(t);
             break;
+        case 4:
+            while (1)
+            {
+                gr_ctx_init_mpn_mod_randtest(ctx, state);
+                if (MPN_MOD_CTX_NLIMBS(ctx) <= 5)
+                    break;
+                gr_ctx_clear(ctx);
+            }
+            break;
+
     }
 }
 

src/gr_mat/minpoly_field.c

@@ -47,7 +47,6 @@ gr_mat_minpoly_field(gr_poly_t p, const gr_mat_t X, gr_ctx_t ctx)
 
     GR_TMP_INIT2(t, h, ctx);
 
-    gr_init(h, ctx);
     gr_poly_init(b, ctx);
     gr_poly_init(g, ctx);
     gr_poly_init(r, ctx);

src/gr_poly/test/t-gcd.c

@@ -13,6 +13,8 @@
 #include "ulong_extras.h"
 #include "gr_poly.h"
 
+FLINT_DLL extern gr_static_method_table _ca_methods;
+
 TEST_FUNCTION_START(gr_poly_gcd, state)
 {
     slong iter;
@@ -35,7 +37,7 @@ TEST_FUNCTION_START(gr_poly_gcd, state)
         gr_poly_init(G, ctx);
 
         n = 6;
-        if (ctx->which_ring == GR_CTX_CC_CA || ctx->which_ring == GR_CTX_RR_CA)
+        if (ctx->methods == _ca_methods)
             n = 3;
 
         status = gr_poly_randtest(A, state, n, ctx);

src/gr_poly/test/t-gcd_euclidean.c

@@ -13,6 +13,8 @@
 #include "ulong_extras.h"
 #include "gr_poly.h"
 
+FLINT_DLL extern gr_static_method_table _ca_methods;
+
 TEST_FUNCTION_START(gr_poly_gcd_euclidean, state)
 {
     slong iter;
@@ -35,7 +37,7 @@ TEST_FUNCTION_START(gr_poly_gcd_euclidean, state)
         gr_poly_init(G, ctx);
 
         n = 6;
-        if (ctx->which_ring == GR_CTX_CC_CA || ctx->which_ring == GR_CTX_RR_CA)
+        if (ctx->methods == _ca_methods)
             n = 3;
 
         status = gr_poly_randtest(A, state, n, ctx);

src/gr_poly/test/t-resultant_hgcd.c

@@ -57,7 +57,7 @@ TEST_FUNCTION_START(gr_poly_resultant_hgcd, state)
             cutoff1 = n_randint(state, 100);
             cutoff2 = n_randint(state, 100);
         }
-        else if (ctx->which_ring == GR_CTX_CC_CA || ctx->which_ring == GR_CTX_RR_CA)
+        else if (ctx->methods == _ca_methods)
         {
             n = n_randint(state, 3);
             cutoff1 = n_randint(state, 3);

src/gr_poly/test/t-xgcd.c

@@ -16,6 +16,8 @@
 #include "ulong_extras.h"
 #include "gr_poly.h"
 
+FLINT_DLL extern gr_static_method_table _ca_methods;
+
 TEST_FUNCTION_START(gr_poly_xgcd, state)
 {
     int i;
@@ -33,7 +35,7 @@ TEST_FUNCTION_START(gr_poly_xgcd, state)
 
         if (gr_ctx_is_finite(ctx) == T_TRUE && n_randint(state, 2) == 0)
             n = 20;
-        else if (ctx->which_ring == GR_CTX_CC_CA || ctx->which_ring == GR_CTX_RR_CA)
+        else if (ctx->methods == _ca_methods)
             n = 4;
         else
             n = 6;

src/gr_poly/test/t-xgcd_euclidean.c

@@ -16,6 +16,8 @@
 #include "ulong_extras.h"
 #include "gr_poly.h"
 
+FLINT_DLL extern gr_static_method_table _ca_methods;
+
 TEST_FUNCTION_START(gr_poly_xgcd_euclidean, state)
 {
     int i;
@@ -33,7 +35,7 @@ TEST_FUNCTION_START(gr_poly_xgcd_euclidean, state)
 
         if (gr_ctx_is_finite(ctx) == T_TRUE && n_randint(state, 2) == 0)
             n = 20;
-        else if (ctx->which_ring == GR_CTX_CC_CA || ctx->which_ring == GR_CTX_RR_CA)
+        else if (ctx->methods == _ca_methods)
             n = 4;
         else
             n = 6;

src/gr_poly/test/t-xgcd_hgcd.c

@@ -16,6 +16,8 @@
 #include "ulong_extras.h"
 #include "gr_poly.h"
 
+FLINT_DLL extern gr_static_method_table _ca_methods;
+
 TEST_FUNCTION_START(gr_poly_xgcd_hgcd, state)
 {
     int i;
@@ -40,7 +42,7 @@ TEST_FUNCTION_START(gr_poly_xgcd_hgcd, state)
             cutoff1 = n_randint(state, 100);
             cutoff2 = n_randint(state, 100);
         }
-        else if (ctx->which_ring == GR_CTX_CC_CA || ctx->which_ring == GR_CTX_RR_CA)
+        else if (ctx->methods == _ca_methods)
             n = 4;
         else
             n = 6;

src/mpn_mod.h

@@ -212,6 +212,7 @@ int mpn_mod_mat_lu_classical_delayed(slong * res_rank, slong * P, gr_mat_t A, co
 int mpn_mod_mat_lu(slong * rank, slong * P, gr_mat_t LU, const gr_mat_t A, int rank_check, gr_ctx_t ctx);
 int mpn_mod_mat_det(nn_ptr res, const gr_mat_t A, gr_ctx_t ctx);
 int _mpn_mod_mat_charpoly(nn_ptr res, const gr_mat_t mat, gr_ctx_t ctx);
+int mpn_mod_mat_reduce_row(slong * column, gr_mat_t A, slong * P, slong * L, slong m, gr_ctx_t ctx);
 
 /* Polynomial algorithms */
 

src/mpn_mod/ctx.c

@@ -134,6 +134,7 @@ gr_method_tab_input _mpn_mod_methods_input[] =
     {GR_METHOD_MAT_LU,          (gr_funcptr) mpn_mod_mat_lu},
     {GR_METHOD_MAT_DET,         (gr_funcptr) mpn_mod_mat_det},
     {GR_METHOD_MAT_CHARPOLY,    (gr_funcptr) _mpn_mod_mat_charpoly},
+    {GR_METHOD_MAT_REDUCE_ROW,  (gr_funcptr) mpn_mod_mat_reduce_row},
     {0,                         (gr_funcptr) NULL},
 };
 #if defined(__GNUC__)

src/mpn_mod/mat_reduce_row.c

@@ -0,0 +1,112 @@
+/*
+    Copyright (C) 2025 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 3 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include "mpn_mod.h"
+
+/* todo: optimize for when 2n or 2n-1 rather than 2n+1 limbs suffice */
+int mpn_mod_mat_reduce_row(slong * column, gr_mat_t A, slong * P, slong * L,
+                                         slong m, gr_ctx_t ctx)
+{
+    slong n = A->c, i, j, r, res = -WORD(1);
+    int status = GR_SUCCESS;
+    slong Astride = A->stride;
+    slong nlimbs = MPN_MOD_CTX_NLIMBS(ctx);
+    slong tnlimbs = 2 * nlimbs + 1;
+    nn_ptr aa = A->entries;
+    nn_ptr tmp, h, u;
+    slong alloc;
+
+    TMP_INIT;
+    TMP_START;
+
+    alloc = (n * tnlimbs) + (2 * nlimbs) + (nlimbs);
+
+    tmp = TMP_ALLOC(alloc * sizeof(ulong));
+    u = tmp + (n * tnlimbs);
+    h = u + (2 * nlimbs);
+
+#define ENTRY(ii, jj) (aa + (((ii) * Astride + (jj)) * nlimbs))
+#define TMP(ii) (tmp + tnlimbs * (ii))
+
+    for (i = 0; i < n; i++)
+    {
+        flint_mpn_copyi(TMP(i), ENTRY(m, i), nlimbs);
+        flint_mpn_zero(TMP(i) + nlimbs, tnlimbs - nlimbs);
+    }
+
+    for (i = 0; i < n; i++)
+    {
+        if (i != 0)
+        {
+            mpn_mod_set_mpn(ENTRY(m, i), TMP(i), tnlimbs, ctx);
+        }
+
+        if (mpn_mod_is_zero(ENTRY(m, i), ctx) == T_FALSE)
+        {
+            r = P[i];
+
+            if (r != -WORD(1))
+            {
+                mpn_mod_neg(h, ENTRY(m, i), ctx);
+                mpn_mod_zero(ENTRY(m, i), ctx);
+
+                if (nlimbs == 2)
+                {
+                    for (j = i + 1; j < L[r]; j++)
+                    {
+                        ulong t[4];
+                        nn_srcptr ap = ENTRY(r, j);
+                        nn_ptr s = TMP(j);
+
+                        FLINT_MPN_MUL_2X2(t[3], t[2], t[1], t[0], ap[1], ap[0], h[1], h[0]);
+                        add_sssssaaaaaaaaaa(s[4], s[3], s[2], s[1], s[0],
+                                            s[4], s[3], s[2], s[1], s[0],
+                                            0, t[3], t[2], t[1], t[0]);
+                    }
+                }
+                else
+                {
+                    for (j = i + 1; j < L[r]; j++)
+                    {
+                        flint_mpn_mul_n(u, ENTRY(r, j), h, nlimbs);
+                        TMP(j)[tnlimbs - 1] += mpn_add_n(TMP(j), TMP(j), u, 2 * nlimbs);
+                    }
+                }
+            }
+            else
+            {
+                status = mpn_mod_inv(h, ENTRY(m, i), ctx);
+                if (status != GR_SUCCESS)
+                    goto cleanup;
+
+                mpn_mod_one(ENTRY(m, i), ctx);
+
+                for (j = i + 1; j < L[m]; j++)
+                {
+                    mpn_mod_set_mpn(ENTRY(m, j), TMP(j), tnlimbs, ctx);
+                    mpn_mod_mul(ENTRY(m, j), ENTRY(m, j), h, ctx);
+                }
+
+                P[i] = m;
+                res = i;
+                break;
+
+            }
+        }
+    }
+
+cleanup:
+
+    TMP_END;
+
+    *column = res;
+    return status;
+}