diff --git a/flint-extras/src/nmod_poly_mat_extra.h b/flint-extras/src/nmod_poly_mat_extra.h index 70f5d99c..b6d675ae 100644 --- a/flint-extras/src/nmod_poly_mat_extra.h +++ b/flint-extras/src/nmod_poly_mat_extra.h @@ -63,6 +63,7 @@ #include "nmod_poly_mat_kernel.h" + // TODO remove once using flint's comp instead NMOD_POLY_MAT_INLINE void apply_perm_to_vector(slong *res, const slong *initial_vect, @@ -72,9 +73,7 @@ apply_perm_to_vector(slong *res, const slong *initial_vect, res[perm[i]] = initial_vect[i]; } -/** \todo doc - * \todo move in suitable header - **/ +/* TODO move in suitable header */ void nmod_poly_mat_det_iter(nmod_poly_t det, nmod_poly_mat_t mat); // TODO implem + doc @@ -86,6 +85,24 @@ void nmod_poly_mat_det_iter(nmod_poly_t det, nmod_poly_mat_t mat); * \todo triangularization / Hermite */ + +/***************************** +* 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 -*- */ diff --git a/flint-extras/src/nmod_poly_mat_extra/approximant_basis.c b/flint-extras/src/nmod_poly_mat_extra/approximant_basis.c new file mode 100644 index 00000000..9c216269 --- /dev/null +++ b/flint-extras/src/nmod_poly_mat_extra/approximant_basis.c @@ -0,0 +1,61 @@ +/* + 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 "nmod_poly_mat_multiply.h" // for middle product +#include "nmod_poly_mat_approximant.h" +#include "nmod_mat_poly.h" +#include "nmod_poly_mat_utils.h" + +void nmod_poly_mat_mbasis(nmod_poly_mat_t appbas, + slong * shift, + const nmod_poly_mat_t pmat, + ulong order) +{ + nmod_mat_poly_t app, matp; + nmod_mat_poly_init(matp, pmat->r, pmat->c, pmat->modulus); + nmod_mat_poly_set_trunc_from_poly_mat(matp, pmat, order); + nmod_mat_poly_init(app, pmat->r, pmat->r, pmat->modulus); + nmod_mat_poly_mbasis(app, shift, matp, order); + nmod_poly_mat_set_from_mat_poly(appbas, app); + nmod_mat_poly_clear(matp); + nmod_mat_poly_clear(app); +} + +void nmod_poly_mat_pmbasis(nmod_poly_mat_t appbas, + slong * shift, + const nmod_poly_mat_t pmat, + slong order) +{ + if (order <= PMBASIS_THRES) + { + nmod_poly_mat_mbasis(appbas, shift, pmat, order); + return; + } + + const long order1 = order>>1; + const long order2 = order - order1; + nmod_poly_mat_t appbas2, residual; + + nmod_poly_mat_init(appbas2, pmat->r, pmat->r, pmat->modulus); + nmod_poly_mat_init(residual, pmat->r, pmat->c, pmat->modulus); + + nmod_poly_mat_pmbasis(appbas, shift, pmat, order1); + + nmod_poly_mat_middle_product_naive(residual, appbas, pmat, order1, order2-1); + + nmod_poly_mat_pmbasis(appbas2, shift, residual, order2); + + nmod_poly_mat_mul(appbas, appbas2, appbas); + + nmod_poly_mat_clear(appbas2); + nmod_poly_mat_clear(residual); +} diff --git a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_pivots.c b/flint-extras/src/nmod_poly_mat_extra/degrees_pivots_leadingmatrix.c similarity index 55% rename from flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_pivots.c rename to flint-extras/src/nmod_poly_mat_extra/degrees_pivots_leadingmatrix.c index 7fd8a858..f9621b15 100644 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_pivots.c +++ b/flint-extras/src/nmod_poly_mat_extra/degrees_pivots_leadingmatrix.c @@ -1,11 +1,142 @@ +/* + 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 "nmod_poly_mat_forms.h" +/*------------------------------------------------------------*/ +/* row degree - column degree - degree matrix */ +/*------------------------------------------------------------*/ + +void nmod_poly_mat_row_degree(slong *rdeg, const nmod_poly_mat_t mat, const slong *shift) +{ + if (shift == NULL) + { + slong max, d; + for (slong i = 0; i < mat->r; i++) + { + max = -1; + for (slong j = 0; j < mat->c; j++) + { + d = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)); + if (max < d) + max = d; + } + rdeg[i] = max; + } + } + else + { + slong max, d; + slong min_shift = (mat->c > 0) ? shift[0] : 0; + + // find minimum of shift + for (slong j = 0; j < mat->c; ++j) + if (shift[j] < min_shift) + min_shift = shift[j]; + + for (slong i = 0; i < mat->r; i++) + { + max = min_shift-1; // zero rows will have this as rdeg + for (slong j = 0; j < mat->c; j++) + { + d = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)) + shift[j]; + // if new maximum at a nonzero entry, update + if (shift[j] <= d && max < d) + max = d; + } + rdeg[i] = max; + } + } +} + +void nmod_poly_mat_column_degree(slong *cdeg, const nmod_poly_mat_t mat, const slong *shift) +{ + if (shift == NULL) + { + slong max, d; + for (slong j = 0; j < mat->c; j++) + { + max = -1; + for (slong i = 0; i < mat->r; i++) + { + d = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)); + if (max < d) + max = d; + } + cdeg[j] = max; + } + } + else + { + slong max, d; + slong min_shift = (mat->r > 0) ? shift[0] : 0; -/********************************************************************** -* shifted pivot profile, vector * -**********************************************************************/ + // find minimum of shift + for (slong i = 0; i < mat->r; ++i) + if (shift[i] < min_shift) + min_shift = shift[i]; + + for (slong j = 0; j < mat->c; j++) + { + max = min_shift-1; + for (slong i = 0; i < mat->r; i++) + { + d = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)) + shift[i]; + // if new maximum at a nonzero entry, update + if (shift[i] <= d && max < d) + max = d; + } + cdeg[j] = max; + } + } +} + +void nmod_poly_mat_degree_matrix(fmpz_mat_t dmat, const nmod_poly_mat_t mat) +{ + for(slong i = 0; i < mat->r; i++) + for(slong j = 0; j < mat->c; j++) + *fmpz_mat_entry(dmat, i, j) = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)); +} + +void nmod_poly_mat_degree_matrix_shifted(fmpz_mat_t dmat, + const nmod_poly_mat_t mat, + const slong * shift, + orientation_t orient) +{ + if (shift) + { + if (orient == ROW_LOWER || orient == ROW_UPPER) + for(slong i = 0; i < mat->r; i++) + for(slong j = 0; j < mat->c; j++) + *fmpz_mat_entry(dmat, i, j) = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)) + shift[j]; + else // orient == COL_* + for(slong i = 0; i < mat->r; i++) + for(slong j = 0; j < mat->c; j++) + *fmpz_mat_entry(dmat, i, j) = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)) + shift[i]; + } + else + for(slong i = 0; i < mat->r; i++) + for(slong j = 0; j < mat->c; j++) + *fmpz_mat_entry(dmat, i, j) = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)); +} + + +/*------------------------------------------------------------*/ +/* shifted pivot profile (index + degree) */ +/*------------------------------------------------------------*/ void _nmod_poly_vec_pivot_profile(slong * pivind, slong * pivdeg, @@ -98,10 +229,6 @@ void _nmod_poly_vec_pivot_profile(slong * pivind, } } -/********************************************************************** -* shifted pivot index * -**********************************************************************/ - void nmod_poly_mat_pivot_index(slong *pivind, const nmod_poly_mat_t mat, const slong * shift, @@ -128,10 +255,6 @@ void nmod_poly_mat_pivot_index(slong *pivind, } } -/********************************************************************** -* shifted pivot profile * -**********************************************************************/ - void nmod_poly_mat_pivot_profile(slong * pivind, slong * pivdeg, const nmod_poly_mat_t mat, @@ -155,11 +278,6 @@ void nmod_poly_mat_pivot_profile(slong * pivind, } } - -/********************************************************************** -* echelon pivot index * -**********************************************************************/ - void nmod_poly_mat_echelon_pivot_index(slong * pivind, const nmod_poly_mat_t mat, orientation_t orient) { switch (orient) @@ -211,10 +329,6 @@ void nmod_poly_mat_echelon_pivot_index(slong * pivind, const nmod_poly_mat_t mat } } -/********************************************************************** -* echelon pivot profile * -**********************************************************************/ - void nmod_poly_mat_echelon_pivot_profile(slong * pivind, slong * pivdeg, const nmod_poly_mat_t mat, orientation_t orient) { switch (orient) @@ -271,5 +385,76 @@ void nmod_poly_mat_echelon_pivot_profile(slong * pivind, slong * pivdeg, const n } } -/* -*- 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 + +/*------------------------------------------------------------*/ +/* leading matrix */ +/*------------------------------------------------------------*/ + +void nmod_poly_mat_leading_matrix(nmod_mat_t lmat, + const nmod_poly_mat_t mat, + const slong * shift, + orientation_t orient) +{ + if (orient == ROW_LOWER || orient == ROW_UPPER) + { + // find row degrees + slong rdeg[mat->r]; + nmod_poly_mat_row_degree(rdeg, mat, shift); + + // deduce leading matrix + if (shift == NULL) + { + for (slong i = 0; i < mat->r; i++) + if (rdeg[i] >= 0) + for (slong j = 0; j < mat->c; j++) + nmod_mat_set_entry(lmat, i, j, + nmod_poly_get_coeff_ui(nmod_poly_mat_entry(mat, i, j), + rdeg[i])); + else + for (slong j = 0; j < mat->c; j++) + nmod_mat_set_entry(lmat, i, j, 0); + } + else + { + for (slong i = 0; i < mat->r; i++) + for (slong j = 0; j < mat->c; j++) + if (rdeg[i] >= shift[j]) + nmod_mat_set_entry(lmat, i, j, + nmod_poly_get_coeff_ui(nmod_poly_mat_entry(mat, i, j), + rdeg[i] - shift[j])); + else + nmod_mat_set_entry(lmat, i, j, 0); + } + } + else // orient == COL_* + { + // find column degrees + slong cdeg[mat->c]; + nmod_poly_mat_column_degree(cdeg, mat, shift); + + // deduce leading matrix + if (shift == NULL) + { + for (slong j = 0; j < mat->c; j++) + if (cdeg[j] >= 0) + for (slong i = 0; i < mat->r; i++) + nmod_mat_set_entry(lmat, i, j, + nmod_poly_get_coeff_ui(nmod_poly_mat_entry(mat, i, j), + cdeg[j])); + else + for (slong i = 0; i < mat->r; i++) + nmod_mat_set_entry(lmat, i, j, 0); + } + else + { + for (slong j = 0; j < mat->c; j++) + for (slong i = 0; i < mat->r; i++) + if (cdeg[j] >= shift[i]) + nmod_mat_set_entry(lmat, i, j, + nmod_poly_get_coeff_ui(nmod_poly_mat_entry(mat, i, j), + cdeg[j] - shift[i])); + else + nmod_mat_set_entry(lmat, i, j, 0); + } + } +} diff --git a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_kernel.c b/flint-extras/src/nmod_poly_mat_extra/kernel.c similarity index 77% rename from flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_kernel.c rename to flint-extras/src/nmod_poly_mat_extra/kernel.c index 084aee45..f46208d3 100644 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/kernel.c @@ -1,36 +1,47 @@ -#include +/* + 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 "nmod_poly_mat_extra.h" #include "nmod_poly_mat_kernel.h" /** - * + * * The content of the file has not yet been tested very extensively - * - */ + * + */ /** - * - * In place, M is m x n, shifted column degrees - * permute the columns in nondecreasing order - * - * Todo/to see: zero is set to degree 0 (not -1), for specific use + * + * In place, M is m x n, shifted column degrees + * permute the columns in nondecreasing order + * + * Todo/to see: zero is set to degree 0 (not -1), for specific use * in ZLS algorithm for the kernel - * + * * Input: * M m x n * ishift[m] - * + * * Output: - * M is modified in place - * perm[n],initialized outside, the carried out permutation + * M is modified in place + * perm[n],initialized outside, the carried out permutation * sdeg[n] initialized outside - * + * */ - -// Dummy function for qsort +// Dummy function for qsort int cmp(const void *a, const void *b) { const slong *x = a; @@ -41,7 +52,6 @@ return 0; } - void _nmod_poly_mat_sort_permute_columns_zls(nmod_poly_mat_t M, slong *sdeg, \ slong *perm, const slong *ishift) { @@ -51,7 +61,7 @@ void _nmod_poly_mat_sort_permute_columns_zls(nmod_poly_mat_t M, slong *sdeg, \ nmod_poly_mat_column_degree(sdeg, M, ishift); - for (j=0; j= 2, for the order of the order bases - * kappa * s instead of 3 *s in ZLS - * + * "with entries arranged in non-decreasing order and bounding the + * corresponding column degrees of A." + * kappa, a double >= 2, for the order of the order bases + * kappa * s instead of 3 *s in ZLS + * * Output: - * returns the dimension w of the kernel, which may be zero - * N, n x w, is initialized here only when w > 0, + * returns the dimension w of the kernel, which may be zero + * N, n x w, is initialized here only when w > 0, * should be freed outside in that case only, - * gives a minimal basis of the right kernel + * gives a minimal basis of the right kernel * degN[n], initialized outside, its first w entries are concerned, - * they are the ishift shift degrees of the kernel basis - * + * they are the ishift shift degrees of the kernel basis + * */ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat_t A, \ @@ -104,13 +114,13 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat slong min_mn; slong rho=0; - // Tuning the order of the subsequent approximant + // Tuning the order of the subsequent approximant // ---------------------------------------------- - long shift[n]; // temporary variable + long shift[n]; // temporary variable - // In order to sort for computing the order of the approximant, for the test m=1 to be ok + // In order to sort for computing the order of the approximant, for the test m=1 to be ok for (i=0; i n %ld %ld",m,n); // To see, rectangular case - } + for (i=0; i n %ld %ld",m,n); // To see, rectangular case + } slong s; - s = ceil((double) rho/min_mn); + s = ceil((double) rho/min_mn); - // Transposed A for the approximant PT on the left - // PT will then be use without transposing + // Transposed A for the approximant PT on the left + // PT will then be use without transposing nmod_poly_mat_t AT; nmod_poly_mat_init(AT, n, m, A->modulus); @@ -147,33 +157,33 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat nmod_poly_mat_t PT; nmod_poly_mat_init(PT, n, n, A->modulus); - // Approximant PT - // shift is modified in place + // Approximant PT + // shift is modified in place // -------------------------- for (i=0; imodulus); nmod_poly_mat_mul(RT,PT,AT); - + slong cdeg[n]; nmod_poly_mat_row_degree(cdeg, RT, NULL); nmod_poly_mat_clear(AT); - // Zero residue matrix P1 n x n1 + // Zero residue matrix P1 n x n1 // nonzero residues matrix P2 n x n2 //----------------------------------- @@ -181,7 +191,7 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat slong n2; for (j=0; j 0 and m > 1 - // one can proceed to the divide and conquer from the rows + // one can proceed to the divide and conquer from the rows // of the nonzero residue P2, and of A.P2 // -------------------------------------------------------- @@ -259,15 +269,15 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat slong degP2[n2]; - // perm will be used also later for G from RT + // perm will be used also later for G from RT _nmod_poly_mat_sort_permute_columns_zls(P2,degP2,perm,ishift); for (i = 0; i < n2; i++) { - degP2[i]=degP2[i]-ks; // used below for the recursive calls + degP2[i]=degP2[i]-ks; // used below for the recursive calls } - //+++++++++ OLD + //+++++++++ OLD // nmod_poly_mat_t G; // nmod_poly_mat_init(G, m, n2, A->modulus); @@ -279,7 +289,7 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat // nmod_poly_mat_shift_right(TT,G,ks); - + // slong new_m=floor((double) m/2); @@ -302,7 +312,7 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat // } - //++++++++ END OLD + //++++++++ END OLD //++++++++++ NEW ++++++++++++++ @@ -310,7 +320,7 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat nmod_poly_mat_init(TT, m, n2, A->modulus); // We extract G (temporary TT) from RT as we have been extracting P2 from PT - // and apply above ordering of P2 + // and apply above ordering of P2 k=0; for (j = 0; j < n; j++) { @@ -330,7 +340,7 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat nmod_poly_mat_init(G, m, n2, A->modulus); nmod_poly_mat_shift_right(G,TT,ks); - + // We split G for the recursive call // --------------------------------- @@ -359,7 +369,7 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat nmod_poly_mat_clear(G); //++++++++++++++ END NEW ++++++++++++++++++ - // Recursive calls + // Recursive calls // --------------- nmod_poly_mat_t N1; @@ -369,11 +379,11 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat slong c2=0; - c1=nmod_poly_mat_zls_sorted(N1, degN, G1, degP2, kappa); + c1=nmod_poly_mat_zls_sorted(N1, degN, G1, degP2, kappa); if (c1 != 0) { - + nmod_poly_mat_t G3; nmod_poly_mat_init(G3, m-new_m, c1, A->modulus); @@ -383,7 +393,7 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat shift[i]=degN[i]; } - c2=nmod_poly_mat_zls_sorted(N2, degN, G3, shift, kappa); + c2=nmod_poly_mat_zls_sorted(N2, degN, G3, shift, kappa); nmod_poly_mat_clear(G3); } @@ -392,12 +402,12 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat nmod_poly_mat_clear(G1); nmod_poly_mat_clear(G2); - // the recursive calls did not provide anything more + // the recursive calls did not provide anything more if ((c1==0) || (c2==0)){ if (n1==0) { - return 0; + return 0; } else { nmod_poly_mat_init_set(N,P1); @@ -425,15 +435,15 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat nmod_poly_mat_clear(N1); nmod_poly_mat_clear(N2); - nmod_poly_mat_clear(Q1); + nmod_poly_mat_clear(Q1); nmod_poly_mat_clear(P2); if (n1 ==0) { - nmod_poly_mat_init_set(N,Q); // We should not need to copy + nmod_poly_mat_init_set(N,Q); // We should not need to copy nmod_poly_mat_column_degree(degN, Q, ishift); - nmod_poly_mat_clear(Q); + nmod_poly_mat_clear(Q); return c2; } @@ -452,7 +462,7 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat } } - slong odeg[n1+c2]; + slong odeg[n1+c2]; nmod_poly_mat_column_degree(odeg, P1, ishift); @@ -466,54 +476,54 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat degN[i+n1]=odeg[i]; } - nmod_poly_mat_clear(P1); - nmod_poly_mat_clear(Q); - + nmod_poly_mat_clear(P1); + nmod_poly_mat_clear(Q); + return n1+c2; } - return 0; + return 0; } /** - * - * Right shifted kernel of a polynomial matrix, assuming that the columns has been - * sorted by shifted degree - * + * + * Right shifted kernel of a polynomial matrix, assuming that the columns has been + * sorted by shifted degree + * * Algorithm of Wei Zhou, George Labahn, and Arne Storjohann * "Computing Minimal Nullspace Bases" * ISSAC 2012, https://dl.acm.org/doi/abs/10.1145/2442829.2442881 - * - * Calls nmod_poly_mat_zls_sorted after an initial sorting - * - * TODO/TO SEE: - * - * Input: - * iA in m x n - * ishift[n], NULL (the degrees are computed) or initialized outside, - * the shift for the kernel + * + * Calls nmod_poly_mat_zls_sorted after an initial sorting + * + * TODO/TO SEE: + * + * Input: + * iA in m x n + * ishift[n], NULL (the degrees are computed) or initialized outside, + * the shift for the kernel * values should be at least 0 (even for zero columns in A - * "with entries arranged in non-decreasing order and bounding the - * corresponding column degrees of A." - * kappa, a double >= 2, for the order of the order bases - * kappa * s instead of 3 *s in ZLS + * "with entries arranged in non-decreasing order and bounding the + * corresponding column degrees of A." + * kappa, a double >= 2, for the order of the order bases + * kappa * s instead of 3 *s in ZLS * * Output: - * returns the dimension w of the kernel, which may be zero - * N, is initialized n x n outside - * its first w columns give a minimal basis of the kernel + * returns the dimension w of the kernel, which may be zero + * N, is initialized n x n outside + * its first w columns give a minimal basis of the kernel * degN[n], initialized outside, its first w entries are concerned, - * they are the ishift shifted degrees of the kernel basis - * + * they are the ishift shifted degrees of the kernel basis + * */ /** - * Calls nmod_poly_mat_zls_sorted after an initial sorting + * Calls nmod_poly_mat_zls_sorted after an initial sorting */ -int nmod_poly_mat_kernel(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat_t iA, \ +int nmod_poly_mat_kernel_zls(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat_t iA, \ const slong *ishift, const double kappa) { @@ -528,26 +538,26 @@ int nmod_poly_mat_kernel(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat_t i slong * perm = flint_malloc(n * sizeof(slong)); slong sdeg[n]; - // No input shift, simply the column degrees, then ordered + // No input shift, simply the column degrees, then ordered // a permutation is carried out in place // ------------------------------------------------------- if (ishift == NULL) { nmod_poly_mat_column_degree(sdeg, A, NULL); - for (j=0; j 0 - return nz; + return nz; } flint_free(perm); nmod_poly_mat_clear(A); - return 0; - + return 0; + } /** - * Experimental, should not be really considered + * Experimental, should not be really considered * */ @@ -605,11 +615,11 @@ int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong *degN, const nmod_ slong degmax=deg[0]; for (i=1; im) + if (n >m) min_nm=m; nmod_poly_mat_t AT; @@ -620,8 +630,8 @@ int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong *degN, const nmod_ nmod_poly_mat_init(PT, n, n, A->modulus); - // Approximant PT - // shift is modified in place + // Approximant PT + // shift is modified in place // -------------------------- slong shift[n]; @@ -629,18 +639,18 @@ int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong *degN, const nmod_ { shift[i]=ishift[i]; } - + nmod_poly_mat_pmbasis(PT, shift, AT, degmax*(min_nm +1)+1); // Looking for zero residues and non zero residues - // the global residue is reused later + // the global residue is reused later // ----------------------------------------------- nmod_poly_mat_t RT; nmod_poly_mat_init(RT, n, m, A->modulus); nmod_poly_mat_mul(RT,PT,AT); - + slong cdeg[n]; nmod_poly_mat_row_degree(cdeg, RT, NULL); @@ -648,7 +658,7 @@ int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong *degN, const nmod_ nmod_poly_mat_clear(RT); - // Zero residue matrix P1 + // Zero residue matrix P1 //----------------------- slong n1=0; @@ -663,7 +673,7 @@ int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong *degN, const nmod_ if (n1 > 0) { k=0; - + for (j=0; js,f0,{0,g0,(0,\:0,t0,+0,=s diff --git a/flint-extras/src/nmod_poly_mat_extra/middle_product_geometric.c b/flint-extras/src/nmod_poly_mat_extra/middle_product.c similarity index 86% rename from flint-extras/src/nmod_poly_mat_extra/middle_product_geometric.c rename to flint-extras/src/nmod_poly_mat_extra/middle_product.c index 6d1abe27..ec7a460a 100644 --- a/flint-extras/src/nmod_poly_mat_extra/middle_product_geometric.c +++ b/flint-extras/src/nmod_poly_mat_extra/middle_product.c @@ -1,3 +1,15 @@ +/* + Copyright (C) 2025 Vincent Neiger, Éric Schost + + 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 @@ -5,6 +17,20 @@ #include "nmod_extra.h" #include "nmod_poly_mat_multiply.h" +/** Middle product for polynomial matrices + * sets C = ((A * B) div x^dA) mod x^(dB+1), assuming deg(A) <= dA and deg(B) <= dA + dB + * output can alias input + * naive implementation (multiply, shift, truncate) + */ +void nmod_poly_mat_middle_product_naive(nmod_poly_mat_t C, const nmod_poly_mat_t A, const nmod_poly_mat_t B, + const ulong dA, const ulong dB) +{ + nmod_poly_mat_mul(C, A, B); + nmod_poly_mat_shift_right(C, C, dA); + nmod_poly_mat_truncate(C, dB + 1); +} + + /** Middle product for polynomial matrices * sets C = ((A * B) div x^dA) mod x^(dB+1) * output can alias input diff --git a/flint-extras/src/nmod_poly_mat_extra/mul_geometric.c b/flint-extras/src/nmod_poly_mat_extra/mul_geometric.c index 31f73c9a..607cd515 100644 --- a/flint-extras/src/nmod_poly_mat_extra/mul_geometric.c +++ b/flint-extras/src/nmod_poly_mat_extra/mul_geometric.c @@ -1,3 +1,15 @@ +/* + Copyright (C) 2025 Vincent Neiger, Éric Schost + + 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 diff --git a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_mul_waksman.c b/flint-extras/src/nmod_poly_mat_extra/mul_waksman.c similarity index 91% rename from flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_mul_waksman.c rename to flint-extras/src/nmod_poly_mat_extra/mul_waksman.c index da5ebdc9..4c608c3a 100644 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_mul_waksman.c +++ b/flint-extras/src/nmod_poly_mat_extra/mul_waksman.c @@ -1,3 +1,15 @@ +/* + Copyright (C) 2025 Éric Schost + + 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 @@ -12,7 +24,7 @@ void nmod_poly_mat_mul_waksman(nmod_poly_mat_t C, const nmod_poly_mat_t A, cons n = B->r; p = B->c; mod = A->modulus; - + if (m < 1 || n < 1 || p < 1) { nmod_poly_mat_zero(C); @@ -29,7 +41,7 @@ void nmod_poly_mat_mul_waksman(nmod_poly_mat_t C, const nmod_poly_mat_t A, cons return; } - half = (mod+1) >> 1; // 1/2 + half = (mod+1) >> 1; // 1/2 nmod_poly_init(val0, mod); nmod_poly_init(val1, mod); @@ -54,14 +66,14 @@ void nmod_poly_mat_mul_waksman(nmod_poly_mat_t C, const nmod_poly_mat_t A, cons for (j = 1; j <= np; j++) { j2 = (j << 1) - 1; - + for (k = 0; k < p; k++) { nmod_poly_add(val1, nmod_poly_mat_entry(A, 0, j2-1), nmod_poly_mat_entry(B, j2, k)); nmod_poly_add(val2, nmod_poly_mat_entry(A, 0, j2), nmod_poly_mat_entry(B, j2-1, k)); nmod_poly_mul(val1, val1, val2); nmod_poly_add(nmod_poly_mat_entry(C, 0, k), nmod_poly_mat_entry(C, 0, k), val1); - + nmod_poly_sub(val1, nmod_poly_mat_entry(A, 0, j2-1), nmod_poly_mat_entry(B, j2, k)); nmod_poly_sub(val2, nmod_poly_mat_entry(A, 0, j2), nmod_poly_mat_entry(B, j2-1, k)); nmod_poly_mul(val1, val1, val2); @@ -74,13 +86,13 @@ void nmod_poly_mat_mul_waksman(nmod_poly_mat_t C, const nmod_poly_mat_t A, cons nmod_poly_add(val2, nmod_poly_mat_entry(A, l, j2), nmod_poly_mat_entry(B, j2-1, 0)); nmod_poly_mul(val1, val1, val2); nmod_poly_add(nmod_poly_mat_entry(C, l, 0), nmod_poly_mat_entry(C, l, 0), val1); - + nmod_poly_sub(val1, nmod_poly_mat_entry(A, l, j2-1), nmod_poly_mat_entry(B, j2, 0)); nmod_poly_sub(val2, nmod_poly_mat_entry(A, l, j2), nmod_poly_mat_entry(B, j2-1, 0)); nmod_poly_mul(val1, val1, val2); nmod_poly_add(Ccol + l, Ccol + l, val1); } - + for (k = 1; k < p; k++) { for (l = 1; l < m; l++) @@ -92,7 +104,7 @@ void nmod_poly_mat_mul_waksman(nmod_poly_mat_t C, const nmod_poly_mat_t A, cons } } } - + for (l=1; lr, pmat->c, pmat->modulus); - // TODO improve: set init - nmod_mat_poly_set_trunc_from_poly_mat(matp, pmat, order); - nmod_mat_poly_init(app, pmat->r, pmat->r, pmat->modulus); - nmod_mat_poly_mbasis(app, shift, matp, order); - // TODO improve: set init - nmod_poly_mat_set_from_mat_poly(appbas, app); - nmod_mat_poly_clear(matp); - nmod_mat_poly_clear(app); -} - - -/* -*- 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/nmod_poly_mat_approximant_mbasis1.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_approximant_mbasis1.c deleted file mode 100644 index 4020843e..00000000 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_approximant_mbasis1.c +++ /dev/null @@ -1,199 +0,0 @@ -#include -#include "nmod_mat_extra.h" // for left_nullspace, permutation -#include "nmod_poly_mat_extra.h" // for apply_perm -#include "nmod_poly_mat_utils.h" -#include "nmod_poly_mat_approximant.h" - -/* type for sorting permutation */ -typedef struct -{ - slong value; - slong ord; -} slong_pair; - -/* comparison function for quicksort */ -static inline int compare(const void * a, const void * b) -{ - slong_pair aa = * (const slong_pair *) a; - slong_pair bb = * (const slong_pair *) b; - return aa.value - bb.value; -} - -/** Creates a permutation from the sorting of a shift - * After running this, perm is the list of integers which sorts - * shift increasingly, i.e. - * shift[perm[0]] <= shift[perm[1]] < ... < shift[perm[n-1]] - * - * \param perm permutation (list of integers) - * \param shift list of integer to be sorted increasingly - * \param n length of both perm and shift - * - * TODO make stable - */ -static inline void sort_and_create_perm(slong * perm, const slong * shift, slong n) -{ - slong_pair temp[n]; - for (slong i = 0; i < n; i++) - { - temp[i].value = shift[i]; - temp[i].ord = i; - } - - qsort(temp, n, sizeof(slong_pair), compare); - - for (slong i = 0; i < n; i++) - perm[i] = temp[i].ord; -} - - -void mbasis1(nmod_poly_mat_t appbas, - slong *res_shifts, - const nmod_mat_t mat, - const slong *shifts) -{ - slong rdim = mat->r; - slong i, j, alloc, temp[rdim], - rank_mat, rank_kernel, - *P, *P_inv, *comp, *comp_inv, *perm; - nmod_mat_t K, mat_cp; - nmod_poly_t One, constant; - ulong prime = mat->mod.n; - - perm = _perm_init(rdim); - sort_and_create_perm(perm, shifts, rdim); - - /* LU operation on pi*mat and extraction L = [ [Lr, 0], [G, I] ] */ - - nmod_mat_init_set(mat_cp, mat); - - nmod_mat_permute_rows(mat_cp, perm, NULL); - - P = _perm_init(rdim); - - rank_kernel = nmod_mat_left_nullspace_compact(K, P, mat_cp); - rank_mat = rdim - rank_kernel; - - P_inv = _perm_init(rdim); - for (i = 0; i < rank_kernel; i++) - P_inv[rank_mat + i] = P[i]; - for (i = 0; i < rank_mat; i++) - P_inv[i] = P[rank_kernel + i]; - - _perm_inv(P_inv, P_inv, rdim); - _perm_inv(perm, perm, rdim); - - /* Computation of the block matrix [ [xIr, 0], [K, I_{rdim - rank}] ] */ - nmod_poly_init(One, prime); - nmod_poly_set_coeff_ui(One, 0, 1); // 1 - - for (i = rank_mat; i < rdim; i++) - nmod_poly_set(nmod_poly_mat_entry(appbas, i, i), One); - - nmod_poly_shift_left(One, One, 1); // x - for (i = 0; i < rank_mat; i++) - nmod_poly_set(nmod_poly_mat_entry(appbas, i, i), One); - - nmod_poly_init(constant, prime); - for (i = 0; i < rank_kernel; i++) - for (j = 0; j < rank_mat; j++) - { - alloc = (slong) nmod_mat_get_entry(K, i, j); - if (alloc != 0) - { - nmod_poly_set_coeff_ui(constant, 0, alloc); - nmod_poly_set(nmod_poly_mat_entry(appbas, i + rank_mat, j), constant); - } - } - /* Multiply by the permutations */ - comp = _perm_init(rdim); - comp_inv = _perm_init(rdim); - - _perm_compose(comp, P_inv, perm, rdim); - - _perm_inv(comp_inv, comp, rdim); - - nmod_poly_mat_permute_columns(appbas, comp, NULL); - nmod_poly_mat_permute_rows(appbas, comp_inv, NULL); - - /* Compute the new shift */ - apply_perm_to_vector(temp, shifts, comp, rdim); - - for (i = 0; i < rank_mat; i++) - temp[i] += 1; - - apply_perm_to_vector(res_shifts, temp, comp_inv, rdim); - - nmod_mat_clear(mat_cp); - nmod_mat_clear(K); - nmod_poly_clear(constant); - nmod_poly_clear(One); - - _perm_clear(perm); - _perm_clear(P); - _perm_clear(P_inv); - _perm_clear(comp); - _perm_clear(comp_inv); -} - -slong mbasis1_for_mbasis(nmod_mat_t appbas, - slong *res_shifts, - slong *res_perm, - const nmod_mat_t mat, - const slong *shifts) -{ - slong rdim = mat->r; - slong i, rank_mat, rank_kernel, temp[rdim]; - - slong *P, *P_inv, *comp_inv, *perm = _perm_init(rdim); - nmod_mat_t K, mat_cp; - - sort_and_create_perm(perm, shifts, rdim); - - /* compute left kernel of perm*mat with rank profile */ - nmod_mat_init_set(mat_cp, mat); - - nmod_mat_permute_rows(mat_cp, perm, NULL); - - P = _perm_init(rdim); - - rank_kernel = nmod_mat_left_nullspace_compact(K, P, mat_cp); - rank_mat = rdim - rank_kernel; - - P_inv = _perm_init(rdim); - for (i = 0; i < rank_kernel; i++) - P_inv[rank_mat + i] = P[i]; - for (i = 0; i < rank_mat; i++) - P_inv[i] = P[rank_kernel + i]; - - _perm_inv(perm, perm, rdim); - _perm_inv(perm, perm, rdim); - - nmod_mat_init_set(appbas, K); - - /* Permutations */ - _perm_compose(res_perm, P_inv, perm, rdim); - - /* Compute the new shift */ - apply_perm_to_vector(temp, shifts, res_perm, rdim); - - for (i = 0; i < rank_mat; i++) - temp[i] += 1; - - comp_inv = _perm_init(rdim); - - _perm_inv(comp_inv, res_perm, rdim); - apply_perm_to_vector(res_shifts, temp, comp_inv, rdim); - - nmod_mat_clear(mat_cp); - nmod_mat_clear(K); - - _perm_clear(P); - _perm_clear(perm); - _perm_clear(P_inv); - _perm_clear(comp_inv); - - return rank_mat; -} - -/* -*- 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/nmod_poly_mat_approximant_pmbasis.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_approximant_pmbasis.c deleted file mode 100644 index dffe78ea..00000000 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_approximant_pmbasis.c +++ /dev/null @@ -1,35 +0,0 @@ -#include "nmod_poly_mat_multiply.h" // for middle product -#include "nmod_poly_mat_approximant.h" - -void nmod_poly_mat_pmbasis(nmod_poly_mat_t appbas, - slong * shift, - const nmod_poly_mat_t pmat, - slong order) -{ - if (order <= PMBASIS_THRES) - { - nmod_poly_mat_mbasis(appbas, shift, pmat, order); - return; - } - - const long order1 = order>>1; - const long order2 = order - order1; - nmod_poly_mat_t appbas2, residual; - - nmod_poly_mat_init(appbas2, pmat->r, pmat->r, pmat->modulus); - nmod_poly_mat_init(residual, pmat->r, pmat->c, pmat->modulus); - - nmod_poly_mat_pmbasis(appbas, shift, pmat, order1); - - nmod_poly_mat_middle_product_naive(residual, appbas, pmat, order1, order2-1); - - nmod_poly_mat_pmbasis(appbas2, shift, residual, order2); - - nmod_poly_mat_mul(appbas, appbas2, appbas); - - nmod_poly_mat_clear(appbas2); - nmod_poly_mat_clear(residual); -} - -/* -*- 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/nmod_poly_mat_approximant_verification.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_approximant_verification.c deleted file mode 100644 index 769d6a37..00000000 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_approximant_verification.c +++ /dev/null @@ -1,102 +0,0 @@ -#include "nmod_poly_mat_utils.h" -#include "nmod_poly_mat_forms.h" -#include -#include -#include - -/* TODO currently specialized to ROW_LOWER (or at least ROW_stuff) */ -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) -{ - // context - const slong rdim = pmat->r; - const slong cdim = pmat->c; - const ulong prime = pmat->modulus; - - nmod_poly_mat_t residual; - nmod_poly_t pol; - nmod_mat_t CP0; - - nmod_poly_init(pol, prime); - nmod_poly_mat_init(residual, rdim, cdim, prime); - nmod_mat_init(CP0, rdim, rdim+cdim, prime); - - int success = 1; - - // check appbas is square with the right dimension - if (appbas->r != rdim || appbas->c != rdim) - { - printf("basis has wrong row dimension or column dimension\n"); - success = 0; - } - - // check appbas has form at least "form" - if (!nmod_poly_mat_is_ordered_weak_popov(appbas, shift, orient)) - { - printf("basis is not shifted-reduced\n"); - success = 0; - } - - // compute residual, check rows of appbas are approximants - nmod_poly_mat_mul(residual, appbas, pmat); - - for(slong i = 0; i < rdim; i++) - { - for(slong j = 0; j < cdim; j++) - { - nmod_poly_set_trunc(pol, nmod_poly_mat_entry(residual, i, j), order); - if (!nmod_poly_is_zero(pol)) - { - printf("entry %ld, %ld of residual has valuation less than the order\n",i,j); - success = 0; - } - } - } - - // check generation: follows ideas from Algorithm 1 in Giorgi-Neiger, ISSAC 2018 - - // generation, test 1: check determinant of appbas is lambda * x**D - // since ordered weak Popov, deg-det is sum of diagonal degrees - slong D = 0; - for (slong i = 0; i < rdim; i++) - D += nmod_poly_degree(nmod_poly_mat_entry(appbas, i, i)); - nmod_poly_mat_det(pol, appbas); - if (nmod_poly_degree(pol) != D) - { - printf("determinant is not lambda * x**(sum(diag-deg))"); - success = 0; - } - - // generation, test 2: check that [P(0) C] has full rank - // where C = (appbas * pmat * X^{-order}) mod X - // (coefficient "C" of degree order of the residual) - for (slong i = 0; i < rdim; i++) - { - for (slong j = 0; j < rdim; j++) - { - ulong c = nmod_poly_get_coeff_ui(nmod_poly_mat_entry(appbas, i, j), 0); - nmod_mat_set_entry(CP0, i, j, c); - } - for (slong j = 0; j < cdim; j++) - { - ulong c = nmod_poly_get_coeff_ui(nmod_poly_mat_entry(residual, i, j), order); - nmod_mat_set_entry(CP0, i, rdim+j, c); - } - } - - slong rank = nmod_mat_rank(CP0); - if (rank != rdim) - { - printf("generation test (step: [C P(0)] has full rank) failed"); - success = 0; - } - - nmod_poly_clear(pol); - nmod_poly_mat_clear(residual); - nmod_mat_clear(CP0); - - return success; -} diff --git a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_degree_matrix.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_degree_matrix.c deleted file mode 100644 index 9fb88fad..00000000 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_degree_matrix.c +++ /dev/null @@ -1,37 +0,0 @@ -#include -#include -#include "nmod_poly_mat_forms.h" - -void nmod_poly_mat_degree_matrix(fmpz_mat_t dmat, - const nmod_poly_mat_t mat) -{ - for(slong i = 0; i < mat->r; i++) - for(slong j = 0; j < mat->c; j++) - *fmpz_mat_entry(dmat, i, j) = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)); -} - -void nmod_poly_mat_degree_matrix_shifted(fmpz_mat_t dmat, - const nmod_poly_mat_t mat, - const slong * shift, - orientation_t orient) -{ - if (shift) - { - if (orient == ROW_LOWER || orient == ROW_UPPER) - for(slong i = 0; i < mat->r; i++) - for(slong j = 0; j < mat->c; j++) - *fmpz_mat_entry(dmat, i, j) = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)) + shift[j]; - else // orient == COL_* - for(slong i = 0; i < mat->r; i++) - for(slong j = 0; j < mat->c; j++) - *fmpz_mat_entry(dmat, i, j) = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)) + shift[i]; - } - else - for(slong i = 0; i < mat->r; i++) - for(slong j = 0; j < mat->c; j++) - *fmpz_mat_entry(dmat, i, j) = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)); -} - - -/* -*- 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/nmod_poly_mat_description.bak b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_description.bak index b95e9013..e37ed8c3 100644 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_description.bak +++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_description.bak @@ -1,3 +1,15 @@ +/* + 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 @@ -289,12 +301,3 @@ int nmod_poly_mat_kernel(nmod_poly_mat_t N, const nmod_poly_mat_t M, slong delta return nbnull; } - - - - - - - -/* -*- 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/nmod_poly_mat_det.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_det.c index fc87f257..d6f859ef 100644 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_det.c +++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_det.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 #include "nmod_poly_mat_utils.h" // for permute_rows_by_sorting_vec #include "nmod_poly_mat_forms.h" @@ -67,6 +79,3 @@ void nmod_poly_mat_det_iter(nmod_poly_t det, nmod_poly_mat_t mat) nmod_poly_mul(det, det, nmod_poly_mat_entry(view, i, i)); nmod_poly_mat_window_clear(view); } - -/* -*- 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/nmod_poly_mat_dixon.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_dixon.c index 785c8be9..6cee6dd1 100644 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_dixon.c +++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_dixon.c @@ -1,3 +1,15 @@ +/* + 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 @@ -186,5 +198,3 @@ void nmod_poly_mat_dixon(nmod_poly_mat_t X, } -/* -*- 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/nmod_poly_mat_hermite.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_hermite.c index c5afcea8..a57cb892 100644 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_hermite.c +++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_hermite.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 #include #include @@ -1021,6 +1033,3 @@ slong nmod_poly_mat_uref_matrixgcd_iter(nmod_poly_mat_t mat, flint_free(perm); return rk; } - -/* -*- 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/nmod_poly_mat_io.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_io.c index c2449e1a..1e6f8e87 100644 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_io.c +++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_io.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 #include #include @@ -79,6 +91,3 @@ void nmod_poly_mat_leading_matrix_print_pretty(const nmod_poly_mat_t mat, nmod_mat_print_pretty(lmat); nmod_mat_clear(lmat); } - -/* -*- 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/nmod_poly_mat_is_form.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_is_form.c index 81821ebf..5ea873f9 100644 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_is_form.c +++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_is_form.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 // qsort #include #include @@ -300,6 +312,3 @@ int nmod_poly_mat_is_hermite(const nmod_poly_mat_t mat, orientation_t orient) flint_free(pivdeg); return 1; } - -/* -*- 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/nmod_poly_mat_is_unimodular.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_is_unimodular.c index 4d5aa934..2f5819f3 100644 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_is_unimodular.c +++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_is_unimodular.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 "nmod_poly_mat_utils.h" #include #include @@ -48,6 +60,3 @@ int nmod_poly_mat_is_unimodular_randomized(const nmod_poly_mat_t pmat, flint_ran return (det1 != 0 && det1 == det2); } - -/* -*- 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/nmod_poly_mat_leading_matrix.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_leading_matrix.c deleted file mode 100644 index 251ddfe3..00000000 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_leading_matrix.c +++ /dev/null @@ -1,75 +0,0 @@ -#include -#include -#include "nmod_poly_mat_forms.h" - -void nmod_poly_mat_leading_matrix(nmod_mat_t lmat, - const nmod_poly_mat_t mat, - const slong * shift, - orientation_t orient) -{ - if (orient == ROW_LOWER || orient == ROW_UPPER) - { - // find row degrees - slong rdeg[mat->r]; - nmod_poly_mat_row_degree(rdeg, mat, shift); - - // deduce leading matrix - if (shift == NULL) - { - for (slong i = 0; i < mat->r; i++) - if (rdeg[i] >= 0) - for (slong j = 0; j < mat->c; j++) - nmod_mat_set_entry(lmat, i, j, - nmod_poly_get_coeff_ui(nmod_poly_mat_entry(mat, i, j), - rdeg[i])); - else - for (slong j = 0; j < mat->c; j++) - nmod_mat_set_entry(lmat, i, j, 0); - } - else - { - for (slong i = 0; i < mat->r; i++) - for (slong j = 0; j < mat->c; j++) - if (rdeg[i] >= shift[j]) - nmod_mat_set_entry(lmat, i, j, - nmod_poly_get_coeff_ui(nmod_poly_mat_entry(mat, i, j), - rdeg[i] - shift[j])); - else - nmod_mat_set_entry(lmat, i, j, 0); - } - } - else // orient == COL_* - { - // find column degrees - slong cdeg[mat->c]; - nmod_poly_mat_column_degree(cdeg, mat, shift); - - // deduce leading matrix - if (shift == NULL) - { - for (slong j = 0; j < mat->c; j++) - if (cdeg[j] >= 0) - for (slong i = 0; i < mat->r; i++) - nmod_mat_set_entry(lmat, i, j, - nmod_poly_get_coeff_ui(nmod_poly_mat_entry(mat, i, j), - cdeg[j])); - else - for (slong i = 0; i < mat->r; i++) - nmod_mat_set_entry(lmat, i, j, 0); - } - else - { - for (slong j = 0; j < mat->c; j++) - for (slong i = 0; i < mat->r; i++) - if (cdeg[j] >= shift[i]) - nmod_mat_set_entry(lmat, i, j, - nmod_poly_get_coeff_ui(nmod_poly_mat_entry(mat, i, j), - cdeg[j] - shift[i])); - else - nmod_mat_set_entry(lmat, i, j, 0); - } - } -} - -/* -*- 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/nmod_poly_mat_middle_product_naive.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_middle_product_naive.c deleted file mode 100644 index c24d01f8..00000000 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_middle_product_naive.c +++ /dev/null @@ -1,14 +0,0 @@ -#include - -/** Middle product for polynomial matrices - * sets C = ((A * B) div x^dA) mod x^(dB+1), assuming deg(A) <= dA and deg(B) <= dA + dB - * output can alias input - * naive implementation (multiply, shift, truncate) - */ -void nmod_poly_mat_middle_product_naive(nmod_poly_mat_t C, const nmod_poly_mat_t A, const nmod_poly_mat_t B, - const ulong dA, const ulong dB) -{ - nmod_poly_mat_mul(C, A, B); - nmod_poly_mat_shift_right(C, C, dA); - nmod_poly_mat_truncate(C, dB + 1); -} diff --git a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_rdeg_cdeg.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_rdeg_cdeg.c deleted file mode 100644 index 8d3580df..00000000 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_rdeg_cdeg.c +++ /dev/null @@ -1,90 +0,0 @@ -#include - -#include "nmod_poly_mat_forms.h" - -void nmod_poly_mat_row_degree(slong *rdeg, const nmod_poly_mat_t mat, const slong *shift) -{ - if (shift == NULL) - { - slong max, d; - for (slong i = 0; i < mat->r; i++) - { - max = -1; - for (slong j = 0; j < mat->c; j++) - { - d = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)); - if (max < d) - max = d; - } - rdeg[i] = max; - } - } - else - { - slong max, d; - slong min_shift = (mat->c > 0) ? shift[0] : 0; - - // find minimum of shift - for (slong j = 0; j < mat->c; ++j) - if (shift[j] < min_shift) - min_shift = shift[j]; - - for (slong i = 0; i < mat->r; i++) - { - max = min_shift-1; // zero rows will have this as rdeg - for (slong j = 0; j < mat->c; j++) - { - d = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)) + shift[j]; - // if new maximum at a nonzero entry, update - if (shift[j] <= d && max < d) - max = d; - } - rdeg[i] = max; - } - } -} - -void nmod_poly_mat_column_degree(slong *cdeg, const nmod_poly_mat_t mat, const slong *shift) -{ - if (shift == NULL) - { - slong max, d; - for (slong j = 0; j < mat->c; j++) - { - max = -1; - for (slong i = 0; i < mat->r; i++) - { - d = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)); - if (max < d) - max = d; - } - cdeg[j] = max; - } - } - else - { - slong max, d; - slong min_shift = (mat->r > 0) ? shift[0] : 0; - - // find minimum of shift - for (slong i = 0; i < mat->r; ++i) - if (shift[i] < min_shift) - min_shift = shift[i]; - - for (slong j = 0; j < mat->c; j++) - { - max = min_shift-1; - for (slong i = 0; i < mat->r; i++) - { - d = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)) + shift[i]; - // if new maximum at a nonzero entry, update - if (shift[i] <= d && max < d) - max = d; - } - cdeg[j] = max; - } - } -} - -/* -*- 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/nmod_poly_mat_utils.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_utils.c index b2eef20a..7520d17a 100644 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_utils.c +++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_utils.c @@ -1,9 +1,20 @@ +/* + 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 qsort #include #include "nmod_poly_mat_utils.h" #include "nmod_mat_poly.h" - /********************************************************************** * ROW ROTATION DOWNWARD/UPWARD * **********************************************************************/ @@ -165,7 +176,6 @@ void _nmod_poly_mat_permute_columns_by_sorting_vec(nmod_poly_mat_t mat, } - /*------------------------------------------------------------*/ /*------------------------------------------------------------*/ /* SET FROM MATRIX POLYNOMIAL */ @@ -198,6 +208,3 @@ void nmod_poly_mat_set_trunc_from_mat_poly(nmod_poly_mat_t pmat, _nmod_poly_normalise(nmod_poly_mat_entry(pmat, i, j)); } } - -/* -*- 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/nmod_poly_mat_weak_popov.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_weak_popov.c index 4149fa55..3929f2b9 100644 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_weak_popov.c +++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_weak_popov.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 #include #include @@ -274,6 +286,3 @@ slong nmod_poly_mat_ordered_weak_popov_iter(nmod_poly_mat_t mat, flint_free(perm); return rk; } - -/* -*- 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/profile/p-kernel.c b/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c new file mode 100644 index 00000000..7587bfa8 --- /dev/null +++ b/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c @@ -0,0 +1,198 @@ +#include // for atoi + +#include +#include +#include +#include +#include + +#include "nmod_poly_mat_utils.h" +#include "nmod_poly_mat_kernel.h" + +#define MEASURE_SAMPLE 0 + +/* Shift types are: + * -1 -> [TODO not implemented yet] negative shift that represents known degree bound on a kernel basis + * 0 -> uniform, shift == (0,...,0) + * 1 -> input degrees: shift == row degrees of input matrix + * 100 -> shift that yields the kernel basis Hermite form + * */ + +typedef struct +{ + slong rdim; /* row dimension */ + slong cdim; /* column dimension */ + slong deg; /* degree */ + slong rank; /* rank */ + slong stype; /* shift type */ + slong modn; /* modulus */ +} +time_args; + +#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 F; \ + nmod_poly_mat_init(F, rdim, cdim, n); \ + nmod_poly_mat_rand(F, state, deg); \ + \ + /* TMP */ \ + nmod_poly_mat_t Ft; \ + nmod_poly_mat_init(Ft, cdim, rdim, n); \ + nmod_poly_mat_transpose(Ft, F); \ + /* END TMP */ \ + \ + slong degN[rdim]; \ + nmod_poly_mat_t K; \ + nmod_poly_mat_init(K, rdim, rdim, n); \ + \ + double FLINT_SET_BUT_UNUSED(tcpu), twall; \ + \ + TIMEIT_START; \ + nmod_poly_mat_##fun(K, degN, Ft, NULL, 3.); \ + TIMEIT_STOP_VALUES(tcpu, twall); \ + \ + printf("%.2e", twall); \ + \ + nmod_poly_mat_clear(F); \ + nmod_poly_mat_clear(Ft); /* TMP */ \ +} + +TIME_KER(kernel_zls) + +/*-------------------------*/ +/* main */ +/*-------------------------*/ + +int main(int argc, char ** argv) +{ + flint_rand_t state; + 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}; + + /* TODO rank */ + /* TODO shift type */ + + // bench functions + const slong nfuns = 1; + typedef void (*timefun) (time_args, flint_rand_t); + const timefun funs[] = { + time_kernel_zls, // 0 + }; + + // TODO + //typedef void (*samplefun) (void*, ulong); + //const samplefun sfuns[] = { + // sample_mul, // 0 + // sample_mul_geometric, // 1 + //}; + //#if MEASURE_SAMPLE + // const samplefun sfun = sfuns[ifun]; + // double min, max; + // prof_repeat(&min, &max, sfun, (void*) &targs); + // printf("%.2e", min/1000000); + //#else + // const timefun tfun = funs[ifun]; + // tfun(targs, state); + //#endif + + const char * description[] = { + "#0 --> kernel (general interface) ", + "#1 --> .... TODO ", + }; + + if (argc < 4 || argc > 6) // show usage + { + printf("Usage: `%s [nbits] [rdim] [cdim] [deg] [rank] [shift] [fun]`\n", argv[0]); + printf(" No argument shows this help.\n"); + printf(" First 3 arguments are mandatory.\n"); + printf(" [rank] and [shift] not supported yet.\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]); + + return 0; + } + + printf("#warmup...\n"); + for (slong i = 0; i < 3; i++) + { + /* 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); + 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 + { + const slong rdim = atoi(argv[2]); + const slong cdim = atoi(argv[3]); + const slong deg = atoi(argv[4]); + 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]; + 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}; + tfun(targs, state); + printf(" "); + printf("\n"); + } + } + + flint_rand_clear(state); + return 0; +} 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 8851a374..913204c0 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 @@ -108,7 +108,7 @@ int main(int argc, char ** argv) printf(" (nbits == -1 launches full suite)\n"); printf(" - dim: matrices are dim x dim\n"); printf(" - deg: matrices are random of degree < deg\n"); - printf(" - deg: id number of the timed function (see below),\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]); diff --git a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_rand.c b/flint-extras/src/nmod_poly_mat_extra/randomisation.c similarity index 94% rename from flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_rand.c rename to flint-extras/src/nmod_poly_mat_extra/randomisation.c index cb8234ed..e692ef48 100644 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_rand.c +++ b/flint-extras/src/nmod_poly_mat_extra/randomisation.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 #include #include @@ -5,7 +17,6 @@ #include "nmod_poly_mat_utils.h" #include "nmod_poly_mat_forms.h" - /*------------------------------------------------------------*/ /*------------------------------------------------------------*/ /* RANDOM - UNIFORM, DENSE */ @@ -161,7 +172,3 @@ void nmod_poly_mat_rand_popov(nmod_poly_mat_t mat, flint_free(uniform_shift); flint_free(iota); } - - -/* -*- 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/test/main.c b/flint-extras/src/nmod_poly_mat_extra/test/main.c index 165b48de..e89828cd 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_zls.c" #include "t-middle_product_geometric.c" #include "t-mul_geometric.c" #include "t-mbasis.c" @@ -24,23 +25,22 @@ #include "t-mul_waksman.c" #include "t-rand.c" #include "t-weak_popov_form.c" -#include "t-kernel.c" /* Array of test functions ***************************************************/ test_struct tests[] = { - TEST_FUNCTION(nmod_poly_mat_middle_product_geometric), - TEST_FUNCTION(nmod_poly_mat_mul_geometric), TEST_FUNCTION(nmod_poly_mat_det), TEST_FUNCTION(nmod_poly_mat_dixon), TEST_FUNCTION(nmod_poly_mat_hnf), + 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), TEST_FUNCTION(nmod_poly_mat_pmbasis), /* TEST_FUNCTION(nmod_poly_mat_mul_waksman), */ /* TODO */ TEST_FUNCTION(nmod_poly_mat_rand), TEST_FUNCTION(nmod_poly_mat_weak_popov_form), - TEST_FUNCTION(nmod_poly_mat_kernel) }; /* main function *************************************************************/ diff --git a/flint-extras/src/nmod_poly_mat_extra/test/t-dixon.c b/flint-extras/src/nmod_poly_mat_extra/test/t-dixon.c index ababa3ba..801fc059 100644 --- a/flint-extras/src/nmod_poly_mat_extra/test/t-dixon.c +++ b/flint-extras/src/nmod_poly_mat_extra/test/t-dixon.c @@ -246,6 +246,3 @@ TEST_FUNCTION_START(nmod_poly_mat_dixon, state) TEST_FUNCTION_END(state); } } - -/* -*- 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/test/t-kernel.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls.c similarity index 56% rename from flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c rename to flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls.c index 99e54e67..94a6edbe 100644 --- a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel_zls.c @@ -10,131 +10,109 @@ . */ +#include #include #include #include #include -#include "nmod_poly_mat_utils.h" -#include "nmod_poly_mat_extra.h" +#include "nmod_poly_mat_extra.h" +#include "nmod_poly_mat_forms.h" // test one given input -int core_test_kernel(const nmod_poly_mat_t mat) +/* 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) { + slong m = mat->r; + slong n = mat->c; - slong m=mat->r; - slong n=mat->c; - - // Flint rank - // ---------- - - slong rkflint; + /* 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); - rkflint=nmod_poly_mat_rank(mat); - - - // PML nullspace - // ------------- - - nmod_poly_mat_t N; + nmod_poly_mat_t N; nmod_poly_mat_init(N, n, n, mat->modulus); - slong nz; - - int i,j; - - // slong iz[m]; - - // for (i = 0; i < m; i++) - // { - // iz[i]=0; - // } - slong degN[n]; - nz=nmod_poly_mat_kernel(N, degN, mat, NULL, 2); - - int verif; - verif = (n-rkflint == nz); - - if (nz !=0) { - - nmod_poly_mat_t NN; - nmod_poly_mat_init(NN, n, nz, mat->modulus); + slong nz = nmod_poly_mat_kernel_zls(N, degN, mat, NULL, 2.); - for (i = 0; i < n; i++) - for (j = 0; j < nz; j++) { - nmod_poly_set(nmod_poly_mat_entry(NN, i, j), nmod_poly_mat_entry(N, i, j)); - } + 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 Z; - nmod_poly_mat_init(Z, m, nz, mat->modulus); + nmod_poly_mat_t Mt; + nmod_poly_mat_init(Mt, n, m, mat->modulus); + nmod_poly_mat_transpose(Mt, mat); - nmod_poly_mat_mul(Z, mat, NN); + int verif = nmod_poly_mat_is_kernel(Nt, Mt, rdeg, ROW_LOWER); - verif = verif && nmod_poly_mat_is_zero(Z); + nmod_poly_mat_clear(N); + nmod_poly_mat_clear(Nt); + nmod_poly_mat_clear(Mt); + flint_free(rdeg); - nmod_poly_mat_clear(Z); - nmod_poly_mat_clear(NN); - - } - - return verif; + return verif; } - -TEST_FUNCTION_START(nmod_poly_mat_kernel, state) +TEST_FUNCTION_START(nmod_poly_mat_kernel_zls, state) { - flint_rand_t state; - flint_rand_init(state); - srand(time(NULL)); - flint_rand_set_seed(state, rand(), rand()); - int i,result; for (i = 0; i < 16 * flint_test_multiplier(); i++) { - - 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, 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 prime = n_randprime(state, nbits, 1); - nmod_poly_mat_t A; - if (i < 4) { + if (i < 4) + { cdim=rdim; nmod_poly_mat_init(A, rdim, cdim, prime); nmod_poly_mat_randtest_sparse(A, state, deg+1, 1.0); } - else if (i < 8) { + else if (i < 8) + { cdim=rdim+1; nmod_poly_mat_init(A, rdim, cdim, prime); nmod_poly_mat_randtest_sparse(A, state, deg+1, 0.8); } - else if (i < 12) { + else if (i < 12) + { nmod_poly_mat_init(A, rdim, cdim, prime); nmod_poly_mat_randtest_sparse(A, state, deg+1, 0.2); } - else { + else + { nmod_poly_mat_init(A, rdim, cdim, prime); nmod_poly_mat_randtest_sparse(A, state, deg+1, 0.84); - } + } - result = core_test_kernel(A); + result = core_test_kernel_zls(A); nmod_poly_mat_clear(A); - if (!result) { + if (!result) + { TEST_FUNCTION_FAIL("rdim = %wu, cdim = %wu, degree = %wu, p = %wu\n", \ rdim, cdim, deg, prime); } } - // Degree zero + // Degree zero // ----------- for (i = 0; i < 2; i++) @@ -143,18 +121,19 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state) ulong cdim = rdim + n_randint(state, 8); ulong deg = 0; - ulong prime = 7; + ulong prime = 7; nmod_poly_mat_t A; nmod_poly_mat_init(A, rdim, cdim, prime); nmod_poly_mat_randtest_sparse(A, state, deg+1, 0.8); - result = core_test_kernel(A); + result = core_test_kernel_zls(A); nmod_poly_mat_clear(A); - if (!result) { + if (!result) + { TEST_FUNCTION_FAIL("rdim = %wu, cdim = %wu, degree = %wu, p = %wu\n", \ rdim, cdim, deg, prime); } @@ -169,25 +148,23 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state) ulong cdim = rdim -8; ulong deg = 1; - ulong prime = 2; + ulong prime = 2; nmod_poly_mat_t A; nmod_poly_mat_init(A, rdim, cdim, prime); nmod_poly_mat_randtest_sparse(A, state, deg+1, 0.2); - result = core_test_kernel(A); + result = core_test_kernel_zls(A); nmod_poly_mat_clear(A); - if (!result) { + if (!result) + { TEST_FUNCTION_FAIL("rdim = %wu, cdim = %wu, degree = %wu, p = %wu\n", \ rdim, cdim, deg, prime); } } -TEST_FUNCTION_END(state); + TEST_FUNCTION_END(state); } - -/* -*- 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/timings/time_hnf.c b/flint-extras/src/nmod_poly_mat_extra/timings/time_hnf.c index 6bb58a30..13d9367d 100644 --- a/flint-extras/src/nmod_poly_mat_extra/timings/time_hnf.c +++ b/flint-extras/src/nmod_poly_mat_extra/timings/time_hnf.c @@ -452,6 +452,3 @@ int main(int argc, char *argv[]) return 0; } - -/* -*- 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/timings/time_pmbasis.c b/flint-extras/src/nmod_poly_mat_extra/timings/time_pmbasis.c index 8a57cdb7..a527a31a 100644 --- a/flint-extras/src/nmod_poly_mat_extra/timings/time_pmbasis.c +++ b/flint-extras/src/nmod_poly_mat_extra/timings/time_pmbasis.c @@ -147,6 +147,3 @@ int main(int argc, char *argv[]) return 0; } - -/* -*- 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/timings/time_polmatmul.c b/flint-extras/src/nmod_poly_mat_extra/timings/time_polmatmul.c index 66d51024..7088ce01 100644 --- a/flint-extras/src/nmod_poly_mat_extra/timings/time_polmatmul.c +++ b/flint-extras/src/nmod_poly_mat_extra/timings/time_polmatmul.c @@ -329,6 +329,3 @@ int main(int argc, char *argv[]) return 0; } - -/* -*- 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/verification.c b/flint-extras/src/nmod_poly_mat_extra/verification.c new file mode 100644 index 00000000..1ce2fe37 --- /dev/null +++ b/flint-extras/src/nmod_poly_mat_extra/verification.c @@ -0,0 +1,174 @@ +/* + 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 "nmod_poly_mat_forms.h" +#include +#include +#include + +/* TODO currently specialized to ROW_LOWER (or at least ROW_stuff) */ +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) +{ + const slong rdim = pmat->r; + const slong cdim = pmat->c; + const ulong prime = pmat->modulus; + + nmod_poly_mat_t residual; + nmod_poly_t pol; + nmod_mat_t CP0; + + nmod_poly_init(pol, prime); + nmod_poly_mat_init(residual, rdim, cdim, prime); + nmod_mat_init(CP0, rdim, rdim+cdim, prime); + + int success = 1; + + /* check appbas is square with the right dimension */ + if (appbas->r != rdim || appbas->c != rdim) + { + printf("basis has wrong row dimension or column dimension\n"); + success = 0; + } + + /* check appbas is shifted reduced */ + if (!nmod_poly_mat_is_ordered_weak_popov(appbas, shift, orient)) + { + printf("basis is not shifted-weak Popov\n"); + success = 0; + } + + /* compute residual, check rows of appbas are approximants */ + nmod_poly_mat_mul(residual, appbas, pmat); + + for (slong i = 0; i < rdim; i++) + { + for (slong j = 0; j < cdim; j++) + { + nmod_poly_set_trunc(pol, nmod_poly_mat_entry(residual, i, j), order); + if (!nmod_poly_is_zero(pol)) + { + printf("entry %ld, %ld of residual has valuation less than the order\n",i,j); + success = 0; + } + } + } + + /* check generation: follows ideas from Algorithm 1 in Giorgi-Neiger, ISSAC 2018 */ + + /* generation, test 1: check determinant of appbas is lambda * x**D * + * since ordered weak Popov, deg(det(appbas)) is sum of diagonal degrees * + */ + slong D = 0; + for (slong i = 0; i < rdim; i++) + D += nmod_poly_degree(nmod_poly_mat_entry(appbas, i, i)); + nmod_poly_mat_det(pol, appbas); + if (nmod_poly_degree(pol) != D) + { + printf("determinant is not lambda * x**(sum(diag-deg))"); + success = 0; + } + + /* generation, test 2: check that [P(0) C] has full rank * + * where C = (appbas * pmat * X^{-order}) mod X * + * (coefficient "C" of degree order of the residual) * + **/ + for (slong i = 0; i < rdim; i++) + { + for (slong j = 0; j < rdim; j++) + { + ulong c = nmod_poly_get_coeff_ui(nmod_poly_mat_entry(appbas, i, j), 0); + nmod_mat_set_entry(CP0, i, j, c); + } + for (slong j = 0; j < cdim; j++) + { + ulong c = nmod_poly_get_coeff_ui(nmod_poly_mat_entry(residual, i, j), order); + nmod_mat_set_entry(CP0, i, rdim+j, c); + } + } + + slong rank = nmod_mat_rank(CP0); + if (rank != rdim) + { + printf("generation test (step: [C P(0)] has full rank) failed"); + success = 0; + } + + nmod_poly_clear(pol); + nmod_poly_mat_clear(residual); + nmod_mat_clear(CP0); + + 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 */ +int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker, + const nmod_poly_mat_t pmat, + const slong * shift, + orientation_t orient) +{ + 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); + + 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 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)) + { + printf("basis is not shifted-reduced\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 rank */ + slong rk = nmod_poly_mat_rank(pmat); + if (kdim != 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_kernel.h b/flint-extras/src/nmod_poly_mat_kernel.h index 81b4d0eb..20cb54b6 100644 --- a/flint-extras/src/nmod_poly_mat_kernel.h +++ b/flint-extras/src/nmod_poly_mat_kernel.h @@ -20,44 +20,46 @@ extern "C" { #endif /** - * - * Right shifted kernel of a polynomial matrix, assuming that the columns has been - * sorted by shifted degree - * + * + * Right shifted kernel of a polynomial matrix, assuming that the columns has been + * sorted by shifted degree + * * Algorithm of Wei Zhou, George Labahn, and Arne Storjohann * "Computing Minimal Nullspace Bases" * ISSAC 2012, https://dl.acm.org/doi/abs/10.1145/2442829.2442881 - * - * Calls nmod_poly_mat_zls_sorted after an initial sorting - * - * TODO/TO SEE: - * - * Input: - * iA in m x n - * ishift[n], NULL (the degrees are computed) or initialized outside, - * the shift for the kernel + * + * Calls nmod_poly_mat_zls_sorted after an initial sorting + * + * TODO/TO SEE: + * + * Input: + * iA in m x n + * ishift[n], NULL (the degrees are computed) or initialized outside, + * the shift for the kernel * values should be at least 0 (even for zero columns in A - * "with entries arranged in non-decreasing order and bounding the - * corresponding column degrees of A." - * kappa, a double >= 2, for the order of the order bases - * kappa * s instead of 3 *s in ZLS + * "with entries arranged in non-decreasing order and bounding the + * corresponding column degrees of A." + * kappa, a double >= 2, for the order of the order bases + * kappa * s instead of 3 *s in ZLS * * Output: - * returns the dimension w of the kernel, which may be zero - * N, is initialized n x n outside - * its first w columns give a minimal basis of the kernel + * returns the dimension w of the kernel, which may be zero + * N, is initialized n x n outside + * its first w columns give a minimal basis of the kernel * degN[n], initialized outside, its first w entries are concerned, - * they are the ishift shifted degrees of the kernel basis - * + * they are the ishift shifted degrees of the kernel basis + * */ -int nmod_poly_mat_kernel(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat_t A, \ - const slong *ishift, const double kappa); +/* int nmod_poly_mat_kernel(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat_t A, \ */ +/* const slong *ishift, const double kappa); */ +int nmod_poly_mat_kernel_zls(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat_t A, \ + const slong *ishift, const double kappa); /** - * Experimental, should not be really considered + * Experimental, should not be really considered * */ @@ -70,8 +72,3 @@ int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong *degN, const nmod_ #endif #endif - -/* -*- 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 - -