From 83bc9fee2b5d403fee68199d82121d530677420a Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sat, 3 Jan 2026 11:51:23 +0100 Subject: [PATCH 01/12] finalize the removal of mbasis1 --- flint-extras/src/nmod_mat_poly.h | 11 ++- flint-extras/src/nmod_poly_mat_approximant.h | 68 +-------------- .../nmod_poly_mat_extra/test/test_mbasis1.c | 85 ------------------- 3 files changed, 6 insertions(+), 158 deletions(-) delete mode 100644 flint-extras/src/nmod_poly_mat_extra/test/test_mbasis1.c diff --git a/flint-extras/src/nmod_mat_poly.h b/flint-extras/src/nmod_mat_poly.h index e351a35a..08065d1a 100644 --- a/flint-extras/src/nmod_mat_poly.h +++ b/flint-extras/src/nmod_mat_poly.h @@ -587,7 +587,7 @@ nmod_mat_poly_set_from_poly_mat(nmod_mat_poly_t matp, const nmod_poly_mat_t pmat * The functions here compute a `shift`-minimal ordered weak Popov approximant * basis for `(pmat,orders)` in the case where `orders` is given by a single * integer `orders = (order,...order)`. They iterate from `1` to `order`, - * computing at each step a basis at order `1` (see @ref mbasis1) and using it + * computing at each step a basis at order `1` and using it * to update the output `appbas`, the so-called _residual matrix_, and the * considered shift. At step `d`, we have `appbas*pmat = 0 mod x^{d-1}`, and we * want to update `appbas` so that this becomes zero modulo `x^d`. @@ -614,8 +614,8 @@ nmod_mat_poly_set_from_poly_mat(nmod_mat_poly_t matp, const nmod_poly_mat_t pmat * \todo integrate */ // Complexity: pmat is m x n -// - 'order' calls to mbasis1 with dimension m x n, each one gives a -// constant matrix K which is generically m-n x m (may have more rows in +// - 'order' calls to constant nullspace with dimension m x n, each one gives +// a constant matrix K which is generically m-n x m (may have more rows in // exceptional cases) // - order products (X Id + K ) * appbas to update the approximant basis // - order computations of "coeff k of appbas*pmat" to find residuals @@ -639,8 +639,8 @@ nmod_mat_poly_set_from_poly_mat(nmod_mat_poly_t matp, const nmod_poly_mat_t pmat // Residual (X^-d appbas*pmat mod X^(order-d)) is continuously updated along // the iterations // Complexity: pmat is m x n -// - 'order' calls to mbasis1 with dimension m x n, each one gives a -// constant matrix K which is generically m-n x m (may have more rows in +// - 'order' calls to constant nullspace with dimension m x n, each one gives +// a constant matrix K which is generically m-n x m (may have more rows in // exceptional cases) // - order products (X Id + K ) * appbas to update the approximant basis // - order-1 products (X Id + K ) * (matrix of degree order-ord) to update @@ -675,7 +675,6 @@ nmod_mat_poly_set_from_poly_mat(nmod_mat_poly_t matp, const nmod_poly_mat_t pmat // // To understand the threshold (cdim > rdim/2 + 1), see the complexities // // mentioned above for these two variants of mbasis //} -// TODO DOC (see @ref mbasis1) // TODO resupdate version // TODO general version with choice void nmod_mat_poly_mbasis(nmod_mat_poly_t appbas, diff --git a/flint-extras/src/nmod_poly_mat_approximant.h b/flint-extras/src/nmod_poly_mat_approximant.h index 6f99d5d7..387bc7a6 100644 --- a/flint-extras/src/nmod_poly_mat_approximant.h +++ b/flint-extras/src/nmod_poly_mat_approximant.h @@ -243,82 +243,16 @@ int nmod_poly_mat_is_approximant_basis(const nmod_poly_mat_t appbas, /** @name Approximant basis via linear algebra - * \anchor mbasis1 + * \anchor approx_linalg * * These functions compute a shifted minimal or shifted Popov approximant * basis using fast linear algebra on constant matrices. * - * Currently, this is only implemented for the uniform order `(1,...,1)`, following - * the algorithm _mbasis_ described in - * - P. Giorgi, C.-P. Jeannerod, G. Villard. Proceeding ISSAC 2003, - * - P. Giorgi, R. Lebreton. Proceedings ISSAC 2014, - * - C.-P. Jeannerod, V. Neiger, G. Villard. Preprint 2018. - * - * The latter reference explicitly shows how to ensure that we obtain the - * canonical s-Popov approximant basis. - * * \todo implement the general Krylov-based approach from JNSV17, which * efficiently solves the problem when the sum of the orders is small */ //@{ -/** Computes a shifted Popov approximant basis at order `(1,...,1)` using fast - * linear algebra on constant matrices. Thus, in this context, the input matrix - * `pmat` is a constant matrix. This approximant basis consists of two sets of - * rows: - * - rows forming a left kernel basis in reduced row echelon form of some - * row-permutation of `pmat` (the permutation depends on the shift) - * - the other rows are coordinate vectors multiplied by the variable `x`; - * there is one such row at each index `i` which is not among the pivot - * indices of the left kernel basis above. - * - * Here the whole approximant basis is stored in the OUT parameter `appbas`. - * - * \param[out] appbas polynomial matrix - * \param[in] pmat constant matrix - * \param[in] shift shift (vector of integers) - * - * \todo modify shift in place - * \todo does it return Popov? if yes name it popov_mbasis1 and write mbasis1 - * if that makes sense; otherwise write other version popov_mbasis1 - * \todo safety correctness checks - * \todo is output Popov by this method? check - * \todo seems that it does not return owP ? (see test_mbasis1.c) - */ -void mbasis1(nmod_poly_mat_t appbas, - slong * res_shift, - const nmod_mat_t mat, - const slong * shift); - -/** Computes a shifted Popov approximant basis at order `(1,...,1)` using fast - * linear algebra on constant matrices. Thus, in this context, the input matrix - * `pmat` is a constant matrix. This approximant basis consists of two sets of - * rows: - * - rows forming a left kernel basis in reduced row echelon form of some - * row-permutation of `pmat` (the permutation depends on the shift) - * - the other rows are coordinate vectors multiplied by the variable `x`; - * there is one such row at each index `i` which is not among the pivot - * indices of the left kernel basis above. - * - * Here only the first set of rows is stored in the OUT parameter `kerbas`, - * along with permutation... - * \todo improve doc, to say exactly what this computes (and what is the return - * value?) - * - * \param[out] kerbas constant matrix - * \param[in] pmat constant matrix - * \param[in] shift shift (vector of integers) - * - * \todo modify shift in place - * \todo safety correctness checks - * \todo is output Popov by this method? check - */ -slong mbasis1_for_mbasis(nmod_mat_t kerbas, - slong * res_shift, - slong * res_perm, - const nmod_mat_t mat, - const slong * shift); - //@} // doxygen group: Approximant basis via linear algebra /** @name M-Basis algorithm (uniform approximant order) diff --git a/flint-extras/src/nmod_poly_mat_extra/test/test_mbasis1.c b/flint-extras/src/nmod_poly_mat_extra/test/test_mbasis1.c deleted file mode 100644 index d17cd011..00000000 --- a/flint-extras/src/nmod_poly_mat_extra/test/test_mbasis1.c +++ /dev/null @@ -1,85 +0,0 @@ -/* - Copyright (C) 2025 Vincent Neiger, Kevin Tran - - This file is part of PML. - - PML is free software: you can redistribute it and/or modify it under - the terms of the GNU General Public License version 2.0 (GPL-2.0-or-later) - as published by the Free Software Foundation; either version 2 of the - License, or (at your option) any later version. See - . -*/ - -#include - -#include "nmod_poly_mat_approximant.h" -#include "nmod_poly_mat_io.h" // for print_pretty - -#define PRIME_30_BITS 536870923 -#define PRIME_60_BITS 576460752303423619 - -int test_mbasis1(void) -{ - nmod_mat_t mat; - slong rdim = 9, cdim = 3, prime = 7; - - nmod_mat_init(mat, rdim, cdim, prime); - - - flint_rand_t state; - flint_rand_init(state); - - nmod_mat_randtest(mat, state); - - - nmod_poly_mat_t res; - nmod_poly_mat_init(res, rdim, rdim, prime); - nmod_poly_mat_zero(res); - - slong res_shift[rdim]; - - slong shift[rdim]; - - for (slong i = 0; i < rdim; i++) - shift[i] = n_randint(state, 10) - 5; - - mbasis1(res, res_shift, mat, shift); - - - printf("Matrix\n"); - nmod_mat_print_pretty(mat); - printf("Result mbasis1\n"); - nmod_poly_mat_print_pretty(res,"x"); - nmod_poly_mat_clear(res); - flint_rand_clear(state); - - nmod_mat_t res2; - slong res_shift2[rdim]; - slong res_perm2[rdim]; - slong rank = mbasis1_for_mbasis(res2, res_shift2, res_perm2, mat, shift); - - nmod_mat_clear(mat); - - printf("\nMatrix A from basis = [[x,0],[A,1]]\n"); - nmod_mat_print_pretty(res2); - - printf("\nshifts = "); - for(slong i = 0; i < rdim; i++) - printf("%ld ", res_shift2[i]); - - printf("\npermutation = "); - for(slong i = 0; i < rdim; i++) - printf("%ld ", res_perm2[i]); - - - printf("\nrank = %ld\n", rank); - - nmod_mat_clear(res2); - return 0; -} - -int main(void) -{ - test_mbasis1(); - return 0; -} From 8dcc89011e68c04f601a5eca2df0dc744cc0bd12 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sat, 3 Jan 2026 12:05:15 +0100 Subject: [PATCH 02/12] kernel: general conventions/definitions --- flint-extras/src/nmod_poly_mat_kernel.h | 45 +++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/flint-extras/src/nmod_poly_mat_kernel.h b/flint-extras/src/nmod_poly_mat_kernel.h index 20cb54b6..f2f0a083 100644 --- a/flint-extras/src/nmod_poly_mat_kernel.h +++ b/flint-extras/src/nmod_poly_mat_kernel.h @@ -19,6 +19,51 @@ extern "C" { #endif +/** \file nmod_poly_mat_kernel.h + * Definition (shifted minimal kernel basis). + * ------------------------------------------ + * Consider: + * - an m x n matrix of univariate polynomials F, + * - a degree shift s (a list of m integers). + * + * A (left) kernel basis for F is a matrix over the univariate polynomials + * whose rows form a basis for the module + * { p in K[X]^{1 x m} | p * F == 0 }, + * whose rank is nz = m - rank(F), called the nullity of F. Such a basis + * matrix has dimensions nz x m and has full row rank. + * + * A kernel basis for (F,d) is said to be a shift-minimal (resp. + * a shift-ordered weak Popov, resp. the shift-Popov) + * kernel basis if it is in shift-reduced form (resp. in shift-ordered weak + * Popov form, resp. in shift-Popov form). See nmod_poly_mat_forms.h for + * definitions of these forms. + */ + +/** \file nmod_poly_mat_kernel.h + * Conventions. + * ------------ + * Apart from the general interfaces (TODO) which offer several choices of + * orientation, all other functions compute left kernel bases and use the + * following parameters: + * + * \param[out] ker the output kernel basis (cannot alias `pmat`) + * \param[in] pmat the input polynomial matrix (no restriction) + * \param[in,out] shift in: the input shift; and out: the output shifted row + * degree of `ker` (list of integers, length must be the number of rows of + * `pmat`) + * + * The computed `ker` is in shifted ordered weak Popov form, or in the canonical + * shifted Popov form when the name of the function indicates so. + */ + + + + + + + + + /** * * Right shifted kernel of a polynomial matrix, assuming that the columns has been From 82103b6e5c1fcb6c831fc5a9607a3e63ce04df74 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sat, 3 Jan 2026 20:26:27 +0100 Subject: [PATCH 03/12] add kernel via approx, and test --- .../nmod_mat_poly_mbasis.c | 3 - flint-extras/src/nmod_poly_mat_approximant.h | 5 +- flint-extras/src/nmod_poly_mat_extra.h | 20 --- flint-extras/src/nmod_poly_mat_extra/kernel.c | 150 ++++++++++++++++++ .../src/nmod_poly_mat_extra/test/main.c | 2 + .../test/t-kernel_via_approx.c | 70 ++++++++ .../nmod_poly_mat_extra/test/t-kernel_zls.c | 14 +- .../src/nmod_poly_mat_extra/verification.c | 17 +- flint-extras/src/nmod_poly_mat_kernel.h | 83 +++++++++- 9 files changed, 323 insertions(+), 41 deletions(-) create mode 100644 flint-extras/src/nmod_poly_mat_extra/test/t-kernel_via_approx.c diff --git a/flint-extras/src/nmod_mat_poly_extra/nmod_mat_poly_mbasis.c b/flint-extras/src/nmod_mat_poly_extra/nmod_mat_poly_mbasis.c index 1c737577..e5e90b01 100644 --- a/flint-extras/src/nmod_mat_poly_extra/nmod_mat_poly_mbasis.c +++ b/flint-extras/src/nmod_mat_poly_extra/nmod_mat_poly_mbasis.c @@ -376,6 +376,3 @@ _PROFILER_REGION_STOP(t_others) _MBASIS_PROFILER_OUTPUT return; } - -/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ -// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s diff --git a/flint-extras/src/nmod_poly_mat_approximant.h b/flint-extras/src/nmod_poly_mat_approximant.h index 387bc7a6..9bb8feff 100644 --- a/flint-extras/src/nmod_poly_mat_approximant.h +++ b/flint-extras/src/nmod_poly_mat_approximant.h @@ -182,14 +182,13 @@ * \param[in] order order * \param[in] shift shift * \param[in] form required form for `appbas` (see #poly_mat_form_t) - * \param[in] row_wise indicates whether to compute left approximants (working row-wise) or right approximants (working column-wise) + * \param[in] orient indicates the orientation (left/right approximants) and the definition of pivots * \param[in] randomized if `true`, the algorithm may use a Monte Carlo or Las Vegas verification algorithm * * \return boolean, result of the verification * - * \todo add parameter row_wise + * \todo update documentation * \todo support all options, make doc more clear concerning Las Vegas / Monte Carlo - * \todo WARNING! for the moment, does not really check generation! * \todo WARNING! for the moment, hardcoded to check for ordered weak Popov */ int nmod_poly_mat_is_approximant_basis(const nmod_poly_mat_t appbas, diff --git a/flint-extras/src/nmod_poly_mat_extra.h b/flint-extras/src/nmod_poly_mat_extra.h index b6d675ae..c4bfa8ef 100644 --- a/flint-extras/src/nmod_poly_mat_extra.h +++ b/flint-extras/src/nmod_poly_mat_extra.h @@ -86,24 +86,4 @@ void nmod_poly_mat_det_iter(nmod_poly_t det, nmod_poly_mat_t mat); */ -/***************************** -* Verification algorithms * -*****************************/ - -int nmod_poly_mat_is_approximant_basis(const nmod_poly_mat_t appbas, - const nmod_poly_mat_t pmat, - slong order, - const slong * shift, - orientation_t orient); - -int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker, - const nmod_poly_mat_t pmat, - const slong * shift, - orientation_t orient); - - - #endif // NMOD_POLY_MAT_EXTRA_H - -/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ -// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s diff --git a/flint-extras/src/nmod_poly_mat_extra/kernel.c b/flint-extras/src/nmod_poly_mat_extra/kernel.c index f46208d3..f285152f 100644 --- a/flint-extras/src/nmod_poly_mat_extra/kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/kernel.c @@ -1,5 +1,6 @@ /* Copyright (C) 2025 Gilles Villard + Copyright (C) 2026 Vincent Neiger This file is part of PML. @@ -10,12 +11,125 @@ . */ +#include +#include +#include +#include #include #include #include "nmod_poly_mat_extra.h" #include "nmod_poly_mat_kernel.h" +/* /1* requires len >= 1 *1/ */ +/* FLINT_INLINE slong _slong_vec_amplitude(slong * vec, slong len) */ +/* { */ +/* slong min = vec[0]; */ +/* slong max = vec[0]; */ +/* for (slong i = 1; i < len; i++) */ +/* { */ +/* slong v = vec[i]; */ +/* if (v < min) */ +/* min = v; */ +/* if (v > max) */ +/* max = v; */ +/* } */ +/* return max - min; */ +/* } */ + +slong nmod_poly_mat_kernel_via_approx(nmod_poly_mat_t ker, + slong * pivind, + slong * shift, + const nmod_poly_mat_t pmat) +{ + const slong m = pmat->r; + const slong n = pmat->c; + const slong d = nmod_poly_mat_degree(pmat); + + /* pmat == 0 => kernel is identity*/ + if (d == -1) + { + nmod_poly_mat_one(ker); + for (slong i = 0; i < m; i++) + pivind[i] = i; + return m; + } + + /* pmat full row rank => empty kernel */ + /* full row rank implied by m == 1 or (m <= n && constant coefficient has full rank) */ + if (m == 1) + { + nmod_poly_mat_zero(ker); + return 0; + } + + if (m <= n) + { + nmod_mat_t coeff0; + nmod_mat_init(coeff0, m, n, pmat->modulus); + nmod_poly_mat_get_coeff_mat(coeff0, pmat, 0); + /* could use: + * const slong r = nmod_mat_rank(coeff0); + * but let's call more directly its core computation + * (this uses pivind as temporary of length m for row permutation) + */ + const slong r = nmod_mat_lu(pivind, coeff0, 1); + if (r == m) + { + nmod_poly_mat_zero(ker); + return 0; + } + } + + /* compute amplitude of `shift`, make `pivind` a temporary copy of `shift` */ + slong min = shift[0]; + slong amp = shift[0]; + for (slong i = 0; i < m; i++) + { + slong s = shift[i]; + pivind[i] = s; + if (s < min) + min = s; + else if (s > amp) + amp = s; + } + amp = amp - min; + + /* order for approximation: bound on rdeg_s(ker) + deg(pmat) + 1 */ + /* (see notes at end of file for rdeg_s(ker)) */ + const slong order = FLINT_MIN(m, n) * d + amp + d + 1; + + /* compute approximant basis and shifted row degree */ + nmod_poly_mat_t appbas; + nmod_poly_mat_init(appbas, m, m, pmat->modulus); + nmod_poly_mat_pmbasis(appbas, pivind, pmat, order); + + /* find rows which belong to the kernel and record their pivot information */ + /* note: at this stage, pivind == rdeg_s(appbas) */ + slong nz = 0; + for (slong i = 0; i < m; i++) + { + if (pivind[i] - shift[i] + amp < order - d) + { + shift[nz] = pivind[i]; + pivind[nz] = i; + nz += 1; + } + } + + /* swap the kernel rows */ + for (slong i = 0; i < nz; i++) + for (slong j = 0; j < m; j++) + FLINT_SWAP(nmod_poly_struct, + *nmod_poly_mat_entry(ker, i, j), + *nmod_poly_mat_entry(appbas, pivind[i], j)); + + nmod_poly_mat_clear(appbas); + + return nz; +} + + /** * * The content of the file has not yet been tested very extensively @@ -690,3 +804,39 @@ int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong *degN, const nmod_ nmod_poly_mat_clear(PT); return 0; } + + + + +/*------------------------------------------------------------*/ +/* NOTES */ +/*------------------------------------------------------------*/ + +/** Degree bounds for kernel bases. + * + * Assume pmat is m x n and has rank r. Let P be some shifted ordered weak + * Popov kernel basis P of `pmat` (for a given, arbitrary shift). + * + * From Theorem 3.3 in + * Labahn, Neiger, Vu, Zhou + * Proceedings ISSAC 2022 (https://arxiv.org/pdf/2202.09329) + * we get non-strict upper bounds on the sum of pivot degree of P: + * - it is <= sum of degrees of the r largest-degree rows of pmat + * - in particular, <= r * deg(pmat) <= min(m, n) * deg(pmat) + * - it is <= sum of degrees of the r largest-degree columns of pmat + * - in particular, <= sum of cdeg(pmat) + * The third item above is not found directly in the mentioned Theorem 3.3, + * but is deduced from its second item (part with deg(det(S))) + * + * The maximum of pivot degrees can be equal to their sum (although this is not + * expected generically). This gives bounds on the non-pivot part of P: pick + * one of the bounds above, and add the amplitude `max(shift) - min(shift)`. + * Note that this might be pessimistic for shifts of large amplitude, since we + * also have bounds coming from formula with adjugate of some submatrix of `pmat`. + * For example, one can prove that the non-pivot part of the kernel basis in + * shifted Popov form has degree at most n * deg(pmat) (cf Lemma 12 of a research + * internship report by Vu Thi Xuan). + * + * The maximum pivot degree also gives a bound on max(rdeg_s(ker)): pick one of + * the bounds above for sum of pivot degrees, and add max(s) - min(s). + */ diff --git a/flint-extras/src/nmod_poly_mat_extra/test/main.c b/flint-extras/src/nmod_poly_mat_extra/test/main.c index e89828cd..ec0cec9d 100644 --- a/flint-extras/src/nmod_poly_mat_extra/test/main.c +++ b/flint-extras/src/nmod_poly_mat_extra/test/main.c @@ -17,6 +17,7 @@ #include "t-det.c" #include "t-dixon.c" #include "t-hermite_normal_form.c" +#include "t-kernel_via_approx.c" #include "t-kernel_zls.c" #include "t-middle_product_geometric.c" #include "t-mul_geometric.c" @@ -33,6 +34,7 @@ test_struct tests[] = TEST_FUNCTION(nmod_poly_mat_det), TEST_FUNCTION(nmod_poly_mat_dixon), TEST_FUNCTION(nmod_poly_mat_hnf), + TEST_FUNCTION(nmod_poly_mat_kernel_via_approx), TEST_FUNCTION(nmod_poly_mat_kernel_zls), TEST_FUNCTION(nmod_poly_mat_mbasis), TEST_FUNCTION(nmod_poly_mat_middle_product_geometric), diff --git a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_via_approx.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_via_approx.c new file mode 100644 index 00000000..d6ef8acb --- /dev/null +++ b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_via_approx.c @@ -0,0 +1,70 @@ +/* + Copyright (C) 2025 Gilles Villard + + This file is part of PML. + + PML is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License version 2.0 (GPL-2.0-or-later) + as published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. See + . +*/ + +#include +#include +#include +#include +#include +#include + +#include "nmod_poly_mat_kernel.h" + +TEST_FUNCTION_START(nmod_poly_mat_kernel_via_approx, state) +{ + int i, result; + + for (i = 0; i < 16 * flint_test_multiplier(); i++) + { + ulong nbits = 2 + n_randint(state, 63); + slong rdim = n_randint(state, 30); + slong cdim = n_randint(state, 30); + ulong len = 1 + n_randint(state, 100); + flint_printf("TEST iter %ld -- %ld, %ld, %ld\n", i, rdim, cdim, len); + + slong * shift = FLINT_ARRAY_ALLOC(rdim, slong); + for (slong i = 0; i < rdim; i++) + shift[i] = 0; + + slong * pivind = FLINT_ARRAY_ALLOC(rdim, slong); + slong * rdeg = FLINT_ARRAY_ALLOC(rdim, slong); + for (slong i = 0; i < rdim; i++) + rdeg[i] = shift[i]; + + ulong prime = n_randprime(state, nbits, 1); + + nmod_poly_mat_t pmat; + nmod_poly_mat_init(pmat, rdim, cdim, prime); + nmod_poly_mat_randtest(pmat, state, len); + + nmod_poly_mat_t ker; + nmod_poly_mat_init(ker, rdim, rdim, prime); + /* TODO test weak Popov + pivind */ + /* TODO currently does not test generation */ + slong nz = nmod_poly_mat_kernel_via_approx(ker, pivind, rdeg, pmat); + flint_printf("kernel computed, now testing...\n"); + result = nmod_poly_mat_is_kernel(ker, nz, shift, pmat, ROW_LOWER); + + nmod_poly_mat_clear(pmat); + flint_free(shift); + flint_free(rdeg); + flint_free(pivind); + + if (!result) + { + TEST_FUNCTION_FAIL("rdim = %wu, cdim = %wu, length = %wu, p = %wu\n", \ + rdim, cdim, len, prime); + } + } + + TEST_FUNCTION_END(state); +} diff --git a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls.c index 94a6edbe..0e7358df 100644 --- a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls.c +++ b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls.c @@ -49,7 +49,7 @@ int core_test_kernel_zls(const nmod_poly_mat_t mat) nmod_poly_mat_init(Mt, n, m, mat->modulus); nmod_poly_mat_transpose(Mt, mat); - int verif = nmod_poly_mat_is_kernel(Nt, Mt, rdeg, ROW_LOWER); + int verif = nmod_poly_mat_is_kernel(Nt, nz, rdeg, Mt, ROW_LOWER); nmod_poly_mat_clear(N); nmod_poly_mat_clear(Nt); @@ -65,14 +65,10 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel_zls, state) for (i = 0; i < 16 * flint_test_multiplier(); i++) { - ulong nbits = 2 + n_randint(state, 4); - ulong rdim = 1 + n_randint(state, 4); - ulong cdim = rdim + 1 + n_randint(state, 2); - ulong deg = n_randint(state, 5); - /* ulong nbits = 2 + n_randint(state, 30); */ - /* ulong rdim = 1 + n_randint(state, 60); */ - /* ulong cdim = rdim + 1 + n_randint(state, 20); */ - /* ulong deg = n_randint(state, 20); */ + ulong nbits = 2 + n_randint(state, 63); + ulong rdim = 1 + n_randint(state, 60); + ulong cdim = rdim + 1 + n_randint(state, 20); + ulong deg = n_randint(state, 20); ulong prime = n_randprime(state, nbits, 1); diff --git a/flint-extras/src/nmod_poly_mat_extra/verification.c b/flint-extras/src/nmod_poly_mat_extra/verification.c index 1ce2fe37..997f4223 100644 --- a/flint-extras/src/nmod_poly_mat_extra/verification.c +++ b/flint-extras/src/nmod_poly_mat_extra/verification.c @@ -116,8 +116,9 @@ int nmod_poly_mat_is_approximant_basis(const nmod_poly_mat_t appbas, /* TODO currently specialized to ROW_LOWER (or at least ROW_stuff) */ /* -> checks whether ker is a shift-ordered weak Popov left kernel basis of pmat */ int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker, - const nmod_poly_mat_t pmat, + slong nz, const slong * shift, + const nmod_poly_mat_t pmat, orientation_t orient) { const slong rdim = pmat->r; @@ -137,6 +138,15 @@ int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker, success = 0; } + if (ker->r < nz) + { + printf("basis has wrong row dimension\n"); + success = 0; + } + + /* for the rest of tests, pretend ker->r is nz (DIRT: bypass const..) */ + ((nmod_poly_mat_struct *) ker)->r = nz; + /* check kernel is shifted reduced */ if (!nmod_poly_mat_is_ordered_weak_popov(ker, shift, orient)) { @@ -160,7 +170,7 @@ int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker, /* check rank */ slong rk = nmod_poly_mat_rank(pmat); - if (kdim != rdim - rk) + if (nz != rdim - rk) { printf("number of rows does not equal nullity\n"); success = 0; @@ -168,6 +178,9 @@ int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker, /* !! TODO check generation !! */ + /* revert row dimension of ker (DIRT: bypass const..) */ + ((nmod_poly_mat_struct *) ker)->r = kdim; + nmod_poly_mat_clear(residual); return success; diff --git a/flint-extras/src/nmod_poly_mat_kernel.h b/flint-extras/src/nmod_poly_mat_kernel.h index f2f0a083..caa3125f 100644 --- a/flint-extras/src/nmod_poly_mat_kernel.h +++ b/flint-extras/src/nmod_poly_mat_kernel.h @@ -13,7 +13,9 @@ #ifndef NMOD_POLY_MAT_KERNEL_H #define NMOD_POLY_MAT_KERNEL_H -#include "pml.h" +#include + +#include "nmod_poly_mat_forms.h" #ifdef __cplusplus extern "C" { @@ -44,25 +46,98 @@ extern "C" { * ------------ * Apart from the general interfaces (TODO) which offer several choices of * orientation, all other functions compute left kernel bases and use the - * following parameters: + * following parameters, where `nz` is the a priori unknown nullity of `pmat`: * - * \param[out] ker the output kernel basis (cannot alias `pmat`) - * \param[in] pmat the input polynomial matrix (no restriction) + * \param[out] ker the output kernel basis (cannot alias `pmat`, must be + * initialized with at least `nz` rows) + * \param[out] pivind the pivot index of `ker` (list of integers, length must + * be the number of rows of `pmat`, only the first `nz` entries matter in output) * \param[in,out] shift in: the input shift; and out: the output shifted row * degree of `ker` (list of integers, length must be the number of rows of * `pmat`) + * \param[in] pmat the input polynomial matrix (no restriction) + * + * \return slong, the nullity of `pmat` * * The computed `ker` is in shifted ordered weak Popov form, or in the canonical * shifted Popov form when the name of the function indicates so. + * + * FIXME what output dimensions for `ker`? + */ + + +/* general interface */ +/* `form` not implemented yet */ +slong nmod_poly_mat_kernel(nmod_poly_mat_t ker, + slong * pivind, + slong * shift, + const nmod_poly_mat_t pmat, + poly_mat_form_t form, + orientation_t orient); +// ROW_LOWER -> direct call +// COL_UPPER -> transpose input, call, transpose output +// ROW_UPPER -> mirror rows of input, call, mirror columns+rows of output +// COL_LOWER -> + +/* using approximant basis at sufficiently large order */ +/** Computes a `shift`-ordered weak Popov left kernel basis `ker` for `pmat`, + * recovered as a subset of the rows of a minimal approximant basis at + * sufficiently large order. Computes the + * `shift`-pivot index `pivind` of `kerbas`, and `shift` becomes the shifted + * row degree of `kerbas` (for the input shift). */ +slong nmod_poly_mat_kernel_via_approx(nmod_poly_mat_t ker, + slong * pivind, + slong * shift, + const nmod_poly_mat_t pmat); +/* FIXME not implemented yet: using interpolant basis at sufficiently many points */ +/* slong nmod_poly_mat_kernel_via_interp(nmod_poly_mat_t ker, */ +/* slong * shift, */ +/* const nmod_poly_mat_t pmat); */ +/* TODO */ +slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, + slong * pivind, + slong * shift, + const nmod_poly_mat_t pmat); +/* TODO */ +/* slong nmod_poly_mat_kernel_zls_interp(nmod_poly_mat_t ker, */ +/* slong * shift, */ +/* const nmod_poly_mat_t pmat); */ +/** Verifying if a matrix is a minimal kernel basis. + * + * This checks whether the matrix `ker` is a `shift`-minimal kernel basis for + * `pmat` for the required form `form`, with orientation described by `orient`. + * + * \param[in] ker kernel basis + * \param[in] nz nullity associated to ker + * \param[in] shift shift + * \param[in] pmat polynomial matrix + * \param[in] form [not implemented yet] required form for `kerbas` (see #poly_mat_form_t) + * \param[in] orient indicates the orientation (left/right kernel) and the definition of pivots + * \param[in] randomized [not implemented yet] if `true`, the algorithm may use a Monte Carlo or Las Vegas verification algorithm + * + * \return int, result of the verification + * + * \todo support all options, make doc more clear concerning Las Vegas / Monte Carlo + * \todo WARNING! for the moment, does not really check generation! + * \todo WARNING! for the moment, hardcoded to check for ordered weak Popov + */ +int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker, + slong nz, + /* const slong * pivind, */ /* TODO */ + const slong * shift, + const nmod_poly_mat_t pmat, + orientation_t orient); + + /** * From 29f1ab4edddde277ba29616639f7045bf31936b3 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sat, 3 Jan 2026 21:20:09 +0100 Subject: [PATCH 04/12] profile updated --- flint-extras/src/nmod_poly_mat_extra/kernel.c | 19 +++---- .../nmod_poly_mat_extra/profile/p-kernel.c | 49 +++++++++++++++++-- 2 files changed, 50 insertions(+), 18 deletions(-) diff --git a/flint-extras/src/nmod_poly_mat_extra/kernel.c b/flint-extras/src/nmod_poly_mat_extra/kernel.c index f285152f..2bc073c9 100644 --- a/flint-extras/src/nmod_poly_mat_extra/kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/kernel.c @@ -100,11 +100,9 @@ slong nmod_poly_mat_kernel_via_approx(nmod_poly_mat_t ker, const slong order = FLINT_MIN(m, n) * d + amp + d + 1; /* compute approximant basis and shifted row degree */ - nmod_poly_mat_t appbas; - nmod_poly_mat_init(appbas, m, m, pmat->modulus); - nmod_poly_mat_pmbasis(appbas, pivind, pmat, order); + nmod_poly_mat_pmbasis(ker, pivind, pmat, order); - /* find rows which belong to the kernel and record their pivot information */ + /* gather information for rows which belong to the kernel */ /* note: at this stage, pivind == rdeg_s(appbas) */ slong nz = 0; for (slong i = 0; i < m; i++) @@ -113,19 +111,14 @@ slong nmod_poly_mat_kernel_via_approx(nmod_poly_mat_t ker, { shift[nz] = pivind[i]; pivind[nz] = i; + for (slong j = 0; j < m; j++) + FLINT_SWAP(nmod_poly_struct, + *nmod_poly_mat_entry(ker, nz, j), + *nmod_poly_mat_entry(ker, i, j)); nz += 1; } } - /* swap the kernel rows */ - for (slong i = 0; i < nz; i++) - for (slong j = 0; j < m; j++) - FLINT_SWAP(nmod_poly_struct, - *nmod_poly_mat_entry(ker, i, j), - *nmod_poly_mat_entry(appbas, pivind[i], j)); - - nmod_poly_mat_clear(appbas); - return nz; } diff --git a/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c b/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c index 7587bfa8..3c52caea 100644 --- a/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c @@ -29,7 +29,7 @@ typedef struct } time_args; -#define TIME_KER(fun) \ +#define TIME_KER_ZLS(fun) \ void time_##fun(time_args targs, flint_rand_t state) \ { \ const slong rdim = targs.rdim; \ @@ -66,9 +66,47 @@ void time_##fun(time_args targs, flint_rand_t state) \ \ nmod_poly_mat_clear(F); \ nmod_poly_mat_clear(Ft); /* TMP */ \ + nmod_poly_mat_clear(K); \ +} + +#define TIME_KER(fun) \ +void time_##fun(time_args targs, flint_rand_t state) \ +{ \ + const slong rdim = targs.rdim; \ + const slong cdim = targs.cdim; \ + const slong deg = targs.deg; \ + /* const slong rank = targs.rank; */ /* TODO */ \ + /* const slong stype = targs.stype; */ /* TODO */ \ + const slong n = targs.modn; \ + \ + nmod_t mod; \ + nmod_init(&mod, n); \ + \ + nmod_poly_mat_t pmat; \ + nmod_poly_mat_init(pmat, rdim, cdim, n); \ + nmod_poly_mat_rand(pmat, state, deg); \ + \ + slong pivind[rdim]; \ + slong shift[rdim]; \ + for (slong i = 0; i < rdim; i++) \ + shift[i] = 0; \ + nmod_poly_mat_t ker; \ + nmod_poly_mat_init(ker, rdim, rdim, n); \ + \ + double FLINT_SET_BUT_UNUSED(tcpu), twall; \ + \ + TIMEIT_START; \ + nmod_poly_mat_##fun(ker, pivind, shift, pmat); \ + TIMEIT_STOP_VALUES(tcpu, twall); \ + \ + printf("%.2e", twall); \ + \ + nmod_poly_mat_clear(pmat); \ + nmod_poly_mat_clear(ker); \ } -TIME_KER(kernel_zls) +TIME_KER_ZLS(kernel_zls) +TIME_KER(kernel_via_approx) /*-------------------------*/ /* main */ @@ -97,10 +135,11 @@ int main(int argc, char ** argv) /* TODO shift type */ // bench functions - const slong nfuns = 1; + const slong nfuns = 2; typedef void (*timefun) (time_args, flint_rand_t); const timefun funs[] = { - time_kernel_zls, // 0 + time_kernel_via_approx, // 0 + time_kernel_zls, // 1 }; // TODO @@ -148,7 +187,7 @@ int main(int argc, char ** argv) { /* rdim; cdim; deg; rank; stype; modn; */ time_args targs = {8, 4, 1000, 4, 0, n_nextprime(UWORD(1) << 20, 0)}; - time_kernel_zls(targs, state); + time_kernel_via_approx(targs, state); printf(" "); } printf("\n\n"); From c3edea389d3af5898202c739980adaaa41c92fcc Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sun, 4 Jan 2026 10:55:17 +0100 Subject: [PATCH 05/12] zls-approx : splitting columns is ok --- flint-extras/src/nmod_poly_mat_extra/kernel.c | 346 ++++++++++++++++-- .../src/nmod_poly_mat_extra/test/main.c | 2 + .../test/t-kernel_via_approx.c | 3 +- .../test/t-kernel_zls_approx.c | 74 ++++ .../src/nmod_poly_mat_extra/verification.c | 20 +- 5 files changed, 411 insertions(+), 34 deletions(-) create mode 100644 flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls_approx.c diff --git a/flint-extras/src/nmod_poly_mat_extra/kernel.c b/flint-extras/src/nmod_poly_mat_extra/kernel.c index 2bc073c9..c703f58b 100644 --- a/flint-extras/src/nmod_poly_mat_extra/kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/kernel.c @@ -21,22 +21,6 @@ #include "nmod_poly_mat_extra.h" #include "nmod_poly_mat_kernel.h" -/* /1* requires len >= 1 *1/ */ -/* FLINT_INLINE slong _slong_vec_amplitude(slong * vec, slong len) */ -/* { */ -/* slong min = vec[0]; */ -/* slong max = vec[0]; */ -/* for (slong i = 1; i < len; i++) */ -/* { */ -/* slong v = vec[i]; */ -/* if (v < min) */ -/* min = v; */ -/* if (v > max) */ -/* max = v; */ -/* } */ -/* return max - min; */ -/* } */ - slong nmod_poly_mat_kernel_via_approx(nmod_poly_mat_t ker, slong * pivind, slong * shift, @@ -55,13 +39,12 @@ slong nmod_poly_mat_kernel_via_approx(nmod_poly_mat_t ker, return m; } + /* FIXME may be reasonable to also have a case for d == 0 */ + /* pmat full row rank => empty kernel */ /* full row rank implied by m == 1 or (m <= n && constant coefficient has full rank) */ if (m == 1) - { - nmod_poly_mat_zero(ker); return 0; - } if (m <= n) { @@ -75,10 +58,7 @@ slong nmod_poly_mat_kernel_via_approx(nmod_poly_mat_t ker, */ const slong r = nmod_mat_lu(pivind, coeff0, 1); if (r == m) - { - nmod_poly_mat_zero(ker); return 0; - } } /* compute amplitude of `shift`, make `pivind` a temporary copy of `shift` */ @@ -122,6 +102,328 @@ slong nmod_poly_mat_kernel_via_approx(nmod_poly_mat_t ker, return nz; } +/* TODO make pmat non const */ +/* TODO handle empty matrices? */ +/* NOTE requires shift >= rdeg(pmat) entrywise (?) */ +slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, + slong * pivind, + slong * shift, + const nmod_poly_mat_t pmat) +{ + const slong m = pmat->r; + const slong n = pmat->c; + + /* pmat == 0 => kernel is identity*/ + if (nmod_poly_mat_is_zero(pmat)) + { + nmod_poly_mat_one(ker); + for (slong i = 0; i < m; i++) + pivind[i] = i; + return m; + } + + /* FIXME may be reasonable to also have a case for constant matrices */ + + /* pmat full row rank => empty kernel */ + /* full row rank implied by m == 1 or (m <= n && constant coefficient has full rank) */ + if (m == 1) + return 0; + +// if (m <= n) +// { +// nmod_mat_t coeff0; +// nmod_mat_init(coeff0, m, n, pmat->modulus); +// nmod_poly_mat_get_coeff_mat(coeff0, pmat, 0); +// /* could use: +// * const slong r = nmod_mat_rank(coeff0); +// * but let's call more directly its core computation +// * (this uses pivind as temporary of length m for row permutation) +// */ +// const slong r = nmod_mat_lu(pivind, coeff0, 1); +// if (r == m) +// return 0; +// } + + /* if n > m/2 (with here m >= 2), straightforward recursion by splitting */ + /* the matrix into two submatrices with n/2 columns */ + if (2*n > m) + { + /* we split pmat == [ n - n/2 cols | n/2 cols] */ + /* (we want more columns in left part, for the residual computation) */ + slong nz; + + nmod_poly_mat_t submat; /* window only */ + nmod_poly_mat_t residual; /* window only */ + nmod_poly_mat_t ker1; /* window only */ + nmod_poly_mat_t ker2; /* actual init */ + nmod_poly_mat_t ker3; /* actual init */ + + /* first recursive call on left submatrix */ + nmod_poly_mat_window_init(submat, pmat, 0, 0, m, n - n/2); + /* nz = nmod_poly_mat_kernel_zls_approx(ker, pivind, shift, submat); */ + nz = nmod_poly_mat_kernel_zls_approx(ker, pivind, shift, submat); + nmod_poly_mat_window_clear(submat); + + /* if first kernel is empty, early exit */ + if (nz == 0) + return 0; + + /* actual kernel: nz first rows */ + nmod_poly_mat_window_init(ker1, ker, 0, 0, nz, m); + + /* residual: (first ker) * (right submatrix) */ + /* overwrite left columns of pmat */ + nmod_poly_mat_window_init(submat, pmat, 0, n - n/2, m, n); + nmod_poly_mat_window_init(residual, pmat, 0, 0, nz, n/2); + nmod_poly_mat_mul_classical(residual, ker1, submat); + nmod_poly_mat_window_clear(submat); + + /* recursive call 2, on residual */ + nmod_poly_mat_init(ker2, nz, nz, pmat->modulus); + slong * pivind2 = FLINT_ARRAY_ALLOC(nz, slong); + /* nz = nmod_poly_mat_kernel_zls_approx(ker2, pivind2, shift, residual); */ + nz = nmod_poly_mat_kernel_zls_approx(ker2, pivind2, shift, residual); + nmod_poly_mat_window_clear(residual); + + /* multiply bases and update pivind */ + nmod_poly_mat_init(ker3, nz, m, pmat->modulus); + nmod_poly_mat_window_init(submat, ker2, 0, 0, nz, ker2->c); + nmod_poly_mat_mul(ker3, submat, ker1); + nmod_poly_mat_window_clear(submat); + nmod_poly_mat_window_clear(ker1); + + for (slong i = 0; i < nz; i++) + { + for (slong j = 0; j < m; j++) + { + FLINT_SWAP(nmod_poly_struct, + *nmod_poly_mat_entry(ker, i, j), + *nmod_poly_mat_entry(ker3, i, j)); + } + } + + for (slong i = 0; i < nz; i++) + pivind[i] = pivind[pivind2[i]]; + + nmod_poly_mat_clear(ker2); + nmod_poly_mat_clear(ker3); + flint_free(pivind2); + + return nz; + } + + /* here, we are in the case m >= 2, n <= m/2, pmat nonzero */ + return nmod_poly_mat_kernel_via_approx(ker, pivind, shift, pmat); + + /* FIXME see if useful: */ + // add a constant to all the shift entries to ensure that the difference + // diff_shift = min(shift - rdeg(pmat)) is zero + // Recall that currently rdeg is rdeg(pmat) + // long diff_shift = shift[0] - rdeg[0]; + // for (long i = 1; i < m; ++i) + // { + // const long diff = shift[i]-rdeg[i]; + // if (diff_shift > diff) + // diff_shift = diff; + // } + // std::transform(shift.begin(), shift.end(), shift.begin(), + // [diff_shift](long x) -> long { return x - diff_shift;}); + +// // find parameter rho: sum of the n largest entries of (reduced) shift +// VecLong rdeg1(shift); // rdeg1 will be re-used later for another purpose, hence the name +// std::sort(rdeg1.begin(), rdeg1.end()); +// const long rho = std::accumulate(rdeg1.begin()+m-n, rdeg1.end(), 0); +// +// // order for call to approximation +// // choosing this order (with factor 2) is sufficient to ensure that the +// // approximant basis captures the whole kernel when n <= m/2 and pmat is +// // sufficiently generic +// const long order = 2 * ceil((double)rho / n) + 1; +// +// // compute approximant basis, along with shift-row degree and pivot degree +// // --> it does not necessarily capture the whole kernel; but does in two +// // notable cases: if n=1, or in the generic situation mentioned above +// +// // copy of input shift, will be needed at the end to compute shifted pivot index +// VecLong copy_shift(shift); +// Mat appbas; +// pmbasis(appbas, pmat, order, shift); // does not modify pmat +// +// // Identify submatrix of some rows of appbas which are in the kernel +// // Note the criterion: since before the call we have rdeg(pmat) <= shift, +// // the after the call we have rdeg(appbas*pmat) <= shift and therefore rows +// // with shift[i] < order are such that appbas[i] * pmat = 0. +// // Note that this may miss rows in the kernel, if there are some with shift >= order +// // which are in this appbas (this is usually not the case) +// rdeg1.clear(); +// VecLong rows_app_ker; +// VecLong rdeg2, rows_app_notker; +// for (long i=0; i this is often the case (e.g. in the generic situation mentioned +// // above), and testing this avoids going through a few multiplications and +// // the two recursive calls (which only lead to the conclusion that there is +// // no new kernel row to be found) +// long sum_matdeg = 0; +// for (long i = 0; i < m2; ++i) +// sum_matdeg += rdeg[rows_app_notker[i]]; +// long sum_kerdeg = 0; +// for (long i = 0; i < m1; ++i) +// sum_kerdeg += deg(appbas[rows_app_ker[i]][rows_app_ker[i]]); +// +// const bool early_exit = (sum_kerdeg == sum_matdeg); +// +// // if the whole kernel was captured according to the above test, or if +// // there was just one column, or if the kernel is full (matrix was zero), +// // then we have the whole kernel: just copy and return +// if (early_exit || n == 1 || m1 == m) +// { +// kerbas.SetDims(m1, m); +// for (long i = 0; i < m1; ++i) +// kerbas[i].swap(appbas[rows_app_ker[i]]); +// shift.swap(rdeg1); +// return; +// } +// +// // Note: for most input matrices/shifts, the following code will not be +// // executed because we will have taken early exit +// +// // retrieve the non-kernel part of the approximant +// Mat approx(INIT_SIZE, m2, m); +// for (long i = 0; i < m2; ++i) +// approx[i].swap(appbas[rows_app_notker[i]]); +// +// // compute residual: +// // pmat = trunc(trunc(approx, dA+1)*pmat div x^order, deg(pmat)+deg(approx)-order+1) +// middle_product(pmat,approx,pmat,order,deg(pmat)+deg(approx)-order); +// +// // pmat will be column-splitted into two submatrices of column dimension ~ n/2 +// const long n1 = n/2; +// const long n2 = n-n1; +// +// // recursive call 1, with left submatrix of the residual pmat +// Mat pmat_sub(INIT_SIZE, m2, n1); +// for (long i = 0; i < m2; ++i) +// for (long j = 0; j < n1; ++j) +// pmat_sub[i][j] = pmat[i][j]; +// +// Mat kerbas1; +// kernel_basis_zls_via_approximation(kerbas1, pmat_sub, rdeg2); +// +// // recursive call 2, with right submatrix of the residual pmat +// pmat_sub.SetDims(m2, n2); +// for (long i = 0; i < m2; ++i) +// for (long j = 0; j < n2; ++j) +// pmat_sub[i][j] = pmat[i][n1+j]; +// multiply(pmat, kerbas1, pmat_sub); +// pmat_sub.kill(); +// +// kernel_basis_zls_via_approximation(kerbas, pmat, rdeg2); +// +// // if kerbas is empty: (i.e. the approximant basis already captured the +// // whole kernel, although we had not guessed it with early_exit) +// // --> just copy and return +// if (kerbas.NumRows() == 0) +// { +// kerbas.SetDims(m1, m); +// for (long i = 0; i < m1; ++i) +// kerbas[i].swap(appbas[rows_app_ker[i]]); +// shift.swap(rdeg1); +// return; +// } +// +// // kerbas is non-empty: we have found new kernel rows +// // I/ complete the computation of these rows via +// // kerbas = kerbas * kerbas1 * approx +// // II/ merge this with rows from the kernel part of appbas +// +// // I/ complete the computation of new kernel rows +// // We use pmat as a temp, and store the result in approx +// // v1: +// //multiply(pmat, kerbas, kerbas1); +// //multiply(approx, pmat, approx); +// // v2: +// multiply(pmat, kerbas1, approx); +// multiply(approx, kerbas, pmat); +// // TODO when FFT / eval is used, this kind of product may be faster by +// // avoiding interpolation in the middle (?) +// +// // II/ merge with previously obtained kernel rows, ensuring that +// // the resulting basis is in shifted ordered weak Popov form +// +// // Note concerning shift: to ensure that shift contains the correct shifted +// // row degree of kerbas, we add the constant diff_shift that had been +// // removed from the input shift +// +// // Note that `rows_app_ker` is the shifted pivot index of the kernel part +// // of `appbas` +// // --> we need to compute the shifted pivot index of the newly found rows +// // of the kernel, as follows (note that copy_shift is a copy of the input +// // shift, while here `shift` is used as a temporary variable which will +// // store pivot degrees, which we discard): +// VecLong pivots_approx; +// row_pivots(pivots_approx, shift, approx, copy_shift); +// +// const long ker_dim = m1 + approx.NumRows(); +// kerbas.SetDims(ker_dim, m); +// shift.resize(ker_dim); +// long i=0, i_appbas=0, i_approx=0; +// while (i_appbas < m1 && i_approx < approx.NumRows()) +// { +// if (rows_app_ker[i_appbas] < pivots_approx[i_approx]) +// { +// // row already captured in preliminary appbas computation +// kerbas[i].swap(appbas[rows_app_ker[i_appbas]]); +// shift[i] = rdeg1[i_appbas]+diff_shift; +// ++i_appbas; ++i; +// } +// else +// { +// // row computed via recursive calls and stored in `approx` +// kerbas[i].swap(approx[i_approx]); +// shift[i] = rdeg2[i_approx]+diff_shift; +// ++i_approx; ++i; +// } +// } +// while (i_appbas < m1) +// { +// // row already captured in preliminary appbas computation +// kerbas[i].swap(appbas[rows_app_ker[i_appbas]]); +// shift[i] = rdeg1[i_appbas]+diff_shift; +// ++i_appbas; ++i; +// } +// while (i_approx < approx.NumRows()) +// { +// // row computed via recursive calls and stored in `approx` +// kerbas[i].swap(approx[i_approx]); +// shift[i] = rdeg2[i_approx]+diff_shift; +// ++i_approx; ++i; +// } +} + + + + /** * diff --git a/flint-extras/src/nmod_poly_mat_extra/test/main.c b/flint-extras/src/nmod_poly_mat_extra/test/main.c index ec0cec9d..e94da225 100644 --- a/flint-extras/src/nmod_poly_mat_extra/test/main.c +++ b/flint-extras/src/nmod_poly_mat_extra/test/main.c @@ -18,6 +18,7 @@ #include "t-dixon.c" #include "t-hermite_normal_form.c" #include "t-kernel_via_approx.c" +#include "t-kernel_zls_approx.c" #include "t-kernel_zls.c" #include "t-middle_product_geometric.c" #include "t-mul_geometric.c" @@ -35,6 +36,7 @@ test_struct tests[] = TEST_FUNCTION(nmod_poly_mat_dixon), TEST_FUNCTION(nmod_poly_mat_hnf), TEST_FUNCTION(nmod_poly_mat_kernel_via_approx), + TEST_FUNCTION(nmod_poly_mat_kernel_zls_approx), TEST_FUNCTION(nmod_poly_mat_kernel_zls), TEST_FUNCTION(nmod_poly_mat_mbasis), TEST_FUNCTION(nmod_poly_mat_middle_product_geometric), diff --git a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_via_approx.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_via_approx.c index d6ef8acb..480e0403 100644 --- a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_via_approx.c +++ b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_via_approx.c @@ -1,5 +1,5 @@ /* - Copyright (C) 2025 Gilles Villard + Copyright (C) 2025 Vincent Neiger This file is part of PML. @@ -55,6 +55,7 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel_via_approx, state) result = nmod_poly_mat_is_kernel(ker, nz, shift, pmat, ROW_LOWER); nmod_poly_mat_clear(pmat); + nmod_poly_mat_clear(ker); flint_free(shift); flint_free(rdeg); flint_free(pivind); diff --git a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls_approx.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls_approx.c new file mode 100644 index 00000000..75eb37aa --- /dev/null +++ b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls_approx.c @@ -0,0 +1,74 @@ +/* + Copyright (C) 2025 Vincent Neiger + + This file is part of PML. + + PML is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License version 2.0 (GPL-2.0-or-later) + as published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. See + . +*/ + +#include +#include +#include +#include +#include +#include + +#include "nmod_poly_mat_kernel.h" + +TEST_FUNCTION_START(nmod_poly_mat_kernel_zls_approx, state) +{ + int i, result; + + for (i = 0; i < 16 * flint_test_multiplier(); i++) + { + ulong nbits = 2 + n_randint(state, 63); + slong rdim = n_randint(state, 20); + slong cdim = n_randint(state, 20); + ulong len = 1 + n_randint(state, 40); + flint_printf("TEST iter %ld -- %ld, %ld, %ld\n", i, rdim, cdim, len); + + slong * shift = FLINT_ARRAY_ALLOC(rdim, slong); + for (slong i = 0; i < rdim; i++) + shift[i] = 0; + + slong * pivind = FLINT_ARRAY_ALLOC(rdim, slong); + slong * rdeg = FLINT_ARRAY_ALLOC(rdim, slong); + for (slong i = 0; i < rdim; i++) + rdeg[i] = shift[i]; + + ulong prime = n_randprime(state, nbits, 1); + + nmod_poly_mat_t pmat; + nmod_poly_mat_init(pmat, rdim, cdim, prime); + nmod_poly_mat_randtest(pmat, state, len); + nmod_poly_mat_t copy_pmat; + nmod_poly_mat_init_set(copy_pmat, pmat); + + nmod_poly_mat_t ker; + nmod_poly_mat_init(ker, rdim, rdim, prime); + /* TODO test weak Popov + pivind */ + /* TODO currently does not test generation */ + slong nz = nmod_poly_mat_kernel_zls_approx(ker, pivind, rdeg, pmat); + flint_printf("kernel computed, now testing...\n"); + result = nmod_poly_mat_is_kernel(ker, nz, shift, copy_pmat, ROW_LOWER); + + nmod_poly_mat_clear(ker); + nmod_poly_mat_clear(pmat); + nmod_poly_mat_clear(copy_pmat); + flint_free(shift); + flint_free(rdeg); + flint_free(pivind); + + if (!result) + { + TEST_FUNCTION_FAIL("rdim = %wu, cdim = %wu, length = %wu, p = %wu\n", \ + rdim, cdim, len, prime); + } + } + + TEST_FUNCTION_END(state); +} diff --git a/flint-extras/src/nmod_poly_mat_extra/verification.c b/flint-extras/src/nmod_poly_mat_extra/verification.c index 997f4223..9ba51042 100644 --- a/flint-extras/src/nmod_poly_mat_extra/verification.c +++ b/flint-extras/src/nmod_poly_mat_extra/verification.c @@ -11,6 +11,7 @@ */ #include "nmod_poly_mat_forms.h" +#include "nmod_poly_mat_io.h" #include #include #include @@ -123,11 +124,13 @@ int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker, { const slong rdim = pmat->r; const slong cdim = pmat->c; - const slong kdim = ker->r; const ulong prime = pmat->modulus; nmod_poly_mat_t residual; - nmod_poly_mat_init(residual, kdim, cdim, prime); + nmod_poly_mat_init(residual, nz, cdim, prime); + + nmod_poly_mat_t kernz; + nmod_poly_mat_window_init(kernz, ker, 0, 0, nz, rdim); int success = 1; @@ -144,15 +147,12 @@ int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker, success = 0; } - /* for the rest of tests, pretend ker->r is nz (DIRT: bypass const..) */ - ((nmod_poly_mat_struct *) ker)->r = nz; - /* check kernel is shifted reduced */ - if (!nmod_poly_mat_is_ordered_weak_popov(ker, shift, orient)) + if (!nmod_poly_mat_is_ordered_weak_popov(kernz, shift, orient)) { /* printf("basis is not shifted-weak Popov\n"); */ /* TODO temporarily allow reduced but not s-weak Popov */ - if (!nmod_poly_mat_is_reduced(ker, shift, orient)) + if (!nmod_poly_mat_is_reduced(kernz, shift, orient)) { printf("basis is not shifted-reduced\n"); success = 0; @@ -161,7 +161,8 @@ int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker, /* compute residual, check rows of ker are in the kernel */ /* TODO enhancement: offer randomized option, using Freivalds-like randomized strategy (see ntl-extras) */ - nmod_poly_mat_mul(residual, ker, pmat); + nmod_poly_mat_mul(residual, kernz, pmat); + /* nmod_poly_mat_degree_matrix_print_pretty(residual); */ if (!nmod_poly_mat_is_zero(residual)) { printf("not all rows are in the kernel\n"); @@ -178,9 +179,6 @@ int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker, /* !! TODO check generation !! */ - /* revert row dimension of ker (DIRT: bypass const..) */ - ((nmod_poly_mat_struct *) ker)->r = kdim; - nmod_poly_mat_clear(residual); return success; From a01c0bd5e307c51a04f047640d9341cde96b81d8 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sun, 4 Jan 2026 22:26:17 +0100 Subject: [PATCH 06/12] ZLS works, missing one early exit and some cleaning --- flint-extras/src/nmod_poly_mat_extra/kernel.c | 402 ++++++++---------- .../nmod_poly_mat_extra/profile/p-kernel.c | 94 ++-- .../src/nmod_poly_mat_extra/test/main.c | 2 +- .../test/t-kernel_zls_approx.c | 6 +- flint-extras/src/nmod_poly_mat_utils.h | 3 - 5 files changed, 240 insertions(+), 267 deletions(-) diff --git a/flint-extras/src/nmod_poly_mat_extra/kernel.c b/flint-extras/src/nmod_poly_mat_extra/kernel.c index c703f58b..d18cfe0b 100644 --- a/flint-extras/src/nmod_poly_mat_extra/kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/kernel.c @@ -19,7 +19,17 @@ #include #include "nmod_poly_mat_extra.h" +#include "nmod_poly_mat_io.h" #include "nmod_poly_mat_kernel.h" +#include "nmod_poly_mat_multiply.h" + +/* general interface: + * A/ compute rdeg and cdeg + * ->suppress zero columns + * ->trivialize zero rows + * B/ handle case of constant matrices + * C/ go for zls + * */ slong nmod_poly_mat_kernel_via_approx(nmod_poly_mat_t ker, slong * pivind, @@ -51,12 +61,11 @@ slong nmod_poly_mat_kernel_via_approx(nmod_poly_mat_t ker, nmod_mat_t coeff0; nmod_mat_init(coeff0, m, n, pmat->modulus); nmod_poly_mat_get_coeff_mat(coeff0, pmat, 0); - /* could use: - * const slong r = nmod_mat_rank(coeff0); - * but let's call more directly its core computation - * (this uses pivind as temporary of length m for row permutation) - */ + /* could use: const slong r = nmod_mat_rank(coeff0); */ + /* but let's call more directly its core computation */ + /* (this uses pivind as temporary of length m for row permutation) */ const slong r = nmod_mat_lu(pivind, coeff0, 1); + nmod_mat_clear(coeff0); if (r == m) return 0; } @@ -103,46 +112,37 @@ slong nmod_poly_mat_kernel_via_approx(nmod_poly_mat_t ker, } /* TODO make pmat non const */ -/* TODO handle empty matrices? */ -/* NOTE requires shift >= rdeg(pmat) entrywise (?) */ +/* Follows the description of Zhou-Labahn-Storjohann algorithm in [LNVZ22, Algo.1] */ +/* [LNVZ22] Labahn-Neiger-Vu-Zhou, Proceedings ISSAC 2022, arxiv.org/pdf/2202.09329 */ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, slong * pivind, slong * shift, const nmod_poly_mat_t pmat) { + /* flint_printf("call %ld, %ld\n", pmat->r, pmat->c); */ + /* nmod_poly_mat_degree_matrix_print_pretty(pmat); */ const slong m = pmat->r; const slong n = pmat->c; - /* pmat == 0 => kernel is identity*/ - if (nmod_poly_mat_is_zero(pmat)) + /* empty matrices, or single row/column vectors that are zero */ + if (m == 0 || n == 0 + || ((m == 1 || n == 1) && (nmod_poly_mat_max_length(pmat) == 0))) { + /* flint_printf("base case empty||zero!\n"); */ nmod_poly_mat_one(ker); for (slong i = 0; i < m; i++) pivind[i] = i; return m; } - /* FIXME may be reasonable to also have a case for constant matrices */ - - /* pmat full row rank => empty kernel */ - /* full row rank implied by m == 1 or (m <= n && constant coefficient has full rank) */ + /* one (nonzero) row vector */ if (m == 1) return 0; -// if (m <= n) -// { -// nmod_mat_t coeff0; -// nmod_mat_init(coeff0, m, n, pmat->modulus); -// nmod_poly_mat_get_coeff_mat(coeff0, pmat, 0); -// /* could use: -// * const slong r = nmod_mat_rank(coeff0); -// * but let's call more directly its core computation -// * (this uses pivind as temporary of length m for row permutation) -// */ -// const slong r = nmod_mat_lu(pivind, coeff0, 1); -// if (r == m) -// return 0; -// } + /* one (nonzero) column vector */ + /* TODO */ + /* if (n == 1) */ + /* return 0; */ /* if n > m/2 (with here m >= 2), straightforward recursion by splitting */ /* the matrix into two submatrices with n/2 columns */ @@ -160,7 +160,6 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, /* first recursive call on left submatrix */ nmod_poly_mat_window_init(submat, pmat, 0, 0, m, n - n/2); - /* nz = nmod_poly_mat_kernel_zls_approx(ker, pivind, shift, submat); */ nz = nmod_poly_mat_kernel_zls_approx(ker, pivind, shift, submat); nmod_poly_mat_window_clear(submat); @@ -181,7 +180,6 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, /* recursive call 2, on residual */ nmod_poly_mat_init(ker2, nz, nz, pmat->modulus); slong * pivind2 = FLINT_ARRAY_ALLOC(nz, slong); - /* nz = nmod_poly_mat_kernel_zls_approx(ker2, pivind2, shift, residual); */ nz = nmod_poly_mat_kernel_zls_approx(ker2, pivind2, shift, residual); nmod_poly_mat_window_clear(residual); @@ -193,14 +191,10 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, nmod_poly_mat_window_clear(ker1); for (slong i = 0; i < nz; i++) - { for (slong j = 0; j < m; j++) - { FLINT_SWAP(nmod_poly_struct, *nmod_poly_mat_entry(ker, i, j), *nmod_poly_mat_entry(ker3, i, j)); - } - } for (slong i = 0; i < nz; i++) pivind[i] = pivind[pivind2[i]]; @@ -212,71 +206,70 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, return nz; } - /* here, we are in the case m >= 2, n <= m/2, pmat nonzero */ - return nmod_poly_mat_kernel_via_approx(ker, pivind, shift, pmat); - - /* FIXME see if useful: */ - // add a constant to all the shift entries to ensure that the difference - // diff_shift = min(shift - rdeg(pmat)) is zero - // Recall that currently rdeg is rdeg(pmat) - // long diff_shift = shift[0] - rdeg[0]; - // for (long i = 1; i < m; ++i) - // { - // const long diff = shift[i]-rdeg[i]; - // if (diff_shift > diff) - // diff_shift = diff; - // } - // std::transform(shift.begin(), shift.end(), shift.begin(), - // [diff_shift](long x) -> long { return x - diff_shift;}); - -// // find parameter rho: sum of the n largest entries of (reduced) shift -// VecLong rdeg1(shift); // rdeg1 will be re-used later for another purpose, hence the name -// std::sort(rdeg1.begin(), rdeg1.end()); -// const long rho = std::accumulate(rdeg1.begin()+m-n, rdeg1.end(), 0); -// -// // order for call to approximation -// // choosing this order (with factor 2) is sufficient to ensure that the -// // approximant basis captures the whole kernel when n <= m/2 and pmat is -// // sufficiently generic -// const long order = 2 * ceil((double)rho / n) + 1; -// -// // compute approximant basis, along with shift-row degree and pivot degree -// // --> it does not necessarily capture the whole kernel; but does in two -// // notable cases: if n=1, or in the generic situation mentioned above -// -// // copy of input shift, will be needed at the end to compute shifted pivot index -// VecLong copy_shift(shift); -// Mat appbas; -// pmbasis(appbas, pmat, order, shift); // does not modify pmat -// -// // Identify submatrix of some rows of appbas which are in the kernel -// // Note the criterion: since before the call we have rdeg(pmat) <= shift, -// // the after the call we have rdeg(appbas*pmat) <= shift and therefore rows -// // with shift[i] < order are such that appbas[i] * pmat = 0. -// // Note that this may miss rows in the kernel, if there are some with shift >= order -// // which are in this appbas (this is usually not the case) -// rdeg1.clear(); -// VecLong rows_app_ker; -// VecLong rdeg2, rows_app_notker; -// for (long i=0; i maxdeg) + maxdeg = d; + if (diff_shift > diff) + diff_shift = diff; + } + + /* early exit: pmat == 0 => kernel is identity */ + if (maxdeg == -1) + { + /* flint_printf("second base case zero!\n"); */ + nmod_poly_mat_one(ker); + for (slong i = 0; i < m; i++) + pivind[i] = i; + flint_free(buf); + return m; + } + + /* build order for approximation (this choice diverges from [LNVZ22, Algo.1]): */ + /* order = 1 + max(maxdeg, 1 + floor((sum(shift) - 1) / (m - n))) */ + /* -> approx basis "generically" captures the whole kernel (see end of file) */ + /* FIXME could be improved in case of unbalanced cdeg(pmat), */ + /* based on [LNVZ22, Theorem 3.3, Proof of Item (iii)] */ + slong order = - m * diff_shift - 1; + for (slong i = 0; i < m; i++) + order += shift[i]; + order = 1 + order / (m - n); + order = 1 + FLINT_MAX(maxdeg, order); + + nmod_poly_mat_pmbasis(ker, shift, pmat, order); + + /* run an easy degree-based detection of rows in the kernel */ + /* permute ker into [kernel rows \\ other rows], preserving increasing pivots */ + slong nz = 0; + for (slong i = 0; i < m; i++) + { + if (shift[i] < order + diff_shift) + { + pivind[nz] = i; + nz += 1; + } + else + buf[i - nz] = i; + } + for (slong i = nz; i < m; i++) + pivind[i] = buf[i - nz]; + + nmod_poly_mat_permute_rows(ker, pivind, shift); + /* flint_printf("pmbasis ok, found %ld rows\n", nz); */ + /* flint_printf("pivind: %{slong*}\n", pivind, m); */ + /* nmod_poly_mat_degree_matrix_print_pretty(ker); */ + + /* early exit */ + /* TODO */ // // if the sum of pivot degrees in the part of kernel basis computed is // // equal to one of the upper bounds (here, sum of row degrees of pmat for // // non-pivot rows), then we know that we have captured the whole kernel @@ -293,10 +286,9 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, // // const bool early_exit = (sum_kerdeg == sum_matdeg); // -// // if the whole kernel was captured according to the above test, or if -// // there was just one column, or if the kernel is full (matrix was zero), +// // if the whole kernel was captured according to the above test, or if the kernel is full (matrix was zero), // // then we have the whole kernel: just copy and return -// if (early_exit || n == 1 || m1 == m) +// if (early_exit) // { // kerbas.SetDims(m1, m); // for (long i = 0; i < m1; ++i) @@ -304,121 +296,99 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, // shift.swap(rdeg1); // return; // } -// -// // Note: for most input matrices/shifts, the following code will not be -// // executed because we will have taken early exit -// -// // retrieve the non-kernel part of the approximant -// Mat approx(INIT_SIZE, m2, m); -// for (long i = 0; i < m2; ++i) -// approx[i].swap(appbas[rows_app_notker[i]]); -// -// // compute residual: -// // pmat = trunc(trunc(approx, dA+1)*pmat div x^order, deg(pmat)+deg(approx)-order+1) -// middle_product(pmat,approx,pmat,order,deg(pmat)+deg(approx)-order); -// -// // pmat will be column-splitted into two submatrices of column dimension ~ n/2 -// const long n1 = n/2; -// const long n2 = n-n1; -// -// // recursive call 1, with left submatrix of the residual pmat -// Mat pmat_sub(INIT_SIZE, m2, n1); -// for (long i = 0; i < m2; ++i) -// for (long j = 0; j < n1; ++j) -// pmat_sub[i][j] = pmat[i][j]; -// -// Mat kerbas1; -// kernel_basis_zls_via_approximation(kerbas1, pmat_sub, rdeg2); -// -// // recursive call 2, with right submatrix of the residual pmat -// pmat_sub.SetDims(m2, n2); -// for (long i = 0; i < m2; ++i) -// for (long j = 0; j < n2; ++j) -// pmat_sub[i][j] = pmat[i][n1+j]; -// multiply(pmat, kerbas1, pmat_sub); -// pmat_sub.kill(); -// -// kernel_basis_zls_via_approximation(kerbas, pmat, rdeg2); -// -// // if kerbas is empty: (i.e. the approximant basis already captured the -// // whole kernel, although we had not guessed it with early_exit) -// // --> just copy and return -// if (kerbas.NumRows() == 0) -// { -// kerbas.SetDims(m1, m); -// for (long i = 0; i < m1; ++i) -// kerbas[i].swap(appbas[rows_app_ker[i]]); -// shift.swap(rdeg1); -// return; -// } -// -// // kerbas is non-empty: we have found new kernel rows -// // I/ complete the computation of these rows via -// // kerbas = kerbas * kerbas1 * approx -// // II/ merge this with rows from the kernel part of appbas -// -// // I/ complete the computation of new kernel rows -// // We use pmat as a temp, and store the result in approx -// // v1: -// //multiply(pmat, kerbas, kerbas1); -// //multiply(approx, pmat, approx); -// // v2: -// multiply(pmat, kerbas1, approx); -// multiply(approx, kerbas, pmat); -// // TODO when FFT / eval is used, this kind of product may be faster by -// // avoiding interpolation in the middle (?) -// -// // II/ merge with previously obtained kernel rows, ensuring that -// // the resulting basis is in shifted ordered weak Popov form -// -// // Note concerning shift: to ensure that shift contains the correct shifted -// // row degree of kerbas, we add the constant diff_shift that had been -// // removed from the input shift -// -// // Note that `rows_app_ker` is the shifted pivot index of the kernel part -// // of `appbas` -// // --> we need to compute the shifted pivot index of the newly found rows -// // of the kernel, as follows (note that copy_shift is a copy of the input -// // shift, while here `shift` is used as a temporary variable which will -// // store pivot degrees, which we discard): -// VecLong pivots_approx; -// row_pivots(pivots_approx, shift, approx, copy_shift); -// -// const long ker_dim = m1 + approx.NumRows(); -// kerbas.SetDims(ker_dim, m); -// shift.resize(ker_dim); -// long i=0, i_appbas=0, i_approx=0; -// while (i_appbas < m1 && i_approx < approx.NumRows()) -// { -// if (rows_app_ker[i_appbas] < pivots_approx[i_approx]) -// { -// // row already captured in preliminary appbas computation -// kerbas[i].swap(appbas[rows_app_ker[i_appbas]]); -// shift[i] = rdeg1[i_appbas]+diff_shift; -// ++i_appbas; ++i; -// } -// else -// { -// // row computed via recursive calls and stored in `approx` -// kerbas[i].swap(approx[i_approx]); -// shift[i] = rdeg2[i_approx]+diff_shift; -// ++i_approx; ++i; -// } -// } -// while (i_appbas < m1) -// { -// // row already captured in preliminary appbas computation -// kerbas[i].swap(appbas[rows_app_ker[i_appbas]]); -// shift[i] = rdeg1[i_appbas]+diff_shift; -// ++i_appbas; ++i; -// } -// while (i_approx < approx.NumRows()) -// { -// // row computed via recursive calls and stored in `approx` -// kerbas[i].swap(approx[i_approx]); -// shift[i] = rdeg2[i_approx]+diff_shift; -// ++i_approx; ++i; -// } + + /* proceed with remaining rows */ + nmod_poly_mat_t residual; /* actual init */ + nmod_poly_mat_t ker2; /* actual init */ + nmod_poly_mat_t prod; /* actual init */ + nmod_poly_mat_t approx; /* window only */ + nmod_poly_mat_t ker2nz; /* window only */ + + nmod_poly_mat_init(residual, m - nz, n, pmat->modulus); + nmod_poly_mat_window_init(approx, ker, nz, 0, m, m); + nmod_poly_mat_middle_product_naive(residual, approx, pmat, order, order + maxdeg + 1); + /* FIXME provide (and use) general middle_product interface */ + + /* TODO handle zero rows in G? */ + + /* kernel of residual */ + nmod_poly_mat_init(ker2, m - nz, m - nz, pmat->modulus); + const slong nz2 = nmod_poly_mat_kernel_zls_approx(ker2, buf, shift + nz, residual); + /* flint_printf("second kernel ok, found %ld rows\n", nz2); */ + /* flint_printf("pivind: %{slong*}\n", buf, m - nz); */ + + /* if nz2 == 0, whole kernel already computed, just return */ + if (nz2 == 0) + { + nmod_poly_mat_clear(residual); + nmod_poly_mat_window_clear(approx); + nmod_poly_mat_clear(ker2); + flint_free(buf); + return nz; + } + + /* multiply to get missing part of kernel */ + /* flint_printf("nz2 == %ld\n", nz2); */ + /* nmod_poly_mat_degree_matrix_print_pretty(ker2); */ + /* nmod_poly_mat_degree_matrix_print_pretty(approx); */ + nmod_poly_mat_window_init(ker2nz, ker2, 0, 0, nz2, ker2->c); + nmod_poly_mat_init(prod, nz2, m, pmat->modulus); + nmod_poly_mat_mul(prod, ker2nz, approx); + /* flint_printf("before swap::\n"); */ + /* nmod_poly_mat_degree_matrix_print_pretty(prod); */ + /* nmod_poly_mat_degree_matrix_print_pretty(ker); */ + for (slong i = 0; i < nz2; i++) + for (slong j = 0; j < m; j++) + FLINT_SWAP(nmod_poly_struct, + *nmod_poly_mat_entry(prod, i, j), + *nmod_poly_mat_entry(approx, i, j)); + /* flint_printf("after swap::\n"); */ + /* nmod_poly_mat_degree_matrix_print_pretty(ker); */ + + /* update pivind */ + for (slong i = 0; i < nz2; i++) + pivind[nz + i] = pivind[nz + buf[i]]; + + /* permute rows to get increasing pivot indices */ + nmod_poly_mat_window_clear(approx); + nmod_poly_mat_window_init(approx, ker, 0, 0, nz + nz2, m); + slong i1 = 0; + slong i2 = 0; + for (slong i = 0; i < nz + nz2; i++) + { + if (i2 == nz2 || (i1 < nz && pivind[i1] < pivind[nz + i2])) + { + buf[i] = i1; + i1++; + } + else + { + buf[i] = nz + i2; + i2++; + } + } + /* flint_printf("permutation buf: %{slong*}\n", buf, nz + nz2); */ + nmod_poly_mat_permute_rows(approx, buf, NULL); + slong * tmp = FLINT_ARRAY_ALLOC(nz + nz2, slong); /* use TMP_ALLOC */ + for (slong i = 0; i < nz + nz2; i++) + tmp[i] = pivind[i]; + for (slong i = 0; i < nz + nz2; i++) + pivind[i] = tmp[buf[i]]; + for (slong i = 0; i < nz + nz2; i++) + tmp[i] = shift[i]; + for (slong i = 0; i < nz + nz2; i++) + shift[i] = tmp[buf[i]]; + flint_free(tmp); + + nmod_poly_mat_window_clear(ker2nz); + nmod_poly_mat_window_clear(approx); + nmod_poly_mat_clear(residual); + nmod_poly_mat_clear(ker2); + nmod_poly_mat_clear(prod); + flint_free(buf); + + /* nmod_poly_mat_degree_matrix_print_pretty(ker); */ + return nz + nz2; } diff --git a/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c b/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c index 3c52caea..a2c5c90b 100644 --- a/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c @@ -52,7 +52,7 @@ void time_##fun(time_args targs, flint_rand_t state) \ nmod_poly_mat_transpose(Ft, F); \ /* END TMP */ \ \ - slong degN[rdim]; \ + slong * degN = FLINT_ARRAY_ALLOC(rdim, slong); \ nmod_poly_mat_t K; \ nmod_poly_mat_init(K, rdim, rdim, n); \ \ @@ -62,8 +62,9 @@ void time_##fun(time_args targs, flint_rand_t state) \ nmod_poly_mat_##fun(K, degN, Ft, NULL, 3.); \ TIMEIT_STOP_VALUES(tcpu, twall); \ \ - printf("%.2e", twall); \ + flint_printf("%.2e", twall); \ \ + flint_free(degN); \ nmod_poly_mat_clear(F); \ nmod_poly_mat_clear(Ft); /* TMP */ \ nmod_poly_mat_clear(K); \ @@ -86,8 +87,8 @@ void time_##fun(time_args targs, flint_rand_t state) \ nmod_poly_mat_init(pmat, rdim, cdim, n); \ nmod_poly_mat_rand(pmat, state, deg); \ \ - slong pivind[rdim]; \ - slong shift[rdim]; \ + slong * pivind = FLINT_ARRAY_ALLOC(rdim, slong); \ + slong * shift = FLINT_ARRAY_ALLOC(rdim, slong); \ for (slong i = 0; i < rdim; i++) \ shift[i] = 0; \ nmod_poly_mat_t ker; \ @@ -96,17 +97,24 @@ void time_##fun(time_args targs, flint_rand_t state) \ double FLINT_SET_BUT_UNUSED(tcpu), twall; \ \ TIMEIT_START; \ - nmod_poly_mat_##fun(ker, pivind, shift, pmat); \ + nmod_poly_mat_t copy_pmat; \ + nmod_poly_mat_init_set(copy_pmat, pmat); \ + for (slong i = 0; i < rdim; i++) \ + shift[i] = 0; \ + nmod_poly_mat_##fun(ker, pivind, shift, copy_pmat); \ TIMEIT_STOP_VALUES(tcpu, twall); \ \ - printf("%.2e", twall); \ + flint_printf("%.2e", twall); \ \ + flint_free(pivind); \ + flint_free(shift); \ nmod_poly_mat_clear(pmat); \ nmod_poly_mat_clear(ker); \ } TIME_KER_ZLS(kernel_zls) TIME_KER(kernel_via_approx) +TIME_KER(kernel_zls_approx) /*-------------------------*/ /* main */ @@ -135,11 +143,12 @@ int main(int argc, char ** argv) /* TODO shift type */ // bench functions - const slong nfuns = 2; + const slong nfuns = 3; typedef void (*timefun) (time_args, flint_rand_t); const timefun funs[] = { time_kernel_via_approx, // 0 - time_kernel_zls, // 1 + time_kernel_zls_approx, // 1 + time_kernel_zls, // 2 }; // TODO @@ -159,22 +168,23 @@ int main(int argc, char ** argv) //#endif const char * description[] = { - "#0 --> kernel (general interface) ", - "#1 --> .... TODO ", + "#0 --> via approx ", + "#1 --> ZLS via approx ", + "#2 --> .... TODO ", }; if (argc < 4 || argc > 6) // show usage { - printf("Usage: `%s [nbits] [rdim] [cdim] [deg] [rank] [shift] [fun]`\n", argv[0]); + printf("Usage: `%s [fun] [nbits] [rdim] [cdim] [deg] [rank] [shift]`\n", argv[0]); printf(" No argument shows this help.\n"); - printf(" First 3 arguments are mandatory.\n"); + printf(" First 4 arguments are mandatory.\n"); printf(" [rank] and [shift] not supported yet.\n"); + printf(" - fun: id number of the timed function (see below, -1 times all),\n"); printf(" - nbits: number of bits in [2..64] for the modulus, chosen as nextprime(2**(nbits-1))\n"); printf(" - rdim, cdim: input matrix is rdim x cdim\n"); printf(" - deg: matrix is random of degree < deg (default: predefined list)\n"); printf(" - rank: [unsupported] matrix is random of this rank\n"); printf(" - shift: [unsupported] type of shift\n"); - printf(" - fun: id number of the timed function (see below),\n"); printf("\nAvailable functions:\n"); for (slong j = 0; j < nfuns; j++) printf(" %s\n", description[j]); @@ -192,45 +202,41 @@ int main(int argc, char ** argv) } printf("\n\n"); - if (argc == 4) // nbits + rdim + cdim given - { - const slong rdim = atoi(argv[2]); - const slong cdim = atoi(argv[3]); - const slong b = atoi(argv[1]); - printf("bits fun rdim cdim deg\n"); - ulong n = n_nextprime(UWORD(1) << (b-1), 0); - for (slong ifun = 0; ifun < nfuns; ifun++) - { - const timefun tfun = funs[ifun]; - for (slong d = 0; d < ndegs; d++) - { - printf("%-5ld#%-3ld%-5ld%-5ld%-8ld", b, ifun, rdim, cdim, degs[d]); - time_args targs = {rdim, cdim, degs[d], FLINT_MIN(rdim, cdim), 0, n}; - tfun(targs, state); - printf(" "); - printf("\n"); - } - } - } - else if (argc == 5) // nbits + rdim + cdim + deg given + if (argc == 5) // fun + nbits + rdim + cdim given { - const slong rdim = atoi(argv[2]); - const slong cdim = atoi(argv[3]); - const slong deg = atoi(argv[4]); - const slong b = atoi(argv[1]); + const slong ifun = atoi(argv[1]); + const slong bits = atoi(argv[2]); + const slong rdim = atoi(argv[3]); + const slong cdim = atoi(argv[4]); + const timefun tfun = funs[ifun]; + const ulong n = n_nextprime(UWORD(1) << (bits-1), 0); printf("bits fun rdim cdim deg\n"); - ulong n = n_nextprime(UWORD(1) << (b-1), 0); - for (slong ifun = 0; ifun < nfuns; ifun++) + for (slong d = 0; d < ndegs; d++) { - const timefun tfun = funs[ifun]; - printf("%-5ld#%-3ld%-5ld%-5ld%-8ld", b, ifun, rdim, cdim, deg); - /* rdim; cdim; deg; rank; stype; modn; */ - time_args targs = {rdim, cdim, deg, FLINT_MIN(rdim, cdim), 0, n}; + printf("%-5ld#%-3ld%-5ld%-5ld%-8ld", bits, ifun, rdim, cdim, degs[d]); + time_args targs = {rdim, cdim, degs[d], FLINT_MIN(rdim, cdim), 0, n}; tfun(targs, state); printf(" "); printf("\n"); } } + else if (argc == 6) // fun + nbits + rdim + cdim + deg given + { + const slong ifun = atoi(argv[1]); + const slong bits = atoi(argv[2]); + const slong rdim = atoi(argv[3]); + const slong cdim = atoi(argv[4]); + const slong deg = atoi(argv[5]); + const timefun tfun = funs[ifun]; + const ulong n = n_nextprime(UWORD(1) << (bits-1), 0); + printf("bits fun rdim cdim deg\n"); + printf("%-5ld#%-3ld%-5ld%-5ld%-8ld", bits, ifun, rdim, cdim, deg); + /* rdim; cdim; deg; rank; stype; modn; */ + time_args targs = {rdim, cdim, deg, FLINT_MIN(rdim, cdim), 0, n}; + tfun(targs, state); + printf(" "); + printf("\n"); + } flint_rand_clear(state); return 0; diff --git a/flint-extras/src/nmod_poly_mat_extra/test/main.c b/flint-extras/src/nmod_poly_mat_extra/test/main.c index e94da225..f0c10bd6 100644 --- a/flint-extras/src/nmod_poly_mat_extra/test/main.c +++ b/flint-extras/src/nmod_poly_mat_extra/test/main.c @@ -32,11 +32,11 @@ test_struct tests[] = { + TEST_FUNCTION(nmod_poly_mat_kernel_zls_approx), TEST_FUNCTION(nmod_poly_mat_det), TEST_FUNCTION(nmod_poly_mat_dixon), TEST_FUNCTION(nmod_poly_mat_hnf), TEST_FUNCTION(nmod_poly_mat_kernel_via_approx), - TEST_FUNCTION(nmod_poly_mat_kernel_zls_approx), TEST_FUNCTION(nmod_poly_mat_kernel_zls), TEST_FUNCTION(nmod_poly_mat_mbasis), TEST_FUNCTION(nmod_poly_mat_middle_product_geometric), diff --git a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls_approx.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls_approx.c index 75eb37aa..5fa8f39c 100644 --- a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls_approx.c +++ b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls_approx.c @@ -26,9 +26,9 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel_zls_approx, state) for (i = 0; i < 16 * flint_test_multiplier(); i++) { ulong nbits = 2 + n_randint(state, 63); - slong rdim = n_randint(state, 20); - slong cdim = n_randint(state, 20); - ulong len = 1 + n_randint(state, 40); + slong rdim = n_randint(state, 30); + slong cdim = n_randint(state, 30); + ulong len = 1 + n_randint(state, 100); flint_printf("TEST iter %ld -- %ld, %ld, %ld\n", i, rdim, cdim, len); slong * shift = FLINT_ARRAY_ALLOC(rdim, slong); diff --git a/flint-extras/src/nmod_poly_mat_utils.h b/flint-extras/src/nmod_poly_mat_utils.h index ad96d2b0..8e7156ef 100644 --- a/flint-extras/src/nmod_poly_mat_utils.h +++ b/flint-extras/src/nmod_poly_mat_utils.h @@ -639,6 +639,3 @@ nmod_poly_mat_set_from_mat_poly(nmod_poly_mat_t pmat, #endif #endif // NMOD_POLY_MAT_UTILS_H - -/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ -// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s From ae90690dfb1d72e1f137e520eec3f6a59ee11db1 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sun, 4 Jan 2026 23:37:17 +0100 Subject: [PATCH 07/12] make test more complete --- flint-extras/src/nmod_poly_mat_approximant.h | 8 -- flint-extras/src/nmod_poly_mat_extra/kernel.c | 6 +- .../src/nmod_poly_mat_extra/test/main.c | 12 +- .../src/nmod_poly_mat_extra/test/t-kernel.c | 127 ++++++++++++++++++ .../test/t-kernel_via_approx.c | 71 ---------- .../src/nmod_poly_mat_extra/verification.c | 73 +++++----- flint-extras/src/nmod_poly_mat_forms.h | 36 +++-- flint-extras/src/nmod_poly_mat_kernel.h | 63 +++++---- 8 files changed, 231 insertions(+), 165 deletions(-) create mode 100644 flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c delete mode 100644 flint-extras/src/nmod_poly_mat_extra/test/t-kernel_via_approx.c diff --git a/flint-extras/src/nmod_poly_mat_approximant.h b/flint-extras/src/nmod_poly_mat_approximant.h index 9bb8feff..400367d3 100644 --- a/flint-extras/src/nmod_poly_mat_approximant.h +++ b/flint-extras/src/nmod_poly_mat_approximant.h @@ -311,14 +311,6 @@ void nmod_poly_mat_mbasis(nmod_poly_mat_t appbas, */ //@{ -/** - * TODO old doc, to be cleaned (new doc below, maybe not complete enough) - * F \in K[x]^{nxm} <-> F_prime K^{mxn}[x] FIX and - * Compute P_{k-1} \in K[x]^{mxm} with the list structured_multiplication_blocks - * Compute x^{-k} P_{k-1} F mod x iteratively with a naive polynomial multiplication - * - * Use naive polynomial multiplication and list_structured_multiplication_blocks - */ /** Computes a `shift`-ordered weak Popov approximant basis for `(pmat,order)` * using the algorithm PM-Basis (see @ref pmbasis) */ diff --git a/flint-extras/src/nmod_poly_mat_extra/kernel.c b/flint-extras/src/nmod_poly_mat_extra/kernel.c index d18cfe0b..27fc2e2a 100644 --- a/flint-extras/src/nmod_poly_mat_extra/kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/kernel.c @@ -19,7 +19,6 @@ #include #include "nmod_poly_mat_extra.h" -#include "nmod_poly_mat_io.h" #include "nmod_poly_mat_kernel.h" #include "nmod_poly_mat_multiply.h" @@ -111,13 +110,12 @@ slong nmod_poly_mat_kernel_via_approx(nmod_poly_mat_t ker, return nz; } -/* TODO make pmat non const */ /* Follows the description of Zhou-Labahn-Storjohann algorithm in [LNVZ22, Algo.1] */ /* [LNVZ22] Labahn-Neiger-Vu-Zhou, Proceedings ISSAC 2022, arxiv.org/pdf/2202.09329 */ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, slong * pivind, slong * shift, - const nmod_poly_mat_t pmat) + nmod_poly_mat_t pmat) { /* flint_printf("call %ld, %ld\n", pmat->r, pmat->c); */ /* nmod_poly_mat_degree_matrix_print_pretty(pmat); */ @@ -979,7 +977,7 @@ int nmod_poly_mat_kernel_zls(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat * */ -int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat_t A, \ +int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong * FLINT_UNUSED(degN), const nmod_poly_mat_t A, \ const slong *ishift) { diff --git a/flint-extras/src/nmod_poly_mat_extra/test/main.c b/flint-extras/src/nmod_poly_mat_extra/test/main.c index f0c10bd6..682ecf9c 100644 --- a/flint-extras/src/nmod_poly_mat_extra/test/main.c +++ b/flint-extras/src/nmod_poly_mat_extra/test/main.c @@ -17,9 +17,9 @@ #include "t-det.c" #include "t-dixon.c" #include "t-hermite_normal_form.c" -#include "t-kernel_via_approx.c" -#include "t-kernel_zls_approx.c" -#include "t-kernel_zls.c" +#include "t-kernel.c" +/* #include "t-kernel_zls_approx.c" */ +/* #include "t-kernel_zls.c" */ #include "t-middle_product_geometric.c" #include "t-mul_geometric.c" #include "t-mbasis.c" @@ -32,12 +32,12 @@ test_struct tests[] = { - TEST_FUNCTION(nmod_poly_mat_kernel_zls_approx), TEST_FUNCTION(nmod_poly_mat_det), TEST_FUNCTION(nmod_poly_mat_dixon), TEST_FUNCTION(nmod_poly_mat_hnf), - TEST_FUNCTION(nmod_poly_mat_kernel_via_approx), - TEST_FUNCTION(nmod_poly_mat_kernel_zls), + TEST_FUNCTION(nmod_poly_mat_kernel), + /* TEST_FUNCTION(nmod_poly_mat_kernel_zls_approx), */ + /* TEST_FUNCTION(nmod_poly_mat_kernel_zls), */ TEST_FUNCTION(nmod_poly_mat_mbasis), TEST_FUNCTION(nmod_poly_mat_middle_product_geometric), TEST_FUNCTION(nmod_poly_mat_mul_geometric), diff --git a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c new file mode 100644 index 00000000..636d0762 --- /dev/null +++ b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c @@ -0,0 +1,127 @@ +/* + Copyright (C) 2025 Vincent Neiger + + This file is part of PML. + + PML is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License version 2.0 (GPL-2.0-or-later) + as published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. See + . +*/ + +#include +#include +#include +#include +#include +#include +#include + +#include "nmod_poly_mat_forms.h" +#include "nmod_poly_mat_kernel.h" + +TEST_FUNCTION_START(nmod_poly_mat_kernel, state) +{ + int i, result; + + for (i = 0; i < 100 * flint_test_multiplier(); i++) + { + ulong nbits = 2 + n_randint(state, 63); + slong rdim = n_randint(state, 12); + slong cdim = n_randint(state, 12); + ulong len = 1 + n_randint(state, 150); + + slong * shift = FLINT_ARRAY_ALLOC(rdim, slong); + for (slong i = 0; i < rdim; i++) + shift[i] = n_randint(state, len); + + slong * pivind = FLINT_ARRAY_ALLOC(rdim, slong); + slong * pivind_check = FLINT_ARRAY_ALLOC(rdim, slong); + slong * rdeg = FLINT_ARRAY_ALLOC(rdim, slong); + slong * rdeg_check = FLINT_ARRAY_ALLOC(rdim, slong); + + ulong prime = n_randprime(state, nbits, 1); + + nmod_poly_mat_t pmat; + nmod_poly_mat_init(pmat, rdim, cdim, prime); + nmod_poly_mat_randtest(pmat, state, len); + + nmod_poly_mat_t ker; + nmod_poly_mat_t kernz; + nmod_poly_mat_init(ker, rdim, rdim, prime); + + { + for (slong i = 0; i < rdim; i++) + rdeg[i] = shift[i]; + slong nz = nmod_poly_mat_kernel_via_approx(ker, pivind, rdeg, pmat); + + nmod_poly_mat_window_init(kernz, ker, 0, 0, nz, rdim); + result = nmod_poly_mat_is_kernel(kernz, shift, pmat, ORD_WEAK_POPOV, ROW_LOWER); + + if (!result) + { + TEST_FUNCTION_FAIL("(via_approx, ker) -- rdim = %wu, cdim = %wu, length = %wu, p = %wu\n", \ + rdim, cdim, len, prime); + } + + nmod_poly_mat_row_degree(rdeg_check, kernz, shift); + nmod_poly_mat_pivot_index(pivind_check, kernz, shift, ROW_LOWER); + result = 1; + for (slong k = 0; k < nz; k++) + if (rdeg_check[k] != rdeg[k] || pivind_check[k] != pivind[k]) + result = 0; + + if (!result) + { + TEST_FUNCTION_FAIL("(via_approx, rdeg/pivind) -- rdim = %wu, cdim = %wu, length = %wu, p = %wu\n", \ + rdim, cdim, len, prime); + } + nmod_poly_mat_window_clear(kernz); + } + + { + for (slong i = 0; i < rdim; i++) + rdeg[i] = shift[i]; + nmod_poly_mat_t pmat_c; + nmod_poly_mat_init_set(pmat_c, pmat); + slong nz = nmod_poly_mat_kernel_zls_approx(ker, pivind, rdeg, pmat_c); + + nmod_poly_mat_window_init(kernz, ker, 0, 0, nz, rdim); + result = nmod_poly_mat_is_kernel(kernz, shift, pmat, ORD_WEAK_POPOV, ROW_LOWER); + + if (!result) + { + TEST_FUNCTION_FAIL("(zls_approx, ker) -- rdim = %wu, cdim = %wu, length = %wu, p = %wu\n", \ + rdim, cdim, len, prime); + } + + nmod_poly_mat_row_degree(rdeg_check, kernz, shift); + nmod_poly_mat_pivot_index(pivind_check, kernz, shift, ROW_LOWER); + result = 1; + for (slong k = 0; k < nz; k++) + if (rdeg_check[k] != rdeg[k] || pivind_check[k] != pivind[k]) + result = 0; + + if (!result) + { + TEST_FUNCTION_FAIL("(zls_approx, rdeg/pivind) -- rdim = %wu, cdim = %wu, length = %wu, p = %wu\n", \ + rdim, cdim, len, prime); + } + + nmod_poly_mat_clear(pmat_c); + nmod_poly_mat_window_clear(kernz); + } + + nmod_poly_mat_window_clear(kernz); + nmod_poly_mat_clear(pmat); + nmod_poly_mat_clear(ker); + flint_free(shift); + flint_free(rdeg); + flint_free(pivind); + flint_free(rdeg_check); + flint_free(pivind_check); + } + + TEST_FUNCTION_END(state); +} diff --git a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_via_approx.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_via_approx.c deleted file mode 100644 index 480e0403..00000000 --- a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_via_approx.c +++ /dev/null @@ -1,71 +0,0 @@ -/* - Copyright (C) 2025 Vincent Neiger - - This file is part of PML. - - PML is free software: you can redistribute it and/or modify it under - the terms of the GNU General Public License version 2.0 (GPL-2.0-or-later) - as published by the Free Software Foundation; either version 2 of the - License, or (at your option) any later version. See - . -*/ - -#include -#include -#include -#include -#include -#include - -#include "nmod_poly_mat_kernel.h" - -TEST_FUNCTION_START(nmod_poly_mat_kernel_via_approx, state) -{ - int i, result; - - for (i = 0; i < 16 * flint_test_multiplier(); i++) - { - ulong nbits = 2 + n_randint(state, 63); - slong rdim = n_randint(state, 30); - slong cdim = n_randint(state, 30); - ulong len = 1 + n_randint(state, 100); - flint_printf("TEST iter %ld -- %ld, %ld, %ld\n", i, rdim, cdim, len); - - slong * shift = FLINT_ARRAY_ALLOC(rdim, slong); - for (slong i = 0; i < rdim; i++) - shift[i] = 0; - - slong * pivind = FLINT_ARRAY_ALLOC(rdim, slong); - slong * rdeg = FLINT_ARRAY_ALLOC(rdim, slong); - for (slong i = 0; i < rdim; i++) - rdeg[i] = shift[i]; - - ulong prime = n_randprime(state, nbits, 1); - - nmod_poly_mat_t pmat; - nmod_poly_mat_init(pmat, rdim, cdim, prime); - nmod_poly_mat_randtest(pmat, state, len); - - nmod_poly_mat_t ker; - nmod_poly_mat_init(ker, rdim, rdim, prime); - /* TODO test weak Popov + pivind */ - /* TODO currently does not test generation */ - slong nz = nmod_poly_mat_kernel_via_approx(ker, pivind, rdeg, pmat); - flint_printf("kernel computed, now testing...\n"); - result = nmod_poly_mat_is_kernel(ker, nz, shift, pmat, ROW_LOWER); - - nmod_poly_mat_clear(pmat); - nmod_poly_mat_clear(ker); - flint_free(shift); - flint_free(rdeg); - flint_free(pivind); - - if (!result) - { - TEST_FUNCTION_FAIL("rdim = %wu, cdim = %wu, length = %wu, p = %wu\n", \ - rdim, cdim, len, prime); - } - } - - TEST_FUNCTION_END(state); -} diff --git a/flint-extras/src/nmod_poly_mat_extra/verification.c b/flint-extras/src/nmod_poly_mat_extra/verification.c index 9ba51042..79c5c80c 100644 --- a/flint-extras/src/nmod_poly_mat_extra/verification.c +++ b/flint-extras/src/nmod_poly_mat_extra/verification.c @@ -11,7 +11,6 @@ */ #include "nmod_poly_mat_forms.h" -#include "nmod_poly_mat_io.h" #include #include #include @@ -114,23 +113,40 @@ int nmod_poly_mat_is_approximant_basis(const nmod_poly_mat_t appbas, return success; } -/* TODO currently specialized to ROW_LOWER (or at least ROW_stuff) */ -/* -> checks whether ker is a shift-ordered weak Popov left kernel basis of pmat */ +/* TODO currently does not check generation */ int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker, - slong nz, const slong * shift, - const nmod_poly_mat_t pmat, + const nmod_poly_mat_t mat, + poly_mat_form_t form, orientation_t orient) { - const slong rdim = pmat->r; - const slong cdim = pmat->c; - const ulong prime = pmat->modulus; + if (orient == COL_UPPER || orient == COL_LOWER) + { + nmod_poly_mat_t ker_t; + nmod_poly_mat_t mat_t; + nmod_poly_mat_init(ker_t, ker->c, ker->r, ker->modulus); + nmod_poly_mat_transpose(ker_t, ker); + nmod_poly_mat_init(mat_t, mat->c, mat->r, mat->modulus); + nmod_poly_mat_transpose(mat_t, mat); + + int result; + if (orient == COL_UPPER) + result = nmod_poly_mat_is_kernel(ker_t, shift, mat_t, form, ROW_LOWER); + else /* orient == COL_LOWER */ + result = nmod_poly_mat_is_kernel(ker_t, shift, mat_t, form, ROW_UPPER); + + nmod_poly_mat_clear(ker_t); + nmod_poly_mat_clear(mat_t); + + return result; + } - nmod_poly_mat_t residual; - nmod_poly_mat_init(residual, nz, cdim, prime); + const slong rdim = mat->r; + const slong cdim = mat->c; + const ulong prime = mat->modulus; - nmod_poly_mat_t kernz; - nmod_poly_mat_window_init(kernz, ker, 0, 0, nz, rdim); + nmod_poly_mat_t residual; + nmod_poly_mat_init(residual, rdim, cdim, prime); int success = 1; @@ -141,44 +157,29 @@ int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker, success = 0; } - if (ker->r < nz) + /* check rank */ + slong rk = nmod_poly_mat_rank(mat); + if (ker->r != rdim - rk) { - printf("basis has wrong row dimension\n"); + printf("number of rows does not equal nullity\n"); success = 0; } - /* check kernel is shifted reduced */ - if (!nmod_poly_mat_is_ordered_weak_popov(kernz, shift, orient)) + /* check kernel has the right form */ + if (!nmod_poly_mat_is_form(ker, form, shift, orient)) { - /* printf("basis is not shifted-weak Popov\n"); */ - /* TODO temporarily allow reduced but not s-weak Popov */ - if (!nmod_poly_mat_is_reduced(kernz, shift, orient)) - { - printf("basis is not shifted-reduced\n"); - success = 0; - } + printf("basis does not have the required form\n"); + success = 0; } /* compute residual, check rows of ker are in the kernel */ - /* TODO enhancement: offer randomized option, using Freivalds-like randomized strategy (see ntl-extras) */ - nmod_poly_mat_mul(residual, kernz, pmat); - /* nmod_poly_mat_degree_matrix_print_pretty(residual); */ + nmod_poly_mat_mul(residual, ker, mat); if (!nmod_poly_mat_is_zero(residual)) { printf("not all rows are in the kernel\n"); success = 0; } - /* check rank */ - slong rk = nmod_poly_mat_rank(pmat); - if (nz != rdim - rk) - { - printf("number of rows does not equal nullity\n"); - success = 0; - } - - /* !! TODO check generation !! */ - nmod_poly_mat_clear(residual); return success; diff --git a/flint-extras/src/nmod_poly_mat_forms.h b/flint-extras/src/nmod_poly_mat_forms.h index 113edf7c..899f46b2 100644 --- a/flint-extras/src/nmod_poly_mat_forms.h +++ b/flint-extras/src/nmod_poly_mat_forms.h @@ -52,10 +52,9 @@ extern "C" { * \enum poly_mat_form_t * * Note that the assigned numbers follow the "strength" of these forms: - * - 0<-1<-2<-3<-4: Popov implies ordered weak Popov, which implies weak Popov, + * - 0<1<2<3<4: Popov implies ordered weak Popov, which implies weak Popov, * which implies reduced. - * - 0<-5<-6: Hermite implies echelon. - * + * - 0<5<6: Hermite implies echelon. */ typedef enum { @@ -468,13 +467,37 @@ int nmod_poly_mat_is_popov(const nmod_poly_mat_t mat, orientation_t orient); /** Tests whether `mat` is in echelon form (see @ref MatrixForms) */ -int nmod_poly_mat_is_echelon(const nmod_poly_mat_t pmat, +int nmod_poly_mat_is_echelon(const nmod_poly_mat_t mat, orientation_t orient); /** Tests whether `mat` is in Hermite form (see @ref MatrixForms) */ -int nmod_poly_mat_is_hermite(const nmod_poly_mat_t pmat, +int nmod_poly_mat_is_hermite(const nmod_poly_mat_t mat, orientation_t orient); +/** Tests whether `mat` has the requested form (see @ref MatrixForms) */ +NMOD_POLY_MAT_INLINE +int nmod_poly_mat_is_form(const nmod_poly_mat_t mat, + poly_mat_form_t form, + const slong * shift, + orientation_t orient) +{ + if (form == NONE) + return 1; + if (form == REDUCED) + return nmod_poly_mat_is_reduced(mat, shift, orient); + if (form == WEAK_POPOV) + return nmod_poly_mat_is_weak_popov(mat, shift, orient); + if (form == ORD_WEAK_POPOV) + return nmod_poly_mat_is_ordered_weak_popov(mat, shift, orient); + if (form == POPOV) + return nmod_poly_mat_is_popov(mat, shift, orient); + if (form == ECHELON) + return nmod_poly_mat_is_echelon(mat, orient); + if (form == HERMITE) + return nmod_poly_mat_is_echelon(mat, orient); + return -1; +} + //@} // doxygen group: Testing polynomial matrix forms @@ -831,6 +854,3 @@ slong nmod_poly_mat_ordered_weak_popov_iter(nmod_poly_mat_t mat, #endif #endif // NMOD_POLY_MAT_FORMS_H - -/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ -// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s diff --git a/flint-extras/src/nmod_poly_mat_kernel.h b/flint-extras/src/nmod_poly_mat_kernel.h index caa3125f..2534e38a 100644 --- a/flint-extras/src/nmod_poly_mat_kernel.h +++ b/flint-extras/src/nmod_poly_mat_kernel.h @@ -42,32 +42,35 @@ extern "C" { */ /** \file nmod_poly_mat_kernel.h - * Conventions. - * ------------ - * Apart from the general interfaces (TODO) which offer several choices of + * \anchor kernel_basis + * + * Apart from the general interfaces which offer several choices of * orientation, all other functions compute left kernel bases and use the - * following parameters, where `nz` is the a priori unknown nullity of `pmat`: + * following parameters: + * + * \return slong, the nullity of `pmat` (called `nz` below) * * \param[out] ker the output kernel basis (cannot alias `pmat`, must be * initialized with at least `nz` rows) * \param[out] pivind the pivot index of `ker` (list of integers, length must - * be the number of rows of `pmat`, only the first `nz` entries matter in output) + * be the number of rows of `pmat`) * \param[in,out] shift in: the input shift; and out: the output shifted row * degree of `ker` (list of integers, length must be the number of rows of * `pmat`) * \param[in] pmat the input polynomial matrix (no restriction) * - * \return slong, the nullity of `pmat` + * The computed `ker` is in shifted ordered weak Popov form, or in the + * canonical shifted Popov form when the name of the function indicates so. * - * The computed `ker` is in shifted ordered weak Popov form, or in the canonical - * shifted Popov form when the name of the function indicates so. + * Some functions may overwrite data in pmat (check the presence of `const`). * - * FIXME what output dimensions for `ker`? + * `ker` is not resized/reallocated: it keeps its original shape, but only the + * first `nz` rows contain useful entries. Similarly, for `pivind` and `shift` + * only the first `nz` entries matter in output. */ - /* general interface */ -/* `form` not implemented yet */ +/* `form` not really implemented yet, only shifted weak Popov */ slong nmod_poly_mat_kernel(nmod_poly_mat_t ker, slong * pivind, slong * shift, @@ -79,31 +82,32 @@ slong nmod_poly_mat_kernel(nmod_poly_mat_t ker, // ROW_UPPER -> mirror rows of input, call, mirror columns+rows of output // COL_LOWER -> -/* using approximant basis at sufficiently large order */ /** Computes a `shift`-ordered weak Popov left kernel basis `ker` for `pmat`, - * recovered as a subset of the rows of a minimal approximant basis at - * sufficiently large order. Computes the - * `shift`-pivot index `pivind` of `kerbas`, and `shift` becomes the shifted - * row degree of `kerbas` (for the input shift). + * through a minimal approximant basis at sufficiently large order. */ - slong nmod_poly_mat_kernel_via_approx(nmod_poly_mat_t ker, slong * pivind, slong * shift, const nmod_poly_mat_t pmat); +/** Computes a `shift`-ordered weak Popov left kernel basis `ker` for `pmat`, + * using the Zhou-Labahn-Storjohann algorithm: + * Wei Zhou, George Labahn, Arne Storjohann -- "Computing Minimal Nullspace Bases" + * Proceedings ISSAC 2012 -- https://doi.org/10.1145/2442829.2442881 + */ +/* TODO middle_product currently naive */ +slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, + slong * pivind, + slong * shift, + nmod_poly_mat_t pmat); + + /* FIXME not implemented yet: using interpolant basis at sufficiently many points */ /* slong nmod_poly_mat_kernel_via_interp(nmod_poly_mat_t ker, */ /* slong * shift, */ /* const nmod_poly_mat_t pmat); */ -/* TODO */ -slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, - slong * pivind, - slong * shift, - const nmod_poly_mat_t pmat); - -/* TODO */ +/* FIXME not implemented yet: ZLS algo using interpolant basis */ /* slong nmod_poly_mat_kernel_zls_interp(nmod_poly_mat_t ker, */ /* slong * shift, */ /* const nmod_poly_mat_t pmat); */ @@ -117,28 +121,23 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, * `pmat` for the required form `form`, with orientation described by `orient`. * * \param[in] ker kernel basis - * \param[in] nz nullity associated to ker * \param[in] shift shift * \param[in] pmat polynomial matrix - * \param[in] form [not implemented yet] required form for `kerbas` (see #poly_mat_form_t) + * \param[in] form required form for `kerbas` (see #poly_mat_form_t) * \param[in] orient indicates the orientation (left/right kernel) and the definition of pivots * \param[in] randomized [not implemented yet] if `true`, the algorithm may use a Monte Carlo or Las Vegas verification algorithm * * \return int, result of the verification * - * \todo support all options, make doc more clear concerning Las Vegas / Monte Carlo - * \todo WARNING! for the moment, does not really check generation! - * \todo WARNING! for the moment, hardcoded to check for ordered weak Popov */ +/* TODO WARNING! for the moment, does not really check generation! */ int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker, - slong nz, - /* const slong * pivind, */ /* TODO */ const slong * shift, const nmod_poly_mat_t pmat, + poly_mat_form_t form, orientation_t orient); - /** * * Right shifted kernel of a polynomial matrix, assuming that the columns has been From 0bb1805f99c66d182b19ff1bb6bc9356f35a9117 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sun, 4 Jan 2026 23:50:00 +0100 Subject: [PATCH 08/12] make test more complete for ZLS --- .../src/nmod_poly_mat_extra/test/main.c | 6 +- .../nmod_poly_mat_extra/test/t-kernel_zls.c | 55 ++++++-------- .../test/t-kernel_zls_approx.c | 74 ------------------- 3 files changed, 25 insertions(+), 110 deletions(-) delete mode 100644 flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls_approx.c diff --git a/flint-extras/src/nmod_poly_mat_extra/test/main.c b/flint-extras/src/nmod_poly_mat_extra/test/main.c index 682ecf9c..0c8fe9d4 100644 --- a/flint-extras/src/nmod_poly_mat_extra/test/main.c +++ b/flint-extras/src/nmod_poly_mat_extra/test/main.c @@ -18,8 +18,7 @@ #include "t-dixon.c" #include "t-hermite_normal_form.c" #include "t-kernel.c" -/* #include "t-kernel_zls_approx.c" */ -/* #include "t-kernel_zls.c" */ +#include "t-kernel_zls.c" #include "t-middle_product_geometric.c" #include "t-mul_geometric.c" #include "t-mbasis.c" @@ -36,8 +35,7 @@ test_struct tests[] = TEST_FUNCTION(nmod_poly_mat_dixon), TEST_FUNCTION(nmod_poly_mat_hnf), TEST_FUNCTION(nmod_poly_mat_kernel), - /* TEST_FUNCTION(nmod_poly_mat_kernel_zls_approx), */ - /* TEST_FUNCTION(nmod_poly_mat_kernel_zls), */ + TEST_FUNCTION(nmod_poly_mat_kernel_zls), TEST_FUNCTION(nmod_poly_mat_mbasis), TEST_FUNCTION(nmod_poly_mat_middle_product_geometric), TEST_FUNCTION(nmod_poly_mat_mul_geometric), diff --git a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls.c index 0e7358df..0cbce91a 100644 --- a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls.c +++ b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls.c @@ -11,64 +11,55 @@ */ #include +#include +#include #include #include #include #include -#include "nmod_poly_mat_extra.h" #include "nmod_poly_mat_forms.h" +#include "nmod_poly_mat_kernel.h" // test one given input -/* TODO does not test generation */ -/* TODO tests reducedness instead of weak Popov */ -/* TODO does not test shifts */ -int core_test_kernel_zls(const nmod_poly_mat_t mat) +int core_test_kernel_zls(const nmod_poly_mat_t mat, flint_rand_t state) { - slong m = mat->r; slong n = mat->c; - /* careful: this is currently the column degree due to working on the right */ - slong * rdeg = FLINT_ARRAY_ALLOC(n, slong); - nmod_poly_mat_column_degree(rdeg, mat, NULL); + /* pick shift with entries >= cdeg */ + slong * cdeg = FLINT_ARRAY_ALLOC(n, slong); + nmod_poly_mat_column_degree(cdeg, mat, NULL); + slong * shift = FLINT_ARRAY_ALLOC(n, slong); + for (slong j = 0; j < n; j++) + shift[j] = FLINT_MAX(0, cdeg[j]) + n_randint(state, 30); nmod_poly_mat_t N; nmod_poly_mat_init(N, n, n, mat->modulus); - slong degN[n]; + slong nz = nmod_poly_mat_kernel_zls(N, degN, mat, shift, 2.); - slong nz = nmod_poly_mat_kernel_zls(N, degN, mat, NULL, 2.); - - nmod_poly_mat_t Nt; - nmod_poly_mat_init(Nt, nz, n, mat->modulus); - for (long i = 0; i < n; i++) - for (long j = 0; j < nz; j++) - nmod_poly_set(nmod_poly_mat_entry(Nt, j, i), nmod_poly_mat_entry(N, i, j)); - - nmod_poly_mat_t Mt; - nmod_poly_mat_init(Mt, n, m, mat->modulus); - nmod_poly_mat_transpose(Mt, mat); - - int verif = nmod_poly_mat_is_kernel(Nt, nz, rdeg, Mt, ROW_LOWER); + nmod_poly_mat_t Nnz; + nmod_poly_mat_window_init(Nnz, N, 0, 0, n, nz); + int verif = nmod_poly_mat_is_kernel(Nnz, shift, mat, REDUCED, COL_UPPER); nmod_poly_mat_clear(N); - nmod_poly_mat_clear(Nt); - nmod_poly_mat_clear(Mt); - flint_free(rdeg); + nmod_poly_mat_window_clear(Nnz); + flint_free(cdeg); + flint_free(shift); return verif; } TEST_FUNCTION_START(nmod_poly_mat_kernel_zls, state) { - int i,result; + int i, result; for (i = 0; i < 16 * flint_test_multiplier(); i++) { ulong nbits = 2 + n_randint(state, 63); - ulong rdim = 1 + n_randint(state, 60); + ulong rdim = 1 + n_randint(state, 20); ulong cdim = rdim + 1 + n_randint(state, 20); - ulong deg = n_randint(state, 20); + ulong deg = n_randint(state, 100); ulong prime = n_randprime(state, nbits, 1); @@ -97,7 +88,7 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel_zls, state) nmod_poly_mat_randtest_sparse(A, state, deg+1, 0.84); } - result = core_test_kernel_zls(A); + result = core_test_kernel_zls(A, state); nmod_poly_mat_clear(A); @@ -124,7 +115,7 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel_zls, state) nmod_poly_mat_init(A, rdim, cdim, prime); nmod_poly_mat_randtest_sparse(A, state, deg+1, 0.8); - result = core_test_kernel_zls(A); + result = core_test_kernel_zls(A, state); nmod_poly_mat_clear(A); @@ -151,7 +142,7 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel_zls, state) nmod_poly_mat_init(A, rdim, cdim, prime); nmod_poly_mat_randtest_sparse(A, state, deg+1, 0.2); - result = core_test_kernel_zls(A); + result = core_test_kernel_zls(A, state); nmod_poly_mat_clear(A); diff --git a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls_approx.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls_approx.c deleted file mode 100644 index 5fa8f39c..00000000 --- a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls_approx.c +++ /dev/null @@ -1,74 +0,0 @@ -/* - Copyright (C) 2025 Vincent Neiger - - This file is part of PML. - - PML is free software: you can redistribute it and/or modify it under - the terms of the GNU General Public License version 2.0 (GPL-2.0-or-later) - as published by the Free Software Foundation; either version 2 of the - License, or (at your option) any later version. See - . -*/ - -#include -#include -#include -#include -#include -#include - -#include "nmod_poly_mat_kernel.h" - -TEST_FUNCTION_START(nmod_poly_mat_kernel_zls_approx, state) -{ - int i, result; - - for (i = 0; i < 16 * flint_test_multiplier(); i++) - { - ulong nbits = 2 + n_randint(state, 63); - slong rdim = n_randint(state, 30); - slong cdim = n_randint(state, 30); - ulong len = 1 + n_randint(state, 100); - flint_printf("TEST iter %ld -- %ld, %ld, %ld\n", i, rdim, cdim, len); - - slong * shift = FLINT_ARRAY_ALLOC(rdim, slong); - for (slong i = 0; i < rdim; i++) - shift[i] = 0; - - slong * pivind = FLINT_ARRAY_ALLOC(rdim, slong); - slong * rdeg = FLINT_ARRAY_ALLOC(rdim, slong); - for (slong i = 0; i < rdim; i++) - rdeg[i] = shift[i]; - - ulong prime = n_randprime(state, nbits, 1); - - nmod_poly_mat_t pmat; - nmod_poly_mat_init(pmat, rdim, cdim, prime); - nmod_poly_mat_randtest(pmat, state, len); - nmod_poly_mat_t copy_pmat; - nmod_poly_mat_init_set(copy_pmat, pmat); - - nmod_poly_mat_t ker; - nmod_poly_mat_init(ker, rdim, rdim, prime); - /* TODO test weak Popov + pivind */ - /* TODO currently does not test generation */ - slong nz = nmod_poly_mat_kernel_zls_approx(ker, pivind, rdeg, pmat); - flint_printf("kernel computed, now testing...\n"); - result = nmod_poly_mat_is_kernel(ker, nz, shift, copy_pmat, ROW_LOWER); - - nmod_poly_mat_clear(ker); - nmod_poly_mat_clear(pmat); - nmod_poly_mat_clear(copy_pmat); - flint_free(shift); - flint_free(rdeg); - flint_free(pivind); - - if (!result) - { - TEST_FUNCTION_FAIL("rdim = %wu, cdim = %wu, length = %wu, p = %wu\n", \ - rdim, cdim, len, prime); - } - } - - TEST_FUNCTION_END(state); -} From 9b023a2fc075de850eede4abb1e7ae18061d1333 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Mon, 5 Jan 2026 12:11:52 +0100 Subject: [PATCH 09/12] early exit --- flint-extras/src/nmod_poly_mat_extra/kernel.c | 75 ++++++------------- .../nmod_poly_mat_extra/profile/p-kernel.c | 63 ++++++++++++---- .../profile/p-middle_product.c | 12 +++ .../src/nmod_poly_mat_extra/profile/p-mul.c | 12 +++ .../src/nmod_poly_mat_extra/test/t-kernel.c | 2 +- 5 files changed, 95 insertions(+), 69 deletions(-) diff --git a/flint-extras/src/nmod_poly_mat_extra/kernel.c b/flint-extras/src/nmod_poly_mat_extra/kernel.c index 27fc2e2a..697a6b93 100644 --- a/flint-extras/src/nmod_poly_mat_extra/kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/kernel.c @@ -13,12 +13,14 @@ #include #include +#include #include #include #include #include #include "nmod_poly_mat_extra.h" +#include "nmod_poly_mat_forms.h" #include "nmod_poly_mat_kernel.h" #include "nmod_poly_mat_multiply.h" @@ -117,8 +119,6 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, slong * shift, nmod_poly_mat_t pmat) { - /* flint_printf("call %ld, %ld\n", pmat->r, pmat->c); */ - /* nmod_poly_mat_degree_matrix_print_pretty(pmat); */ const slong m = pmat->r; const slong n = pmat->c; @@ -126,7 +126,6 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, if (m == 0 || n == 0 || ((m == 1 || n == 1) && (nmod_poly_mat_max_length(pmat) == 0))) { - /* flint_printf("base case empty||zero!\n"); */ nmod_poly_mat_one(ker); for (slong i = 0; i < m; i++) pivind[i] = i; @@ -137,11 +136,6 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, if (m == 1) return 0; - /* one (nonzero) column vector */ - /* TODO */ - /* if (n == 1) */ - /* return 0; */ - /* if n > m/2 (with here m >= 2), straightforward recursion by splitting */ /* the matrix into two submatrices with n/2 columns */ if (2*n > m) @@ -224,7 +218,6 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, /* early exit: pmat == 0 => kernel is identity */ if (maxdeg == -1) { - /* flint_printf("second base case zero!\n"); */ nmod_poly_mat_one(ker); for (slong i = 0; i < m; i++) pivind[i] = i; @@ -248,6 +241,7 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, /* run an easy degree-based detection of rows in the kernel */ /* permute ker into [kernel rows \\ other rows], preserving increasing pivots */ slong nz = 0; + slong sum_pmatdeg = 0; /* sum of degrees of rows of pmat in complement of pivind */ for (slong i = 0; i < m; i++) { if (shift[i] < order + diff_shift) @@ -256,44 +250,29 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, nz += 1; } else + { + sum_pmatdeg += FLINT_MAX(0, buf[i]); buf[i - nz] = i; + } } for (slong i = nz; i < m; i++) pivind[i] = buf[i - nz]; nmod_poly_mat_permute_rows(ker, pivind, shift); - /* flint_printf("pmbasis ok, found %ld rows\n", nz); */ - /* flint_printf("pivind: %{slong*}\n", pivind, m); */ - /* nmod_poly_mat_degree_matrix_print_pretty(ker); */ /* early exit */ - /* TODO */ -// // if the sum of pivot degrees in the part of kernel basis computed is -// // equal to one of the upper bounds (here, sum of row degrees of pmat for -// // non-pivot rows), then we know that we have captured the whole kernel -// // --> this is often the case (e.g. in the generic situation mentioned -// // above), and testing this avoids going through a few multiplications and -// // the two recursive calls (which only lead to the conclusion that there is -// // no new kernel row to be found) -// long sum_matdeg = 0; -// for (long i = 0; i < m2; ++i) -// sum_matdeg += rdeg[rows_app_notker[i]]; -// long sum_kerdeg = 0; -// for (long i = 0; i < m1; ++i) -// sum_kerdeg += deg(appbas[rows_app_ker[i]][rows_app_ker[i]]); -// -// const bool early_exit = (sum_kerdeg == sum_matdeg); -// -// // if the whole kernel was captured according to the above test, or if the kernel is full (matrix was zero), -// // then we have the whole kernel: just copy and return -// if (early_exit) -// { -// kerbas.SetDims(m1, m); -// for (long i = 0; i < m1; ++i) -// kerbas[i].swap(appbas[rows_app_ker[i]]); -// shift.swap(rdeg1); -// return; -// } + if (nz >= m - n) /* otherwise, there must be some kernel rows not in ker yet */ + { + slong sum_pivdeg = 0; + for (slong i = 0; i < nz; i++) + sum_pivdeg += nmod_poly_degree(nmod_poly_mat_entry(ker, i, pivind[i])); + + if (sum_pivdeg == sum_pmatdeg) /* whole kernel already found */ + { + flint_free(buf); + return nz; + } + } /* proceed with remaining rows */ nmod_poly_mat_t residual; /* actual init */ @@ -306,14 +285,11 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, nmod_poly_mat_window_init(approx, ker, nz, 0, m, m); nmod_poly_mat_middle_product_naive(residual, approx, pmat, order, order + maxdeg + 1); /* FIXME provide (and use) general middle_product interface */ - - /* TODO handle zero rows in G? */ + /* note: on non-generic input, one might want to check for zero rows in residual */ /* kernel of residual */ nmod_poly_mat_init(ker2, m - nz, m - nz, pmat->modulus); const slong nz2 = nmod_poly_mat_kernel_zls_approx(ker2, buf, shift + nz, residual); - /* flint_printf("second kernel ok, found %ld rows\n", nz2); */ - /* flint_printf("pivind: %{slong*}\n", buf, m - nz); */ /* if nz2 == 0, whole kernel already computed, just return */ if (nz2 == 0) @@ -326,22 +302,14 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, } /* multiply to get missing part of kernel */ - /* flint_printf("nz2 == %ld\n", nz2); */ - /* nmod_poly_mat_degree_matrix_print_pretty(ker2); */ - /* nmod_poly_mat_degree_matrix_print_pretty(approx); */ nmod_poly_mat_window_init(ker2nz, ker2, 0, 0, nz2, ker2->c); nmod_poly_mat_init(prod, nz2, m, pmat->modulus); nmod_poly_mat_mul(prod, ker2nz, approx); - /* flint_printf("before swap::\n"); */ - /* nmod_poly_mat_degree_matrix_print_pretty(prod); */ - /* nmod_poly_mat_degree_matrix_print_pretty(ker); */ for (slong i = 0; i < nz2; i++) for (slong j = 0; j < m; j++) FLINT_SWAP(nmod_poly_struct, *nmod_poly_mat_entry(prod, i, j), *nmod_poly_mat_entry(approx, i, j)); - /* flint_printf("after swap::\n"); */ - /* nmod_poly_mat_degree_matrix_print_pretty(ker); */ /* update pivind */ for (slong i = 0; i < nz2; i++) @@ -365,9 +333,9 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, i2++; } } - /* flint_printf("permutation buf: %{slong*}\n", buf, nz + nz2); */ + nmod_poly_mat_permute_rows(approx, buf, NULL); - slong * tmp = FLINT_ARRAY_ALLOC(nz + nz2, slong); /* use TMP_ALLOC */ + slong * tmp = FLINT_ARRAY_ALLOC(nz + nz2, slong); for (slong i = 0; i < nz + nz2; i++) tmp[i] = pivind[i]; for (slong i = 0; i < nz + nz2; i++) @@ -385,7 +353,6 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, nmod_poly_mat_clear(prod); flint_free(buf); - /* nmod_poly_mat_degree_matrix_print_pretty(ker); */ return nz + nz2; } diff --git a/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c b/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c index a2c5c90b..f25ab657 100644 --- a/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c @@ -1,3 +1,15 @@ +/* + Copyright (C) 2026 Vincent Neiger + + This file is part of PML. + + PML is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License version 2.0 (GPL-2.0-or-later) + as published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. See + . +*/ + #include // for atoi #include @@ -112,6 +124,36 @@ void time_##fun(time_args targs, flint_rand_t state) \ nmod_poly_mat_clear(ker); \ } +void time_nullspace(time_args targs, flint_rand_t state) +{ + const slong rdim = targs.rdim; + const slong cdim = targs.cdim; + const slong deg = targs.deg; + /* const slong rank = targs.rank; */ /* TODO */ + const slong n = targs.modn; + + nmod_t mod; + nmod_init(&mod, n); + + nmod_poly_mat_t pmat; + nmod_poly_mat_init(pmat, cdim, rdim, n); + nmod_poly_mat_rand(pmat, state, deg); + + nmod_poly_mat_t ker; + nmod_poly_mat_init(ker, rdim, rdim, n); + + double FLINT_SET_BUT_UNUSED(tcpu), twall; + + TIMEIT_START; + nmod_poly_mat_nullspace(ker, pmat); + TIMEIT_STOP_VALUES(tcpu, twall); + + flint_printf("%.2e", twall); + + nmod_poly_mat_clear(pmat); + nmod_poly_mat_clear(ker); +} + TIME_KER_ZLS(kernel_zls) TIME_KER(kernel_via_approx) TIME_KER(kernel_zls_approx) @@ -126,15 +168,6 @@ int main(int argc, char ** argv) flint_rand_init(state); flint_rand_set_seed(state, time(NULL), time(NULL)+129384125L); - // modulus bitsize - /* const slong nbits = 7; */ - /* const ulong bits[] = {12, 24, 30, 40, 50, 60, 63}; */ - - // matrix dimensions - /* const slong ndims = 10; */ - /* const ulong rdims[] = {2, 4, 6, 8, 11, 15, 20, 30, 50, 100}; */ - /* const ulong cdims[] = {1, 3, 5, 7, 9, 13, 17, 25, 40, 75}; */ - // matrix degrees const slong ndegs = 13; const ulong degs[] = {2, 5, 10, 20, 40, 80, 160, 320, 640, 1280, 2560, 5120, 10240}; @@ -143,12 +176,13 @@ int main(int argc, char ** argv) /* TODO shift type */ // bench functions - const slong nfuns = 3; + const slong nfuns = 4; typedef void (*timefun) (time_args, flint_rand_t); const timefun funs[] = { time_kernel_via_approx, // 0 time_kernel_zls_approx, // 1 - time_kernel_zls, // 2 + time_nullspace, // 2 + time_kernel_zls, // 3 }; // TODO @@ -168,9 +202,10 @@ int main(int argc, char ** argv) //#endif const char * description[] = { - "#0 --> via approx ", - "#1 --> ZLS via approx ", - "#2 --> .... TODO ", + "#0 --> via approx ", + "#1 --> ZLS via approx ", + "#2 --> FLINT's nullspace ", + "#3 --> ZLS via approx(bis) ", }; if (argc < 4 || argc > 6) // show usage diff --git a/flint-extras/src/nmod_poly_mat_extra/profile/p-middle_product.c b/flint-extras/src/nmod_poly_mat_extra/profile/p-middle_product.c index 93ed4fe7..af579728 100644 --- a/flint-extras/src/nmod_poly_mat_extra/profile/p-middle_product.c +++ b/flint-extras/src/nmod_poly_mat_extra/profile/p-middle_product.c @@ -1,3 +1,15 @@ +/* + Copyright (C) 2025 Vincent Neiger + + This file is part of PML. + + PML is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License version 2.0 (GPL-2.0-or-later) + as published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. See + . +*/ + #include // for atoi #include diff --git a/flint-extras/src/nmod_poly_mat_extra/profile/p-mul.c b/flint-extras/src/nmod_poly_mat_extra/profile/p-mul.c index 913204c0..d8cf659c 100644 --- a/flint-extras/src/nmod_poly_mat_extra/profile/p-mul.c +++ b/flint-extras/src/nmod_poly_mat_extra/profile/p-mul.c @@ -1,3 +1,15 @@ +/* + Copyright (C) 2025 Vincent Neiger + + This file is part of PML. + + PML is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License version 2.0 (GPL-2.0-or-later) + as published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. See + . +*/ + #include // for atoi #include diff --git a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c index 636d0762..cdf6753e 100644 --- a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c @@ -28,8 +28,8 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state) for (i = 0; i < 100 * flint_test_multiplier(); i++) { ulong nbits = 2 + n_randint(state, 63); - slong rdim = n_randint(state, 12); slong cdim = n_randint(state, 12); + slong rdim = n_randint(state, 12); ulong len = 1 + n_randint(state, 150); slong * shift = FLINT_ARRAY_ALLOC(rdim, slong); From 1791de25f3d5d61c6add9c8725e7f949d83c5ab6 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Mon, 5 Jan 2026 13:38:44 +0100 Subject: [PATCH 10/12] general interface --- flint-extras/src/nmod_poly_mat_extra/kernel.c | 101 +++++++++++++- .../src/nmod_poly_mat_extra/test/t-kernel.c | 131 ++++++++++++++++- .../src/nmod_poly_mat_extra/verification.c | 132 +++++++++++------- flint-extras/src/nmod_poly_mat_forms.h | 6 +- flint-extras/src/nmod_poly_mat_kernel.h | 15 +- 5 files changed, 310 insertions(+), 75 deletions(-) diff --git a/flint-extras/src/nmod_poly_mat_extra/kernel.c b/flint-extras/src/nmod_poly_mat_extra/kernel.c index 697a6b93..3ae3977c 100644 --- a/flint-extras/src/nmod_poly_mat_extra/kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/kernel.c @@ -23,14 +23,101 @@ #include "nmod_poly_mat_forms.h" #include "nmod_poly_mat_kernel.h" #include "nmod_poly_mat_multiply.h" +#include "nmod_poly_mat_utils.h" + +slong nmod_poly_mat_kernel(nmod_poly_mat_t ker, + slong * pivind, + slong * shift, + const nmod_poly_mat_t pmat, + poly_mat_form_t form, + orientation_t orient) +{ + if (form > ORD_WEAK_POPOV) + flint_throw(FLINT_ERROR, "Exception (nmod_poly_mat_kernel). form > ORD_WEAK_POPOV not implemented."); + + if (orient == ROW_LOWER) + { + nmod_poly_mat_t mat; + nmod_poly_mat_init_set(mat, pmat); + slong nz = nmod_poly_mat_kernel_zls_approx(ker, pivind, shift, mat); + nmod_poly_mat_clear(mat); + return nz; + } + + if (orient == COL_UPPER) + { + /* transpose input, call ROW_LOWER, transpose output */ + nmod_poly_mat_t mat_t; + nmod_poly_mat_t ker_t; + nmod_poly_mat_init(mat_t, pmat->c, pmat->r, pmat->modulus); + nmod_poly_mat_init(ker_t, ker->c, ker->r, ker->modulus); + nmod_poly_mat_transpose(mat_t, pmat); + slong nz = nmod_poly_mat_kernel_zls_approx(ker_t, pivind, shift, mat_t); + nmod_poly_mat_transpose(ker, ker_t); + nmod_poly_mat_clear(ker_t); + nmod_poly_mat_clear(mat_t); + return nz; + } + + if (orient == ROW_UPPER) + { + /* mirror rows of input, call ROW_LOWER, mirror columns+rows of output */ + nmod_poly_mat_t mat_i; + nmod_poly_mat_init_set(mat_i, pmat); + nmod_poly_mat_invert_rows(mat_i, shift); + slong nz = nmod_poly_mat_kernel_zls_approx(ker, pivind, shift, mat_i); + nmod_poly_mat_t kernz; + nmod_poly_mat_window_init(kernz, ker, 0, 0, nz, ker->c); + nmod_poly_mat_invert_rows(kernz, shift); + nmod_poly_mat_invert_columns(kernz, NULL); + for (slong i = 0; i < nz/2; i++) + { + slong tmp = ker->c - 1 - pivind[i]; + pivind[i] = ker->c - 1 - pivind[nz - 1 - i]; + pivind[nz - 1 - i] = tmp; + } + if (nz % 2) + pivind[nz/2] = ker->c - 1 - pivind[nz/2]; + nmod_poly_mat_window_clear(kernz); + nmod_poly_mat_clear(mat_i); + return nz; + } + + if (orient == COL_LOWER) + { + /* transpose input, call ROW_UPPER, transpose output */ + nmod_poly_mat_t mat_it; + nmod_poly_mat_t ker_it; + nmod_poly_mat_init(mat_it, pmat->c, pmat->r, pmat->modulus); + nmod_poly_mat_init(ker_it, ker->c, ker->r, ker->modulus); + nmod_poly_mat_transpose(mat_it, pmat); + + nmod_poly_mat_invert_rows(mat_it, shift); + slong nz = nmod_poly_mat_kernel_zls_approx(ker_it, pivind, shift, mat_it); + nmod_poly_mat_t kernz; + nmod_poly_mat_window_init(kernz, ker_it, 0, 0, nz, ker_it->c); + nmod_poly_mat_invert_rows(kernz, shift); + nmod_poly_mat_invert_columns(kernz, NULL); + for (slong i = 0; i < nz/2; i++) + { + slong tmp = ker_it->c - 1 - pivind[i]; + pivind[i] = ker_it->c - 1 - pivind[nz - 1 - i]; + pivind[nz - 1 - i] = tmp; + } + if (nz % 2) + pivind[nz/2] = ker_it->c - 1 - pivind[nz/2]; + nmod_poly_mat_window_clear(kernz); + + nmod_poly_mat_transpose(ker, ker_it); + nmod_poly_mat_clear(ker_it); + nmod_poly_mat_clear(mat_it); + return nz; + } + + flint_throw(FLINT_ERROR, "Exception (nmod_poly_mat_kernel). Requested orientation not implemented."); +} + -/* general interface: - * A/ compute rdeg and cdeg - * ->suppress zero columns - * ->trivialize zero rows - * B/ handle case of constant matrices - * C/ go for zls - * */ slong nmod_poly_mat_kernel_via_approx(nmod_poly_mat_t ker, slong * pivind, diff --git a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c index cdf6753e..5565c7df 100644 --- a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c @@ -23,14 +23,15 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state) { - int i, result; + int k, result; - for (i = 0; i < 100 * flint_test_multiplier(); i++) + /* left kernel */ + for (k = 0; k < 50 * flint_test_multiplier(); k++) { ulong nbits = 2 + n_randint(state, 63); slong cdim = n_randint(state, 12); slong rdim = n_randint(state, 12); - ulong len = 1 + n_randint(state, 150); + ulong len = 1 + n_randint(state, 75); slong * shift = FLINT_ARRAY_ALLOC(rdim, slong); for (slong i = 0; i < rdim; i++) @@ -51,6 +52,7 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state) nmod_poly_mat_t kernz; nmod_poly_mat_init(ker, rdim, rdim, prime); + /* kernel_via_approx */ { for (slong i = 0; i < rdim; i++) rdeg[i] = shift[i]; @@ -68,8 +70,8 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state) nmod_poly_mat_row_degree(rdeg_check, kernz, shift); nmod_poly_mat_pivot_index(pivind_check, kernz, shift, ROW_LOWER); result = 1; - for (slong k = 0; k < nz; k++) - if (rdeg_check[k] != rdeg[k] || pivind_check[k] != pivind[k]) + for (slong i = 0; i < nz; i++) + if (rdeg_check[i] != rdeg[i] || pivind_check[i] != pivind[i]) result = 0; if (!result) @@ -80,6 +82,7 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state) nmod_poly_mat_window_clear(kernz); } + /* kernel_zls_approx */ { for (slong i = 0; i < rdim; i++) rdeg[i] = shift[i]; @@ -99,8 +102,8 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state) nmod_poly_mat_row_degree(rdeg_check, kernz, shift); nmod_poly_mat_pivot_index(pivind_check, kernz, shift, ROW_LOWER); result = 1; - for (slong k = 0; k < nz; k++) - if (rdeg_check[k] != rdeg[k] || pivind_check[k] != pivind[k]) + for (slong i = 0; i < nz; i++) + if (rdeg_check[i] != rdeg[i] || pivind_check[i] != pivind[i]) result = 0; if (!result) @@ -113,6 +116,44 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state) nmod_poly_mat_window_clear(kernz); } + /* kernel interface */ + { + poly_mat_form_t form = ORD_WEAK_POPOV; + orientation_t orient = (n_randint(state, 2)) ? ROW_LOWER : ROW_UPPER; + + for (slong i = 0; i < rdim; i++) + rdeg[i] = shift[i]; + slong nz = nmod_poly_mat_kernel(ker, pivind, rdeg, pmat, form, orient); + + nmod_poly_mat_window_init(kernz, ker, 0, 0, nz, rdim); + result = nmod_poly_mat_is_kernel(kernz, shift, pmat, form, orient); + + if (!result) + { + TEST_FUNCTION_FAIL("(kernel, ker) -- rdim = %wu, cdim = %wu, length = %wu, p = %wu, orient = %wu\n", \ + rdim, cdim, len, prime, orient); + } + + nmod_poly_mat_row_degree(rdeg_check, kernz, shift); + nmod_poly_mat_pivot_index(pivind_check, kernz, shift, orient); + result = 1; + for (slong i = 0; i < nz; i++) + if (rdeg_check[i] != rdeg[i] || pivind_check[i] != pivind[i]) + result = 0; + + if (!result) + { + flint_printf("true rdeg == %{slong*}\n", rdeg, nz); + flint_printf("got rdeg == %{slong*}\n", rdeg_check, nz); + flint_printf("true pivind == %{slong*}\n", pivind, nz); + flint_printf("got pivind == %{slong*}\n", pivind_check, nz); + TEST_FUNCTION_FAIL("(kernel, rdeg/pivind) -- rdim = %wu, cdim = %wu, length = %wu, p = %wu, orient = %wu\n", \ + rdim, cdim, len, prime, orient); + } + + nmod_poly_mat_window_clear(kernz); + } + nmod_poly_mat_window_clear(kernz); nmod_poly_mat_clear(pmat); nmod_poly_mat_clear(ker); @@ -123,5 +164,81 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state) flint_free(pivind_check); } + /* right kernel */ + for (k = 0; k < 50 * flint_test_multiplier(); k++) + { + ulong nbits = 2 + n_randint(state, 63); + slong cdim = n_randint(state, 12); + slong rdim = n_randint(state, 12); + ulong len = 1 + n_randint(state, 75); + + slong * shift = FLINT_ARRAY_ALLOC(cdim, slong); + for (slong j = 0; j < cdim; j++) + shift[j] = n_randint(state, len); + + slong * pivind = FLINT_ARRAY_ALLOC(cdim, slong); + slong * pivind_check = FLINT_ARRAY_ALLOC(cdim, slong); + slong * cdeg = FLINT_ARRAY_ALLOC(cdim, slong); + slong * cdeg_check = FLINT_ARRAY_ALLOC(cdim, slong); + + ulong prime = n_randprime(state, nbits, 1); + + nmod_poly_mat_t pmat; + nmod_poly_mat_init(pmat, rdim, cdim, prime); + nmod_poly_mat_randtest(pmat, state, len); + + nmod_poly_mat_t ker; + nmod_poly_mat_t kernz; + nmod_poly_mat_init(ker, cdim, cdim, prime); + + /* kernel interface */ + { + poly_mat_form_t form = ORD_WEAK_POPOV; + /* orientation_t orient = (n_randint(state, 2)) ? COL_LOWER : COL_UPPER; */ + orientation_t orient = COL_UPPER; + + for (slong j = 0; j < cdim; j++) + cdeg[j] = shift[j]; + slong nz = nmod_poly_mat_kernel(ker, pivind, cdeg, pmat, form, orient); + + nmod_poly_mat_window_init(kernz, ker, 0, 0, cdim, nz); + result = nmod_poly_mat_is_kernel(kernz, shift, pmat, form, orient); + + if (!result) + { + TEST_FUNCTION_FAIL("(kernel, ker) -- rdim = %wu, cdim = %wu, length = %wu, p = %wu, orient = %wu\n", \ + rdim, cdim, len, prime, orient); + } + + nmod_poly_mat_column_degree(cdeg_check, kernz, shift); + nmod_poly_mat_pivot_index(pivind_check, kernz, shift, orient); + result = 1; + for (slong k = 0; k < nz; k++) + if (cdeg_check[k] != cdeg[k] || pivind_check[k] != pivind[k]) + result = 0; + + if (!result) + { + flint_printf("true cdeg == %{slong*}\n", cdeg, nz); + flint_printf("got cdeg == %{slong*}\n", cdeg_check, nz); + flint_printf("true pivind == %{slong*}\n", pivind, nz); + flint_printf("got pivind == %{slong*}\n", pivind_check, nz); + TEST_FUNCTION_FAIL("(kernel, rdeg/pivind) -- rdim = %wu, cdim = %wu, length = %wu, p = %wu, orient = %wu\n", \ + rdim, cdim, len, prime, orient); + } + + nmod_poly_mat_window_clear(kernz); + } + + nmod_poly_mat_window_clear(kernz); + nmod_poly_mat_clear(pmat); + nmod_poly_mat_clear(ker); + flint_free(shift); + flint_free(cdeg); + flint_free(pivind); + flint_free(cdeg_check); + flint_free(pivind_check); + } + TEST_FUNCTION_END(state); } diff --git a/flint-extras/src/nmod_poly_mat_extra/verification.c b/flint-extras/src/nmod_poly_mat_extra/verification.c index 79c5c80c..35a20dd1 100644 --- a/flint-extras/src/nmod_poly_mat_extra/verification.c +++ b/flint-extras/src/nmod_poly_mat_extra/verification.c @@ -120,67 +120,97 @@ int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker, poly_mat_form_t form, orientation_t orient) { - if (orient == COL_UPPER || orient == COL_LOWER) + if (orient == ROW_UPPER || orient == ROW_LOWER) { - nmod_poly_mat_t ker_t; - nmod_poly_mat_t mat_t; - nmod_poly_mat_init(ker_t, ker->c, ker->r, ker->modulus); - nmod_poly_mat_transpose(ker_t, ker); - nmod_poly_mat_init(mat_t, mat->c, mat->r, mat->modulus); - nmod_poly_mat_transpose(mat_t, mat); - - int result; - if (orient == COL_UPPER) - result = nmod_poly_mat_is_kernel(ker_t, shift, mat_t, form, ROW_LOWER); - else /* orient == COL_LOWER */ - result = nmod_poly_mat_is_kernel(ker_t, shift, mat_t, form, ROW_UPPER); - - nmod_poly_mat_clear(ker_t); - nmod_poly_mat_clear(mat_t); - - return result; - } + const slong rdim = mat->r; + const slong cdim = mat->c; + const ulong prime = mat->modulus; - const slong rdim = mat->r; - const slong cdim = mat->c; - const ulong prime = mat->modulus; + nmod_poly_mat_t residual; + nmod_poly_mat_init(residual, rdim, cdim, prime); - nmod_poly_mat_t residual; - nmod_poly_mat_init(residual, rdim, cdim, prime); + int success = 1; - int success = 1; + /* check kernel has the right column dimension */ + if (ker->c != rdim) + { + printf("basis has wrong column dimension\n"); + success = 0; + } - /* check kernel has the right column dimension */ - if (ker->c != rdim) - { - printf("basis has wrong column dimension\n"); - success = 0; - } + /* check rank */ + slong rk = nmod_poly_mat_rank(mat); + if (ker->r != rdim - rk) + { + printf("number of rows does not equal nullity\n"); + success = 0; + } - /* check rank */ - slong rk = nmod_poly_mat_rank(mat); - if (ker->r != rdim - rk) - { - printf("number of rows does not equal nullity\n"); - success = 0; - } + /* check kernel has the right form */ + if (!nmod_poly_mat_is_form(ker, form, shift, orient)) + { + printf("basis does not have the required form\n"); + success = 0; + } - /* check kernel has the right form */ - if (!nmod_poly_mat_is_form(ker, form, shift, orient)) - { - printf("basis does not have the required form\n"); - success = 0; + /* compute residual, check rows of ker are in the kernel */ + nmod_poly_mat_mul(residual, ker, mat); + if (!nmod_poly_mat_is_zero(residual)) + { + printf("not all rows are in the kernel\n"); + success = 0; + } + + nmod_poly_mat_clear(residual); + + return success; } - /* compute residual, check rows of ker are in the kernel */ - nmod_poly_mat_mul(residual, ker, mat); - if (!nmod_poly_mat_is_zero(residual)) + if (orient == COL_UPPER || orient == COL_LOWER) { - printf("not all rows are in the kernel\n"); - success = 0; - } + const slong rdim = mat->r; + const slong cdim = mat->c; + const ulong prime = mat->modulus; - nmod_poly_mat_clear(residual); + nmod_poly_mat_t residual; + nmod_poly_mat_init(residual, rdim, cdim, prime); - return success; + int success = 1; + + /* check kernel has the right column dimension */ + if (ker->r != cdim) + { + printf("basis has wrong row dimension\n"); + success = 0; + } + + /* check rank */ + slong rk = nmod_poly_mat_rank(mat); + if (ker->c != cdim - rk) + { + printf("number of rows does not equal nullity\n"); + success = 0; + } + + /* check kernel has the right form */ + if (!nmod_poly_mat_is_form(ker, form, shift, orient)) + { + printf("basis does not have the required form\n"); + success = 0; + } + + /* compute residual, check rows of ker are in the kernel */ + nmod_poly_mat_mul(residual, mat, ker); + if (!nmod_poly_mat_is_zero(residual)) + { + printf("not all rows are in the kernel\n"); + success = 0; + } + + nmod_poly_mat_clear(residual); + + return success; + } + + flint_throw(FLINT_ERROR, "Exception (nmod_poly_mat_is_kernel). Requested orientation not implemented."); } diff --git a/flint-extras/src/nmod_poly_mat_forms.h b/flint-extras/src/nmod_poly_mat_forms.h index 899f46b2..3182a5e3 100644 --- a/flint-extras/src/nmod_poly_mat_forms.h +++ b/flint-extras/src/nmod_poly_mat_forms.h @@ -30,12 +30,10 @@ * */ +#include "pml.h" #include // for fmpz_mat (degree matrix) -//#include #include // fmpz_types not enough, need NMOD_POLY_MAT_INLINE -#include "pml.h" - #ifdef __cplusplus extern "C" { #endif @@ -495,7 +493,7 @@ int nmod_poly_mat_is_form(const nmod_poly_mat_t mat, return nmod_poly_mat_is_echelon(mat, orient); if (form == HERMITE) return nmod_poly_mat_is_echelon(mat, orient); - return -1; + flint_throw(FLINT_ERROR, "Exception (nmod_poly_mat_is_form). Requested form not implemented."); } //@} // doxygen group: Testing polynomial matrix forms diff --git a/flint-extras/src/nmod_poly_mat_kernel.h b/flint-extras/src/nmod_poly_mat_kernel.h index 2534e38a..30d0fdf4 100644 --- a/flint-extras/src/nmod_poly_mat_kernel.h +++ b/flint-extras/src/nmod_poly_mat_kernel.h @@ -69,18 +69,14 @@ extern "C" { * only the first `nz` entries matter in output. */ -/* general interface */ -/* `form` not really implemented yet, only shifted weak Popov */ +/** general interface */ +/* TODO `form` not really implemented yet, only shifted weak Popov */ slong nmod_poly_mat_kernel(nmod_poly_mat_t ker, slong * pivind, slong * shift, const nmod_poly_mat_t pmat, poly_mat_form_t form, orientation_t orient); -// ROW_LOWER -> direct call -// COL_UPPER -> transpose input, call, transpose output -// ROW_UPPER -> mirror rows of input, call, mirror columns+rows of output -// COL_LOWER -> /** Computes a `shift`-ordered weak Popov left kernel basis `ker` for `pmat`, * through a minimal approximant basis at sufficiently large order. @@ -112,6 +108,13 @@ slong nmod_poly_mat_kernel_zls_approx(nmod_poly_mat_t ker, /* slong * shift, */ /* const nmod_poly_mat_t pmat); */ +/* FIXME could add a "prepare pmat" function: + * A/ compute rdeg and cdeg + * ->suppress zero columns + * ->trivialize zero rows + * B/ handle case of constant matrices + * C/ go for actual work + * */ From 9ec768582737aa1a7b6091e7e494f0012d634d69 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Mon, 5 Jan 2026 14:14:45 +0100 Subject: [PATCH 11/12] add notes about choice of approximation order --- flint-extras/src/nmod_poly_mat_extra/kernel.c | 33 +++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/flint-extras/src/nmod_poly_mat_extra/kernel.c b/flint-extras/src/nmod_poly_mat_extra/kernel.c index 3ae3977c..dcc1b3a1 100644 --- a/flint-extras/src/nmod_poly_mat_extra/kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/kernel.c @@ -1157,3 +1157,36 @@ int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong * FLINT_UNUSED(deg * The maximum pivot degree also gives a bound on max(rdeg_s(ker)): pick one of * the bounds above for sum of pivot degrees, and add max(s) - min(s). */ + +/** Expected degree / order for approximation in ZLS algorithm. + * + * For the moment, the code in zls_approx, when n <= m/2, uses the approximation order + * order = 1 + max(maxdeg, 1 + floor((sum(shift) - 1) / (m - n))) + * where maxdeg is deg(pmat) and shift >= rdeg(pmat) entry-wise + * (more precisely, the code implicitly ensures min(shift - rdeg(pmat)) == 0) + * + * Note that any order > 0 is enough to guarantee that the algorithm terminates, + * but order too small may impact performance. We pick order that is large + * enough to provide the whole kernel in some generic/frequent situations. + * + * -> `1 + maxdeg` is because any lower order is too low: it means ignoring top + * degree coefficients of `pmat` + * + * -> the other term, `1 + dbound` where + * dbound = 1 + floor((sum(shift) - 1) / (m - n)) + * is because this is the expected maximum entry of the shift-row degree of a + * shift-weak Popov kernel basis P, if degrees are balanced in a generic way. + * + * Indeed: let dbound = max(rdeg_s(P)) + * . sum(rdeg_s(P)) <= sum(s) is known (Zhou-Labahn-Storjohann, 2012) + * . if we assume degrees are well balanced (which happens generically): + * then rdeg_s(P) consists of nz values all equal to dbound or dbound - 1 + * -> dbound + (nz - 1) * (dbound - 1) <= sum(s) + * -> dbound <= floor((sum(s) + (nz - 1)) / nz) + * = 1 + floor((sum(s) - 1) / nz) + * <= 1 + floor((sum(s) - 1) / (m - n)). + * + * Note also that since s >= rdeg(pmat), we have rdeg(p*pmat) <= rdeg_s(p) < order + * for any row vector p such that rdeg_s(p) <= dbound. So, approximants of pmat + * whose s-degree is <= dbound are necessarily actual kernel rows p*pmat == 0. + */ From 3422c0ea2fce2630481bb737dc38bdcc334c2894 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Mon, 5 Jan 2026 14:44:42 +0100 Subject: [PATCH 12/12] general interface: support NULL shift and pivind --- flint-extras/src/nmod_poly_mat_extra/kernel.c | 76 +++++++++++++------ .../src/nmod_poly_mat_extra/test/t-kernel.c | 11 +-- flint-extras/src/nmod_poly_mat_kernel.h | 11 ++- 3 files changed, 62 insertions(+), 36 deletions(-) diff --git a/flint-extras/src/nmod_poly_mat_extra/kernel.c b/flint-extras/src/nmod_poly_mat_extra/kernel.c index dcc1b3a1..774cd250 100644 --- a/flint-extras/src/nmod_poly_mat_extra/kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/kernel.c @@ -35,16 +35,42 @@ slong nmod_poly_mat_kernel(nmod_poly_mat_t ker, if (form > ORD_WEAK_POPOV) flint_throw(FLINT_ERROR, "Exception (nmod_poly_mat_kernel). form > ORD_WEAK_POPOV not implemented."); + slong nz = -1; + + slong * _pivind = pivind; + slong * _shift = shift; + + if (_pivind == NULL) + { + if (orient == ROW_LOWER || orient == ROW_UPPER) /* left kernel */ + _pivind = FLINT_ARRAY_ALLOC(pmat->r, slong); + else /* right kernel */ + _pivind = FLINT_ARRAY_ALLOC(pmat->c, slong); + } + + if (_shift == NULL) + { + if (orient == ROW_LOWER || orient == ROW_UPPER) /* left kernel */ + { + _shift = FLINT_ARRAY_ALLOC(pmat->r, slong); + nmod_poly_mat_row_degree(_shift, pmat, NULL); + } + else /* right kernel */ + { + _shift = FLINT_ARRAY_ALLOC(pmat->c, slong); + nmod_poly_mat_column_degree(_shift, pmat, NULL); + } + } + if (orient == ROW_LOWER) { nmod_poly_mat_t mat; nmod_poly_mat_init_set(mat, pmat); - slong nz = nmod_poly_mat_kernel_zls_approx(ker, pivind, shift, mat); + nz = nmod_poly_mat_kernel_zls_approx(ker, _pivind, _shift, mat); nmod_poly_mat_clear(mat); - return nz; } - if (orient == COL_UPPER) + else if (orient == COL_UPPER) { /* transpose input, call ROW_LOWER, transpose output */ nmod_poly_mat_t mat_t; @@ -52,38 +78,36 @@ slong nmod_poly_mat_kernel(nmod_poly_mat_t ker, nmod_poly_mat_init(mat_t, pmat->c, pmat->r, pmat->modulus); nmod_poly_mat_init(ker_t, ker->c, ker->r, ker->modulus); nmod_poly_mat_transpose(mat_t, pmat); - slong nz = nmod_poly_mat_kernel_zls_approx(ker_t, pivind, shift, mat_t); + nz = nmod_poly_mat_kernel_zls_approx(ker_t, _pivind, _shift, mat_t); nmod_poly_mat_transpose(ker, ker_t); nmod_poly_mat_clear(ker_t); nmod_poly_mat_clear(mat_t); - return nz; } - if (orient == ROW_UPPER) + else if (orient == ROW_UPPER) { /* mirror rows of input, call ROW_LOWER, mirror columns+rows of output */ nmod_poly_mat_t mat_i; nmod_poly_mat_init_set(mat_i, pmat); - nmod_poly_mat_invert_rows(mat_i, shift); - slong nz = nmod_poly_mat_kernel_zls_approx(ker, pivind, shift, mat_i); + nmod_poly_mat_invert_rows(mat_i, _shift); + nz = nmod_poly_mat_kernel_zls_approx(ker, _pivind, _shift, mat_i); nmod_poly_mat_t kernz; nmod_poly_mat_window_init(kernz, ker, 0, 0, nz, ker->c); - nmod_poly_mat_invert_rows(kernz, shift); + nmod_poly_mat_invert_rows(kernz, _shift); nmod_poly_mat_invert_columns(kernz, NULL); for (slong i = 0; i < nz/2; i++) { - slong tmp = ker->c - 1 - pivind[i]; - pivind[i] = ker->c - 1 - pivind[nz - 1 - i]; - pivind[nz - 1 - i] = tmp; + slong tmp = ker->c - 1 - _pivind[i]; + _pivind[i] = ker->c - 1 - _pivind[nz - 1 - i]; + _pivind[nz - 1 - i] = tmp; } if (nz % 2) - pivind[nz/2] = ker->c - 1 - pivind[nz/2]; + _pivind[nz/2] = ker->c - 1 - _pivind[nz/2]; nmod_poly_mat_window_clear(kernz); nmod_poly_mat_clear(mat_i); - return nz; } - if (orient == COL_LOWER) + else if (orient == COL_LOWER) { /* transpose input, call ROW_UPPER, transpose output */ nmod_poly_mat_t mat_it; @@ -92,29 +116,33 @@ slong nmod_poly_mat_kernel(nmod_poly_mat_t ker, nmod_poly_mat_init(ker_it, ker->c, ker->r, ker->modulus); nmod_poly_mat_transpose(mat_it, pmat); - nmod_poly_mat_invert_rows(mat_it, shift); - slong nz = nmod_poly_mat_kernel_zls_approx(ker_it, pivind, shift, mat_it); + nmod_poly_mat_invert_rows(mat_it, _shift); + nz = nmod_poly_mat_kernel_zls_approx(ker_it, _pivind, _shift, mat_it); nmod_poly_mat_t kernz; nmod_poly_mat_window_init(kernz, ker_it, 0, 0, nz, ker_it->c); - nmod_poly_mat_invert_rows(kernz, shift); + nmod_poly_mat_invert_rows(kernz, _shift); nmod_poly_mat_invert_columns(kernz, NULL); for (slong i = 0; i < nz/2; i++) { - slong tmp = ker_it->c - 1 - pivind[i]; - pivind[i] = ker_it->c - 1 - pivind[nz - 1 - i]; - pivind[nz - 1 - i] = tmp; + slong tmp = ker_it->c - 1 - _pivind[i]; + _pivind[i] = ker_it->c - 1 - _pivind[nz - 1 - i]; + _pivind[nz - 1 - i] = tmp; } if (nz % 2) - pivind[nz/2] = ker_it->c - 1 - pivind[nz/2]; + _pivind[nz/2] = ker_it->c - 1 - _pivind[nz/2]; nmod_poly_mat_window_clear(kernz); nmod_poly_mat_transpose(ker, ker_it); nmod_poly_mat_clear(ker_it); nmod_poly_mat_clear(mat_it); - return nz; } - flint_throw(FLINT_ERROR, "Exception (nmod_poly_mat_kernel). Requested orientation not implemented."); + if (pivind == NULL) + flint_free(_pivind); + if (shift == NULL) + flint_free(_shift); + + return nz; } diff --git a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c index 5565c7df..04628375 100644 --- a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c @@ -143,10 +143,6 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state) if (!result) { - flint_printf("true rdeg == %{slong*}\n", rdeg, nz); - flint_printf("got rdeg == %{slong*}\n", rdeg_check, nz); - flint_printf("true pivind == %{slong*}\n", pivind, nz); - flint_printf("got pivind == %{slong*}\n", pivind_check, nz); TEST_FUNCTION_FAIL("(kernel, rdeg/pivind) -- rdim = %wu, cdim = %wu, length = %wu, p = %wu, orient = %wu\n", \ rdim, cdim, len, prime, orient); } @@ -194,8 +190,7 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state) /* kernel interface */ { poly_mat_form_t form = ORD_WEAK_POPOV; - /* orientation_t orient = (n_randint(state, 2)) ? COL_LOWER : COL_UPPER; */ - orientation_t orient = COL_UPPER; + orientation_t orient = (n_randint(state, 2)) ? COL_LOWER : COL_UPPER; for (slong j = 0; j < cdim; j++) cdeg[j] = shift[j]; @@ -219,10 +214,6 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state) if (!result) { - flint_printf("true cdeg == %{slong*}\n", cdeg, nz); - flint_printf("got cdeg == %{slong*}\n", cdeg_check, nz); - flint_printf("true pivind == %{slong*}\n", pivind, nz); - flint_printf("got pivind == %{slong*}\n", pivind_check, nz); TEST_FUNCTION_FAIL("(kernel, rdeg/pivind) -- rdim = %wu, cdim = %wu, length = %wu, p = %wu, orient = %wu\n", \ rdim, cdim, len, prime, orient); } diff --git a/flint-extras/src/nmod_poly_mat_kernel.h b/flint-extras/src/nmod_poly_mat_kernel.h index 30d0fdf4..2598dd2e 100644 --- a/flint-extras/src/nmod_poly_mat_kernel.h +++ b/flint-extras/src/nmod_poly_mat_kernel.h @@ -69,8 +69,15 @@ extern "C" { * only the first `nz` entries matter in output. */ -/** general interface */ -/* TODO `form` not really implemented yet, only shifted weak Popov */ +/** general interface + * -> all four orientations are supported (ROW_LOWER|ROW_UPPER yield a left + * kernel basis while COL_LOWER|COL_UPPER yield a right kernel basis) + * -> pivind may be NULL if that information is not wanted + * -> shift may be NULL in input, in which case it is left NULL in output, and + * the basis is computed in weak Popov form with respect to row degrees of + * `pmat` (if left kernel) or column degrees of `pmat` (if right kernel) + */ +/* TODO `form` not implemented yet, only shifted weak Popov */ slong nmod_poly_mat_kernel(nmod_poly_mat_t ker, slong * pivind, slong * shift,