From 2fe41af761d3cda34b1279654d7cfbc7b25e62af Mon Sep 17 00:00:00 2001 From: Gilles Villard Date: Wed, 29 Apr 2026 10:26:46 +0200 Subject: [PATCH 1/5] initial truncation in pmbasis --- flint-extras/src/nmod_poly_mat_approximant.h | 7 ++++ .../nmod_poly_mat_extra/approximant_basis.c | 40 ++++++++++++++++++- 2 files changed, 46 insertions(+), 1 deletion(-) diff --git a/flint-extras/src/nmod_poly_mat_approximant.h b/flint-extras/src/nmod_poly_mat_approximant.h index 400367d3..9cf241cb 100644 --- a/flint-extras/src/nmod_poly_mat_approximant.h +++ b/flint-extras/src/nmod_poly_mat_approximant.h @@ -316,6 +316,13 @@ void nmod_poly_mat_mbasis(nmod_poly_mat_t appbas, // TODO doc // TODO middle_product currently naive + + +void nmod_poly_mat_pmbasis_old(nmod_poly_mat_t appbas, + slong * shift, + const nmod_poly_mat_t pmat, + slong order); + void nmod_poly_mat_pmbasis(nmod_poly_mat_t appbas, slong * shift, const nmod_poly_mat_t pmat, diff --git a/flint-extras/src/nmod_poly_mat_extra/approximant_basis.c b/flint-extras/src/nmod_poly_mat_extra/approximant_basis.c index 9c216269..428f1b60 100644 --- a/flint-extras/src/nmod_poly_mat_extra/approximant_basis.c +++ b/flint-extras/src/nmod_poly_mat_extra/approximant_basis.c @@ -30,7 +30,7 @@ void nmod_poly_mat_mbasis(nmod_poly_mat_t appbas, nmod_mat_poly_clear(app); } -void nmod_poly_mat_pmbasis(nmod_poly_mat_t appbas, +void nmod_poly_mat_pmbasis_old(nmod_poly_mat_t appbas, slong * shift, const nmod_poly_mat_t pmat, slong order) @@ -59,3 +59,41 @@ void nmod_poly_mat_pmbasis(nmod_poly_mat_t appbas, nmod_poly_mat_clear(appbas2); nmod_poly_mat_clear(residual); } + +void nmod_poly_mat_pmbasis(nmod_poly_mat_t appbas, + slong * shift, + const nmod_poly_mat_t ipmat, + slong order) +{ + + nmod_poly_mat_t pmat; + + nmod_poly_mat_init(pmat, ipmat->r, ipmat->c, ipmat->modulus); + nmod_poly_mat_set_trunc(pmat,ipmat,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); +} + + From eeea60624b82b0391cf9a07554cf88c51f9ad5fe Mon Sep 17 00:00:00 2001 From: Gilles Villard Date: Wed, 29 Apr 2026 10:43:47 +0200 Subject: [PATCH 2/5] new middle_product_naive --- .../src/nmod_poly_mat_extra/middle_product.c | 27 ++++++++++++++++++- flint-extras/src/nmod_poly_mat_multiply.h | 10 +++++++ 2 files changed, 36 insertions(+), 1 deletion(-) diff --git a/flint-extras/src/nmod_poly_mat_extra/middle_product.c b/flint-extras/src/nmod_poly_mat_extra/middle_product.c index ec7a460a..a38805c2 100644 --- a/flint-extras/src/nmod_poly_mat_extra/middle_product.c +++ b/flint-extras/src/nmod_poly_mat_extra/middle_product.c @@ -22,7 +22,7 @@ * 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, +void nmod_poly_mat_middle_product_naive_old(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); @@ -31,6 +31,31 @@ void nmod_poly_mat_middle_product_naive(nmod_poly_mat_t C, const nmod_poly_mat_t } +/** Middle product for polynomial matrices + * sets C = ((A * B) div x^d1) mod x^(d2+1), assuming deg(A) <= d1 and deg(B) <= d1 + d2 + * 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 d1, const ulong d2) +{ + slong degA; + + degA=nmod_poly_mat_degree(A); + + nmod_poly_mat_t BT; + nmod_poly_mat_init(BT, A->c, B->c, B->modulus); + nmod_poly_mat_shift_right(BT, B, d1-degA); + + nmod_poly_mat_mul(C,A,BT); + + nmod_poly_mat_shift_right(C, C, degA); + + nmod_poly_mat_truncate(C, d2+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_multiply.h b/flint-extras/src/nmod_poly_mat_multiply.h index b2839d39..4d99a25a 100644 --- a/flint-extras/src/nmod_poly_mat_multiply.h +++ b/flint-extras/src/nmod_poly_mat_multiply.h @@ -78,9 +78,19 @@ void nmod_poly_mat_mul_waksman(nmod_poly_mat_t C, const nmod_poly_mat_t A, cons * output can alias input * naive implementation (multiply, shift, truncate) */ +void nmod_poly_mat_middle_product_naive_old(nmod_poly_mat_t C, const nmod_poly_mat_t A, const nmod_poly_mat_t B, + const ulong dA, const ulong dB); + + +/** Middle product for polynomial matrices + * sets C = ((A * B) div x^d1) mod x^(d2+1), assuming deg(A) <= d1 and deg(B) <= d1 + d2 + * 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); + /** Middle product for polynomial matrices * sets C = ((A * B) div x^dA) mod x^(dB+1) * output can alias input From a8862b5c40ea952990ad5fc79e4ac209b520a25f Mon Sep 17 00:00:00 2001 From: Gilles Villard Date: Wed, 29 Apr 2026 12:47:57 +0200 Subject: [PATCH 3/5] new calls in old --- flint-extras/src/nmod_poly_mat_extra/approximant_basis.c | 6 +++--- flint-extras/src/nmod_poly_mat_extra/middle_product.c | 1 - 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/flint-extras/src/nmod_poly_mat_extra/approximant_basis.c b/flint-extras/src/nmod_poly_mat_extra/approximant_basis.c index 428f1b60..a93e642e 100644 --- a/flint-extras/src/nmod_poly_mat_extra/approximant_basis.c +++ b/flint-extras/src/nmod_poly_mat_extra/approximant_basis.c @@ -48,11 +48,11 @@ void nmod_poly_mat_pmbasis_old(nmod_poly_mat_t appbas, 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_pmbasis_old(appbas, shift, pmat, order1); - nmod_poly_mat_middle_product_naive(residual, appbas, pmat, order1, order2-1); + nmod_poly_mat_middle_product_naive_old(residual, appbas, pmat, order1, order2-1); - nmod_poly_mat_pmbasis(appbas2, shift, residual, order2); + nmod_poly_mat_pmbasis_old(appbas2, shift, residual, order2); nmod_poly_mat_mul(appbas, appbas2, appbas); diff --git a/flint-extras/src/nmod_poly_mat_extra/middle_product.c b/flint-extras/src/nmod_poly_mat_extra/middle_product.c index a38805c2..3843122e 100644 --- a/flint-extras/src/nmod_poly_mat_extra/middle_product.c +++ b/flint-extras/src/nmod_poly_mat_extra/middle_product.c @@ -41,7 +41,6 @@ void nmod_poly_mat_middle_product_naive(nmod_poly_mat_t C, const nmod_poly_mat_t const nmod_poly_mat_t B, const ulong d1, const ulong d2) { slong degA; - degA=nmod_poly_mat_degree(A); nmod_poly_mat_t BT; From b40e9b21a10ca8a9d28302938f25609394150e83 Mon Sep 17 00:00:00 2001 From: Gilles Villard Date: Wed, 29 Apr 2026 17:53:15 +0200 Subject: [PATCH 4/5] clear, no aliasing specified for nmod_poly_mat_shift_right? --- flint-extras/src/nmod_poly_mat_extra/middle_product.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/flint-extras/src/nmod_poly_mat_extra/middle_product.c b/flint-extras/src/nmod_poly_mat_extra/middle_product.c index 3843122e..5b9b45ae 100644 --- a/flint-extras/src/nmod_poly_mat_extra/middle_product.c +++ b/flint-extras/src/nmod_poly_mat_extra/middle_product.c @@ -45,13 +45,13 @@ void nmod_poly_mat_middle_product_naive(nmod_poly_mat_t C, const nmod_poly_mat_t nmod_poly_mat_t BT; nmod_poly_mat_init(BT, A->c, B->c, B->modulus); - nmod_poly_mat_shift_right(BT, B, d1-degA); + nmod_poly_mat_shift_right(BT, B, d1-degA); // No aliasing specified for that? nmod_poly_mat_mul(C,A,BT); - - nmod_poly_mat_shift_right(C, C, degA); - + nmod_poly_mat_shift_right(C, C, degA); nmod_poly_mat_truncate(C, d2+1); + + nmod_poly_mat_clear(BT); } From d33f91ac9a4ff89062043bf212745fd0849d2db2 Mon Sep 17 00:00:00 2001 From: Gilles Villard Date: Wed, 29 Apr 2026 18:11:35 +0200 Subject: [PATCH 5/5] starting for the use of the geometric product --- flint-extras/src/nmod_poly_mat_approximant.h | 10 ++++++ .../nmod_poly_mat_extra/approximant_basis.c | 36 +++++++++++++++++++ 2 files changed, 46 insertions(+) diff --git a/flint-extras/src/nmod_poly_mat_approximant.h b/flint-extras/src/nmod_poly_mat_approximant.h index 9cf241cb..f101af23 100644 --- a/flint-extras/src/nmod_poly_mat_approximant.h +++ b/flint-extras/src/nmod_poly_mat_approximant.h @@ -328,6 +328,16 @@ void nmod_poly_mat_pmbasis(nmod_poly_mat_t appbas, const nmod_poly_mat_t pmat, slong order); +/** Using geometric matrix product for the partial bases + * and linearises the product residual x basis + */ + +void nmod_poly_mat_pmbasis_geometric(nmod_poly_mat_t appbas, + slong * shift, + const nmod_poly_mat_t pmat, + slong order); + + /** Computes a `shift`-Popov approximant basis for `(pmat,order)` using the * algorithm PM-Basis (see @ref pmbasis) twice: the first call yields an diff --git a/flint-extras/src/nmod_poly_mat_extra/approximant_basis.c b/flint-extras/src/nmod_poly_mat_extra/approximant_basis.c index a93e642e..579fe6e2 100644 --- a/flint-extras/src/nmod_poly_mat_extra/approximant_basis.c +++ b/flint-extras/src/nmod_poly_mat_extra/approximant_basis.c @@ -96,4 +96,40 @@ void nmod_poly_mat_pmbasis(nmod_poly_mat_t appbas, nmod_poly_mat_clear(residual); } +void nmod_poly_mat_pmbasis_geometric(nmod_poly_mat_t appbas, + slong * shift, + const nmod_poly_mat_t ipmat, + slong order) +{ + + nmod_poly_mat_t pmat; + + nmod_poly_mat_init(pmat, ipmat->r, ipmat->c, ipmat->modulus); + nmod_poly_mat_set_trunc(pmat,ipmat,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_geometric(appbas, shift, pmat, order1); + + nmod_poly_mat_middle_product_naive(residual, appbas, pmat, order1, order2-1); + + nmod_poly_mat_pmbasis_geometric(appbas2, shift, residual, order2); + + nmod_poly_mat_mul_geometric(appbas, appbas2, appbas); + + nmod_poly_mat_clear(appbas2); + nmod_poly_mat_clear(residual); +} +