From b1d7b8632989b62ae26ff22c445a96c66f8a8076 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Thu, 23 Jan 2025 17:28:03 +0100 Subject: [PATCH 1/4] Use delayed canonicalization in fmpz_mod_mat_reduce_row; also handle error for gr --- src/fmpz_mod_mat.h | 2 +- src/fmpz_mod_mat/reduce_row.c | 50 +++++++++++++++++++++++++---------- src/gr/fmpz_mod.c | 6 +++++ 3 files changed, 43 insertions(+), 15 deletions(-) diff --git a/src/fmpz_mod_mat.h b/src/fmpz_mod_mat.h index 96a10d387a..eec69fb672 100644 --- a/src/fmpz_mod_mat.h +++ b/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); diff --git a/src/fmpz_mod_mat/reduce_row.c b/src/fmpz_mod_mat/reduce_row.c index 6d4d625af1..d82d66f157 100644 --- a/src/fmpz_mod_mat/reduce_row.c +++ b/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 . */ +#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; } diff --git a/src/gr/fmpz_mod.c b/src/gr/fmpz_mod.c index d0e3f4aa5b..4db2090ea8 100644 --- a/src/gr/fmpz_mod.c +++ b/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}, }; From 1c1de5c434ed46e7852fbb4cb2931e8a20aa1815 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Thu, 23 Jan 2025 18:47:29 +0100 Subject: [PATCH 2/4] fix memory leak in gr_mat_minpoly_field --- src/gr_mat/minpoly_field.c | 1 - 1 file changed, 1 deletion(-) diff --git a/src/gr_mat/minpoly_field.c b/src/gr_mat/minpoly_field.c index b3692c3826..5f501e2407 100644 --- a/src/gr_mat/minpoly_field.c +++ b/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); From f9764e58942865f99996f966ede9ef7ff28898f4 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Thu, 23 Jan 2025 18:55:45 +0100 Subject: [PATCH 3/4] Delayed canonicalization mpn_mod_mat_reduce_row --- src/gr/init_random.c | 13 +++- src/mpn_mod.h | 1 + src/mpn_mod/ctx.c | 1 + src/mpn_mod/mat_reduce_row.c | 112 +++++++++++++++++++++++++++++++++++ 4 files changed, 126 insertions(+), 1 deletion(-) create mode 100644 src/mpn_mod/mat_reduce_row.c diff --git a/src/gr/init_random.c b/src/gr/init_random.c index f056ba1adb..849c9b5bfd 100644 --- a/src/gr/init_random.c +++ b/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; + } } diff --git a/src/mpn_mod.h b/src/mpn_mod.h index 92de076c46..f6b13a49d2 100644 --- a/src/mpn_mod.h +++ b/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 */ diff --git a/src/mpn_mod/ctx.c b/src/mpn_mod/ctx.c index ea69979ab4..f39e361e8f 100644 --- a/src/mpn_mod/ctx.c +++ b/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__) diff --git a/src/mpn_mod/mat_reduce_row.c b/src/mpn_mod/mat_reduce_row.c new file mode 100644 index 0000000000..ca779e85a8 --- /dev/null +++ b/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 . +*/ + +#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; +} From e4738ce55eb9ebdb6ccf140e04263c5b045c2ba2 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Thu, 23 Jan 2025 23:55:12 +0100 Subject: [PATCH 4/4] safeguard more tests against slow calcium contexts --- src/gr_poly/test/t-gcd.c | 4 +++- src/gr_poly/test/t-gcd_euclidean.c | 4 +++- src/gr_poly/test/t-resultant_hgcd.c | 2 +- src/gr_poly/test/t-xgcd.c | 4 +++- src/gr_poly/test/t-xgcd_euclidean.c | 4 +++- src/gr_poly/test/t-xgcd_hgcd.c | 4 +++- 6 files changed, 16 insertions(+), 6 deletions(-) diff --git a/src/gr_poly/test/t-gcd.c b/src/gr_poly/test/t-gcd.c index 1559b7397f..e90d97d8db 100644 --- a/src/gr_poly/test/t-gcd.c +++ b/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); diff --git a/src/gr_poly/test/t-gcd_euclidean.c b/src/gr_poly/test/t-gcd_euclidean.c index aedf9d8207..4e65470e75 100644 --- a/src/gr_poly/test/t-gcd_euclidean.c +++ b/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); diff --git a/src/gr_poly/test/t-resultant_hgcd.c b/src/gr_poly/test/t-resultant_hgcd.c index ea944d463e..81284b6a95 100644 --- a/src/gr_poly/test/t-resultant_hgcd.c +++ b/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); diff --git a/src/gr_poly/test/t-xgcd.c b/src/gr_poly/test/t-xgcd.c index 92217de49f..4620c4dcba 100644 --- a/src/gr_poly/test/t-xgcd.c +++ b/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; diff --git a/src/gr_poly/test/t-xgcd_euclidean.c b/src/gr_poly/test/t-xgcd_euclidean.c index a6a357c839..7dcb266ad9 100644 --- a/src/gr_poly/test/t-xgcd_euclidean.c +++ b/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; diff --git a/src/gr_poly/test/t-xgcd_hgcd.c b/src/gr_poly/test/t-xgcd_hgcd.c index cdf05bc663..b1f9bd87ca 100644 --- a/src/gr_poly/test/t-xgcd_hgcd.c +++ b/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;