Merge pull request #1759 from flintlib/divexact

nmod_poly_divexact

Author: fredrik-johansson

doc/source/nmod_poly.rst

@@ -1032,6 +1032,12 @@ Division
 
     Computes the remainder `R` on polynomial division of `A` by `B`.
 
+.. function:: void _nmod_poly_divexact(mp_ptr Q, mp_srcptr A, slong lenA, mp_srcptr B, slong lenB, nmod_t mod)
+              void nmod_poly_divexact(nmod_poly_t Q, const nmod_poly_t A, const nmod_poly_t B)
+
+    Computes the quotient `Q` of `A` and `B` assuming that the division
+    is exact.
+
 .. function:: void _nmod_poly_inv_series_basecase(mp_ptr Qinv, mp_srcptr Q, slong Qlen, slong n, nmod_t mod)
 
     Given ``Q`` of length ``Qlen`` whose leading coefficient is invertible

src/gr/nmod.c

@@ -1002,6 +1002,19 @@ _gr_nmod_poly_divrem(mp_ptr Q, mp_ptr R, mp_srcptr A, slong lenA,
     }
 }
 
+int
+_gr_nmod_poly_divexact(mp_ptr Q, mp_srcptr A, slong lenA, mp_srcptr B, slong lenB, gr_ctx_t ctx)
+{
+    slong lenQ = lenA - lenB + 1;
+
+    if (lenB <= 40 || lenQ <= 20)
+        return _gr_poly_divexact_basecase(Q, A, lenA, B, lenB, ctx);
+    else if (lenB <= 60 || lenQ <= 60)
+        return _gr_poly_divexact_basecase_bidirectional(Q, A, lenA, B, lenB, ctx);
+    else
+        return _gr_poly_divexact_bidirectional(Q, A, lenA, B, lenB, ctx);
+}
+
 /* basecase -> Newton cutoffs */
 /* todo: better unbalanced tuning */
 
@@ -1434,6 +1447,7 @@ gr_method_tab_input __gr_nmod_methods_input[] =
     {GR_METHOD_VEC_RECIPROCALS, (gr_funcptr) _gr_nmod_vec_reciprocals},
     {GR_METHOD_POLY_MULLOW,     (gr_funcptr) _gr_nmod_poly_mullow},
     {GR_METHOD_POLY_DIVREM,     (gr_funcptr) _gr_nmod_poly_divrem},
+    {GR_METHOD_POLY_DIVEXACT,   (gr_funcptr) _gr_nmod_poly_divexact},
     {GR_METHOD_POLY_INV_SERIES, (gr_funcptr) _gr_nmod_poly_inv_series},
     {GR_METHOD_POLY_INV_SERIES_BASECASE, (gr_funcptr) _gr_nmod_poly_inv_series_basecase},
     {GR_METHOD_POLY_DIV_SERIES, (gr_funcptr) _gr_nmod_poly_div_series},

src/gr_mat/minpoly_field.c

@@ -165,6 +165,8 @@ gr_mat_minpoly_field(gr_poly_t p, const gr_mat_t X, gr_ctx_t ctx)
         }
         _gr_poly_set_length(b, r1 + 1, ctx);
 
+        /* todo: poly_divexact */
+        /* todo: compute as (p * b) / g or (p / g) * b or p * (g / b) ? */
         status |= gr_poly_gcd(g, p, b, ctx);
         status |= gr_poly_mul(p, p, b, ctx);
         status |= gr_poly_divrem(p, r, p, g, ctx);

src/gr_poly/divexact_basecase.c

@@ -48,6 +48,11 @@ _gr_poly_divexact_basecase(gr_ptr Q,
     slong sz = ctx->sizeof_elem;
     gr_ptr invB;
 
+    /* note: if we wanted to be clever, we could pick a pair of coefficients such
+             that the division is easy */
+    if (Alen == Blen)
+        return gr_divexact(Q, GR_ENTRY(A, Alen - 1, sz), GR_ENTRY(B, Blen - 1, sz), ctx);
+
     GR_TMP_INIT(invB, ctx);
 
     /* todo: we sometimes want to keep dividing, e.g. over RR with small coefficient */

src/n_poly.h

@@ -410,6 +410,22 @@ void _n_poly_mod_div(n_poly_t Q, const n_poly_t A, const n_poly_t B, nmod_t mod)
     Q->length = lenA - lenB + 1;
 }
 
