diff --git a/flint-extras/src/nmod_poly_mat_approximant.h b/flint-extras/src/nmod_poly_mat_approximant.h
index 400367d3..f101af23 100644
--- a/flint-extras/src/nmod_poly_mat_approximant.h
+++ b/flint-extras/src/nmod_poly_mat_approximant.h
@@ -316,11 +316,28 @@ 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,
                            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 9c216269..579fe6e2 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)
@@ -48,6 +48,42 @@ void nmod_poly_mat_pmbasis(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_old(appbas, shift, pmat, order1);
+
+    nmod_poly_mat_middle_product_naive_old(residual, appbas, pmat, order1, order2-1);
+
+    nmod_poly_mat_pmbasis_old(appbas2, shift, residual, order2);
+
+    nmod_poly_mat_mul(appbas, appbas2, 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);
@@ -59,3 +95,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_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);
+}
+
+
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..5b9b45ae 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,30 @@ 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); // No aliasing specified for that?
+
+    nmod_poly_mat_mul(C,A,BT);  
+    nmod_poly_mat_shift_right(C, C, degA);  
+    nmod_poly_mat_truncate(C, d2+1);
+
+    nmod_poly_mat_clear(BT);
+}
+
+
 /** 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