Kernel basis: add profile

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 -*- */

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);
+}

flint-extras/src/nmod_poly_mat_extra/nmod_poly_mat_pivots.c → flint-extras/src/nmod_poly_mat_extra/degrees_pivots_leadingmatrix.c

@@ -1,11 +1,142 @@
+/*
+    Copyright (C) 2025 Vincent Neiger
+
+    This file is part of PML.
+
+    PML is free software: you can redistribute it and/or modify it under
+    the terms of the GNU General Public License version 2.0 (GPL-2.0-or-later)
+    as published by the Free Software Foundation; either version 2 of the
+    License, or (at your option) any later version. See
+    <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"
 
+/*------------------------------------------------------------*/
+/*        row degree - column degree - degree matrix          */
+/*------------------------------------------------------------*/
+
+void nmod_poly_mat_row_degree(slong *rdeg, const nmod_poly_mat_t mat, const slong *shift)
+{
+    if (shift == NULL)
+    {
+        slong max, d;
+        for (slong i = 0; i < mat->r; i++)
+        {
+            max = -1;
+            for (slong j = 0; j < mat->c; j++)
+            {
+                d = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j));
+                if (max < d)
+                    max = d;
+            }
+            rdeg[i] = max;
+        }
+    }
+    else
+    {
+        slong max, d;
+        slong min_shift = (mat->c > 0) ? shift[0] : 0;
+
+        // find minimum of shift
+        for (slong j = 0; j < mat->c; ++j)
+            if (shift[j] < min_shift)
+                min_shift = shift[j];
+
+        for (slong i = 0; i < mat->r; i++)
+        {
+            max = min_shift-1; // zero rows will have this as rdeg
+            for (slong j = 0; j < mat->c; j++)
+            {
+                d = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)) + shift[j];
+                // if new maximum at a nonzero entry, update
+                if (shift[j] <= d && max < d)
+                    max = d;
+            }
+            rdeg[i] = max;
+        }
+    }
+}
+
+void nmod_poly_mat_column_degree(slong *cdeg, const nmod_poly_mat_t mat, const slong *shift)
+{
+    if (shift == NULL)
+    {
+        slong max, d;
+        for (slong j = 0; j < mat->c; j++)
+        {
+            max = -1;
+            for (slong i = 0; i < mat->r; i++)
+            {
+                d = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j));
+                if (max < d)
+                    max = d;
+            }
+            cdeg[j] = max;
+        }
+    }
+    else
+    {
+        slong max, d;
+        slong min_shift = (mat->r > 0) ? shift[0] : 0;
 
-/**********************************************************************
-*                    shifted pivot profile, vector                    *
-**********************************************************************/
+        // find minimum of shift
+        for (slong i = 0; i < mat->r; ++i)
+            if (shift[i] < min_shift)
+                min_shift = shift[i];
+
+        for (slong j = 0; j < mat->c; j++)
+        {
+            max = min_shift-1;
+            for (slong i = 0; i < mat->r; i++)
+            {
+                d = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)) + shift[i];
+                // if new maximum at a nonzero entry, update
+                if (shift[i] <= d && max < d)
+                    max = d;
+            }
+            cdeg[j] = max;
+        }
+    }
+}
+
+void nmod_poly_mat_degree_matrix(fmpz_mat_t dmat, const nmod_poly_mat_t mat)
+{
+    for(slong i = 0; i < mat->r; i++)
+        for(slong j = 0; j < mat->c; j++)
+            *fmpz_mat_entry(dmat, i, j) = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j));
+}
+
+void nmod_poly_mat_degree_matrix_shifted(fmpz_mat_t dmat,
+                                         const nmod_poly_mat_t mat,
+                                         const slong * shift,
+                                         orientation_t orient)
+{
+    if (shift)
+    {
+        if (orient == ROW_LOWER || orient == ROW_UPPER)
+            for(slong i = 0; i < mat->r; i++)
+                for(slong j = 0; j < mat->c; j++)
+                    *fmpz_mat_entry(dmat, i, j) = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)) + shift[j];
+        else // orient == COL_*
+            for(slong i = 0; i < mat->r; i++)
+                for(slong j = 0; j < mat->c; j++)
+                    *fmpz_mat_entry(dmat, i, j) = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j)) + shift[i];
+    }
+    else
+        for(slong i = 0; i < mat->r; i++)
+            for(slong j = 0; j < mat->c; j++)
+                *fmpz_mat_entry(dmat, i, j) = nmod_poly_degree(nmod_poly_mat_entry(mat, i, j));
+}
+
+
+/*------------------------------------------------------------*/
+/*          shifted pivot profile (index + degree)            */
+/*------------------------------------------------------------*/
 
 void _nmod_poly_vec_pivot_profile(slong * pivind,
                                   slong * pivdeg,
@@ -98,10 +229,6 @@ void _nmod_poly_vec_pivot_profile(slong * pivind,
     }
 }
 