+FLINT_FORCE_INLINE
+void _n_poly_mod_divexact(n_poly_t Q, const n_poly_t A, const n_poly_t B, nmod_t mod)
+{
+    const slong lenA = A->length, lenB = B->length;
+    FLINT_ASSERT(lenB > 0);
+    FLINT_ASSERT(Q != A && Q != B);
+    if (lenA < lenB)
+    {
+        n_poly_zero(Q);
+        return;
+    }
+    n_poly_fit_length(Q, lenA - lenB + 1);
+    _nmod_poly_divexact(Q->coeffs, A->coeffs, lenA, B->coeffs, lenB, mod);
+    Q->length = lenA - lenB + 1;
+}
+
 FLINT_FORCE_INLINE
 void _n_poly_mod_rem(n_poly_t R, const n_poly_t A, const n_poly_t B, nmod_t mod)
 {
@@ -472,6 +488,9 @@ void n_poly_mod_rem(n_poly_t R, const n_poly_t A, const n_poly_t B,
 void n_poly_mod_divrem(n_poly_t Q, n_poly_t R, const n_poly_t A,
                                                  const n_poly_t B, nmod_t mod);
 
+void n_poly_mod_divexact(n_poly_t Q, const n_poly_t A, const n_poly_t B,
+                                                                   nmod_t mod);
+
 void n_poly_mod_mulmod(n_poly_t res, const n_poly_t poly1,
                            const n_poly_t poly2, const n_poly_t f, nmod_t mod);
 

src/n_poly/n_bpoly_mod.c

@@ -81,7 +81,7 @@ void n_bpoly_mod_divexact_last(n_bpoly_t A, const n_poly_t b, nmod_t ctx)
         if (A->coeffs[i].length < 1)
             continue;
 
-        n_poly_mod_div(t, A->coeffs + i, b, ctx);
+        n_poly_mod_divexact(t, A->coeffs + i, b, ctx);
         n_poly_swap(A->coeffs + i, t);
     }
 }

src/n_poly/n_bpoly_mod_gcd.c

