diff --git a/src/fmpz_mpoly/test/t-scalar_divexact_fmpz.c b/src/fmpz_mpoly/test/t-scalar_divexact_fmpz.c
index 5e43e46e0b..c68e836351 100644
--- a/src/fmpz_mpoly/test/t-scalar_divexact_fmpz.c
+++ b/src/fmpz_mpoly/test/t-scalar_divexact_fmpz.c
@@ -92,7 +92,7 @@ TEST_FUNCTION_START(fmpz_mpoly_scalar_divexact_fmpz, state)
fmpz_mpoly_randtest_bits(f, state, len, coeff_bits, exp_bits, ctx);
fmpz_mpoly_randtest_bits(h, state, len, coeff_bits, exp_bits, ctx);
- fmpz_randtest(c, state, n_randint(state, 200));
+ fmpz_randtest_not_zero(c, state, n_randint(state, 200) + 1);
fmpz_mpoly_scalar_mul_fmpz(f, f, c, ctx);
diff --git a/src/fmpz_vec/profile/p-divexact_ui.c b/src/fmpz_vec/profile/p-divexact_ui.c
new file mode 100644
index 0000000000..8250965349
--- /dev/null
+++ b/src/fmpz_vec/profile/p-divexact_ui.c
@@ -0,0 +1,60 @@
+#include "ulong_extras.h"
+#include "fmpz.h"
+#include "fmpz_vec.h"
+#include "profiler.h"
+
+void
+_fmpz_vec_scalar_divexact_ui_naive(fmpz * vec1, const fmpz * vec2,
+ slong len2, ulong c)
+{
+ slong i;
+ for (i = 0; i < len2; i++)
+ fmpz_divexact_ui(vec1 + i, vec2 + i, c);
+}
+
+int main()
+{
+ slong len, bits;
+
+ flint_rand_t state;
+ flint_rand_init(state);
+
+ fmpz *A, *Ac, *B;
+ ulong c;
+
+ double t1, t2, FLINT_SET_BUT_UNUSED(__);
+
+ flint_printf(" len bits(A*c) bits(c) old new speedup\n");
+
+ for (len = 1; len <= 10; len++)
+ {
+ for (bits = 5; bits <= 10000; bits *= 2)
+ {
+ A = _fmpz_vec_init(len);
+ Ac = _fmpz_vec_init(len);
+ B = _fmpz_vec_init(len);
+
+ _fmpz_vec_randtest(A, state, len, bits);
+ c = n_randtest_not_zero(state);
+ _fmpz_vec_scalar_mul_ui(Ac, A, len, c);
+
+ TIMEIT_START
+ _fmpz_vec_scalar_divexact_ui_naive(B, Ac, len, c);
+ TIMEIT_STOP_VALUES(__, t1)
+ TIMEIT_START
+ _fmpz_vec_scalar_divexact_ui(B, Ac, len, c);
+ TIMEIT_STOP_VALUES(__, t2)
+
+ flint_printf("%6wd %6wd %6wd %8g %8g %.3f\n",
+ len, _fmpz_vec_max_bits(Ac, len), FLINT_BIT_COUNT(c),
+ t1, t2, t1 / t2);
+
+ _fmpz_vec_clear(A, len);
+ _fmpz_vec_clear(Ac, len);
+ _fmpz_vec_clear(B, len);
+ }
+ }
+
+ flint_rand_clear(state);
+}
+
diff --git a/src/fmpz_vec/scalar.c b/src/fmpz_vec/scalar.c
index 7f8c9db7c0..a56a440a57 100644
--- a/src/fmpz_vec/scalar.c
+++ b/src/fmpz_vec/scalar.c
@@ -126,46 +126,6 @@ _fmpz_vec_scalar_addmul_ui(fmpz * vec1, const fmpz * vec2, slong len2, ulong c)
fmpz_addmul_ui(vec1 + ix, vec2 + ix, c);
}
-void
-_fmpz_vec_scalar_divexact_fmpz(fmpz * vec1, const fmpz * vec2,
- slong len2, const fmpz_t x)
-{
- fmpz c = *x;
-
- if (!COEFF_IS_MPZ(c))
- {
- if (c == 1)
- _fmpz_vec_set(vec1, vec2, len2);
- else if (c == -1)
- _fmpz_vec_neg(vec1, vec2, len2);
- else
- _fmpz_vec_scalar_divexact_si(vec1, vec2, len2, c);
- }
- else
- {
- slong i;
- for (i = 0; i < len2; i++)
- fmpz_divexact(vec1 + i, vec2 + i, x);
- }
-}
-
-void
-_fmpz_vec_scalar_divexact_si(fmpz * vec1, const fmpz * vec2, slong len2, slong c)
-{
- slong i;
- for (i = 0; i < len2; i++)
- fmpz_divexact_si(vec1 + i, vec2 + i, c);
-}
-
-void
-_fmpz_vec_scalar_divexact_ui(fmpz * vec1, const fmpz * vec2,
- slong len2, ulong c)
-{
- slong i;
- for (i = 0; i < len2; i++)
- fmpz_divexact_ui(vec1 + i, vec2 + i, c);
-}
-
void
_fmpz_vec_scalar_fdiv_q_2exp(fmpz * vec1, const fmpz * vec2, slong len2,
ulong exp)
@@ -209,6 +169,8 @@ _fmpz_vec_scalar_fdiv_r_2exp(fmpz * vec1, const fmpz * vec2, slong len2,
fmpz_fdiv_r_2exp(vec1 + i, vec2 + i, exp);
}
+/* todo: preinvert p when appropriate; or better, make sure we don't
+ call this function in places where we have an fmpz_mod */
void _fmpz_vec_scalar_mod_fmpz(fmpz *res, const fmpz *vec, slong len, const fmpz_t p)
{
slong i;
@@ -267,6 +229,7 @@ _fmpz_vec_scalar_mul_ui(fmpz * vec1, const fmpz * vec2, slong len2, ulong c)
fmpz_mul_ui(vec1 + i, vec2 + i, c);
}
+/* todo: preinvert p when appropriate */
void _fmpz_vec_scalar_smod_fmpz(fmpz *res, const fmpz *vec, slong len, const fmpz_t p)
{
slong i;
diff --git a/src/fmpz_vec/scalar_divexact.c b/src/fmpz_vec/scalar_divexact.c
new file mode 100644
index 0000000000..6815e811ab
--- /dev/null
+++ b/src/fmpz_vec/scalar_divexact.c
@@ -0,0 +1,432 @@
+/*
+ Copyright 2000-2003, 2005, 2013 Free Software Foundation, Inc.
+ Copyright 1991-2018, 2021, 2022 Free Software Foundation, Inc.
+ Copyright (C) 2025 Fredrik Johansson
+
+ This file is part of FLINT.
+
+ FLINT is free software: you can redistribute it and/or modify it under
+ the terms of the GNU Lesser General Public License (LGPL) as published
+ by the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version. See .
+*/
+
+/* This file includes code adapted from GMP (gmp-impl.h, mpn/generic/dive_1.c). */
+
+#include
+#include "longlong.h"
+#include "fmpz.h"
+#include "fmpz_vec.h"
+
+/* Helper functions that could go in longlong.h, ulong_extras.h
+ and mpn_extras.h as soon as they find a use elsewhere. */
+
+/* Tuned for Zen 3; with assembly code we would not need GMP. */
+#define DIVEXACT_1_ODD_GMP_CUTOFF 35
+#define DIVEXACT_1_EVEN_GMP_CUTOFF 25
+
+/* The table is slightly faster but requires GMP internals. */
+#define USE_BINVERT_TABLE 0
+
+#if USE_BINVERT_TABLE
+extern const unsigned char __gmp_binvert_limb_table[128];
+
+/* Invert n mod 2^FLINT_BITS */
+static ulong n_binvert(ulong n)
+{
+ ulong r;
+
+ FLINT_ASSERT(n % 2);
+
+ r = __gmp_binvert_limb_table[(n / 2) & 0x7F];
+ r = 2 * r - r * r * n; /* 8 bits */
+ r = 2 * r - r * r * n; /* 16 bits */
+ r = 2 * r - r * r * n; /* 32 bits */
+#if FLINT_BITS == 64
+ r = 2 * r - r * r * n; /* 64 bits */
+#endif
+
+ return r;
+}
+#else
+
+/* Use Hurchalla's algorithm https://arxiv.org/abs/2204.04342 */
+static ulong n_binvert(ulong a)
+{
+ ulong r, y;
+
+ r = (3 * a) ^ 2; /* 5 bits */
+ y = 1 - a * r;
+ r = r * (1 + y); /* 10 bits */
+ y *= y;
+ r = r * (1 + y); /* 20 bits */
+ y *= y;
+ r = r * (1 + y); /* 40 bits */
+#if FLINT_BITS == 64
+ y *= y;
+ r = r * (1 + y); /* 80 bits */
+#endif
+ return r;
+}
+#endif
+
+FLINT_FORCE_INLINE
+ulong n_mulhi(ulong a, ulong b)
+{
+ ulong h, l;
+ umul_ppmm(h, l, a, b);
+ return h;
+}
+
+FLINT_FORCE_INLINE
+ulong n_subc(ulong * c, ulong a, ulong b)
+{
+ ulong d = a - b;
+ *c = d > a; /* the compiler should recognize this idiom */
+ return d;
+}
+
+static
+void nn_divexact_2_1_odd(nn_ptr res, nn_srcptr x, ulong d, ulong dinv)
+{
+ ulong s, l, c;
+
+ s = x[0];
+ l = s * dinv;
+ res[0] = l;
+ c = n_mulhi(l, d);
+ s = x[1];
+ l = s - c;
+ l = l * dinv;
+ res[1] = l;
+}
+
+static
+void nn_divexact_2_1_even(nn_ptr res, nn_srcptr x, ulong d, ulong dinv, unsigned int norm)
+{
+ ulong s, l, c;
+ ulong ls, s_next;
+ ulong h;
+
+ d >>= norm;
+ c = 0;
+ s = x[0];
+ s_next = x[1];
+ ls = (s >> norm) | (s_next << (FLINT_BITS - norm));
+ s = s_next;
+ l = n_subc(&c, ls, c);
+ l = l * dinv;
+ res[0] = l;
+ h = n_mulhi(l, d);
+ c += h;
+ ls = s >> norm;
+ l = ls - c;
+ l = l * dinv;
+ res[1] = l;
+}
+
+static
+void nn_divexact_1_odd(nn_ptr res, nn_srcptr x, slong xn, ulong d, ulong dinv)
+{
+ ulong s, l, c;
+ slong i;
+
+ if (xn >= DIVEXACT_1_ODD_GMP_CUTOFF)
+ {
+ mpn_divexact_1(res, x, xn, d);
+ return;
+ }
+
+ s = x[0];
+ l = s * dinv;
+ res[0] = l;
+ c = 0;
+
+ for (i = 1; i < xn; i++)
+ {
+ c += n_mulhi(l, d);
+ s = x[i];
+ l = n_subc(&c, s, c);
+ l = l * dinv;
+ res[i] = l;
+ }
+}
+
+static
+void nn_divexact_1_even(nn_ptr res, nn_srcptr x, slong xn, ulong d, ulong dinv, unsigned int norm)
+{
+ ulong s, l, c, h;
+ slong i;
+ ulong ls, s_next;
+
+ if (xn >= DIVEXACT_1_EVEN_GMP_CUTOFF)
+ {
+ mpn_divexact_1(res, x, xn, d);
+ return;
+ }
+
+ d >>= norm;
+ c = 0;
+ s = x[0];
+
+ for (i = 1; i < xn; i++)
+ {
+ s_next = x[i];
+ ls = (s >> norm) | (s_next << (FLINT_BITS - norm));
+ s = s_next;
+ l = n_subc(&c, ls, c);
+ l = l * dinv;
+ res[i - 1] = l;
+ h = n_mulhi(l, d);
+ c += h;
+ }
+
+ ls = s >> norm;
+ l = ls - c;
+ l = l * dinv;
+ res[xn - 1] = l;
+}
+
+#define FLINT_MPZ_GET_MPN(zd, zn, zsgnbit, z) \
+ do { \
+ (zsgnbit) = (z)->_mp_size < 0; \
+ (zn) = FLINT_ABS((z)->_mp_size); \
+ (zd) = (z)->_mp_d; \
+ } \
+ while (0)
+
+/* Disable for a small speedup, which however risks creating corrupted fmpzs
+ if the user somehow did not provide exact quotients. */
+#define SAFE_VERSION 1
+
+static void
+_fmpz_vec_divexact_ui(fmpz * res, const fmpz * x, slong len, ulong c, int negc)
+{
+ ulong cinv;
+ slong i;
+ unsigned int norm;
+ mpz_ptr z;
+ slong xn;
+ int negative;
+ nn_srcptr xd;
+ ulong cinv_signed;
+
+ FLINT_ASSERT(c != 0);
+
+ if (c == 1)
+ {
+ if (negc)
+ _fmpz_vec_neg(res, x, len);
+ else
+ _fmpz_vec_set(res, x, len);
+ return;
+ }
+
+ if (c & 1)
+ {
+ cinv = n_binvert(c);
+ cinv_signed = negc ? -cinv : cinv;
+
+ for (i = 0; i < len; i++)
+ {
+ slong v = x[i];
+
+ if (!COEFF_IS_MPZ(v))
+ {
+#if SAFE_VERSION
+ fmpz_set_si(res + i, v * cinv_signed);
+#else
+ _fmpz_demote(res + i);
+ res[i] = v * cinv_signed;
+#endif
+ }
+ else
+ {
+ z = COEFF_TO_PTR(v);
+ FLINT_MPZ_GET_MPN(xd, xn, negative, z);
+ negative ^= negc;
+
+ if (xn == 1)
+ {
+ /* cannot overflow if c >= 2 */
+ fmpz_set_si(res + i, negative ? -xd[0] * cinv : xd[0] * cinv);
+ }
+ else if (xn == 2)
+ {
+ ulong t[2];
+ nn_divexact_2_1_odd(t, xd, c, cinv);
+
+ if (negative)
+ fmpz_neg_uiui(res + i, t[1], t[0]);
+ else
+ fmpz_set_uiui(res + i, t[1], t[0]);
+ }
+ else
+ {
+ mpz_ptr rp = _fmpz_promote(res + i);
+ nn_ptr zd = FLINT_MPZ_REALLOC(rp, xn);
+
+ nn_divexact_1_odd(zd, xd, xn, c, cinv);
+
+ xn -= (zd[xn - 1] == 0);
+ FLINT_ASSERT(zd[xn - 1] != 0);
+ rp->_mp_size = negative ? -xn : xn;
+#if SAFE_VERSION
+ _fmpz_demote_val(res + i);
+#endif
+ }
+ }
+ }
+ }
+ else
+ {
+ norm = flint_ctz(c);
+
+ /* Special case for powers of 2 */
+ if ((c & (c - 1)) == 0)
+ {
+ if (negc)
+ {
+ for (i = 0; i < len; i++)
+ {
+ slong v = x[i];
+
+ if (!COEFF_IS_MPZ(v))
+ {
+ _fmpz_demote(res + i);
+ res[i] = (-v) >> norm;
+ }
+ else
+ {
+ fmpz_tdiv_q_2exp(res + i, x + i, norm);
+ fmpz_neg(res + i, res + i);
+ }
+ }
+ }
+ else
+ {
+ for (i = 0; i < len; i++)
+ {
+ slong v = x[i];
+
+ if (!COEFF_IS_MPZ(v))
+ {
+ _fmpz_demote(res + i);
+ res[i] = v >> norm;
+ }
+ else
+ fmpz_tdiv_q_2exp(res + i, x + i, norm);
+ }
+ }
+
+ return;
+ }
+
+ cinv = n_binvert(c >> norm);
+ cinv_signed = negc ? -cinv : cinv;
+
+ for (i = 0; i < len; i++)
+ {
+ slong v = x[i];
+
+ if (!COEFF_IS_MPZ(v))
+ {
+#if SAFE_VERSION
+ fmpz_set_si(res + i, ((slong) (v * cinv_signed)) >> norm);
+#else
+ _fmpz_demote(res + i);
+ res[i] = ((slong) (v * cinv_signed)) >> norm;
+#endif
+ }
+ else
+ {
+ z = COEFF_TO_PTR(v);
+ FLINT_MPZ_GET_MPN(xd, xn, negative, z);
+ negative ^= negc;
+
+ if (xn == 1)
+ {
+ if (negative)
+ fmpz_neg_ui(res + i, (xd[0] * cinv) >> norm);
+ else
+ fmpz_set_ui(res + i, (xd[0] * cinv) >> norm);
+ }
+ else if (xn == 2)
+ {
+ ulong t[2];
+ nn_divexact_2_1_even(t, xd, c, cinv, norm);
+
+ if (negative)
+ fmpz_neg_uiui(res + i, t[1], t[0]);
+ else
+ fmpz_set_uiui(res + i, t[1], t[0]);
+ }
+ else
+ {
+ mpz_ptr rp = _fmpz_promote(res + i);
+ nn_ptr zd = FLINT_MPZ_REALLOC(rp, xn);
+
+ nn_divexact_1_even(zd, xd, xn, c, cinv, norm);
+ xn -= (zd[xn - 1] == 0);
+ FLINT_ASSERT(zd[xn - 1] != 0);
+ rp->_mp_size = negative ? -xn : xn;
+#if SAFE_VERSION
+ _fmpz_demote_val(res + i);
+#endif
+ }
+ }
+ }
+ }
+}
+
+void
+_fmpz_vec_scalar_divexact_si(fmpz * vec1, const fmpz * vec2, slong len2, slong c)
+{
+ if (len2 == 1)
+ {
+ fmpz_divexact_si(vec1, vec2, c);
+ return;
+ }
+
+ if (c > 0)
+ _fmpz_vec_divexact_ui(vec1, vec2, len2, c, 0);
+ else
+ _fmpz_vec_divexact_ui(vec1, vec2, len2, -(ulong) c, 1);
+}
+
+void
+_fmpz_vec_scalar_divexact_ui(fmpz * vec1, const fmpz * vec2,
+ slong len2, ulong c)
+{
+ if (len2 == 1)
+ {
+ fmpz_divexact_ui(vec1, vec2, c);
+ return;
+ }
+
+ _fmpz_vec_divexact_ui(vec1, vec2, len2, c, 0);
+}
+
+void
+_fmpz_vec_scalar_divexact_fmpz(fmpz * vec1, const fmpz * vec2,
+ slong len2, const fmpz_t x)
+{
+ fmpz c = *x;
+
+ if (!COEFF_IS_MPZ(c))
+ {
+ if (c == 1)
+ _fmpz_vec_set(vec1, vec2, len2);
+ else if (c == -1)
+ _fmpz_vec_neg(vec1, vec2, len2);
+ else
+ _fmpz_vec_scalar_divexact_si(vec1, vec2, len2, c);
+ }
+ else
+ {
+ slong i;
+ for (i = 0; i < len2; i++)
+ fmpz_divexact(vec1 + i, vec2 + i, x);
+ }
+}
+
diff --git a/src/gr/fmpz.c b/src/gr/fmpz.c
index 8f0fee64e0..20c1e99872 100644
--- a/src/gr/fmpz.c
+++ b/src/gr/fmpz.c
@@ -878,6 +878,33 @@ _gr_fmpz_vec_sub(fmpz * res, const fmpz * vec1, const fmpz * vec2, slong len, gr
return GR_SUCCESS;
}
+int
+_gr_fmpz_vec_divexact_scalar_fmpz(fmpz * res, const fmpz * vec1, slong len, const fmpz_t c, gr_ctx_t ctx)
+{
+ if (fmpz_is_zero(c))
+ return GR_DOMAIN;
+ _fmpz_vec_scalar_divexact_fmpz(res, vec1, len, c);
+ return GR_SUCCESS;
+}
+
+int
+_gr_fmpz_vec_divexact_scalar_ui(fmpz * res, const fmpz * vec1, slong len, ulong c, gr_ctx_t ctx)
+{
+ if (c == 0)
+ return GR_DOMAIN;
+ _fmpz_vec_scalar_divexact_ui(res, vec1, len, c);
+ return GR_SUCCESS;
+}
+
+int
+_gr_fmpz_vec_divexact_scalar_si(fmpz * res, const fmpz * vec1, slong len, slong c, gr_ctx_t ctx)
+{
+ if (c == 0)
+ return GR_DOMAIN;
+ _fmpz_vec_scalar_divexact_si(res, vec1, len, c);
+ return GR_SUCCESS;
+}
+
int
_gr_fmpz_vec_sum(fmpz_t res, const fmpz * vec, slong len, gr_ctx_t ctx)
{
@@ -1207,6 +1234,10 @@ gr_method_tab_input _fmpz_methods_input[] =
{GR_METHOD_VEC_EQUAL, (gr_funcptr) _gr_fmpz_vec_equal},
{GR_METHOD_VEC_ADD, (gr_funcptr) _gr_fmpz_vec_add},
{GR_METHOD_VEC_SUB, (gr_funcptr) _gr_fmpz_vec_sub},
+ {GR_METHOD_VEC_DIVEXACT_SCALAR, (gr_funcptr) _gr_fmpz_vec_divexact_scalar_fmpz},
+ {GR_METHOD_VEC_DIVEXACT_SCALAR_FMPZ, (gr_funcptr) _gr_fmpz_vec_divexact_scalar_fmpz},
+ {GR_METHOD_VEC_DIVEXACT_SCALAR_UI, (gr_funcptr) _gr_fmpz_vec_divexact_scalar_ui},
+ {GR_METHOD_VEC_DIVEXACT_SCALAR_SI, (gr_funcptr) _gr_fmpz_vec_divexact_scalar_si},
{GR_METHOD_VEC_SUM, (gr_funcptr) _gr_fmpz_vec_sum},
{GR_METHOD_VEC_DOT, (gr_funcptr) _gr_fmpz_vec_dot},
{GR_METHOD_VEC_DOT_REV, (gr_funcptr) _gr_fmpz_vec_dot_rev},
diff --git a/src/gr/test_ring.c b/src/gr/test_ring.c
index 2341afbe93..52b89186fa 100644
--- a/src/gr/test_ring.c
+++ b/src/gr/test_ring.c
@@ -3458,9 +3458,105 @@ int gr_test_vec_add(gr_ctx_t R, flint_rand_t state, int test_flags) { return gr_
int gr_test_vec_sub(gr_ctx_t R, flint_rand_t state, int test_flags) { return gr_test_vec_binary_op(R, "vec_sub", gr_sub, _gr_vec_sub, state, test_flags); }
int gr_test_vec_mul(gr_ctx_t R, flint_rand_t state, int test_flags) { return gr_test_vec_binary_op(R, "vec_mul", gr_mul, _gr_vec_mul, state, test_flags); }
int gr_test_vec_div(gr_ctx_t R, flint_rand_t state, int test_flags) { return gr_test_vec_binary_op(R, "vec_div", gr_div, _gr_vec_div, state, test_flags); }
-int gr_test_vec_divexact(gr_ctx_t R, flint_rand_t state, int test_flags) { return gr_test_vec_binary_op(R, "vec_divexact", gr_divexact, _gr_vec_divexact, state, test_flags); }
int gr_test_vec_pow(gr_ctx_t R, flint_rand_t state, int test_flags) { return gr_test_vec_binary_op(R, "vec_pow", gr_pow, _gr_vec_pow, state, test_flags); }
+int gr_test_vec_divexact(gr_ctx_t R, flint_rand_t state, int test_flags)
+{
+ int status, aliasing, ref_aliasing;
+ slong i, len;
+ gr_ptr x, y, xy1, xy2;
+
+ len = n_randint(state, 5);
+
+ GR_TMP_INIT_VEC(x, len, R);
+ GR_TMP_INIT_VEC(y, len, R);
+ GR_TMP_INIT_VEC(xy1, len, R);
+ GR_TMP_INIT_VEC(xy2, len, R);
+
+ GR_MUST_SUCCEED(_gr_vec_randtest(x, state, len, R));
+ GR_MUST_SUCCEED(_gr_vec_randtest(y, state, len, R));
+ GR_MUST_SUCCEED(_gr_vec_randtest(xy1, state, len, R));
+ GR_MUST_SUCCEED(_gr_vec_randtest(xy2, state, len, R));
+
+ status = GR_SUCCESS;
+ for (i = 0; i < len && status == GR_SUCCESS; i++)
+ {
+ status |= gr_mul(GR_ENTRY(x, i, R->sizeof_elem),
+ GR_ENTRY(x, i, R->sizeof_elem),
+ GR_ENTRY(y, i, R->sizeof_elem), R);
+
+ /* Check that we can perform division */
+ status |= gr_div(xy1, GR_ENTRY(x, i, R->sizeof_elem),
+ GR_ENTRY(y, i, R->sizeof_elem), R);
+ }
+
+ if (status == GR_SUCCESS)
+ {
+ aliasing = n_randint(state, 4);
+ ref_aliasing = 0;
+
+ switch (aliasing)
+ {
+ case 0:
+ status |= _gr_vec_set(xy1, x, len, R);
+ status |= _gr_vec_divexact(xy1, xy1, y, len, R);
+ break;
+ case 1:
+ status |= _gr_vec_set(xy1, y, len, R);
+ status |= _gr_vec_divexact(xy1, x, xy1, len, R);
+ break;
+ case 2:
+ status |= _gr_vec_set(y, x, len, R);
+ status |= _gr_vec_divexact(xy1, x, x, len, R);
+ break;
+ case 3:
+ status |= _gr_vec_set(y, x, len, R);
+ status |= _gr_vec_set(xy1, x, len, R);
+ status |= _gr_vec_divexact(xy1, xy1, xy1, len, R);
+ break;
+ default:
+ status |= _gr_vec_divexact(xy1, x, y, len, R);
+ }
+
+ for (i = 0; i < len; i++)
+ if (ref_aliasing)
+ status |= gr_divexact(GR_ENTRY(xy2, i, R->sizeof_elem),
+ GR_ENTRY(x, i, R->sizeof_elem),
+ GR_ENTRY(x, i, R->sizeof_elem), R);
+ else
+ status |= gr_divexact(GR_ENTRY(xy2, i, R->sizeof_elem),
+ GR_ENTRY(x, i, R->sizeof_elem),
+ GR_ENTRY(y, i, R->sizeof_elem), R);
+
+ if (status == GR_SUCCESS && _gr_vec_equal(xy1, xy2, len, R) == T_FALSE)
+ {
+ status = GR_TEST_FAIL;
+ }
+
+ if ((test_flags & GR_TEST_ALWAYS_ABLE) && (status & GR_UNABLE))
+ status = GR_TEST_FAIL;
+
+ if ((test_flags & GR_TEST_VERBOSE) || status == GR_TEST_FAIL)
+ {
+ flint_printf("divexact\n");
+ gr_ctx_println(R);
+ flint_printf("aliasing: %d\n", aliasing);
+ _gr_vec_print(x, len, R); flint_printf("\n");
+ _gr_vec_print(y, len, R); flint_printf("\n");
+ _gr_vec_print(xy1, len, R); flint_printf("\n");
+ _gr_vec_print(xy2, len, R); flint_printf("\n");
+ }
+ }
+
+ GR_TMP_CLEAR_VEC(x, len, R);
+ GR_TMP_CLEAR_VEC(y, len, R);
+ GR_TMP_CLEAR_VEC(xy1, len, R);
+ GR_TMP_CLEAR_VEC(xy2, len, R);
+
+ return status;
+}
+
+
int
gr_test_vec_binary_op_scalar(gr_ctx_t R, const char * opname, int (*gr_op)(gr_ptr, gr_srcptr, gr_srcptr, gr_ctx_t),
int (*_gr_vec_op)(gr_ptr, gr_srcptr, slong, gr_srcptr, gr_ctx_t), flint_rand_t state, int test_flags)