From be5d4e59393fb698d589cef74dca669c0923d1a8 Mon Sep 17 00:00:00 2001
From: Vincent Neiger <vneiger@users.noreply.github.com>
Date: Fri, 2 Jan 2026 22:08:47 +0100
Subject: [PATCH 1/8] working on p-kernel; adjust some parts to flint style

---
 .../nmod_poly_mat_kernel.c                    |   3 -
 .../nmod_poly_mat_extra/profile/p-kernel.c    | 272 ++++++++++++++++++
 .../src/nmod_poly_mat_extra/profile/p-mul.c   |   2 +-
 .../src/nmod_poly_mat_extra/test/t-kernel.c   |  87 +++---
 4 files changed, 309 insertions(+), 55 deletions(-)
 create mode 100644 flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c

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 084aee45..67d1c71b 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
@@ -680,6 +680,3 @@ int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong *degN, const nmod_
     nmod_poly_mat_clear(PT);
     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/profile/p-kernel.c b/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c
new file mode 100644
index 00000000..9fb194b3
--- /dev/null
+++ b/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c
@@ -0,0 +1,272 @@
+#include <stdlib.h>  // for atoi
+
+#include <flint/ulong_extras.h>
+#include <flint/profiler.h>
+#include <flint/nmod.h>
+#include <flint/nmod_vec.h>
+#include <flint/nmod_poly_mat.h>
+
+#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 */                                       \
+                                                        \
+    nmod_poly_mat_t K;                                  \
+    nmod_poly_mat_init(K, rdim, rdim, n);               \
+    nmod_poly_mat_kernel(K, degN, Ft, ishift, 3.);      \
+                                                        \
+    double FLINT_SET_BUT_UNUSED(tcpu), twall;           \
+                                                        \
+    TIMEIT_START;                                       \
+    nmod_poly_mat_##fun(C, A, B);                       \
+    TIMEIT_STOP_VALUES(tcpu, twall);                    \
+                                                        \
+    printf("%.2e", twall);                              \
+                                                        \
+    nmod_poly_mat_clear(A);                             \
+    nmod_poly_mat_clear(B);                             \
+    nmod_poly_mat_clear(C);                             \
+}
+
+TIME_KER(mul)
+
+/*-------------------------*/
+/*  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 (all square for the moment)
+    const slong ndims = 10;
+    const ulong dims[] = {2, 4, 6, 8, 11, 15, 20, 30, 50, 100};
+
+    // matrix degrees
+    const slong ndegs = 12;
+    const ulong degs[] = {5, 10, 20, 40, 80, 160, 320, 640, 1280, 2560, 5120, 10240};
+
+    // bench functions
+    const slong nfuns = 2;
+    typedef void (*timefun) (time_args, flint_rand_t);
+    const timefun funs[] = {
+        time_mul,                      // 0
+    };
+
+    // TODO
+    //typedef void (*samplefun) (void*, ulong);
+    //const samplefun sfuns[] = {
+    //    sample_mul,                      // 0
+    //    sample_mul_geometric,            // 1
+    //};
+
+    const char * description[] = {
+        "#0  --> mul                          ",
+        "#1  --> mul_geometric                ",
+    };
+
+    if (argc == 1)  // show usage
+    {
+        printf("Usage: `%s [nbits] [rdim] [cdim] [deg] [rank] [fun]`\n", argv[0]);
+        printf("   Each argument is optional; no argument shows this help.\n");
+        printf("   - nbits: number of bits (in (1..64]) for the modulus, chosen as nextprime(2**(nbits-1))\n");
+        printf("        (nbits == -1 launches full suite)\n");
+        printf("   - rdim, cdim: input matrix is rdim x cdim\n");
+        printf("   - deg: matrix is random of degree < deg\n");
+        printf("   - rank: matrix is random of this rank\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++)
+    {
+        time_args targs = {4, 4, 4, 1000, UWORD(1) << 20};
+        time_mul(targs, state);
+        printf(" ");
+    }
+    printf("\n\n");
+
+    if (argc == 2 && atoi(argv[1]) == -1)  // launching full suite
+    {
+        printf("           dim");
+        for (slong i = 0; i < ndims; i++)
+            printf("%17ld", dims[i]);
+        printf("\n");
+        printf("bits fun deg\n");
+        for (slong j = 0; j < nbits; j++)
+        {
+            const slong b = bits[j];
+            ulong n;
+            n = n_nextprime(UWORD(1) << (b-1), 0);
+            for (slong ifun = 0; ifun < nfuns; ifun++)
+            {
+                for (slong d = 0; d < ndegs; d++)
+                {
+                    printf("%-5ld#%-3ld%-8ld", b, ifun, degs[d]);
+                    for (slong i = 0; i < ndims; i++)
+                    {
+                        time_args targs = {dims[i], dims[i], dims[i], degs[d], n};
+
+#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
+                        printf(" ");
+                    }
+                }
+                printf("\n");
+            }
+        }
+    }
+    else if (argc == 2)  // nbits is given
+    {
+        printf("       dim");
+        for (slong i = 0; i < ndims; i++)
+            printf("%17ld", dims[i]);
+        printf("\n");
+        printf("bits fun deg\n");
+        const slong b = atoi(argv[1]);
+        ulong n;
+        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%-8ld", b, ifun, degs[d]);
+                for (slong i = 0; i < ndims; i++)
+                {
+                    time_args targs = {dims[i], dims[i], dims[i], degs[d], n};
+                    tfun(targs, state);
+                    printf(" ");
+                }
+            }
+            printf("\n");
+        }
+    }
+    else if (argc == 3)  // nbits + dim given
+    {
+        const slong dim = atoi(argv[2]);
+        printf("       dim");
+        printf("%17ld", dim);
+        printf("\n");
+        printf("bits fun deg\n");
+        const slong b = atoi(argv[1]);
+        ulong n;
+        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%-8ld", b, ifun, degs[d]);
+                time_args targs = {dim, dim, dim, degs[d], n};
+                tfun(targs, state);
+                printf(" ");
+                printf("\n");
+            }
+        }
+    }
+    else if (argc == 4)  // nbits + dim + deg given
+    {
+        const slong dim = atoi(argv[2]);
+        const slong deg = atoi(argv[3]);
+        printf("       dim");
+        printf("%17ld", dim);
+        printf("\n");
+        printf("bits fun deg\n");
+        const slong b = atoi(argv[1]);
+        ulong n;
+        n = n_nextprime(UWORD(1) << (b-1), 0);
+        for (slong ifun = 0; ifun < nfuns; ifun++)
+        {
+            const timefun tfun = funs[ifun];
+            printf("%-5ld#%-3ld%-8ld", b, ifun, deg);
+            time_args targs = {dim, dim, dim, deg, n};
+            tfun(targs, state);
+            printf(" ");
+            printf("\n");
+        }
+    }
+    else if (argc == 5)  // nbits + dim + deg + fun given
+    {
+        const slong b = atoi(argv[1]);
+        const slong dim = atoi(argv[2]);
+        const slong deg = atoi(argv[3]);
+        const slong ifun = atoi(argv[4]);
+        const timefun tfun = funs[ifun];
+        printf("       dim");
+        printf("%17ld", dim);
+        printf("\n");
+        printf("bits fun deg\n");
+        ulong n;
+        n = n_nextprime(UWORD(1) << (b-1), 0);
+
+        printf("%-5ld#%-3ld%-8ld", b, ifun, deg);
+        time_args targs = {dim, dim, dim, deg, 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/test/t-kernel.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c
index 99e54e67..db3b519b 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
@@ -16,85 +16,67 @@
 #include <flint/test_helpers.h>
 
 #include "nmod_poly_mat_utils.h"
-#include "nmod_poly_mat_extra.h"  
+#include "nmod_poly_mat_extra.h"
 
 // test one given input
+/* TODO does not test generation */
 int core_test_kernel(const nmod_poly_mat_t mat)
 {
+    slong m = mat->r;
+    slong n = mat->c;
 
-    slong m=mat->r;
-    slong n=mat->c;
-    
-    // Flint rank 
+    // Flint rank
     // ----------
-    
-    slong rkflint;
 
-    rkflint=nmod_poly_mat_rank(mat);
-    
+    slong rkflint;
+    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);
+    nz = nmod_poly_mat_kernel(N, degN, mat, NULL, 2);
 
     int verif;
