From cfe81a2e534ebc7151d3a9246f3df3e4a03a2bd3 Mon Sep 17 00:00:00 2001 From: Gilles Villard Date: Sun, 14 Dec 2025 16:55:41 +0100 Subject: [PATCH 1/8] Modifications of headers --- .../src/nmod_poly_mat_description.bak | 46 +++++++++++ flint-extras/src/nmod_poly_mat_extra.h | 4 +- flint-extras/src/nmod_poly_mat_kernel.h | 77 +++++++++++++++++++ flint-extras/src/nmod_poly_mat_utils.h | 7 ++ 4 files changed, 133 insertions(+), 1 deletion(-) create mode 100644 flint-extras/src/nmod_poly_mat_description.bak create mode 100644 flint-extras/src/nmod_poly_mat_kernel.h diff --git a/flint-extras/src/nmod_poly_mat_description.bak b/flint-extras/src/nmod_poly_mat_description.bak new file mode 100644 index 00000000..a9035e7d --- /dev/null +++ b/flint-extras/src/nmod_poly_mat_description.bak @@ -0,0 +1,46 @@ +/* + 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 + . +*/ + +#ifndef NMOD_POLY_MAT_DESCRIPTION_H +#define NMOD_POLY_MAT_DESCRIPTION_H + +#include +#include + +#ifdef __cplusplus +extern "C" { +#endif + +slong nmod_poly_mat_left_description(nmod_poly_mat_t N, nmod_poly_mat_t D, + const nmod_poly_mat_t H, + slong delta); + + +slong nmod_poly_mat_right_description(nmod_poly_mat_t N, nmod_poly_mat_t D, + const nmod_poly_mat_t H, + slong delta); + + +int nmod_poly_mat_kernel(nmod_poly_mat_t N, const nmod_poly_mat_t M, + slong delta); + + +#ifdef __cplusplus +} +#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 + + diff --git a/flint-extras/src/nmod_poly_mat_extra.h b/flint-extras/src/nmod_poly_mat_extra.h index 690b2e87..70f5d99c 100644 --- a/flint-extras/src/nmod_poly_mat_extra.h +++ b/flint-extras/src/nmod_poly_mat_extra.h @@ -51,7 +51,6 @@ #include "nmod_poly_mat_approximant.h" // #include "nmod_poly_mat_interpolant.h" -// #include "nmod_poly_mat_kernel.h" // #include "nmod_poly_mat_determinant.h" // #include "nmod_poly_mat_inverse.h" @@ -61,6 +60,9 @@ #include "nmod_poly_mat_dixon.h" +#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, diff --git a/flint-extras/src/nmod_poly_mat_kernel.h b/flint-extras/src/nmod_poly_mat_kernel.h new file mode 100644 index 00000000..8e280841 --- /dev/null +++ b/flint-extras/src/nmod_poly_mat_kernel.h @@ -0,0 +1,77 @@ +/* + 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 + . +*/ + +#ifndef NMOD_POLY_MAT_KERNEL_H +#define NMOD_POLY_MAT_KERNEL_H + +#include "pml.h" + +#ifdef __cplusplus +extern "C" { +#endif + +/** + * + * 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 + * 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, an integer >= 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 + * degN[n], initialized outside, its first w entries are concerned, + * 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 slong kappa); + + + +/** + * Experimental, should not be really considered + * + */ + +int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat_t A, \ + const slong *ishift); + + +#ifdef __cplusplus +} +#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 + + diff --git a/flint-extras/src/nmod_poly_mat_utils.h b/flint-extras/src/nmod_poly_mat_utils.h index 88d0d5b9..ad96d2b0 100644 --- a/flint-extras/src/nmod_poly_mat_utils.h +++ b/flint-extras/src/nmod_poly_mat_utils.h @@ -138,6 +138,13 @@ void _nmod_poly_mat_permute_rows_by_sorting_vec(nmod_poly_mat_t mat, slong * perm); + +void _nmod_poly_mat_permute_columns_by_sorting_vec(nmod_poly_mat_t mat, + slong r, + slong * vec, + slong * perm); + + /*------------------------------------------------------------*/ /*------------------------------------------------------------*/ /* SWAP, PERMUTE, TRANSPOSE */ From 60c2bdd9398b98ddcb4a14a32202b3306b0fd449 Mon Sep 17 00:00:00 2001 From: Gilles Villard Date: Sun, 14 Dec 2025 17:00:59 +0100 Subject: [PATCH 2/8] old header --- .../src/nmod_poly_mat_description.bak | 308 ++++++++++++++++-- 1 file changed, 281 insertions(+), 27 deletions(-) diff --git a/flint-extras/src/nmod_poly_mat_description.bak b/flint-extras/src/nmod_poly_mat_description.bak index a9035e7d..b95e9013 100644 --- a/flint-extras/src/nmod_poly_mat_description.bak +++ b/flint-extras/src/nmod_poly_mat_description.bak @@ -1,46 +1,300 @@ -/* - Copyright (C) 2025 Gilles Villard +#include +#include +#include +#include - 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_extra.h" +#include "nmod_poly_mat_io.h" -#ifndef NMOD_POLY_MAT_DESCRIPTION_H -#define NMOD_POLY_MAT_DESCRIPTION_H +#include "nmod_mat_extra.h" +#include "nmod_poly_mat_dixon.h" +#include "nmod_poly_mat_description.h" -#include -#include -#ifdef __cplusplus -extern "C" { -#endif +/** + * Left description computation for H(x) n x m in K(x) (power series) + * with target degree delta + * + * requires enough precision in input: ie at least (m+n)*delta/m +1 + * + * returns nbrows and a partial (or full) description when 0 < nbrows <= n rows, + * or zero if no candidates + * + * N n x m and and D n x n are initialized outside + * hence the result is part of N and D + * + */ -slong nmod_poly_mat_left_description(nmod_poly_mat_t N, nmod_poly_mat_t D, +slong nmod_poly_mat_left_description(nmod_poly_mat_t N, nmod_poly_mat_t D, const nmod_poly_mat_t H, - slong delta); + slong delta) + +{ + slong i,j; + + slong n = H->r; + slong m = H->c; + + slong sigma = ceil((double) (m+n)*delta/m +1); + + nmod_poly_mat_t M; + nmod_poly_mat_init(M, n+m, m, H->modulus); + + nmod_poly_mat_set(M, H); // GV, correct even the dimension of M is larger ? + + nmod_poly_t mone; + nmod_poly_init(mone, H->modulus); + nmod_poly_set_coeff_ui(mone, 0, H->modulus -1); // 1 + + for (i = 0; i < m; i++) + nmod_poly_set(nmod_poly_mat_entry(M, n+i, i), mone); + + // printf("\n"); + // nmod_poly_mat_print_pretty(M, "x"); + // printf("\n"); + + nmod_poly_mat_t B; + nmod_poly_mat_init(B, n+m, n+m, H->modulus); + + // Appropriate shift + slong shift[n+m]; + + for (i = 0; i < n+m; i++) + shift[i]=0; + + // Shifted approximant computation + + nmod_poly_mat_pmbasis(B, shift, M, sigma); + + + slong rows[n+m]; + slong nbrows=0; + + for (i = 0; i < n+m; i++) + { + if (shift[i] <= delta) + { + rows[nbrows]=i; + nbrows +=1; + } + } + + if (nbrows != n) + { + flint_printf("\nA complete description of degree at most %{slong} probably doesn't exist\n", delta); + } + + if (nbrows == 0) + return 0; + else + { + for (i = 0; i < nbrows; i++) + { + for (j = 0; j < m; j++) + nmod_poly_set(nmod_poly_mat_entry(N, i, j), nmod_poly_mat_entry(B, rows[i], n+j)); + } + + for (i = 0; i < nbrows; i++) + { + for (j = 0; j < n; j++) + nmod_poly_set(nmod_poly_mat_entry(D, i, j), nmod_poly_mat_entry(B, rows[i], j)); + } + } -slong nmod_poly_mat_right_description(nmod_poly_mat_t N, nmod_poly_mat_t D, + nmod_poly_mat_clear(M); + nmod_poly_mat_clear(B); + + nmod_poly_clear(mone); + + return nbrows; + +} + + +/** + * Right description computation for H(x) m x n in K(x) (power series) + * with target degree delta + * + * requires enough precision in input: ie at least (m+n)*delta/m +1 + * + * returns nbcols and a partial (or full) description when 0 < nbcols <= n rows, + * or zero if no candidates + * + * N m x n and and D n x n are initialized outside + * hence the result is part of N and D + * + */ + + +slong nmod_poly_mat_right_description(nmod_poly_mat_t N, nmod_poly_mat_t D, const nmod_poly_mat_t H, - slong delta); + slong delta) + +{ + slong nbcols; + + nmod_poly_mat_t NT; + nmod_poly_mat_init(NT, H->c, H->r, H->modulus); + nmod_poly_mat_t DT; + nmod_poly_mat_init(DT, H->c, H->c, H->modulus); -int nmod_poly_mat_kernel(nmod_poly_mat_t N, const nmod_poly_mat_t M, - slong delta); + nmod_poly_mat_t HT; + nmod_poly_mat_init(HT, H->c, H->r, H->modulus); -#ifdef __cplusplus + nmod_poly_mat_transpose(HT,H); + + nbcols = nmod_poly_mat_left_description(NT,DT,HT,delta); + + if (nbcols == 0) + return nbcols; + + nmod_poly_mat_transpose(N,NT); + nmod_poly_mat_transpose(D,DT); + + + // printf("\n"); + // nmod_poly_mat_print_pretty(N, "x"); + // printf("\n"); + // nmod_poly_mat_print_pretty(D, "x"); + // printf("\n"); + + nmod_poly_mat_clear(NT); + nmod_poly_mat_clear(DT); + nmod_poly_mat_clear(HT); + + return nbcols; + } -#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 +/** + * + * Right kernel computation for M(x) n x m in K[x] + * with target degree delta + * + * N is intitialized inside for correct dimensions, n x nbnull + * returns the number nbnull and N, or zero (N not initialized) +* +*/ +int nmod_poly_mat_kernel(nmod_poly_mat_t N, const nmod_poly_mat_t M, slong delta) +{ + + // printf("\n"); + // nmod_poly_mat_print_pretty(M, "x"); + // printf("\n"); + slong i,j; + + slong n = M->r; + slong m = M-> c; + + + nmod_poly_mat_t A; + nmod_poly_mat_init(A, n, n, M->modulus); + + nmod_poly_mat_t B; + nmod_poly_mat_init(B, n, m-n, M->modulus); + + + for (i = 0; i < n; i++) + { + for (j = 0; j < n; j++) + nmod_poly_set(nmod_poly_mat_entry(A, i, j), nmod_poly_mat_entry(M, i, j)); + } + + for (i = 0; i < n; i++) + { + for (j = 0; j < m-n; j++) + nmod_poly_set(nmod_poly_mat_entry(B, i, j), nmod_poly_mat_entry(M, i, n+j)); + } + + slong sigma; + sigma = ceil((double) m*delta/n +1); // sigma >= ceil((double) (rdim + Bcdim)*delta/rdim +1); + + nmod_poly_mat_t X; + nmod_poly_mat_init(X, n, m-n, M->modulus); + + slong order = nmod_poly_mat_degree(A); + + nmod_poly_mat_dixon(X, A, B, order, sigma); + + // printf("\n"); + // nmod_poly_mat_print_pretty(A, "x"); + // printf("\n"); + // nmod_poly_mat_print_pretty(B, "x"); + // printf("sigma %ld\n",sigma); + // printf("\n"); + + + nmod_poly_mat_t P; + nmod_poly_mat_init(P, n, m-n, A->modulus); + + nmod_poly_mat_t Q; + nmod_poly_mat_init(Q, m-n, m-n, A->modulus); + + slong nbnull; + nbnull=nmod_poly_mat_right_description(P, Q, X, delta); + + + if (nbnull != 0) + { + nmod_poly_mat_init(N, m, nbnull, A->modulus); + + + for (i = 0; i < n; i++) + { + for (j = 0; j < nbnull; j++) + nmod_poly_set(nmod_poly_mat_entry(N, i, j), nmod_poly_mat_entry(P, i, j)); + } + + + for (i = 0; i < m-n; i++) + { + for (j = 0; j < nbnull; j++) + { + nmod_poly_set(nmod_poly_mat_entry(N, n+i, j), nmod_poly_mat_entry(Q, i, j)); + nmod_poly_neg(nmod_poly_mat_entry(N, n+i, j), nmod_poly_mat_entry(N, n+i, j)); + } + } + + nmod_poly_mat_t T; + nmod_poly_mat_init(T, n, nbnull, M->modulus); + + nmod_poly_mat_mul(T,M,N); + + // printf("--------- %ld \n",nbnull); + // nmod_poly_mat_print_pretty(N, "x"); + // printf("\n"); + + if (nmod_poly_mat_is_zero(T) == 0) + { + nmod_poly_mat_clear(N); + nbnull = 0; + } + + nmod_poly_mat_clear(T); + } + + + nmod_poly_mat_clear(A); + nmod_poly_mat_clear(B); + nmod_poly_mat_clear(X); + nmod_poly_mat_clear(P); + nmod_poly_mat_clear(Q); + + 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 From 739d7c61edbbc9edef9168c4dce10488a1cf94d8 Mon Sep 17 00:00:00 2001 From: Gilles Villard Date: Sun, 14 Dec 2025 17:02:43 +0100 Subject: [PATCH 3/8] sources .c for kernel --- .../nmod_poly_mat_description.bak | 300 +++++++++ .../nmod_poly_mat_kernel.c | 619 ++++++++++++++++++ 2 files changed, 919 insertions(+) create mode 100644 flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_description.bak create mode 100644 flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_kernel.c 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 new file mode 100644 index 00000000..b95e9013 --- /dev/null +++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_description.bak @@ -0,0 +1,300 @@ +#include +#include +#include +#include + + +#include "nmod_poly_mat_extra.h" +#include "nmod_poly_mat_io.h" + +#include "nmod_mat_extra.h" +#include "nmod_poly_mat_dixon.h" +#include "nmod_poly_mat_description.h" + + +/** + * Left description computation for H(x) n x m in K(x) (power series) + * with target degree delta + * + * requires enough precision in input: ie at least (m+n)*delta/m +1 + * + * returns nbrows and a partial (or full) description when 0 < nbrows <= n rows, + * or zero if no candidates + * + * N n x m and and D n x n are initialized outside + * hence the result is part of N and D + * + */ + +slong nmod_poly_mat_left_description(nmod_poly_mat_t N, nmod_poly_mat_t D, + const nmod_poly_mat_t H, + slong delta) + +{ + slong i,j; + + slong n = H->r; + slong m = H->c; + + slong sigma = ceil((double) (m+n)*delta/m +1); + + nmod_poly_mat_t M; + nmod_poly_mat_init(M, n+m, m, H->modulus); + + nmod_poly_mat_set(M, H); // GV, correct even the dimension of M is larger ? + + nmod_poly_t mone; + nmod_poly_init(mone, H->modulus); + nmod_poly_set_coeff_ui(mone, 0, H->modulus -1); // 1 + + for (i = 0; i < m; i++) + nmod_poly_set(nmod_poly_mat_entry(M, n+i, i), mone); + + // printf("\n"); + // nmod_poly_mat_print_pretty(M, "x"); + // printf("\n"); + + nmod_poly_mat_t B; + nmod_poly_mat_init(B, n+m, n+m, H->modulus); + + // Appropriate shift + slong shift[n+m]; + + for (i = 0; i < n+m; i++) + shift[i]=0; + + // Shifted approximant computation + + nmod_poly_mat_pmbasis(B, shift, M, sigma); + + + slong rows[n+m]; + slong nbrows=0; + + for (i = 0; i < n+m; i++) + { + if (shift[i] <= delta) + { + rows[nbrows]=i; + nbrows +=1; + } + } + + if (nbrows != n) + { + flint_printf("\nA complete description of degree at most %{slong} probably doesn't exist\n", delta); + } + + if (nbrows == 0) + return 0; + else + { + for (i = 0; i < nbrows; i++) + { + for (j = 0; j < m; j++) + nmod_poly_set(nmod_poly_mat_entry(N, i, j), nmod_poly_mat_entry(B, rows[i], n+j)); + } + + + for (i = 0; i < nbrows; i++) + { + for (j = 0; j < n; j++) + nmod_poly_set(nmod_poly_mat_entry(D, i, j), nmod_poly_mat_entry(B, rows[i], j)); + } + } + + nmod_poly_mat_clear(M); + nmod_poly_mat_clear(B); + + nmod_poly_clear(mone); + + return nbrows; + +} + + +/** + * Right description computation for H(x) m x n in K(x) (power series) + * with target degree delta + * + * requires enough precision in input: ie at least (m+n)*delta/m +1 + * + * returns nbcols and a partial (or full) description when 0 < nbcols <= n rows, + * or zero if no candidates + * + * N m x n and and D n x n are initialized outside + * hence the result is part of N and D + * + */ + + +slong nmod_poly_mat_right_description(nmod_poly_mat_t N, nmod_poly_mat_t D, + const nmod_poly_mat_t H, + slong delta) + +{ + slong nbcols; + + nmod_poly_mat_t NT; + nmod_poly_mat_init(NT, H->c, H->r, H->modulus); + + nmod_poly_mat_t DT; + nmod_poly_mat_init(DT, H->c, H->c, H->modulus); + + nmod_poly_mat_t HT; + nmod_poly_mat_init(HT, H->c, H->r, H->modulus); + + + nmod_poly_mat_transpose(HT,H); + + nbcols = nmod_poly_mat_left_description(NT,DT,HT,delta); + + if (nbcols == 0) + return nbcols; + + nmod_poly_mat_transpose(N,NT); + nmod_poly_mat_transpose(D,DT); + + + // printf("\n"); + // nmod_poly_mat_print_pretty(N, "x"); + // printf("\n"); + // nmod_poly_mat_print_pretty(D, "x"); + // printf("\n"); + + nmod_poly_mat_clear(NT); + nmod_poly_mat_clear(DT); + nmod_poly_mat_clear(HT); + + return nbcols; + +} + + +/** + * + * Right kernel computation for M(x) n x m in K[x] + * with target degree delta + * + * N is intitialized inside for correct dimensions, n x nbnull + * returns the number nbnull and N, or zero (N not initialized) +* +*/ + +int nmod_poly_mat_kernel(nmod_poly_mat_t N, const nmod_poly_mat_t M, slong delta) +{ + + // printf("\n"); + // nmod_poly_mat_print_pretty(M, "x"); + // printf("\n"); + + slong i,j; + + slong n = M->r; + slong m = M-> c; + + + nmod_poly_mat_t A; + nmod_poly_mat_init(A, n, n, M->modulus); + + nmod_poly_mat_t B; + nmod_poly_mat_init(B, n, m-n, M->modulus); + + + for (i = 0; i < n; i++) + { + for (j = 0; j < n; j++) + nmod_poly_set(nmod_poly_mat_entry(A, i, j), nmod_poly_mat_entry(M, i, j)); + } + + for (i = 0; i < n; i++) + { + for (j = 0; j < m-n; j++) + nmod_poly_set(nmod_poly_mat_entry(B, i, j), nmod_poly_mat_entry(M, i, n+j)); + } + + slong sigma; + sigma = ceil((double) m*delta/n +1); // sigma >= ceil((double) (rdim + Bcdim)*delta/rdim +1); + + nmod_poly_mat_t X; + nmod_poly_mat_init(X, n, m-n, M->modulus); + + slong order = nmod_poly_mat_degree(A); + + nmod_poly_mat_dixon(X, A, B, order, sigma); + + // printf("\n"); + // nmod_poly_mat_print_pretty(A, "x"); + // printf("\n"); + // nmod_poly_mat_print_pretty(B, "x"); + // printf("sigma %ld\n",sigma); + // printf("\n"); + + + nmod_poly_mat_t P; + nmod_poly_mat_init(P, n, m-n, A->modulus); + + nmod_poly_mat_t Q; + nmod_poly_mat_init(Q, m-n, m-n, A->modulus); + + slong nbnull; + nbnull=nmod_poly_mat_right_description(P, Q, X, delta); + + + if (nbnull != 0) + { + nmod_poly_mat_init(N, m, nbnull, A->modulus); + + + for (i = 0; i < n; i++) + { + for (j = 0; j < nbnull; j++) + nmod_poly_set(nmod_poly_mat_entry(N, i, j), nmod_poly_mat_entry(P, i, j)); + } + + + for (i = 0; i < m-n; i++) + { + for (j = 0; j < nbnull; j++) + { + nmod_poly_set(nmod_poly_mat_entry(N, n+i, j), nmod_poly_mat_entry(Q, i, j)); + nmod_poly_neg(nmod_poly_mat_entry(N, n+i, j), nmod_poly_mat_entry(N, n+i, j)); + } + } + + nmod_poly_mat_t T; + nmod_poly_mat_init(T, n, nbnull, M->modulus); + + nmod_poly_mat_mul(T,M,N); + + // printf("--------- %ld \n",nbnull); + // nmod_poly_mat_print_pretty(N, "x"); + // printf("\n"); + + if (nmod_poly_mat_is_zero(T) == 0) + { + nmod_poly_mat_clear(N); + nbnull = 0; + } + + nmod_poly_mat_clear(T); + } + + + nmod_poly_mat_clear(A); + nmod_poly_mat_clear(B); + nmod_poly_mat_clear(X); + nmod_poly_mat_clear(P); + nmod_poly_mat_clear(Q); + + 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_kernel.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_kernel.c new file mode 100644 index 00000000..bafc01dc --- /dev/null +++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_kernel.c @@ -0,0 +1,619 @@ +#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 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 + * sdeg[n] initialized outside + * + */ + + +void _nmod_poly_mat_sort_permute_columns_zls(nmod_poly_mat_t M, slong *sdeg, \ + slong *perm, const slong *ishift) +{ + + slong n = M->c; + slong j; + + 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 + * + * Output: + * 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 + * degN[n], initialized outside, its first w entries are concerned, + * 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, \ + const slong *ishift, const slong kappa) +{ + + slong i,j,k; + + slong m = A->r; + slong n = A->c; + + + slong min_mn; + slong rho=0; + + // Tuning the order of the subsequent approximant + // ---------------------------------------------- + + if (m <= n) + { + min_mn=m; + for (i=n-m; i n %ld %ld",m,n); // To see, rectangular case + } + + slong s; + s = ceil((double) rho/min_mn); + + + // 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); + nmod_poly_mat_transpose(AT,A); + + nmod_poly_mat_t PT; + nmod_poly_mat_init(PT, n, n, A->modulus); + + // Approximant PT + // shift is modified in place + // -------------------------- + + slong shift[n]; + 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 + // nonzero residues matrix P2 n x n2 + //----------------------------------- + + slong n1=0; + slong n2; + + for (j=0; j0) { + nmod_poly_mat_init(P1, n, n1, A->modulus); + + k=0; + for (j = 0; j < n; j++) + { + if (cdeg[j]<0) { + + for (i = 0; i < n; i++) + nmod_poly_set(nmod_poly_mat_entry(P1, i, k), nmod_poly_mat_entry(PT, j, i)); + k+=1; + } + } + } + + n2=n-n1; + + // the kernel is found + if (n2==0) { + nmod_poly_mat_init_set(N,P1); + nmod_poly_mat_column_degree(degN, P1, ishift); + + nmod_poly_mat_clear(P1); + return n1; + } + + nmod_poly_mat_t P2; + nmod_poly_mat_init(P2, n, n2, A->modulus); + + k=0; + for (j = 0; j < n; j++) + { + if (cdeg[j]>=0) { + + for (i = 0; i < n; i++) + nmod_poly_set(nmod_poly_mat_entry(P2, i, k), nmod_poly_mat_entry(PT, j, i)); + k+=1; + } + + } + + nmod_poly_mat_clear(PT); + + // Special case m=1 + // before the divide and conquer on m + // the kernel is found + // ----------------------------------- + + // !!!!! CHECK + if (m==1){ // Then n2 should be 1, thanks to the choice of kappa + + if (n1==0) + return 0; + else { + nmod_poly_mat_init_set(N,P1); + nmod_poly_mat_column_degree(degN, P1, ishift); + + nmod_poly_mat_clear(P1); + return n1; + } + } + + // Now n2 <> 0 and m > 1 + // one can proceed to the divide and conquer from the rows + // of the nonzero residue P2, and of A.P2 + // -------------------------------------------------------- + + slong * perm = flint_malloc(n2 * sizeof(slong)); + + _nmod_poly_mat_sort_permute_columns_zls(P2,degN,perm,ishift); + + for (i = 0; i < n2; i++) { + degN[i]=degN[i]-kappa*s; // used below for the recursive calls + } + + + nmod_poly_mat_t TT; + 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 + k=0; + for (j = 0; j < n; j++) + { + if (cdeg[j]>=0) { + + for (i = 0; i < m; i++) + nmod_poly_set(nmod_poly_mat_entry(TT, i, k), nmod_poly_mat_entry(RT, j, i)); + k+=1; + } + + } + nmod_poly_mat_clear(RT); + + nmod_poly_mat_permute_columns(TT, perm, NULL); + + nmod_poly_mat_t G; + nmod_poly_mat_init(G, m, n2, A->modulus); + + nmod_poly_mat_shift_right(G,TT,kappa*s); + + + // We split G for the recursive call + // --------------------------------- + + slong new_m=floor((double) m/2); + + nmod_poly_mat_t G1; + nmod_poly_mat_init(G1, new_m, n2, A->modulus); + + for (i = 0; i < new_m; i++){ + for (j = 0; j < n2; j++) { + nmod_poly_set(nmod_poly_mat_entry(G1, i, j), nmod_poly_mat_entry(G, i, j)); + } + } + + nmod_poly_mat_t G2; + nmod_poly_mat_init(G2, m-new_m, n2, A->modulus); + + for (i = 0; i < m-new_m; i++) { + for (j = 0; j < n2; j++) { + nmod_poly_set(nmod_poly_mat_entry(G2, i, j), nmod_poly_mat_entry(G, i+new_m, j)); + } + } + + nmod_poly_mat_clear(TT); + nmod_poly_mat_clear(G); + + + // Recursive calls + // --------------- + + for (i=0; imodulus); + + nmod_poly_mat_mul(G3, G2, N1); + + for (i=0; imodulus); + + nmod_poly_mat_mul(Q1, P2, N1); + + nmod_poly_mat_t Q; + nmod_poly_mat_init(Q, n, c2, A->modulus); + + nmod_poly_mat_mul(Q, Q1, N2); + + nmod_poly_mat_clear(N1); + nmod_poly_mat_clear(N2); + 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_column_degree(degN, Q, ishift); + + nmod_poly_mat_clear(Q); + return c2; + + } + else { + nmod_poly_mat_init(N, n, n1+c2, A->modulus); + + for (i = 0; i < n; i++) { + for (j = 0; j < n1; j++) { + nmod_poly_set(nmod_poly_mat_entry(N, i, j), nmod_poly_mat_entry(P1,i,j)); + } + } + + for (i = 0; i < n; i++) { + for (j = 0; j < c2; j++) { + nmod_poly_set(nmod_poly_mat_entry(N, i, j+n1), nmod_poly_mat_entry(Q, i, j)); + } + } + + slong odeg[n1+c2]; // Too big, to see + + nmod_poly_mat_column_degree(odeg, P1, ishift); + + for (i=0; i= 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 + * degN[n], initialized outside, its first w entries are concerned, + * they are the ishift shifted degrees of the kernel basis + * + */ + +/** + * 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, \ + const slong *ishift, const slong kappa) +{ + + slong i,j,k; + + slong m = iA->r; + slong n = iA->c; + + + nmod_poly_mat_t A; + nmod_poly_mat_init_set(A, iA); + + slong * perm = flint_malloc(n * sizeof(slong)); + slong sdeg[n]; + + // 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; + } + + flint_free(perm); + nmod_poly_mat_clear(A); + + return 0; + +} + + +/** + * Experimental, should not be really considered + * + */ + +int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat_t A, \ + const slong *ishift) +{ + + slong i,j,k; + + slong m = A->r; + slong n = A->c; + + + slong delta; + slong deg[n]; + + nmod_poly_mat_column_degree(deg, A, NULL); + + slong degmax=deg[0]; + for (i=1; im) + min_nm=m; + + nmod_poly_mat_t AT; + nmod_poly_mat_init(AT, n, m, A->modulus); + nmod_poly_mat_transpose(AT,A); + + nmod_poly_mat_t PT; + nmod_poly_mat_init(PT, n, n, A->modulus); + + + // Approximant PT + // shift is modified in place + // -------------------------- + + slong shift[n]; + 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); + nmod_poly_mat_clear(RT); + + + // Zero residue matrix P1 + //----------------------- + + slong n1=0; + + for (j=0; j 0) { + + k=0; + + for (j=0; js,f0,{0,g0,(0,\:0,t0,+0,=s From 2441dc4c90ff27d5ce44af09857a703447c42721 Mon Sep 17 00:00:00 2001 From: Gilles Villard Date: Sun, 14 Dec 2025 17:08:57 +0100 Subject: [PATCH 4/8] .c --- .../src/nmod_poly_mat_extra/nmod_poly_mat_utils.c | 15 +++++++++++++++ 1 file changed, 15 insertions(+) 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 e0fc3a7f..b2eef20a 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 @@ -151,6 +151,21 @@ void _nmod_poly_mat_permute_rows_by_sorting_vec(nmod_poly_mat_t mat, } +void _nmod_poly_mat_permute_columns_by_sorting_vec(nmod_poly_mat_t mat, + slong r, + slong * vec, + slong * perm) +{ + slong_pair * tmp = flint_malloc(r * sizeof(slong_pair)); + _vec_sort_permutation(perm, vec, vec, r, tmp); + for (slong i = r; i < mat->c; i++) + perm[i] = i; + flint_free(tmp); + nmod_poly_mat_permute_columns(mat, perm, NULL); +} + + + /*------------------------------------------------------------*/ /*------------------------------------------------------------*/ /* SET FROM MATRIX POLYNOMIAL */ From 95c0cc197455c14d44f2ec719be7addbfb1a2c07 Mon Sep 17 00:00:00 2001 From: Gilles Villard Date: Sun, 14 Dec 2025 17:14:50 +0100 Subject: [PATCH 5/8] test files --- .../src/nmod_poly_mat_extra/test/main.c | 2 + .../src/nmod_poly_mat_extra/test/t-dixon.c | 7 +- .../src/nmod_poly_mat_extra/test/t-kernel.c | 187 ++++++++++++++++++ 3 files changed, 194 insertions(+), 2 deletions(-) create mode 100644 flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c diff --git a/flint-extras/src/nmod_poly_mat_extra/test/main.c b/flint-extras/src/nmod_poly_mat_extra/test/main.c index 75661fb7..165b48de 100644 --- a/flint-extras/src/nmod_poly_mat_extra/test/main.c +++ b/flint-extras/src/nmod_poly_mat_extra/test/main.c @@ -24,6 +24,7 @@ #include "t-mul_waksman.c" #include "t-rand.c" #include "t-weak_popov_form.c" +#include "t-kernel.c" /* Array of test functions ***************************************************/ @@ -39,6 +40,7 @@ test_struct tests[] = /* 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 691aba09..ababa3ba 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 @@ -144,9 +144,9 @@ slong prime[] = slong Bcdim[] = {1, 10}; - slong degA[] = {1, 20}; + slong degA[] = {1, 12}; - slong degB[] = {1, 2, 20}; + slong degB[] = {1, 8}; slong order; slong sigma; @@ -246,3 +246,6 @@ 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.c new file mode 100644 index 00000000..c773de0a --- /dev/null +++ b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c @@ -0,0 +1,187 @@ +/* + Copyright (C) 2025 Gilles Villard + + This file is part of PML. + + PML is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License version 2.0 (GPL-2.0-or-later) + as published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. See + . +*/ + +#include +#include +#include +#include + +#include "nmod_poly_mat_utils.h" +#include "nmod_poly_mat_extra.h" + +// test one given input +int core_test_kernel(const nmod_poly_mat_t mat) +{ + + slong m=mat->r; + slong n=mat->c; + + // Flint rank + // ---------- + + slong rkflint; + + rkflint=nmod_poly_mat_rank(mat); + + + // PML nullspace + // ------------- + + 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); + + 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 Z; + nmod_poly_mat_init(Z, m, nz, mat->modulus); + + nmod_poly_mat_mul(Z, mat, NN); + + verif = verif && nmod_poly_mat_is_zero(Z); + + nmod_poly_mat_clear(Z); + nmod_poly_mat_clear(NN); + + } + + return verif; +} + + +TEST_FUNCTION_START(nmod_poly_mat_kernel, 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 prime = n_randprime(state, nbits, 1); + + + nmod_poly_mat_t A; + + 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) { + 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) { + nmod_poly_mat_init(A, rdim, cdim, prime); + nmod_poly_mat_randtest_sparse(A, state, deg+1, 0.2); + } + 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); + + nmod_poly_mat_clear(A); + + if (!result) { + TEST_FUNCTION_FAIL("rdim = %wu, cdim = %wu, degree = %wu, p = %wu\n", \ + rdim, cdim, deg, prime); + } + } + + for (i = 0; i < 2; i++) + { + ulong rdim = 40 + n_randint(state, 8); + ulong cdim = rdim + n_randint(state, 8); + ulong deg = 0; + + 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); + + nmod_poly_mat_clear(A); + + if (!result) { + TEST_FUNCTION_FAIL("rdim = %wu, cdim = %wu, degree = %wu, p = %wu\n", \ + rdim, cdim, deg, prime); + } + } + + for (i = 0; i < 2; i++) + { + ulong rdim = 20 + n_randint(state, 8); + ulong cdim = rdim -8; + ulong deg = 1; + + 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); + + nmod_poly_mat_clear(A); + + if (!result) { + TEST_FUNCTION_FAIL("rdim = %wu, cdim = %wu, degree = %wu, p = %wu\n", \ + rdim, cdim, deg, prime); + } + } + +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 From 3a108ccf8dcf3234eedb7d330ca05402dcd39f4e Mon Sep 17 00:00:00 2001 From: Gilles Villard Date: Sun, 14 Dec 2025 18:02:38 +0100 Subject: [PATCH 6/8] unused variables --- flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_kernel.c | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_kernel.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_kernel.c index bafc01dc..49059994 100644 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_kernel.c @@ -448,9 +448,8 @@ int nmod_poly_mat_kernel(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat_t i const slong *ishift, const slong kappa) { - slong i,j,k; + slong j,k; - slong m = iA->r; slong n = iA->c; @@ -531,8 +530,6 @@ int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong *degN, const nmod_ slong m = A->r; slong n = A->c; - - slong delta; slong deg[n]; nmod_poly_mat_column_degree(deg, A, NULL); From f219abee00f2dafdefc8636cd7864fe460792288 Mon Sep 17 00:00:00 2001 From: Gilles Villard Date: Sun, 14 Dec 2025 18:14:05 +0100 Subject: [PATCH 7/8] Check tests --- flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c index c773de0a..99e54e67 100644 --- a/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c @@ -134,6 +134,9 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state) } } + // Degree zero + // ----------- + for (i = 0; i < 2; i++) { ulong rdim = 40 + n_randint(state, 8); @@ -157,6 +160,9 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state) } } + // Column dimension smaller that row dimension + // ------------------------------------------- + for (i = 0; i < 2; i++) { ulong rdim = 20 + n_randint(state, 8); From aedfa16f42e2ad4ec380858d94b6e0d2d927785c Mon Sep 17 00:00:00 2001 From: Gilles Villard Date: Fri, 19 Dec 2025 12:37:23 +0100 Subject: [PATCH 8/8] Bug fixed for correct pull request --- .../nmod_poly_mat_kernel.c | 113 ++++++++++++++---- flint-extras/src/nmod_poly_mat_kernel.h | 4 +- 2 files changed, 93 insertions(+), 24 deletions(-) diff --git a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_kernel.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_kernel.c index 49059994..084aee45 100644 --- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_kernel.c +++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_kernel.c @@ -1,4 +1,5 @@ #include +#include #include "nmod_poly_mat_extra.h" #include "nmod_poly_mat_kernel.h" @@ -29,6 +30,18 @@ */ +// Dummy function for qsort + int cmp(const void *a, const void *b) +{ + const slong *x = a; + const slong *y = b; + + if (*x < *y) return -1; + if (*x > *y) return 1; + return 0; +} + + void _nmod_poly_mat_sort_permute_columns_zls(nmod_poly_mat_t M, slong *sdeg, \ slong *perm, const slong *ishift) { @@ -65,7 +78,7 @@ void _nmod_poly_mat_sort_permute_columns_zls(nmod_poly_mat_t M, slong *sdeg, \ * 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, an integer >= 2, for the order of the order bases + * kappa, a double >= 2, for the order of the order bases * kappa * s instead of 3 *s in ZLS * * Output: @@ -79,7 +92,7 @@ void _nmod_poly_mat_sort_permute_columns_zls(nmod_poly_mat_t M, slong *sdeg, \ */ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat_t A, \ - const slong *ishift, const slong kappa) + const slong *ishift, const double kappa) { slong i,j,k; @@ -94,24 +107,36 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat // Tuning the order of the subsequent approximant // ---------------------------------------------- + long shift[n]; // temporary variable + + + // 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 } + slong s; s = ceil((double) rho/min_mn); - // Transposed A for the approximant PT on the left // PT will then be use without transposing @@ -126,13 +151,14 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat // shift is modified in place // -------------------------- - slong shift[n]; for (i=0; imodulus); + + // nmod_poly_mat_mul(G, A, P2); + + // nmod_poly_mat_t TT; + // nmod_poly_mat_init(TT, m, n2, A->modulus); + + + // nmod_poly_mat_shift_right(TT,G,ks); + + + // slong new_m=floor((double) m/2); + + // nmod_poly_mat_t G1; + // nmod_poly_mat_init(G1, new_m, n2, A->modulus); + + // for (i = 0; i < new_m; i++){ + // for (j = 0; j < n2; j++) { + // nmod_poly_set(nmod_poly_mat_entry(G1, i, j), nmod_poly_mat_entry(TT, i, j)); + // } + // } + + // nmod_poly_mat_t G2; + // nmod_poly_mat_init(G2, m-new_m, n2, A->modulus); + + // for (i = 0; i < m-new_m; i++) { + // for (j = 0; j < n2; j++) { + // nmod_poly_set(nmod_poly_mat_entry(G2, i, j), nmod_poly_mat_entry(TT, i+new_m, j)); + // } + // } + + + //++++++++ END OLD + + + //++++++++++ NEW ++++++++++++++ nmod_poly_mat_t TT; nmod_poly_mat_init(TT, m, n2, A->modulus); @@ -260,7 +329,7 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat nmod_poly_mat_t G; nmod_poly_mat_init(G, m, n2, A->modulus); - nmod_poly_mat_shift_right(G,TT,kappa*s); + nmod_poly_mat_shift_right(G,TT,ks); // We split G for the recursive call @@ -288,22 +357,20 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat nmod_poly_mat_clear(TT); nmod_poly_mat_clear(G); - + //++++++++++++++ END NEW ++++++++++++++++++ // Recursive calls // --------------- - for (i=0; imodulus); @@ -383,7 +452,7 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat } } - slong odeg[n1+c2]; // Too big, to see + slong odeg[n1+c2]; nmod_poly_mat_column_degree(odeg, P1, ishift); @@ -428,7 +497,7 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat * 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, an integer >= 2, for the order of the order bases + * kappa, a double >= 2, for the order of the order bases * kappa * s instead of 3 *s in ZLS * * Output: @@ -445,7 +514,7 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat */ int nmod_poly_mat_kernel(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat_t iA, \ - const slong *ishift, const slong kappa) + const slong *ishift, const double kappa) { slong j,k; diff --git a/flint-extras/src/nmod_poly_mat_kernel.h b/flint-extras/src/nmod_poly_mat_kernel.h index 8e280841..81b4d0eb 100644 --- a/flint-extras/src/nmod_poly_mat_kernel.h +++ b/flint-extras/src/nmod_poly_mat_kernel.h @@ -39,7 +39,7 @@ extern "C" { * 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, an integer >= 2, for the order of the order bases + * kappa, a double >= 2, for the order of the order bases * kappa * s instead of 3 *s in ZLS * * Output: @@ -52,7 +52,7 @@ extern "C" { */ int nmod_poly_mat_kernel(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat_t A, \ - const slong *ishift, const slong kappa); + const slong *ishift, const double kappa);