@@ -225,8 +225,8 @@ int n_bpoly_mod_gcd_brown_smprime(
 
     n_poly_mod_gcd(cG, cA, cB, ctx);
 
-    n_poly_mod_div(cAbar, cA, cG, ctx);
-    n_poly_mod_div(cBbar, cB, cG, ctx);
+    n_poly_mod_divexact(cAbar, cA, cG, ctx);
+    n_poly_mod_divexact(cBbar, cB, cG, ctx);
 
     n_poly_mod_gcd(gamma, A->coeffs + A->length - 1, B->coeffs + B->length - 1, ctx);
 
@@ -311,11 +311,11 @@ int n_bpoly_mod_gcd_brown_smprime(
     else
     {
         n_poly_mod_gcd(Gevalp, Aevalp, Bevalp, ctx);
-        n_poly_mod_div(Abarevalp, Aevalp, Gevalp, ctx);
-        n_poly_mod_div(Bbarevalp, Bevalp, Gevalp, ctx);
+        n_poly_mod_divexact(Abarevalp, Aevalp, Gevalp, ctx);
+        n_poly_mod_divexact(Bbarevalp, Bevalp, Gevalp, ctx);
         n_poly_mod_gcd(Gevalm, Aevalm, Bevalm, ctx);
-        n_poly_mod_div(Abarevalm, Aevalm, Gevalm, ctx);
-        n_poly_mod_div(Bbarevalm, Bevalm, Gevalm, ctx);
+        n_poly_mod_divexact(Abarevalm, Aevalm, Gevalm, ctx);
+        n_poly_mod_divexact(Bbarevalm, Bevalm, Gevalm, ctx);
         gstab = astab = bstab = 0;
     }
 

src/n_poly/n_poly_mod.c

@@ -339,6 +339,57 @@ void n_poly_mod_div(n_poly_t Q, const n_poly_t A, const n_poly_t B, nmod_t ctx)
     Q->length = A_len - B_len + 1;
 }
 
+void n_poly_mod_divexact(n_poly_t Q, const n_poly_t A, const n_poly_t B, nmod_t ctx)
+{
+    n_poly_t tQ;
+    mp_ptr q;
+    slong A_len, B_len;
+
+    B_len = B->length;
+
+    if (B_len == 0)
+    {
+        if (ctx.n == 1)
+        {
+            n_poly_set(Q, A);
+            return;
+        }
+        else
+        {
+            flint_throw(FLINT_ERROR, "Exception (n_poly_mod_divexact). Division by zero.\n");
+        }
+    }
+
+    A_len = A->length;
+
+    if (A_len < B_len)
+    {
+        n_poly_zero(Q);
+        return;
+    }
+
+    if (Q == A || Q == B)
+    {
+        n_poly_init2(tQ, A_len - B_len + 1);
+        q = tQ->coeffs;
+    }
+    else
+    {
+        n_poly_fit_length(Q, A_len - B_len + 1);
+        q = Q->coeffs;
+    }
+
+    _nmod_poly_divexact(q, A->coeffs, A_len, B->coeffs, B_len, ctx);
+
+    if (Q == A || Q == B)
+    {
+        n_poly_swap(tQ, Q);
+        n_poly_clear(tQ);
+    }
+
+    Q->length = A_len - B_len + 1;
+}
+
 void n_poly_mod_rem(n_poly_t R, const n_poly_t A, const n_poly_t B, nmod_t ctx)
 {
     const slong lenA = A->length, lenB = B->length;

src/n_poly/n_polyu1n_gcd.c

@@ -298,8 +298,8 @@ int n_polyu1n_mod_gcd_brown_smprime(
     _n_poly_vec_mod_remove_content(cB, B->coeffs, B->length, ctx);
 
     n_poly_mod_gcd(cG, cA, cB, ctx);
-    n_poly_mod_div(cAbar, cA, cG, ctx);
-    n_poly_mod_div(cBbar, cB, cG, ctx);
+    n_poly_mod_divexact(cAbar, cA, cG, ctx);
+    n_poly_mod_divexact(cBbar, cB, cG, ctx);
 
     n_poly_mod_gcd(gamma, A->coeffs + 0, B->coeffs + 0, ctx);
 
@@ -377,11 +377,11 @@ int n_polyu1n_mod_gcd_brown_smprime(
     else
     {
         n_poly_mod_gcd(Gevalp, Aevalp, Bevalp, ctx);
-        n_poly_mod_div(Abarevalp, Aevalp, Gevalp, ctx);
-        n_poly_mod_div(Bbarevalp, Bevalp, Gevalp, ctx);
+        n_poly_mod_divexact(Abarevalp, Aevalp, Gevalp, ctx);
+        n_poly_mod_divexact(Bbarevalp, Bevalp, Gevalp, ctx);
         n_poly_mod_gcd(Gevalm, Aevalm, Bevalm, ctx);
-        n_poly_mod_div(Abarevalm, Aevalm, Gevalm, ctx);
-        n_poly_mod_div(Bbarevalm, Bevalm, Gevalm, ctx);
+        n_poly_mod_divexact(Abarevalm, Aevalm, Gevalm, ctx);
+        n_poly_mod_divexact(Bbarevalm, Bevalm, Gevalm, ctx);
         gstab = astab = bstab = 0;
     }
 

src/nmod_mat/minpoly.c

@@ -164,8 +164,9 @@ void nmod_mat_minpoly_with_gens(nmod_poly_t p, const nmod_mat_t X, ulong * P)
       _nmod_poly_set_length(b, r1 + 1);
 
       nmod_poly_gcd(g, p, b);
+      /* todo: compute as (p * b) / g or (p / g) * b or p * (g / b) ? */
       nmod_poly_mul(p, p, b);
-      nmod_poly_div(p, p, g);
+      nmod_poly_divexact(p, p, g);
 
       if (first_poly && r2 < n)
       {