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_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 6f99d5d7..400367d3 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,
@@ -243,82 +242,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)
@@ -378,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.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..774cd250 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,11 +11,469 @@
.
*/
+#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"
+#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.");
+
+ 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);
+ nz = nmod_poly_mat_kernel_zls_approx(ker, _pivind, _shift, mat);
+ nmod_poly_mat_clear(mat);
+ }
+
+ else 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);
+ 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);
+ }
+
+ 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);
+ 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);
+ }
+
+ else 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);
+ 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);
+ }
+
+ if (pivind == NULL)
+ flint_free(_pivind);
+ if (shift == NULL)
+ flint_free(_shift);
+
+ return nz;
+}
+
+
+
+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;
+ }
+
+ /* 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)
+ 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);
+ nmod_mat_clear(coeff0);
+ if (r == m)
+ 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_pmbasis(ker, pivind, pmat, order);
+
+ /* 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++)
+ {
+ if (pivind[i] - shift[i] + amp < order - d)
+ {
+ 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;
+ }
+ }
+
+ return nz;
+}
+
+/* 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,
+ nmod_poly_mat_t pmat)
+{
+ const slong m = pmat->r;
+ const slong n = pmat->c;
+
+ /* 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)))
+ {
+ nmod_poly_mat_one(ker);
+ for (slong i = 0; i < m; i++)
+ pivind[i] = i;
+ return m;
+ }
+
+ /* one (nonzero) row vector */
+ if (m == 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)
+ {
+ /* 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);
+ 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);
+ 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 1 <= n <= m/2 */
+
+ /* find rdeg, maxdeg, and diff_shift = min(shift - max(0, rdeg(pmat))) */
+ slong * buf = FLINT_ARRAY_ALLOC(m, slong);
+ nmod_poly_mat_row_degree(buf, pmat, NULL);
+ slong maxdeg = buf[0];
+ slong diff_shift = shift[0] - FLINT_MAX(0, buf[0]);
+ for (slong i = 1; i < m; i++)
+ {
+ slong d = buf[i];
+ slong diff = shift[i] - FLINT_MAX(0, d);
+ if (d > maxdeg)
+ maxdeg = d;
+ if (diff_shift > diff)
+ diff_shift = diff;
+ }
+
+ /* early exit: pmat == 0 => kernel is identity */
+ if (maxdeg == -1)
+ {
+ 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;
+ 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)
+ {
+ pivind[nz] = i;
+ 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);
+
+ /* early exit */
+ 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 */
+ 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 */
+ /* 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);
+
+ /* 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 */
+ 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);
+ 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));
+
+ /* 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++;
+ }
+ }
+
+ nmod_poly_mat_permute_rows(approx, buf, NULL);
+ 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++)
+ 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);
+
+ return nz + nz2;
+}
+
+
+
+
/**
*
@@ -600,7 +1059,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)
{
@@ -690,3 +1149,72 @@ 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).
+ */
+
+/** 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.
+ */
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..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
@@ -29,7 +41,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; \
@@ -52,7 +64,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,13 +74,89 @@ 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); \
}
-TIME_KER(kernel_zls)
+#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 = 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; \
+ nmod_poly_mat_init(ker, rdim, rdim, n); \
+ \
+ double FLINT_SET_BUT_UNUSED(tcpu), twall; \
+ \
+ TIMEIT_START; \
+ 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); \
+ \
+ flint_printf("%.2e", twall); \
+ \
+ flint_free(pivind); \
+ flint_free(shift); \
+ nmod_poly_mat_clear(pmat); \
+ 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)
/*-------------------------*/
/* main */
@@ -80,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};
@@ -97,10 +176,13 @@ int main(int argc, char ** argv)
/* TODO shift type */
// bench functions
- const slong nfuns = 1;
+ const slong nfuns = 4;
typedef void (*timefun) (time_args, flint_rand_t);
const timefun funs[] = {
- time_kernel_zls, // 0
+ time_kernel_via_approx, // 0
+ time_kernel_zls_approx, // 1
+ time_nullspace, // 2
+ time_kernel_zls, // 3
};
// TODO
@@ -120,22 +202,24 @@ int main(int argc, char ** argv)
//#endif
const char * description[] = {
- "#0 --> kernel (general interface) ",
- "#1 --> .... TODO ",
+ "#0 --> via approx ",
+ "#1 --> ZLS via approx ",
+ "#2 --> FLINT's nullspace ",
+ "#3 --> ZLS via approx(bis) ",
};
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]);
@@ -148,50 +232,46 @@ 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");
- 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/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/main.c b/flint-extras/src/nmod_poly_mat_extra/test/main.c
index e89828cd..0c8fe9d4 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.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),
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.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c
new file mode 100644
index 00000000..04628375
--- /dev/null
+++ b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c
@@ -0,0 +1,235 @@
+/*
+ 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 k, result;
+
+ /* 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, 75);
+
+ 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);
+
+ /* kernel_via_approx */
+ {
+ 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 i = 0; i < nz; i++)
+ if (rdeg_check[i] != rdeg[i] || pivind_check[i] != pivind[i])
+ 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);
+ }
+
+ /* kernel_zls_approx */
+ {
+ 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 i = 0; i < nz; i++)
+ if (rdeg_check[i] != rdeg[i] || pivind_check[i] != pivind[i])
+ 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);
+ }
+
+ /* 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)
+ {
+ 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(rdeg);
+ flint_free(pivind);
+ flint_free(rdeg_check);
+ 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;
+
+ 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)
+ {
+ 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/test/t-kernel_zls.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls.c
index 94a6edbe..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,68 +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, Mt, rdeg, 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, 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, 20);
+ ulong cdim = rdim + 1 + n_randint(state, 20);
+ ulong deg = n_randint(state, 100);
ulong prime = n_randprime(state, nbits, 1);
@@ -101,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);
@@ -128,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);
@@ -155,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/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;
-}
diff --git a/flint-extras/src/nmod_poly_mat_extra/verification.c b/flint-extras/src/nmod_poly_mat_extra/verification.c
index 1ce2fe37..35a20dd1 100644
--- a/flint-extras/src/nmod_poly_mat_extra/verification.c
+++ b/flint-extras/src/nmod_poly_mat_extra/verification.c
@@ -113,62 +113,104 @@ 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,
- const nmod_poly_mat_t pmat,
const slong * shift,
+ 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 slong kdim = ker->r;
- const ulong prime = pmat->modulus;
+ if (orient == ROW_UPPER || orient == ROW_LOWER)
+ {
+ 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, kdim, 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 kernel is shifted reduced */
- if (!nmod_poly_mat_is_ordered_weak_popov(ker, 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))
+ /* check rank */
+ slong rk = nmod_poly_mat_rank(mat);
+ if (ker->r != rdim - rk)
{
- printf("basis is not shifted-reduced\n");
+ printf("number of rows does not equal nullity\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, ker, pmat);
- if (!nmod_poly_mat_is_zero(residual))
- {
- printf("not all rows are in the kernel\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;
}
- /* check rank */
- slong rk = nmod_poly_mat_rank(pmat);
- if (kdim != rdim - rk)
+ if (orient == COL_UPPER || orient == COL_LOWER)
{
- printf("number of rows does not equal nullity\n");
- success = 0;
- }
+ const slong rdim = mat->r;
+ const slong cdim = mat->c;
+ const ulong prime = mat->modulus;
- /* !! TODO check generation !! */
+ nmod_poly_mat_t residual;
+ nmod_poly_mat_init(residual, rdim, cdim, prime);
- nmod_poly_mat_clear(residual);
+ int success = 1;
- return success;
+ /* 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 113edf7c..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
@@ -52,10 +50,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 +465,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);
+ flint_throw(FLINT_ERROR, "Exception (nmod_poly_mat_is_form). Requested form not implemented.");
+}
+
//@} // doxygen group: Testing polynomial matrix forms
@@ -831,6 +852,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 20cb54b6..2598dd2e 100644
--- a/flint-extras/src/nmod_poly_mat_kernel.h
+++ b/flint-extras/src/nmod_poly_mat_kernel.h
@@ -13,12 +13,141 @@
#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" {
#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
+ * \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:
+ *
+ * \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`)
+ * \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)
+ *
+ * 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`).
+ *
+ * `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
+ * -> 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,
+ const nmod_poly_mat_t pmat,
+ poly_mat_form_t form,
+ orientation_t orient);
+
+/** Computes a `shift`-ordered weak Popov left kernel basis `ker` for `pmat`,
+ * 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); */
+
+/* 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); */
+
+/* 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
+ * */
+
+
+
+/** 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] shift shift
+ * \param[in] pmat polynomial matrix
+ * \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 WARNING! for the moment, does not really check generation! */
+int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker,
+ 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
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