Merge pull request #1759 from flintlib/divexact

nmod_poly_divexact

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)
       {

src/nmod_mpoly/gcd.c

@@ -960,14 +960,14 @@ int _do_univar(
                                                  I->Gmin_exp, I->Gstride, ctx);
     if (Abar != NULL)
     {
-        nmod_poly_div(t, a, g);
+        nmod_poly_divexact(t, a, g);
         _nmod_mpoly_from_nmod_poly_inflate(Abar, I->Abarbits, t, v_in_both,
                                               I->Abarmin_exp, I->Gstride, ctx);
     }
 
     if (Bbar != NULL)
     {
-        nmod_poly_div(t, b, g);
+        nmod_poly_divexact(t, b, g);
         _nmod_mpoly_from_nmod_poly_inflate(Bbar, I->Bbarbits, t, v_in_both,
                                               I->Bbarmin_exp, I->Gstride, ctx);
     }

src/nmod_mpoly/gcd_brown.c

@@ -155,11 +155,11 @@ static void _splitworker_bivar(void * varg)
         else
         {
             n_poly_mod_gcd(Gevalp, Aevalp, Bevalp, ctx->mod);
-            n_poly_mod_div(Abarevalp, Aevalp, Gevalp, ctx->mod);
-            n_poly_mod_div(Bbarevalp, Bevalp, Gevalp, ctx->mod);
+            n_poly_mod_divexact(Abarevalp, Aevalp, Gevalp, ctx->mod);
+            n_poly_mod_divexact(Bbarevalp, Bevalp, Gevalp, ctx->mod);
             n_poly_mod_gcd(Gevalm, Aevalm, Bevalm, ctx->mod);
-            n_poly_mod_div(Abarevalm, Aevalm, Gevalm, ctx->mod);
-            n_poly_mod_div(Bbarevalm, Bevalm, Gevalm, ctx->mod);
+            n_poly_mod_divexact(Abarevalm, Aevalm, Gevalm, ctx->mod);
+            n_poly_mod_divexact(Bbarevalm, Bevalm, Gevalm, ctx->mod);
         }
 
         FLINT_ASSERT(Gevalp->length > 0);
@@ -874,8 +874,8 @@ int nmod_mpolyn_gcd_brown_smprime_threaded_pool(
 
     n_poly_init(cAbar);
     n_poly_init(cBbar);
-    n_poly_mod_div(cAbar, cA, cG, ctx->mod);
-    n_poly_mod_div(cBbar, cB, cG, ctx->mod);
+    n_poly_mod_divexact(cAbar, cA, cG, ctx->mod);
+    n_poly_mod_divexact(cBbar, cB, cG, ctx->mod);
 
     n_poly_init(gamma);
     n_poly_mod_gcd(gamma, nmod_mpolyn_leadcoeff_poly(A, ctx),

src/nmod_mpoly/mpolyn_divides_threaded.c

@@ -1817,7 +1817,7 @@ int nmod_mpolyn_divides_threaded_pool(
     qexps = (ulong *) TMP_ALLOC(N*sizeof(ulong));
 
     mpoly_monomial_sub(qexps + N*0, A->exps + N*0, B->exps + N*0, N);
-    n_poly_mod_div(qcoeff, A->coeffs + 0, B->coeffs + 0, ctx->mod); /* already checked */
+    n_poly_mod_divexact(qcoeff, A->coeffs + 0, B->coeffs + 0, ctx->mod); /* already checked */
 
     nmod_mpolyn_ts_init(H->polyQ, qcoeff, qexps, 1, H->bits, H->N, ctx);
 

src/nmod_mpoly/mpolyn_gcd_brown.c

@@ -95,8 +95,8 @@ int nmod_mpolyn_gcd_brown_smprime_bivar(
 
     n_poly_mod_gcd(cG, cA, cB, ctx->mod);
 
-    n_poly_mod_div(cAbar, cA, cG, ctx->mod);
-    n_poly_mod_div(cBbar, cB, cG, ctx->mod);
+    n_poly_mod_divexact(cAbar, cA, cG, ctx->mod);
+    n_poly_mod_divexact(cBbar, cB, cG, ctx->mod);
 
     n_poly_mod_gcd(gamma, nmod_mpolyn_leadcoeff_poly(A, ctx),
                           nmod_mpolyn_leadcoeff_poly(B, ctx), ctx->mod);
@@ -186,11 +186,11 @@ int nmod_mpolyn_gcd_brown_smprime_bivar(
     else
     {
         n_poly_mod_gcd(Gevalp, Aevalp, Bevalp, ctx->mod);
-        n_poly_mod_div(Abarevalp, Aevalp, Gevalp, ctx->mod);
-        n_poly_mod_div(Bbarevalp, Bevalp, Gevalp, ctx->mod);
+        n_poly_mod_divexact(Abarevalp, Aevalp, Gevalp, ctx->mod);
+        n_poly_mod_divexact(Bbarevalp, Bevalp, Gevalp, ctx->mod);
         n_poly_mod_gcd(Gevalm, Aevalm, Bevalm, ctx->mod);
-        n_poly_mod_div(Abarevalm, Aevalm, Gevalm, ctx->mod);
-        n_poly_mod_div(Bbarevalm, Bevalm, Gevalm, ctx->mod);
+        n_poly_mod_divexact(Abarevalm, Aevalm, Gevalm, ctx->mod);
+        n_poly_mod_divexact(Bbarevalm, Bevalm, Gevalm, ctx->mod);
         gstab = astab = bstab = 0;
     }
 
@@ -412,8 +412,8 @@ int nmod_mpolyn_gcd_brown_smprime(
 
     n_poly_mod_gcd(cG, cA, cB, ctx->mod);
 
-    n_poly_mod_div(cAbar, cA, cG, ctx->mod);
-    n_poly_mod_div(cBbar, cB, cG, ctx->mod);
+    n_poly_mod_divexact(cAbar, cA, cG, ctx->mod);
+    n_poly_mod_divexact(cBbar, cB, cG, ctx->mod);
 
     n_poly_mod_gcd(gamma, nmod_mpolyn_leadcoeff_poly(A, ctx),
                           nmod_mpolyn_leadcoeff_poly(B, ctx), ctx->mod);
@@ -739,8 +739,8 @@ int nmod_mpolyn_gcd_brown_lgprime_bivar(
 
     n_poly_init(cAbar);
     n_poly_init(cBbar);
-    n_poly_mod_div(cAbar, cA, cG, ctx->mod);
-    n_poly_mod_div(cBbar, cB, cG, ctx->mod);
+    n_poly_mod_divexact(cAbar, cA, cG, ctx->mod);
+    n_poly_mod_divexact(cBbar, cB, cG, ctx->mod);
 
     n_poly_init(gamma);
     n_poly_mod_gcd(gamma, nmod_mpolyn_leadcoeff_poly(A, ctx),
@@ -983,8 +983,8 @@ int nmod_mpolyn_gcd_brown_lgprime(
 
     n_poly_init(cAbar);
     n_poly_init(cBbar);
-    n_poly_mod_div(cAbar, cA, cG, ctx->mod);
-    n_poly_mod_div(cBbar, cB, cG, ctx->mod);
+    n_poly_mod_divexact(cAbar, cA, cG, ctx->mod);
+    n_poly_mod_divexact(cBbar, cB, cG, ctx->mod);
 
     n_poly_init(gamma);
     n_poly_mod_gcd(gamma, nmod_mpolyn_leadcoeff_poly(A, ctx),

src/nmod_mpoly/mpolyu_gcdp_zippel.c

@@ -628,9 +628,9 @@ static int nmod_mpolyu_gcdp_zippel_univar(
     nmod_mpolyu_cvtto_poly(b, B, ctx);
     nmod_poly_gcd(g, a, b);
     nmod_mpolyu_cvtfrom_poly(G, g, ctx);
-    nmod_poly_div(t, a, g);
+    nmod_poly_divexact(t, a, g);
     nmod_mpolyu_cvtfrom_poly(Abar, t, ctx);
-    nmod_poly_div(t, b, g);
+    nmod_poly_divexact(t, b, g);
     nmod_mpolyu_cvtfrom_poly(Bbar, t, ctx);
     nmod_poly_clear(a);
     nmod_poly_clear(b);

src/nmod_mpoly_factor/lcc_wang.c

@@ -88,7 +88,7 @@ int nmod_mpoly_factor_lcc_wang(
             while (n_poly_degree(R) > 0)
             {
                 n_poly_mod_gcd(R, R, Q, ctx->mod);
-                n_poly_mod_div(Q, Q, R, ctx->mod);
+                n_poly_mod_divexact(Q, Q, R, ctx->mod);
                 if (n_poly_degree(Q) < 1)
                 {
                     success = 0;

src/nmod_mpoly_factor/n_bpoly_mod_factor_smprime.c

@@ -437,7 +437,7 @@ static void _n_bpoly_mod_lift_build_steps(n_bpoly_mod_lift_t L, nmod_t ctx)
     for (k = 0; k < r; k++)
     {
         /* s[k] = (prod_{i!=k} B[i].coeffs[0])^-1 (mod B[k].coeffs[0]) */
-        n_poly_mod_div(t, A->coeffs + 0, B[k].coeffs + 0, ctx);
+        n_poly_mod_divexact(t, A->coeffs + 0, B[k].coeffs + 0, ctx);
         if (!n_poly_mod_invmod(s + k, t, B[k].coeffs + 0, ctx))
             flint_throw(FLINT_IMPINV, "n_bpoly_mod_lift: bad inverse");
 

src/nmod_mpoly_factor/n_bpoly_mod_hlift.c

@@ -140,7 +140,7 @@ int n_bpoly_mod_hlift2_cubic(
         n_poly_mod_rem(u, t, B0->coeffs + 0, ctx);
         n_poly_mod_mul(t, u, B1->coeffs + 0, ctx);
         n_poly_mod_sub(c, c, t, ctx);
-        n_poly_mod_div(v, c, B0->coeffs + 0, ctx);
+        n_poly_mod_divexact(v, c, B0->coeffs + 0, ctx);
 
         if (!n_poly_is_zero(u))
         {
@@ -275,7 +275,7 @@ int n_bpoly_mod_hlift2(
         n_poly_mod_rem(u, t, B0->coeffs + 0, ctx);
         n_poly_mod_mul(t, u, B1->coeffs + 0, ctx);
         n_poly_mod_sub(c, c, t, ctx);
-        n_poly_mod_div(v, c, B0->coeffs + 0, ctx);
+        n_poly_mod_divexact(v, c, B0->coeffs + 0, ctx);
 
         if (j < B0->length)
             n_poly_mod_add(B0->coeffs + j, B0->coeffs + j, u, ctx);
@@ -430,7 +430,7 @@ int n_bpoly_mod_hlift_cubic(
     for (k = 0; k < r; k++)
     {
         /* s[k] = (prod_{i!=k} B[i].coeffs[0])^-1 (mod B[k].coeffs[0]) */
-        n_poly_mod_div(t, A->coeffs + 0, B[k].coeffs + 0, ctx);
+        n_poly_mod_divexact(t, A->coeffs + 0, B[k].coeffs + 0, ctx);
         if (!n_poly_mod_invmod(s[k], t, B[k].coeffs + 0, ctx))
         {
             success = -1;