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..b95e9013 --- /dev/null +++ b/flint-extras/src/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.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_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..084aee45 --- /dev/null +++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_kernel.c @@ -0,0 +1,685 @@ +#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 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 + * + */ + + +// 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) +{ + + 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 double 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 + // ---------------------------------------------- + + 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 + + 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 + // -------------------------- + + 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 + // ----------------------------------- + + if (m==1){ + + 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)); + + + slong degP2[n2]; + + // 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 + } + + + //+++++++++ OLD + + // nmod_poly_mat_t G; + // nmod_poly_mat_init(G, m, n2, A->modulus); + + // 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); + + // 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,ks); + + + // 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); + //++++++++++++++ END NEW ++++++++++++++++++ + + // Recursive calls + // --------------- + + nmod_poly_mat_t N1; + nmod_poly_mat_t N2; + + slong c1=0; + slong c2=0; + + + 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); + + 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]; + + 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 double kappa) +{ + + slong j,k; + + 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 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 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 */ 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..99e54e67 --- /dev/null +++ b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c @@ -0,0 +1,193 @@ +/* + 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); + } + } + + // Degree zero + // ----------- + + 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); + } + } + + // Column dimension smaller that row dimension + // ------------------------------------------- + + 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 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..81b4d0eb --- /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, 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 + * 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 double 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 */