-    verif =  (n-rkflint == nz);
+    verif = (n-rkflint == nz);
 
-    if (nz !=0) {
-
-        nmod_poly_mat_t NN; 
+    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++) {
+            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); 
+        verif = verif && nmod_poly_mat_is_zero(Z);
 
         nmod_poly_mat_clear(Z);
         nmod_poly_mat_clear(NN);
-
     }
-    
-    return verif; 
-}
 
+    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);
@@ -102,39 +84,43 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state)
 
         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);
 
         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,7 +129,7 @@ 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;
 
@@ -154,7 +140,8 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state)
 
         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,7 +156,7 @@ 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;
 
@@ -180,14 +167,12 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state)
 
         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

From 32cc1f68bcb2f5bdeb4ca7c8aa9026334fa18419 Mon Sep 17 00:00:00 2001
From: Vincent Neiger <vneiger@users.noreply.github.com>
Date: Fri, 2 Jan 2026 23:48:20 +0100
Subject: [PATCH 2/8] p-kernel: first version

---
 .../nmod_poly_mat_extra/profile/p-kernel.c    | 192 ++++++------------
 .../src/nmod_poly_mat_extra/test/t-kernel.c   |   1 +
 2 files changed, 60 insertions(+), 133 deletions(-)

diff --git a/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c b/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c
index 9fb194b3..658c449a 100644
--- a/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c
+++ b/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c
@@ -35,8 +35,8 @@ 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 rank = targs.rank; */ /* TODO */           \
+    /* const slong stype = targs.stype; */ /* TODO */         \
     const slong n = targs.modn;                         \
                                                         \
     nmod_t mod;                                         \
@@ -52,24 +52,23 @@ void time_##fun(time_args targs, flint_rand_t state)    \
     nmod_poly_mat_transpose(Ft, F);                     \
     /* END TMP */                                       \
                                                         \
+    slong degN[rdim];                                   \
     nmod_poly_mat_t K;                                  \
     nmod_poly_mat_init(K, rdim, rdim, n);               \
-    nmod_poly_mat_kernel(K, degN, Ft, ishift, 3.);      \
                                                         \
     double FLINT_SET_BUT_UNUSED(tcpu), twall;           \
                                                         \
     TIMEIT_START;                                       \
-    nmod_poly_mat_##fun(C, A, B);                       \
+    nmod_poly_mat_##fun(K, degN, Ft, NULL, 3.);         \
     TIMEIT_STOP_VALUES(tcpu, twall);                    \
                                                         \
     printf("%.2e", twall);                              \
                                                         \
-    nmod_poly_mat_clear(A);                             \
-    nmod_poly_mat_clear(B);                             \
-    nmod_poly_mat_clear(C);                             \
+    nmod_poly_mat_clear(F);                             \
+    nmod_poly_mat_clear(Ft); /* TMP */                  \
 }
 
-TIME_KER(mul)
+TIME_KER(kernel)
 
 /*-------------------------*/
 /*  main                   */
@@ -82,22 +81,26 @@ int main(int argc, char ** argv)
     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};
+    /* const slong nbits = 7; */
+    /* const ulong bits[] = {12, 24, 30, 40, 50, 60, 63}; */
 
-    // matrix dimensions (all square for the moment)
-    const slong ndims = 10;
-    const ulong dims[] = {2, 4, 6, 8, 11, 15, 20, 30, 50, 100};
+    // 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 = 12;
-    const ulong degs[] = {5, 10, 20, 40, 80, 160, 320, 640, 1280, 2560, 5120, 10240};
+    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 = 2;
+    const slong nfuns = 1;
     typedef void (*timefun) (time_args, flint_rand_t);
     const timefun funs[] = {
-        time_mul,                      // 0
+        time_kernel,                      // 0
     };
 
     // TODO
@@ -106,21 +109,32 @@ int main(int argc, char ** argv)
     //    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  --> mul                          ",
-        "#1  --> mul_geometric                ",
+        "#0  --> kernel (general interface)   ",
+        "#1  --> .... TODO                ",
     };
 