-/**********************************************************************
-*                        shifted pivot index                         *
-**********************************************************************/
-
 void nmod_poly_mat_pivot_index(slong *pivind,
                                const nmod_poly_mat_t mat,
                                const slong * shift,
@@ -128,10 +255,6 @@ void nmod_poly_mat_pivot_index(slong *pivind,
     }
 }
 
-/**********************************************************************
-*                       shifted pivot profile                        *
-**********************************************************************/
-
 void nmod_poly_mat_pivot_profile(slong * pivind,
                                  slong * pivdeg,
                                  const nmod_poly_mat_t mat,
@@ -155,11 +278,6 @@ void nmod_poly_mat_pivot_profile(slong * pivind,
     }
 }
 
-
-/**********************************************************************
-*                        echelon pivot index                         *
-**********************************************************************/
-
 void nmod_poly_mat_echelon_pivot_index(slong * pivind, const nmod_poly_mat_t mat, orientation_t orient)
 {
     switch (orient)
@@ -211,10 +329,6 @@ void nmod_poly_mat_echelon_pivot_index(slong * pivind, const nmod_poly_mat_t mat
     }
 }
 
-/**********************************************************************
-*                       echelon pivot profile                        *
-**********************************************************************/
-
 void nmod_poly_mat_echelon_pivot_profile(slong * pivind, slong * pivdeg, const nmod_poly_mat_t mat, orientation_t orient)
 {
     switch (orient)
@@ -271,5 +385,76 @@ void nmod_poly_mat_echelon_pivot_profile(slong * pivind, slong * pivdeg, const n
     }
 }
 
-/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
-// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
+
+/*------------------------------------------------------------*/
+/*                      leading matrix                        */
+/*------------------------------------------------------------*/
+
+void nmod_poly_mat_leading_matrix(nmod_mat_t lmat,
+                                  const nmod_poly_mat_t mat,
+                                  const slong * shift,
+                                  orientation_t orient)
+{
+    if (orient == ROW_LOWER || orient == ROW_UPPER)
+    {
+        // find row degrees
+        slong rdeg[mat->r];
+        nmod_poly_mat_row_degree(rdeg, mat, shift);
+
+        // deduce leading matrix
+        if (shift == NULL)
+        {
+            for (slong i = 0; i < mat->r; i++)
+                if (rdeg[i] >= 0)
+                    for (slong j = 0; j < mat->c; j++)
+                        nmod_mat_set_entry(lmat, i, j,
+                                           nmod_poly_get_coeff_ui(nmod_poly_mat_entry(mat, i, j),
+                                                                  rdeg[i]));
+                else
+                    for (slong j = 0; j < mat->c; j++)
+                        nmod_mat_set_entry(lmat, i, j, 0);
+        }
+        else
+        {
+            for (slong i = 0; i < mat->r; i++)
+                for (slong j = 0; j < mat->c; j++)
+                    if (rdeg[i] >= shift[j])
+                        nmod_mat_set_entry(lmat, i, j,
+                                           nmod_poly_get_coeff_ui(nmod_poly_mat_entry(mat, i, j),
+                                                                  rdeg[i] - shift[j]));
+                    else
+                        nmod_mat_set_entry(lmat, i, j, 0);
+        }
+    }
+    else // orient == COL_*
+    {
+        // find column degrees
+        slong cdeg[mat->c];
+        nmod_poly_mat_column_degree(cdeg, mat, shift);
+
+        // deduce leading matrix
+        if (shift == NULL)
+        {
+            for (slong j = 0; j < mat->c; j++)
+                if (cdeg[j] >= 0)
+                    for (slong i = 0; i < mat->r; i++)
+                        nmod_mat_set_entry(lmat, i, j,
+                                           nmod_poly_get_coeff_ui(nmod_poly_mat_entry(mat, i, j),
+                                                                  cdeg[j]));
+                else
+                    for (slong i = 0; i < mat->r; i++)
+                        nmod_mat_set_entry(lmat, i, j, 0);
+        }
+        else
+        {
+            for (slong j = 0; j < mat->c; j++)
+                for (slong i = 0; i < mat->r; i++)
+                    if (cdeg[j] >= shift[i])
+                        nmod_mat_set_entry(lmat, i, j,
+                                           nmod_poly_get_coeff_ui(nmod_poly_mat_entry(mat, i, j),
+                                                                  cdeg[j] - shift[i]));
+                    else
+                        nmod_mat_set_entry(lmat, i, j, 0);
+        }
+    }
+}