-    if (argc == 1)  // show usage
+    if (argc < 4 || argc > 6)  // show usage
     {
-        printf("Usage: `%s [nbits] [rdim] [cdim] [deg] [rank] [fun]`\n", argv[0]);
-        printf("   Each argument is optional; no argument shows this help.\n");
-        printf("   - nbits: number of bits (in (1..64]) for the modulus, chosen as nextprime(2**(nbits-1))\n");
-        printf("        (nbits == -1 launches full suite)\n");
+        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\n");
-        printf("   - rank: matrix is random of this rank\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++)
@@ -132,140 +146,52 @@ int main(int argc, char ** argv)
     printf("#warmup...\n");
     for (slong i = 0; i < 3; i++)
     {
-        time_args targs = {4, 4, 4, 1000, UWORD(1) << 20};
-        time_mul(targs, state);
+        /* rdim; cdim; deg; rank; stype; modn; */
+        time_args targs = {8, 4, 1000, 4, 0, n_nextprime(UWORD(1) << 20, 0)};
+        time_kernel(targs, state);
         printf(" ");
     }
     printf("\n\n");
 
-    if (argc == 2 && atoi(argv[1]) == -1)  // launching full suite
-    {
-        printf("           dim");
-        for (slong i = 0; i < ndims; i++)
-            printf("%17ld", dims[i]);
-        printf("\n");
-        printf("bits fun deg\n");
-        for (slong j = 0; j < nbits; j++)
-        {
-            const slong b = bits[j];
-            ulong n;
-            n = n_nextprime(UWORD(1) << (b-1), 0);
-            for (slong ifun = 0; ifun < nfuns; ifun++)
-            {
-                for (slong d = 0; d < ndegs; d++)
-                {
-                    printf("%-5ld#%-3ld%-8ld", b, ifun, degs[d]);
-                    for (slong i = 0; i < ndims; i++)
-                    {
-                        time_args targs = {dims[i], dims[i], dims[i], degs[d], n};
-
-#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
-                        printf(" ");
-                    }
-                }
-                printf("\n");
-            }
-        }
-    }
-    else if (argc == 2)  // nbits is given
-    {
-        printf("       dim");
-        for (slong i = 0; i < ndims; i++)
-            printf("%17ld", dims[i]);
-        printf("\n");
-        printf("bits fun deg\n");
-        const slong b = atoi(argv[1]);
-        ulong n;
-        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%-8ld", b, ifun, degs[d]);
-                for (slong i = 0; i < ndims; i++)
-                {
-                    time_args targs = {dims[i], dims[i], dims[i], degs[d], n};
-                    tfun(targs, state);
-                    printf(" ");
-                }
-            }
-            printf("\n");
-        }
-    }
-    else if (argc == 3)  // nbits + dim given
+    if (argc == 4)  // nbits + rdim + cdim given
     {
-        const slong dim = atoi(argv[2]);
-        printf("       dim");
-        printf("%17ld", dim);
-        printf("\n");
-        printf("bits fun deg\n");
+        const slong rdim = atoi(argv[2]);
+        const slong cdim = atoi(argv[3]);
         const slong b = atoi(argv[1]);
-        ulong n;
-        n = n_nextprime(UWORD(1) << (b-1), 0);
+        printf("bits fun rdim cdim deg\n");
+        ulong n = n_nextprime(UWORD(1) << (b-1), 0);
         for (slong ifun = 0; ifun < nfuns; ifun++)
         {
             const timefun tfun = funs[ifun];
             for (slong d = 0; d < ndegs; d++)
             {
-                printf("%-5ld#%-3ld%-8ld", b, ifun, degs[d]);
-                time_args targs = {dim, dim, dim, degs[d], n};
+                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 == 4)  // nbits + dim + deg given
+    else if (argc == 5)  // nbits + rdim + cdim + deg given
     {
-        const slong dim = atoi(argv[2]);
-        const slong deg = atoi(argv[3]);
-        printf("       dim");
-        printf("%17ld", dim);
-        printf("\n");
-        printf("bits fun deg\n");
+        const slong rdim = atoi(argv[2]);
+        const slong cdim = atoi(argv[3]);
+        const slong deg  = atoi(argv[4]);
         const slong b = atoi(argv[1]);
-        ulong n;
-        n = n_nextprime(UWORD(1) << (b-1), 0);
+        printf("bits fun rdim cdim deg\n");
+        ulong n = n_nextprime(UWORD(1) << (b-1), 0);
         for (slong ifun = 0; ifun < nfuns; ifun++)
         {
             const timefun tfun = funs[ifun];
-            printf("%-5ld#%-3ld%-8ld", b, ifun, deg);
-            time_args targs = {dim, dim, dim, deg, n};
+            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");
         }
     }
-    else if (argc == 5)  // nbits + dim + deg + fun given
-    {
-        const slong b = atoi(argv[1]);
-        const slong dim = atoi(argv[2]);
-        const slong deg = atoi(argv[3]);
-        const slong ifun = atoi(argv[4]);
-        const timefun tfun = funs[ifun];
-        printf("       dim");
-        printf("%17ld", dim);
-        printf("\n");
-        printf("bits fun deg\n");
-        ulong n;
-        n = n_nextprime(UWORD(1) << (b-1), 0);
-
-        printf("%-5ld#%-3ld%-8ld", b, ifun, deg);
-        time_args targs = {dim, dim, dim, deg, n};
-        tfun(targs, state);
-
-        printf(" ");
-        printf("\n");
-    }
 
     flint_rand_clear(state);
     return 0;
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 db3b519b..9cbc2785 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
@@ -20,6 +20,7 @@
 
 // test one given input
 /* TODO does not test generation */
+/* TODO does not seem tested with shifts */
 int core_test_kernel(const nmod_poly_mat_t mat)
 {
     slong m = mat->r;

From 938e9b9d2285ac99496c5158fa83379730c45a04 Mon Sep 17 00:00:00 2001
From: Vincent Neiger <vneiger@users.noreply.github.com>
Date: Sat, 3 Jan 2026 00:30:56 +0100
Subject: [PATCH 3/8] add is_kernel function

---
 flint-extras/src/nmod_poly_mat_extra.h        | 23 +++++-
 ...ication.c => nmod_poly_mat_verification.c} | 82 ++++++++++++++++---
 .../src/nmod_poly_mat_extra/test/t-kernel.c   |  1 +
 3 files changed, 90 insertions(+), 16 deletions(-)
 rename flint-extras/src/nmod_poly_mat_extra/{nmod_poly_mat_approximant_verification.c => nmod_poly_mat_verification.c} (51%)

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/nmod_poly_mat_approximant_verification.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_verification.c
similarity index 51%
rename from flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_approximant_verification.c
rename to flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_verification.c
index 769d6a37..e1acb0e0 100644
--- 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_verification.c
@@ -11,7 +11,6 @@ int nmod_poly_mat_is_approximant_basis(const nmod_poly_mat_t appbas,
                                        const slong * shift,
                                        orientation_t orient)
 {
-    // context
     const slong rdim = pmat->r;
     const slong cdim = pmat->c;
     const ulong prime = pmat->modulus;
@@ -26,26 +25,26 @@ int nmod_poly_mat_is_approximant_basis(const nmod_poly_mat_t appbas,
 
     int success = 1;
 
-    // check appbas is square with the right dimension
+    /* 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"
+    /* check appbas is shifted reduced */
     if (!nmod_poly_mat_is_ordered_weak_popov(appbas, shift, orient))
     {
-        printf("basis is not shifted-reduced\n");
+        printf("basis is not shifted-weak Popov\n");
         success = 0;
     }
 
-    // compute residual, check rows of appbas are approximants
+    /* compute residual, check rows of appbas are approximants */
     nmod_poly_mat_mul(residual, appbas, pmat);
 
-    for(slong i = 0; i < rdim; i++)
+    for (slong i = 0; i < rdim; i++)
     {
-        for(slong j = 0; j < cdim; j++)
+        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))
@@ -56,10 +55,11 @@ int nmod_poly_mat_is_approximant_basis(const nmod_poly_mat_t appbas,
         }
     }
 
-    // check generation: follows ideas from Algorithm 1 in Giorgi-Neiger, ISSAC 2018
+    /* 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
+    /* 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));
@@ -70,9 +70,10 @@ int nmod_poly_mat_is_approximant_basis(const nmod_poly_mat_t appbas,
         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)
+    /* 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++)
@@ -100,3 +101,58 @@ int nmod_poly_mat_is_approximant_basis(const nmod_poly_mat_t appbas,
 
     return success;
 }
+
+/* TODO currently specialized to ROW_LOWER (or at least ROW_stuff) */
+/* -> checks whether ker is a shift-ordered weak Popov left kernel basis of pmat */
+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");
+        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_extra/test/t-kernel.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c
index 9cbc2785..8f769eb7 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
@@ -20,6 +20,7 @@
 
 // test one given input
 /* TODO does not test generation */
+/* TODO does not test reducedness */
 /* TODO does not seem tested with shifts */
 int core_test_kernel(const nmod_poly_mat_t mat)
 {

From 642b264c778275916f24466a442cd1a34c51badf Mon Sep 17 00:00:00 2001
From: Vincent Neiger <vneiger@users.noreply.github.com>
Date: Sat, 3 Jan 2026 00:59:40 +0100
Subject: [PATCH 4/8] consolidate tests; currently output is not weak Popov

---
 .../nmod_poly_mat_kernel.c                    | 270 +++++++++---------
 .../nmod_poly_mat_verification.c              |   9 +-
 .../src/nmod_poly_mat_extra/test/main.c       |   6 +-
 .../src/nmod_poly_mat_extra/test/t-kernel.c   |  75 +++--
 flint-extras/src/nmod_poly_mat_kernel.h       |  57 ++--
 5 files changed, 205 insertions(+), 212 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 67d1c71b..f54ff663 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,36 +1,36 @@
-#include <math.h> 
+#include <math.h>
 #include <stdlib.h>
 
 #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;
@@ -51,7 +51,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<n; j++) 
+    for (j=0; j<n; j++)
         if (sdeg[j] < 0) sdeg[j]=0;
 
     _nmod_poly_mat_permute_columns_by_sorting_vec(M, n, sdeg, perm);
@@ -59,36 +59,36 @@ void _nmod_poly_mat_sort_permute_columns_zls(nmod_poly_mat_t M, slong *sdeg, \
 }
 
 /**
- *  
- *  Internal function 
- * 
- *  Right shifted kernel of a polynomial matrix, assuming that the columns has been 
- *     sorted by shifted degree 
- * 
+ *
+ *  Internal function
+ *
+ *  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
- * 
- *  TODO/TO SEE: 
- *    
+ *
+ *  TODO/TO SEE:
+ *
  *  Input:
- *    A in m x n 
- *    ishift[n], NULL (the degrees are computed) or initialized outside, 
- *      the shift for the kernel   
+ *    A 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, 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 +104,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; i++)
     {
         shift[i]=ishift[i];
@@ -119,26 +119,26 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat
 
     qsort(shift, n, sizeof(slong), cmp);
 
-    if (m <= n) 
+    if (m <= n)
     {
         min_mn=m;
-        for (i=n-m; i<n; i++)   
-            rho+=shift[i];   
+        for (i=n-m; i<n; i++)
+            rho+=shift[i];
     }
-    else 
-    {   // No sort needed, we consider all degrees 
+    else
+    {   // No sort needed, we consider all degrees
         min_mn=n;
-        for (i=0; i<n; i++)   
-            rho+=ishift[i];  
-        // printf("\n m > n  %ld  %ld",m,n);   // To see, rectangular case 
-    }  
+        for (i=0; i<n; i++)
+            rho+=ishift[i];
+        // printf("\n m > 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 +147,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; i<n; i++)
     {
         shift[i]=ishift[i];
     }
-    
+
     slong ks;
-    ks=ceil((double) kappa*s); 
+    ks=ceil((double) kappa*s);
     nmod_poly_mat_pmbasis(PT, shift, AT, ks+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);
 
     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 +181,7 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat
     slong n2;
 
     for (j=0; j<n; j++) {
-        if (cdeg[j]<0) 
+        if (cdeg[j]<0)
             n1+=1;
     }
 
@@ -201,10 +201,10 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat
             }
         }
     }
-    
+
     n2=n-n1;
 
-    // the kernel is found 
+    // the kernel is found
     if (n2==0) {
         nmod_poly_mat_init_set(N,P1);
         nmod_poly_mat_column_degree(degN, P1, ishift);
@@ -232,13 +232,13 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat
 
     //  Special case m=1
     //  before the divide and conquer on m
-    //  the kernel is found 
+    //  the kernel is found
     // -----------------------------------
 
-    if (m==1){  
+    if (m==1){
 
-        if (n1==0) 
-            return 0; 
+        if (n1==0)
+            return 0;
         else {
 
             nmod_poly_mat_init_set(N,P1);
@@ -250,7 +250,7 @@ int nmod_poly_mat_zls_sorted(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat
     }
 
     //  Now n2 <> 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 +259,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 +279,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 +302,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 +310,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 +330,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 +359,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 +369,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 +383,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 +392,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 +425,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 +452,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 +466,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 +528,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<n; j++) 
+        for (j=0; j<n; j++)
             if (sdeg[j] < 0) sdeg[j]=0;
 
         _nmod_poly_mat_permute_columns_by_sorting_vec(A, n, sdeg, perm);
     }
-    // Input shift, we sort  
+    // Input shift, we sort
     // --------------------
-    else {  
+    else {
 
-        for (j=0; j<n; j++) 
+        for (j=0; j<n; j++)
             sdeg[j]=ishift[j];
 
-        _nmod_poly_mat_permute_columns_by_sorting_vec(A, n, sdeg, perm);        
+        _nmod_poly_mat_permute_columns_by_sorting_vec(A, n, sdeg, perm);
     }
 
     // Call to ZLS
@@ -574,19 +574,19 @@ int nmod_poly_mat_kernel(nmod_poly_mat_t N, slong *degN, const nmod_poly_mat_t i
 
         nmod_poly_mat_clear(NT); // Since nz > 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 +605,11 @@ int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong *degN, const nmod_
 
     slong degmax=deg[0];
     for (i=1; i<n; i++)
-        if (degmax < deg[i]) 
+        if (degmax < deg[i])
             degmax=deg[i];
 
     slong min_nm=n;
-    if (n >m) 
+    if (n >m)
         min_nm=m;
 
     nmod_poly_mat_t AT;
@@ -620,8 +620,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 +629,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 +648,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 +663,7 @@ int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong *degN, const nmod_
     if (n1 > 0) {
 
         k=0;
-        
+
         for (j=0; j<n; j++) {
 
             if (cdeg[j]<0) {
@@ -676,7 +676,7 @@ int nmod_poly_mat_approximant_kernel(nmod_poly_mat_t N, slong *degN, const nmod_
         nmod_poly_mat_clear(PT);
         return n1;
     }
-    
+
     nmod_poly_mat_clear(PT);
-    return 0; 
+    return 0;
 }
diff --git a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_verification.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_verification.c
index e1acb0e0..85a9c91e 100644
--- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_verification.c
+++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_verification.c
@@ -129,8 +129,13 @@ int nmod_poly_mat_is_kernel(const nmod_poly_mat_t ker,
     /* check kernel is shifted reduced */
     if (!nmod_poly_mat_is_ordered_weak_popov(ker, shift, orient))
     {
-        printf("basis is not shifted-weak Popov\n");
-        success = 0;
+        /* 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 */
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..b71e7128 100644
--- a/flint-extras/src/nmod_poly_mat_extra/test/main.c
+++ b/flint-extras/src/nmod_poly_mat_extra/test/main.c
@@ -30,17 +30,17 @@
 
 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-kernel.c b/flint-extras/src/nmod_poly_mat_extra/test/t-kernel.c
index 8f769eb7..3439c5ba 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
@@ -10,79 +10,70 @@
     <https://www.gnu.org/licenses/>.
 */
 
+#include <flint/flint.h>
 #include <flint/nmod_poly_mat.h>
 #include <flint/nmod_types.h>
 #include <flint/profiler.h>
 #include <flint/test_helpers.h>
 
-#include "nmod_poly_mat_utils.h"
 #include "nmod_poly_mat_extra.h"
+#include "nmod_poly_mat_forms.h"
+#include "nmod_poly_mat_io.h"
 
 // test one given input
 /* TODO does not test generation */
 /* TODO does not test reducedness */
 /* TODO does not seem tested with shifts */
-int core_test_kernel(const nmod_poly_mat_t mat)
+int core_test_kernel_zls(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
-    // -------------
+    /* 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);
 
     nmod_poly_mat_t N;
     nmod_poly_mat_init(N, n, n, mat->modulus);
 
-    slong nz;
-
-    int i,j;
-
     slong degN[n];
 
-    nz = nmod_poly_mat_kernel(N, degN, mat, NULL, 2);
+    slong nz = nmod_poly_mat_kernel_zls(N, degN, mat, NULL, 2.);
 
-    int verif;
-    verif = (n-rkflint == nz);
+    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));        
 
-    if (nz !=0)
-    {
-        nmod_poly_mat_t NN;
-        nmod_poly_mat_init(NN, n, nz, mat->modulus);
+    nmod_poly_mat_t Mt;
+    nmod_poly_mat_init(Mt, n, m, mat->modulus);
+    nmod_poly_mat_transpose(Mt, mat);
 
-        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));
+    int verif = nmod_poly_mat_is_kernel(Nt, Mt, rdeg, ROW_LOWER);
 
-        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);
-    }
+    nmod_poly_mat_clear(N);
+    nmod_poly_mat_clear(Nt);
+    nmod_poly_mat_clear(Mt);
+    flint_free(rdeg);
 
     return verif;
 }
 
-TEST_FUNCTION_START(nmod_poly_mat_kernel, state)
+TEST_FUNCTION_START(nmod_poly_mat_kernel_zls, state)
 {
     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);
 
@@ -111,7 +102,7 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state)
             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);
 
@@ -138,7 +129,7 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state)
         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);
 
@@ -165,7 +156,7 @@ TEST_FUNCTION_START(nmod_poly_mat_kernel, state)
         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);
 
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
-
-

From d5490573dfc3d092e7c8aeaec38ba5a2b8388e97 Mon Sep 17 00:00:00 2001
From: Vincent Neiger <vneiger@users.noreply.github.com>
Date: Sat, 3 Jan 2026 09:49:52 +0100
Subject: [PATCH 5/8] fix profile

---
 flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c b/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c
index 658c449a..7587bfa8 100644
--- a/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c
+++ b/flint-extras/src/nmod_poly_mat_extra/profile/p-kernel.c
@@ -68,7 +68,7 @@ void time_##fun(time_args targs, flint_rand_t state)    \
     nmod_poly_mat_clear(Ft); /* TMP */                  \
 }
 
-TIME_KER(kernel)
+TIME_KER(kernel_zls)
 
 /*-------------------------*/
 /*  main                   */
@@ -100,7 +100,7 @@ int main(int argc, char ** argv)
     const slong nfuns = 1;
     typedef void (*timefun) (time_args, flint_rand_t);
     const timefun funs[] = {
-        time_kernel,                      // 0
+        time_kernel_zls,                      // 0
     };
 
     // TODO
@@ -148,7 +148,7 @@ int main(int argc, char ** argv)
     {
         /* rdim; cdim; deg; rank; stype; modn; */
         time_args targs = {8, 4, 1000, 4, 0, n_nextprime(UWORD(1) << 20, 0)};
-        time_kernel(targs, state);
+        time_kernel_zls(targs, state);
         printf(" ");
     }
     printf("\n\n");

From 42b9a1e67cbca5182fd1d68e6174edb4fb67612e Mon Sep 17 00:00:00 2001
From: Vincent Neiger <vneiger@users.noreply.github.com>
Date: Sat, 3 Jan 2026 10:01:26 +0100
Subject: [PATCH 6/8] rename

---
 .../{nmod_poly_mat_kernel.c => kernel.c}                   | 0
 flint-extras/src/nmod_poly_mat_extra/test/main.c           | 2 +-
 .../test/{t-kernel.c => t-kernel_zls.c}                    | 7 +++----
 .../{nmod_poly_mat_verification.c => verification.c}       | 0
 4 files changed, 4 insertions(+), 5 deletions(-)
 rename flint-extras/src/nmod_poly_mat_extra/{nmod_poly_mat_kernel.c => kernel.c} (100%)
 rename flint-extras/src/nmod_poly_mat_extra/test/{t-kernel.c => t-kernel_zls.c} (96%)
 rename flint-extras/src/nmod_poly_mat_extra/{nmod_poly_mat_verification.c => verification.c} (100%)

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 100%
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
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 b71e7128..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,7 +25,6 @@
 #include "t-mul_waksman.c"
 #include "t-rand.c"
 #include "t-weak_popov_form.c"
-#include "t-kernel.c"
 
 /* Array of test functions ***************************************************/
 
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 96%
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 3439c5ba..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
@@ -18,12 +18,11 @@
 
 #include "nmod_poly_mat_extra.h"
 #include "nmod_poly_mat_forms.h"
-#include "nmod_poly_mat_io.h"
 
 // test one given input
 /* TODO does not test generation */
-/* TODO does not test reducedness */
-/* TODO does not seem tested with shifts */
+/* 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;
@@ -44,7 +43,7 @@ int core_test_kernel_zls(const nmod_poly_mat_t mat)
     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_set(nmod_poly_mat_entry(Nt, j, i), nmod_poly_mat_entry(N, i, j));
 
     nmod_poly_mat_t Mt;
     nmod_poly_mat_init(Mt, n, m, mat->modulus);
diff --git a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_verification.c b/flint-extras/src/nmod_poly_mat_extra/verification.c
similarity index 100%
rename from flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_verification.c
rename to flint-extras/src/nmod_poly_mat_extra/verification.c

From 11450fc94199b2d61823ce93955e826e0dad34ba Mon Sep 17 00:00:00 2001
From: Vincent Neiger <vneiger@users.noreply.github.com>
Date: Sat, 3 Jan 2026 10:22:40 +0100
Subject: [PATCH 7/8] maintenance (copyrights, renamings, removing modelines)

---
 .../nmod_poly_mat_extra/approximant_basis.c   |  61 ++++++
 flint-extras/src/nmod_poly_mat_extra/kernel.c |  14 +-
 .../middle_product_geometric.c                |  12 ++
 .../src/nmod_poly_mat_extra/mul_geometric.c   |  12 ++
 .../nmod_poly_mat_approximant_mbasis.c        |  25 ---
 .../nmod_poly_mat_approximant_mbasis1.c       | 199 ------------------
 .../nmod_poly_mat_approximant_pmbasis.c       |  35 ---
 .../nmod_poly_mat_degree_matrix.c             |  16 +-
 .../nmod_poly_mat_description.bak             |  21 +-
 .../nmod_poly_mat_extra/nmod_poly_mat_det.c   |  15 +-
 .../nmod_poly_mat_extra/nmod_poly_mat_dixon.c |  14 +-
 .../nmod_poly_mat_hermite.c                   |  15 +-
 .../nmod_poly_mat_extra/nmod_poly_mat_io.c    |  15 +-
 .../nmod_poly_mat_is_form.c                   |  15 +-
 .../nmod_poly_mat_is_unimodular.c             |  15 +-
 .../nmod_poly_mat_leading_matrix.c            |  15 +-
 .../nmod_poly_mat_middle_product_naive.c      |  12 ++
 .../nmod_poly_mat_mul_waksman.c               |  14 +-
 .../nmod_poly_mat_pivots.c                    |  17 +-
 .../nmod_poly_mat_extra/nmod_poly_mat_rand.c  |  17 +-
 .../nmod_poly_mat_rdeg_cdeg.c                 |  15 +-
 .../nmod_poly_mat_extra/nmod_poly_mat_utils.c |  17 +-
 .../nmod_poly_mat_weak_popov.c                |  15 +-
 .../src/nmod_poly_mat_extra/test/t-dixon.c    |   3 -
 .../nmod_poly_mat_extra/timings/time_hnf.c    |   3 -
 .../timings/time_pmbasis.c                    |   3 -
 .../timings/time_polmatmul.c                  |   3 -
 .../src/nmod_poly_mat_extra/verification.c    |  13 +-
 28 files changed, 301 insertions(+), 330 deletions(-)
 create mode 100644 flint-extras/src/nmod_poly_mat_extra/approximant_basis.c
 delete mode 100644 flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_approximant_mbasis.c
 delete mode 100644 flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_approximant_mbasis1.c
 delete mode 100644 flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_approximant_pmbasis.c

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
+    <https://www.gnu.org/licenses/>.
+*/
+
+#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/kernel.c b/flint-extras/src/nmod_poly_mat_extra/kernel.c
index f54ff663..f46208d3 100644
--- a/flint-extras/src/nmod_poly_mat_extra/kernel.c
+++ b/flint-extras/src/nmod_poly_mat_extra/kernel.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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <math.h>
 #include <stdlib.h>
 
@@ -29,7 +41,6 @@
  *
  */
 
-
 // Dummy function for qsort
  int cmp(const void *a, const void *b)
 {
@@ -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)
 {
diff --git a/flint-extras/src/nmod_poly_mat_extra/middle_product_geometric.c b/flint-extras/src/nmod_poly_mat_extra/middle_product_geometric.c
index 6d1abe27..2eb644dc 100644
--- a/flint-extras/src/nmod_poly_mat_extra/middle_product_geometric.c
+++ b/flint-extras/src/nmod_poly_mat_extra/middle_product_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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <flint/nmod_mat.h>
 #include <flint/nmod_poly.h>
 #include <flint/nmod_poly_mat.h>
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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <flint/nmod.h>
 #include <flint/nmod_vec.h>
 #include <flint/nmod_mat.h>
diff --git a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_approximant_mbasis.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_approximant_mbasis.c
deleted file mode 100644
index 188d146a..00000000
--- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_approximant_mbasis.c
+++ /dev/null
@@ -1,25 +0,0 @@
-#include "nmod_mat_extra.h"
-#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);
-    // 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 <stdlib.h>
-#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_degree_matrix.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_degree_matrix.c
index 9fb88fad..e5e2a098 100644
--- 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
@@ -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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <flint/fmpz_mat.h>
 #include <flint/nmod_poly.h>
 #include "nmod_poly_mat_forms.h"
@@ -31,7 +43,3 @@ void nmod_poly_mat_degree_matrix_shifted(fmpz_mat_t dmat,
             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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <stdlib.h>
 #include <math.h> 
 #include <flint/nmod.h> 
@@ -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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <flint/nmod_vec.h>
 #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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <stdlib.h>
 #include <math.h> 
 
@@ -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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <flint/ulong_extras.h>
 #include <flint/nmod_vec.h>
 #include <flint/nmod.h>
@@ -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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <flint/fmpz_mat.h>
 #include <flint/nmod_mat.h>
 #include <flint/nmod_poly.h>
@@ -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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <stdlib.h> // qsort
 #include <flint/nmod_mat.h>
 #include <flint/nmod_poly.h>
@@ -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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include "nmod_poly_mat_utils.h"
 #include <flint/flint.h>
 #include <flint/nmod_poly_mat.h>
@@ -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
index 251ddfe3..5f641cb1 100644
--- 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
@@ -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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <flint/nmod_mat.h>
 #include <flint/nmod_poly.h>
 #include "nmod_poly_mat_forms.h"
@@ -70,6 +82,3 @@ void nmod_poly_mat_leading_matrix(nmod_mat_t 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_middle_product_naive.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_middle_product_naive.c
index c24d01f8..8968d11c 100644
--- 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
@@ -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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <flint/nmod_poly_mat.h>
 
 /** Middle product for polynomial matrices
diff --git a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_mul_waksman.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_mul_waksman.c
index da5ebdc9..1062c48c 100644
--- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_mul_waksman.c
+++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <flint/nmod_poly.h>
 #include <flint/nmod_poly_mat.h>
 
@@ -142,5 +154,3 @@ void nmod_poly_mat_mul_waksman(nmod_poly_mat_t C, const nmod_poly_mat_t A,  cons
     nmod_poly_clear(val2);
     nmod_poly_clear(crow);
 }
-
-
diff --git a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_pivots.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_pivots.c
index 7fd8a858..6d850604 100644
--- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_pivots.c
+++ b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_pivots.c
@@ -1,7 +1,17 @@
-#include <flint/nmod_poly.h>
+/*
+    Copyright (C) 2025 Vincent Neiger
 
-#include "nmod_poly_mat_forms.h"
+    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
+    <https://www.gnu.org/licenses/>.
+*/
 
+#include <flint/nmod_poly.h>
+#include "nmod_poly_mat_forms.h"
 
 /**********************************************************************
 *                    shifted pivot profile, vector                    *
@@ -270,6 +280,3 @@ void nmod_poly_mat_echelon_pivot_profile(slong * pivind, slong * pivdeg, const n
         break;
     }
 }
-
-/* -*- 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_rand.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_rand.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/nmod_poly_mat_rand.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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <flint/fmpz_mat.h>
 #include <flint/fmpz_vec.h>
 #include <flint/nmod_poly.h>
@@ -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/nmod_poly_mat_rdeg_cdeg.c b/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_rdeg_cdeg.c
index 8d3580df..3ac43899 100644
--- 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
@@ -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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <flint/nmod_poly.h>
 
 #include "nmod_poly_mat_forms.h"
@@ -85,6 +97,3 @@ void nmod_poly_mat_column_degree(slong *cdeg, const nmod_poly_mat_t mat, const s
         }
     }
 }
-
-/* -*- 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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <stdlib.h> // for qsort
 #include <flint/nmod_poly.h>
 #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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include <flint/nmod.h>
 #include <flint/nmod_vec.h>
 #include <flint/nmod_poly.h>
@@ -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/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/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
index 85a9c91e..1ce2fe37 100644
--- a/flint-extras/src/nmod_poly_mat_extra/verification.c
+++ b/flint-extras/src/nmod_poly_mat_extra/verification.c
@@ -1,4 +1,15 @@
-#include "nmod_poly_mat_utils.h"
+/*
+    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
+    <https://www.gnu.org/licenses/>.
+*/
+
 #include "nmod_poly_mat_forms.h"
 #include <flint/nmod_mat.h>
 #include <flint/nmod_poly.h>

From 9b5fbf18f2bd0cd89ec84028e988c39bd4cb3350 Mon Sep 17 00:00:00 2001
From: Vincent Neiger <vneiger@users.noreply.github.com>
Date: Sat, 3 Jan 2026 11:31:53 +0100
Subject: [PATCH 8/8] some more maintenance

---
 ...ivots.c => degrees_pivots_leadingmatrix.c} | 218 ++++++++++++++++--
 ...e_product_geometric.c => middle_product.c} |  14 ++
 ...d_poly_mat_mul_waksman.c => mul_waksman.c} |  26 +--
 .../nmod_poly_mat_degree_matrix.c             |  45 ----
 .../nmod_poly_mat_leading_matrix.c            |  84 -------
 .../nmod_poly_mat_middle_product_naive.c      |  26 ---
 .../nmod_poly_mat_rdeg_cdeg.c                 |  99 --------
 .../{nmod_poly_mat_rand.c => randomisation.c} |   0
 8 files changed, 225 insertions(+), 287 deletions(-)
 rename flint-extras/src/nmod_poly_mat_extra/{nmod_poly_mat_pivots.c => degrees_pivots_leadingmatrix.c} (58%)
 rename flint-extras/src/nmod_poly_mat_extra/{middle_product_geometric.c => middle_product.c} (92%)
 rename flint-extras/src/nmod_poly_mat_extra/{nmod_poly_mat_mul_waksman.c => mul_waksman.c} (97%)
 delete mode 100644 flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_degree_matrix.c
 delete mode 100644 flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_leading_matrix.c
 delete mode 100644 flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_middle_product_naive.c
 delete mode 100644 flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_rdeg_cdeg.c
 rename flint-extras/src/nmod_poly_mat_extra/{nmod_poly_mat_rand.c => randomisation.c} (100%)

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 58%
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 6d850604..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
@@ -10,12 +10,133 @@
     <https://www.gnu.org/licenses/>.
 */
 
+#include <flint/fmpz_mat.h>
+#include <flint/nmod_mat.h>
 #include <flint/nmod_poly.h>
+
 #include "nmod_poly_mat_forms.h"
 
-/**********************************************************************
-*                    shifted pivot profile, vector                    *
-**********************************************************************/
+/*------------------------------------------------------------*/
+/*        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;
+
+        // 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,
@@ -108,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,
@@ -138,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,
@@ -165,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)
@@ -221,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)
@@ -280,3 +384,77 @@ void nmod_poly_mat_echelon_pivot_profile(slong * pivind, slong * pivdeg, const n
         break;
     }
 }
+
+
+/*------------------------------------------------------------*/
+/*                      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/middle_product_geometric.c b/flint-extras/src/nmod_poly_mat_extra/middle_product.c
similarity index 92%
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 2eb644dc..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
@@ -17,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/nmod_poly_mat_mul_waksman.c b/flint-extras/src/nmod_poly_mat_extra/mul_waksman.c
similarity index 97%
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 1062c48c..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
@@ -24,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);
@@ -41,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);
@@ -66,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);
@@ -86,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++)
@@ -104,7 +104,7 @@ void nmod_poly_mat_mul_waksman(nmod_poly_mat_t C, const nmod_poly_mat_t A,  cons
             }
         }
     }
-    
+
     for (l=1; l<m; l++)
     {
         nmod_poly_add(val1, Ccol + l, nmod_poly_mat_entry(C, l, 0));
@@ -115,21 +115,21 @@ void nmod_poly_mat_mul_waksman(nmod_poly_mat_t C, const nmod_poly_mat_t A,  cons
     nmod_poly_add(val1, Crow, nmod_poly_mat_entry(C, 0, 0));
     nmod_poly_scalar_mul_nmod(val0, val1, half);
     nmod_poly_sub(nmod_poly_mat_entry(C, 0, 0), nmod_poly_mat_entry(C, 0, 0), val0);
-    
+
     for (k = 1; k < p; k++)
     {
         nmod_poly_add(crow, Crow + k, nmod_poly_mat_entry(C, 0, k));
         nmod_poly_scalar_mul_nmod(val1, crow, half);
         nmod_poly_sub(nmod_poly_mat_entry(C, 0, k), nmod_poly_mat_entry(C, 0, k), val1);
         nmod_poly_sub(crow, val1, val0);
-        
+
         for (l = 1; l < m; l++)
         {
             nmod_poly_sub(val2, nmod_poly_mat_entry(C, l, k), crow);
             nmod_poly_sub(nmod_poly_mat_entry(C, l, k), val2, Ccol + l);
         }
     }
-    
+
     if (n & 1)
     {
         for (l = 0; l < m; l++)
@@ -142,13 +142,13 @@ void nmod_poly_mat_mul_waksman(nmod_poly_mat_t C, const nmod_poly_mat_t A,  cons
 
     for (i = 0; i < p; i++)
         nmod_poly_clear(Crow+i);
-    
+
     for (i = 0; i < m; i++)
         nmod_poly_clear(Ccol+i);
-    
+
     flint_free(Crow);
     flint_free(Ccol);
-    
+
     nmod_poly_clear(val0);
     nmod_poly_clear(val1);
     nmod_poly_clear(val2);
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 e5e2a098..00000000
--- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_degree_matrix.c
+++ /dev/null
@@ -1,45 +0,0 @@
-/*
-    Copyright (C) 2025 Vincent Neiger
-
-    This file is part of PML.
-
-    PML is free software: you can redistribute it and/or modify it under
-    the terms of the GNU General Public License version 2.0 (GPL-2.0-or-later)
-    as published by the Free Software Foundation; either version 2 of the
-    License, or (at your option) any later version. See
-    <https://www.gnu.org/licenses/>.
-*/
-
-#include <flint/fmpz_mat.h>
-#include <flint/nmod_poly.h>
-#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));
-}
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 5f641cb1..00000000
--- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_leading_matrix.c
+++ /dev/null
@@ -1,84 +0,0 @@
-/*
-    Copyright (C) 2025 Vincent Neiger
-
-    This file is part of PML.
-
-    PML is free software: you can redistribute it and/or modify it under
-    the terms of the GNU General Public License version 2.0 (GPL-2.0-or-later)
-    as published by the Free Software Foundation; either version 2 of the
-    License, or (at your option) any later version. See
-    <https://www.gnu.org/licenses/>.
-*/
-
-#include <flint/nmod_mat.h>
-#include <flint/nmod_poly.h>
-#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);
-        }
-    }
-}
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 8968d11c..00000000
--- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_middle_product_naive.c
+++ /dev/null
@@ -1,26 +0,0 @@
-/*
-    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
-    <https://www.gnu.org/licenses/>.
-*/
-
-#include <flint/nmod_poly_mat.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);
-}
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 3ac43899..00000000
--- a/flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_rdeg_cdeg.c
+++ /dev/null
@@ -1,99 +0,0 @@
-/*
-    Copyright (C) 2025 Vincent Neiger
-
-    This file is part of PML.
-
-    PML is free software: you can redistribute it and/or modify it under
-    the terms of the GNU General Public License version 2.0 (GPL-2.0-or-later)
-    as published by the Free Software Foundation; either version 2 of the
-    License, or (at your option) any later version. See
-    <https://www.gnu.org/licenses/>.
-*/
-
-#include <flint/nmod_poly.h>
-
-#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;
-        }
-    }
-}
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 100%
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