diff --git a/doc/source/fmpz_mod_mpoly.rst b/doc/source/fmpz_mod_mpoly.rst
index 46823749f4..7a40c58829 100644
--- a/doc/source/fmpz_mod_mpoly.rst
+++ b/doc/source/fmpz_mod_mpoly.rst
@@ -436,6 +436,7 @@ Evaluation
     Return `1` for success and `0` for failure.
 
 .. function:: int fmpz_mod_mpoly_compose_fmpz_mod_mpoly_geobucket(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, fmpz_mod_mpoly_struct * const * C, const fmpz_mod_mpoly_ctx_t ctxB, const fmpz_mod_mpoly_ctx_t ctxAC)
+              int fmpz_mod_mpoly_compose_fmpz_mod_mpoly_horner(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, fmpz_mod_mpoly_struct * const * C, const fmpz_mod_mpoly_ctx_t ctxB, const fmpz_mod_mpoly_ctx_t ctxAC)
               int fmpz_mod_mpoly_compose_fmpz_mod_mpoly(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, fmpz_mod_mpoly_struct * const * C, const fmpz_mod_mpoly_ctx_t ctxB, const fmpz_mod_mpoly_ctx_t ctxAC)
 
     Set *A* to the evaluation of *B* where the variables are replaced by the corresponding elements of the array *C*.
diff --git a/src/fmpq_mpoly/test/t-compose_fmpq_mpoly.c b/src/fmpq_mpoly/test/t-compose_fmpq_mpoly.c
index 1ac4529a82..9a33e6c177 100644
--- a/src/fmpq_mpoly/test/t-compose_fmpq_mpoly.c
+++ b/src/fmpq_mpoly/test/t-compose_fmpq_mpoly.c
@@ -17,13 +17,15 @@ TEST_FUNCTION_START(fmpq_mpoly_compose_fmpq_mpoly, state)
     slong i, j, v;
 
     {
-        fmpq_mpoly_t A, A1, A2, B;
+        slong nvarsAC = 2, nvarsB = 3;
+        fmpq_mpoly_t A, A1, A2, B, B1;
         fmpq_mpoly_struct * Cp[3];
         fmpq_mpoly_struct C[3];
         fmpq_mpoly_ctx_t ctxAC, ctxB;
+        slong * c;
 
-        fmpq_mpoly_ctx_init(ctxB, 3, ORD_LEX);
-        fmpq_mpoly_ctx_init(ctxAC, 2, ORD_LEX);
+        fmpq_mpoly_ctx_init(ctxB, nvarsB, ORD_LEX);
+        fmpq_mpoly_ctx_init(ctxAC, nvarsAC, ORD_LEX);
 
         fmpq_mpoly_init(B, ctxB);
         fmpq_mpoly_init(A, ctxAC);
@@ -42,39 +44,61 @@ TEST_FUNCTION_START(fmpq_mpoly_compose_fmpq_mpoly, state)
         fmpq_mpoly_set_str_pretty(C + 1, "x1 - x2", NULL, ctxAC);
         fmpq_mpoly_set_str_pretty(C + 2, "1", NULL, ctxAC);
         if (fmpq_mpoly_compose_fmpq_mpoly(A, B, Cp, ctxB, ctxAC))
-        {
-            printf("FAIL\n");
-            flint_printf("Check non-example 1\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check non-example 1\n");
 
         fmpq_mpoly_set_str_pretty(C + 0, "x1", NULL, ctxAC);
         fmpq_mpoly_set_str_pretty(C + 1, "2*x2", NULL, ctxAC);
         fmpq_mpoly_set_str_pretty(C + 2, "1", NULL, ctxAC);
         if (fmpq_mpoly_compose_fmpq_mpoly(A, B, Cp, ctxB, ctxAC))
-        {
-            printf("FAIL\n");
-            flint_printf("Check non-example 2\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check non-example 2\n");
 
         fmpq_mpoly_set_str_pretty(C + 0, "2*x1", NULL, ctxAC);
         fmpq_mpoly_set_str_pretty(C + 1, "x2", NULL, ctxAC);
         fmpq_mpoly_set_str_pretty(C + 2, "1", NULL, ctxAC);
         if (!fmpq_mpoly_compose_fmpq_mpoly(A, B, Cp, ctxB, ctxAC))
-        {
-            printf("FAIL\n");
-            flint_printf("Check example 3\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check example 3\n");
+
+        /* Aliased generator composition */
+        c = (slong *) flint_malloc(nvarsB*sizeof(slong));
+        fmpq_mpoly_init(B1, ctxB);
+        fmpq_mpoly_set(B1, B, ctxB);
+        for (i = 0; i < nvarsB; i++)
+            c[i] = i;
+
+        fmpq_mpoly_compose_fmpq_mpoly_gen(B1, B1, c, ctxB, ctxB);
+        if (!fmpq_mpoly_equal(B, B1, ctxB))
+            TEST_FUNCTION_FAIL("Check composition with aliased generators\n");
+
+        /* Reverse the generators, twice */
+        for (i = 0; i < nvarsB; i++)
+            c[i] = nvarsB - i - 1;
+
+        fmpq_mpoly_compose_fmpq_mpoly_gen(B1, B1, c, ctxB, ctxB);
+        if (fmpq_mpoly_equal(B, B1, ctxB))
+            TEST_FUNCTION_FAIL("Check composition with reversed aliased generators\n");
+
+        fmpq_mpoly_compose_fmpq_mpoly_gen(B1, B1, c, ctxB, ctxB);
+        if (!fmpq_mpoly_equal(B, B1, ctxB))
+            TEST_FUNCTION_FAIL("Check composition with un-reversed aliased generators\n");
+
+        /* Composition with zero polys */
+        fmpq_mpoly_zero(B1, ctxB);
+
+        fmpq_mpoly_compose_fmpq_mpoly_gen(A, B1, c, ctxB, ctxAC);
+        if (!fmpq_mpoly_is_zero(B1, ctxB))
+            TEST_FUNCTION_FAIL("Check composition with generators of zero poly\n");
+
+        fmpq_mpoly_set_str_pretty(C + 0, "2*x1", NULL, ctxAC);
+        fmpq_mpoly_set_str_pretty(C + 1, "x2", NULL, ctxAC);
+        fmpq_mpoly_set_str_pretty(C + 2, "1", NULL, ctxAC);
+        if (!fmpq_mpoly_compose_fmpq_mpoly(A, B1, Cp, ctxB, ctxAC))
+            TEST_FUNCTION_FAIL("Check example 4\n");
+
+        if (!fmpq_mpoly_is_zero(A, ctxAC))
+            TEST_FUNCTION_FAIL("Check composition with zero poly\n");
 
         fmpq_mpoly_clear(B, ctxB);
         fmpq_mpoly_clear(A, ctxAC);
-        fmpq_mpoly_clear(A1, ctxAC);
-        fmpq_mpoly_clear(A2, ctxAC);
         for (i = 0; i < 3; i++)
             fmpq_mpoly_clear(C + i, ctxAC);
 
@@ -131,22 +155,12 @@ TEST_FUNCTION_START(fmpq_mpoly_compose_fmpq_mpoly, state)
             fmpq_mpoly_compose_fmpq_mpoly_gen(A, B, c, ctxB, ctxAC);
 
             if (!fmpq_mpoly_compose_fmpq_mpoly(A1, B, C, ctxB, ctxAC))
-            {
-                printf("FAIL\n");
-                flint_printf("Check composition success with generators\n"
-                                                     "i: %wd, j: %wd\n", i, j);
-                fflush(stdout);
-                flint_abort();
-            }
+                TEST_FUNCTION_FAIL("Check composition success with generators\n"
+                                   "i: %wd, j: %wd\n", i, j);
 
             if (!fmpq_mpoly_equal(A, A1, ctxAC))
-            {
-                printf("FAIL\n");
-                flint_printf("Check composition with generators\n"
-                                                     "i: %wd, j: %wd\n", i, j);
-                fflush(stdout);
-                flint_abort();
-            }
+                TEST_FUNCTION_FAIL("Check composition with generators\n"
+                                   "i: %wd, j: %wd\n", i, j);
 
             fmpq_mpoly_assert_canonical(A, ctxAC);
             fmpq_mpoly_assert_canonical(A1, ctxAC);
@@ -196,20 +210,10 @@ TEST_FUNCTION_START(fmpq_mpoly_compose_fmpq_mpoly, state)
         coeff_bits = n_randint(state, 100) + 1;
         fmpq_mpoly_randtest_bits(f, state, len1, coeff_bits, exp_bits, ctx);
         if (!fmpq_mpoly_compose_fmpq_mpoly(g, f, vals1, ctx, ctx))
-        {
-            printf("FAIL\n");
-            flint_printf("Check composition success\ni: %wd\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check composition success\ni: %wd\n", i);
 
         if (!fmpq_mpoly_equal(f, g, ctx))
-        {
-            printf("FAIL\n");
-            flint_printf("Check composition with identity\ni: %wd\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check composition with identity\ni: %wd\n", i);
 
         fmpq_mpoly_clear(g, ctx);
         fmpq_mpoly_clear(f, ctx);
@@ -276,41 +280,21 @@ TEST_FUNCTION_START(fmpq_mpoly_compose_fmpq_mpoly, state)
             vals3[v] = (fmpq *) flint_malloc(sizeof(fmpq));
             fmpq_init(vals3[v]);
             if (!fmpq_mpoly_evaluate_all_fmpq(vals3[v], vals1[v], vals2, ctx2))
-            {
-                printf("FAIL\n");
-                flint_printf("Check evaluation success\ni: %wd\n", i);
-                fflush(stdout);
-                flint_abort();
-            }
+                TEST_FUNCTION_FAIL("Check evaluation success\ni: %wd\n", i);
         }
 
         fmpq_mpoly_randtest_bound(f, state, len1, coeff_bits, exp_bound1, ctx1);
 
         if (!fmpq_mpoly_compose_fmpq_mpoly(g, f, vals1, ctx1, ctx2))
-        {
-            printf("FAIL\n");
-            flint_printf("Check composition success\ni: %wd\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check composition success\ni: %wd\n", i);
         fmpq_mpoly_assert_canonical(g, ctx2);
 
         if (!fmpq_mpoly_evaluate_all_fmpq(fe, f, vals3, ctx1) ||
             !fmpq_mpoly_evaluate_all_fmpq(ge, g, vals2, ctx2))
-        {
-            printf("FAIL\n");
-            flint_printf("Check evaluation success\ni: %wd\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check evaluation success\ni: %wd\n", i);
 
         if (!fmpq_equal(fe, ge))
-        {
-            printf("FAIL\n");
-            flint_printf("Check composition and evalall commute\ni: %wd\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check composition and evalall commute\ni: %wd\n", i);
 
         for (v = 0; v < nvars1; v++)
         {
diff --git a/src/fmpz_mod_mpoly.h b/src/fmpz_mod_mpoly.h
index 4c29b71afe..e103095110 100644
--- a/src/fmpz_mod_mpoly.h
+++ b/src/fmpz_mod_mpoly.h
@@ -570,10 +570,19 @@ int fmpz_mod_mpoly_compose_fmpz_mod_mpoly_geobucket(fmpz_mod_mpoly_t A,
             const fmpz_mod_mpoly_t B, fmpz_mod_mpoly_struct * const * C,
             const fmpz_mod_mpoly_ctx_t ctxB, const fmpz_mod_mpoly_ctx_t ctxAC);
 
+int fmpz_mod_mpoly_compose_fmpz_mod_mpoly_horner(fmpz_mod_mpoly_t A,
+            const fmpz_mod_mpoly_t B, fmpz_mod_mpoly_struct * const * C,
+            const fmpz_mod_mpoly_ctx_t ctxB, const fmpz_mod_mpoly_ctx_t ctxAC);
+
 int fmpz_mod_mpoly_compose_fmpz_mod_mpoly(fmpz_mod_mpoly_t A,
             const fmpz_mod_mpoly_t B, fmpz_mod_mpoly_struct * const * C,
             const fmpz_mod_mpoly_ctx_t ctxB, const fmpz_mod_mpoly_ctx_t ctxAC);
 
+void fmpz_mod_mpoly_compose_fmpz_mod_mpoly_gen(fmpz_mod_mpoly_t A,
+            const fmpz_mod_mpoly_t B, const slong * c,
+            const fmpz_mod_mpoly_ctx_t ctxB, const fmpz_mod_mpoly_ctx_t ctxAC);
+
+
 /* Multiplication ************************************************************/
 
 void fmpz_mod_mpoly_mul(fmpz_mod_mpoly_t A,
diff --git a/src/fmpz_mod_mpoly/compose_fmpz_mod_mpoly.c b/src/fmpz_mod_mpoly/compose_fmpz_mod_mpoly.c
index a354e5a693..984cb6e531 100644
--- a/src/fmpz_mod_mpoly/compose_fmpz_mod_mpoly.c
+++ b/src/fmpz_mod_mpoly/compose_fmpz_mod_mpoly.c
@@ -64,14 +64,14 @@ int fmpz_mod_mpoly_compose_fmpz_mod_mpoly(
 matrix_no_good:
 
     fmpz_mat_clear(M);
-/*
+
     for (i = 0; i < ctxB->minfo->nvars; i++)
     {
         if (C[i]->length > 1)
         {
-            return nmod_mpoly_compose_nmod_mpoly_horner(A, B, C, ctxB, ctxAC);
+            return fmpz_mod_mpoly_compose_fmpz_mod_mpoly_horner(A, B, C, ctxB, ctxAC);
         }
     }
-*/
+
     return fmpz_mod_mpoly_compose_fmpz_mod_mpoly_geobucket(A, B, C, ctxB, ctxAC);
 }
diff --git a/src/fmpz_mod_mpoly/compose_fmpz_mod_mpoly_gen.c b/src/fmpz_mod_mpoly/compose_fmpz_mod_mpoly_gen.c
new file mode 100644
index 0000000000..4d16f42566
--- /dev/null
+++ b/src/fmpz_mod_mpoly/compose_fmpz_mod_mpoly_gen.c
@@ -0,0 +1,48 @@
+/*
+    Copyright (C) 2020 Daniel Schultz
+
+    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 <https://www.gnu.org/licenses/>.
+*/
+
+#include "fmpz_mat.h"
+#include "mpoly.h"
+#include "fmpz_mod_mpoly.h"
+
+/* evaluate B(x_1,...,x_n) at x_i = y_c[i], y_j are vars of ctxAC */
+void fmpz_mod_mpoly_compose_fmpz_mod_mpoly_gen(fmpz_mod_mpoly_t A,
+            const fmpz_mod_mpoly_t B, const slong * c,
+            const fmpz_mod_mpoly_ctx_t ctxB, const fmpz_mod_mpoly_ctx_t ctxAC)
+{
+    fmpz_mat_t M;
+
+    if (B->length == 0)
+    {
+        fmpz_mod_mpoly_zero(A, ctxAC);
+        return;
+    }
+
+    fmpz_mat_init(M, ctxAC->minfo->nfields + 1, ctxB->minfo->nfields);
+    mpoly_compose_mat_gen(M, c, ctxB->minfo, ctxAC->minfo);
+
+    if (A == B)
+    {
+        fmpz_mod_mpoly_t T;
+        fmpz_mod_mpoly_init(T, ctxAC);
+        _fmpz_mod_mpoly_compose_mat(T, B, M, ctxB, ctxAC);
+        fmpz_mod_mpoly_swap(A, T, ctxAC);
+        fmpz_mod_mpoly_clear(T, ctxAC);
+    }
+    else
+    {
+        _fmpz_mod_mpoly_compose_mat(A, B, M, ctxB, ctxAC);
+    }
+
+    fmpz_mat_clear(M);
+
+    return;
+}
diff --git a/src/fmpz_mod_mpoly/compose_fmpz_mod_mpoly_horner.c b/src/fmpz_mod_mpoly/compose_fmpz_mod_mpoly_horner.c
new file mode 100644
index 0000000000..24d01119c3
--- /dev/null
+++ b/src/fmpz_mod_mpoly/compose_fmpz_mod_mpoly_horner.c
@@ -0,0 +1,439 @@
+/*
+    Copyright (C) 2018 Daniel Schultz
+
+    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 <https://www.gnu.org/licenses/>.
+*/
+
+#include "fmpz.h"
+#include "fmpz_vec.h"
+#include "mpoly.h"
+#include "fmpz_mod_mpoly.h"
+
+/*
+The conversion to Horner form can be stated as recursive. However, the call
+stack has depth proportial to the length of the input polynomial in the worst
+case. Therefore, we must convert it to an iterative algorithm.
+
+The procedure is
+
+HornerForm(f):
+
+    if f is simple to evaluate
+
+        return eval(f)
+
+    else
+        choose a variable v and the smallest non zero exponent e appearing
+            in the terms of f
+
+        write f = q * v^e + r  where r is independent of the variable v
+
+        return  HornerForm(q) * v^e + HornerForm(r)
+*/
+
+typedef struct
+{
+    slong f;
+    slong r;
+    slong v_var;
+    fmpz_t v_exp;   /* will be managed as stack grows / shrinks */
+    int ret;
+} stack_entry_struct;
+
+typedef stack_entry_struct stack_entry_t[1];
+
+
+/* A = A * X^pow */
+static int _fmpz_mod_mpoly_pmul(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t X,
+                  const fmpz_t pow, fmpz_mod_mpoly_t T, const fmpz_mod_mpoly_ctx_t ctx)
+{
+    ulong p;
+    FLINT_ASSERT(fmpz_sgn(pow) > 0);
+
+    if (!fmpz_fits_si(pow))
+    {
+        if (!fmpz_mod_mpoly_pow_fmpz(T, X, pow, ctx))
+        {
+            fmpz_mod_mpoly_zero(A, ctx);
+            return 0;
+        }
+
+        fmpz_mod_mpoly_mul(A, A, T, ctx);
+        return 1;
+    }
+
+    p = fmpz_get_ui(pow);
+
+    if (X->length <= WORD(2) || (slong) (A->length/p) < X->length)
+    {
+        if (!fmpz_mod_mpoly_pow_ui(T, X, p, ctx))
+        {
+            fmpz_mod_mpoly_zero(A, ctx);
+            return 0;
+        }
+
+        fmpz_mod_mpoly_mul(A, A, T, ctx);
+    }
+    else
+    {
+        while (p >= 1)
+        {
+            fmpz_mod_mpoly_mul(T, A, X, ctx);
+            fmpz_mod_mpoly_swap(A, T, ctx);
+            p--;
+        }
+    }
+
+    return 1;
+}
+
+/* evaluate B(xbar) at xbar = C */
+int fmpz_mod_mpoly_compose_fmpz_mod_mpoly_horner(fmpz_mod_mpoly_t A,
+                  const fmpz_mod_mpoly_t B, fmpz_mod_mpoly_struct * const * C,
+                     const fmpz_mod_mpoly_ctx_t ctxB, const fmpz_mod_mpoly_ctx_t ctxAC)
+{
+    int success = 1;
+    int ret;
+    slong nvars = ctxB->minfo->nvars;
+    slong i, j, k, cur, next, f, r, f_prev, r_prev, v;
+    slong sp, rp;
+    stack_entry_struct * stack;
+    fmpz_mod_mpoly_struct * regs;
+    fmpz_mod_mpoly_t temp;
+    slong * rtypes;
+    ulong totalcounts, maxcounts;
+    ulong * counts;
+    slong Blen = B->length;
+    slong * Blist;
+    const fmpz * Bcoeff = B->coeffs;
+    const ulong * Bexp = B->exps;
+    flint_bitcnt_t Bbits = B->bits;
+    slong BN = mpoly_words_per_exp(Bbits, ctxB->minfo);
+    fmpz * Buexp;
+    fmpz * mdegs;
+    fmpz_t score, tz;
+    TMP_INIT;
+
+    if (Blen < 1)
+    {
+        fmpz_mod_mpoly_zero(A, ctxAC);
+        return 1;
+    }
+
+    if (nvars < 1)
+    {
+       FLINT_ASSERT(Blen == 1);
+       fmpz_mod_mpoly_set_fmpz(A, B->coeffs + 0, ctxAC);
+       return 1;
+    }
+
+    FLINT_ASSERT(A != B);
+    FLINT_ASSERT(Blen > 0);
+
+    TMP_START;
+
+    fmpz_init(score);
+    fmpz_init(tz);
+
+    /* unpack B exponents */
+    Buexp = _fmpz_vec_init(nvars*Blen);
+    for (i = 0; i < Blen; i++)
+        mpoly_get_monomial_ffmpz(Buexp + nvars*i, Bexp + BN*i, Bbits, ctxB->minfo);
+
+    counts = (ulong *) TMP_ALLOC(nvars*sizeof(ulong));
+    mdegs = _fmpz_vec_init(nvars);
+
+    /* stack */
+    sp = -WORD(1); /* start with empty stack */
+    stack = (stack_entry_struct *) TMP_ALLOC(nvars*(Blen + 1)*sizeof(stack_entry_struct));
+    Blist = (slong *) TMP_ALLOC(Blen*sizeof(slong));
+
+    /* registers of polynomials */
+    rp = 0;
+    rtypes = (slong *) TMP_ALLOC((nvars + 1)*sizeof(slong));
+    regs   = (fmpz_mod_mpoly_struct *) TMP_ALLOC(nvars*sizeof(fmpz_mod_mpoly_struct));
+    for (i = 0; i < nvars; i++)
+        fmpz_mod_mpoly_init(regs + i, ctxAC);
+    fmpz_mod_mpoly_init(temp, ctxAC);
+
+    /* polynomials will be stored as link lists */
+    for (i = 0; i + 1 < Blen; i++)
+        Blist[i] = i + 1;
+    Blist[i] = -WORD(1);
+
+    sp++;
+    fmpz_init((stack + sp)->v_exp);
+    (stack + sp)->ret = 0;
+    (stack + sp)->f = 0;
+
+HornerForm:
+
+    f = (stack + sp)->f;
+
+    FLINT_ASSERT(f != -WORD(1)); /* f is not supposed to be zero */
+
+    /* obtain a count of the number of terms containing each variable */
+    for (i = 0; i < nvars; i++)
+    {
+        counts[i] = 0;
+        fmpz_set_si(mdegs + i, -WORD(1));
+    }
+
+    for (j = f; j != -WORD(1); j = Blist[j])
+    {
+        for (i = 0; i < nvars; i++)
+        {
+            if (!fmpz_is_zero(Buexp + nvars*j + i ))
+            {
+                counts[i]++;
+                if (fmpz_sgn(mdegs + i) < 0
+                    || fmpz_cmp(mdegs + i, Buexp + nvars*j + i) > 0)
+                {
+                    fmpz_set(mdegs + i, Buexp + nvars*j + i);
+                }
+            }
+        }
+    }
+
+    totalcounts = 0;
+    maxcounts = 0;
+    v = -WORD(1);
+    for (i = 0; i < nvars; i++)
+    {
+        maxcounts = FLINT_MAX(maxcounts, counts[i]);
+        totalcounts += counts[i];
+        if (counts[i] != 0)
+            v = i;
+    }
+
+    /* handle simple cases */
+    if (totalcounts == 0)
+    {
+        FLINT_ASSERT(Blist[f] == -WORD(1));    /* f should have had only one term */
+        rtypes[rp] = f;
+        goto HornerFormReturn;
+    }
+    else if (totalcounts == 1)
+    {
+        FLINT_ASSERT(!fmpz_is_zero(Buexp + nvars*f + v)); /* this term should not be a scalar */
+        if (!fmpz_mod_mpoly_pow_fmpz(regs + rp, C[v], Buexp + nvars*f + v, ctxAC))
+        {
+            success = 0;
+        }
+        fmpz_mod_mpoly_scalar_mul_fmpz(regs + rp, regs + rp, Bcoeff + f, ctxAC);
+
+        if (Blist[f] != -WORD(1)) /* if f has a second term */
+        {
+            /* this term should be a scalar */
+            FLINT_ASSERT(fmpz_is_zero(Buexp + nvars*Blist[f] + v));
+            fmpz_mod_mpoly_add_fmpz(regs + rp, regs + rp,  Bcoeff + Blist[f], ctxAC);
+        }
+
+        rtypes[rp] = -WORD(1);
+
+        goto HornerFormReturn;
+    }
+
+    /* pick best power to pull out */
+    k = 0;
+    if (maxcounts == 1)
+    {
+        fmpz_set_si(score, -WORD(1));
+        for (i = 0; i < nvars; i++)
+        {
+            if (counts[i] == 1 && (fmpz_sgn(score) < 0
+                                   || fmpz_cmp(mdegs + i, score) < 0))
+            {
+                FLINT_ASSERT(fmpz_sgn(mdegs + i) > 0);
+                fmpz_set(score, mdegs + i);
+                k = i;
+            }
+        }
+    }
+    else
+    {
+        fmpz_zero(score);
+        for (i = 0; i < nvars; i++)
+        {
+            if (counts[i] > 1)
+            {
+                FLINT_ASSERT(fmpz_sgn(mdegs + i) > 0);
+                fmpz_mul_ui(tz, mdegs + i, counts[i] - 1);
+                if (fmpz_cmp(tz, score) > 0)
+                {
+                    fmpz_swap(score, tz);
+                    k = i;
+                }
+            }
+        }
+    }
+
+    /* set variable power v */
+    (stack + sp)->v_var = k;
+    fmpz_set((stack + sp)->v_exp, mdegs + k);
+
+    /* scan f and split into q and v with f = q*v + r then set f = q */
+    r = -WORD(1);
+    cur = f;
+    f_prev = -WORD(1);
+    r_prev = -WORD(1);
+    while (cur != -WORD(1))
+    {
+        next = Blist[cur];
+        if (fmpz_is_zero(Buexp + nvars*cur + k))
+        {
+            if (f_prev == -WORD(1))
+                f = Blist[cur];
+            else
+                Blist[f_prev] = Blist[cur];
+
+            if (r_prev == -WORD(1))
+                r = cur;
+            else
+                Blist[r_prev] = cur;
+
+            Blist[cur] = -WORD(1);
+            r_prev = cur;
+        }
+        else
+        {
+            /* mdegs[k] should be minimum non zero exponent */
+            fmpz_sub(Buexp + nvars*cur + k, Buexp + nvars*cur + k, mdegs + k);
+            FLINT_ASSERT(fmpz_sgn(Buexp + nvars*cur + k) >= 0);
+            f_prev = cur;
+        }
+        cur = next;
+    }
+    (stack + sp)->r = r;
+
+    /* convert the quotient */
+    sp++;
+    fmpz_init((stack + sp)->v_exp);
+    (stack + sp)->ret = 1;
+    (stack + sp)->f = f;
+    goto HornerForm;
+
+HornerForm1:
+
+    /* convert the remainder */
+    r = (stack + sp)->r;
+    if (r != -WORD(1))
+    {
+        /* remainder is non zero */
+        rp++;
+        FLINT_ASSERT(0 <= rp && rp <= nvars);
+        sp++;
+        fmpz_init((stack + sp)->v_exp);
+        (stack + sp)->ret = 2;
+        (stack + sp)->f = r;
+        goto HornerForm;
+
+HornerForm2:
+
+        if (rtypes[rp - 1] == -WORD(1) && rtypes[rp] == -WORD(1))
+        {
+            /* both quotient and remainder are polynomials */
+            if (!_fmpz_mod_mpoly_pmul(regs + rp - 1, C[(stack + sp)->v_var],
+                                             (stack + sp)->v_exp, temp, ctxAC))
+            {
+                success = 0;
+            }
+            fmpz_mod_mpoly_add(temp, regs + rp - 1, regs + rp, ctxAC);
+            fmpz_mod_mpoly_swap(temp, regs + rp - 1, ctxAC);
+        }
+        else if (rtypes[rp - 1] == -WORD(1) && rtypes[rp] != -WORD(1))
+        {
+            /* quotient is a polynomial, remainder is a scalar */
+            if (!_fmpz_mod_mpoly_pmul(regs + rp - 1, C[(stack + sp)->v_var],
+                                             (stack + sp)->v_exp, temp, ctxAC))
+            {
+                success = 0;
+            }
+            fmpz_mod_mpoly_add_fmpz(regs + rp - 1, regs + rp - 1,
+                                                   Bcoeff + rtypes[rp], ctxAC);
+        }
+        else if (rtypes[rp - 1] != -WORD(1) && rtypes[rp] == -WORD(1))
+        {
+            /* quotient is a scalar, remainder is a polynomial */
+            if (!fmpz_mod_mpoly_pow_fmpz(temp, C[(stack + sp)->v_var],
+                                                   (stack + sp)->v_exp, ctxAC))
+            {
+                success = 0;
+            }
+            fmpz_mod_mpoly_scalar_mul_fmpz(temp, temp, Bcoeff + rtypes[rp - 1], ctxAC);
+            fmpz_mod_mpoly_add(regs + rp - 1, temp, regs + rp, ctxAC);
+        }
+        else
+        {
+            /* quotient is a scalar, remainder is a scalar */
+            FLINT_ASSERT(0);    /* this should have been handled by simple case */
+        }
+        rp--;
+        FLINT_ASSERT(0 <= rp && rp <= nvars);
+    }
+    else
+    {
+        /* remainder is zero */
+        FLINT_ASSERT(rtypes[rp] == -WORD(1)); /* quotient is not a scalar */
+
+        /* quotient is a polynomial */
+        if (!_fmpz_mod_mpoly_pmul(regs + rp, C[(stack + sp)->v_var],
+                                            (stack + sp)->v_exp, temp, ctxAC))
+        {
+            success = 0;
+        }
+    }
+
+    rtypes[rp] = -WORD(1);
+
+HornerFormReturn:
+
+    if (!success)
+    {
+        while (sp >= 0)
+        {
+            fmpz_clear((stack + sp)->v_exp);
+            sp--;
+        }
+
+        goto cleanup;
+    }
+
+    ret = (stack + sp)->ret;
+    fmpz_clear((stack + sp)->v_exp);
+    sp--;
+    if (ret == 1) goto HornerForm1;
+    if (ret == 2) goto HornerForm2;
+
+    FLINT_ASSERT(rp == 0);
+    FLINT_ASSERT(sp == -WORD(1));
+
+    if (rtypes[rp] == -WORD(1))
+    {
+        fmpz_mod_mpoly_swap(A, regs + rp, ctxAC);
+    }
+    else
+    {
+        fmpz_mod_mpoly_set_fmpz(A, Bcoeff + rtypes[rp], ctxAC);
+    }
+
+cleanup:
+
+    for (i = 0; i < nvars; i++)
+        fmpz_mod_mpoly_clear(regs + i, ctxAC);
+    fmpz_mod_mpoly_clear(temp, ctxAC);
+
+    fmpz_clear(score);
+    fmpz_clear(tz);
+    _fmpz_vec_clear(mdegs, nvars);
+    _fmpz_vec_clear(Buexp, nvars*Blen);
+
+    TMP_END;
+
+    return success;
+}
diff --git a/src/fmpz_mod_mpoly/test/main.c b/src/fmpz_mod_mpoly/test/main.c
index c29f6419ec..7331db8fd3 100644
--- a/src/fmpz_mod_mpoly/test/main.c
+++ b/src/fmpz_mod_mpoly/test/main.c
@@ -15,6 +15,7 @@
 #include "t-add_sub_fmpz.c"
 #include "t-add_sub_si.c"
 #include "t-cmp.c"
+#include "t-compose_fmpz_mod_mpoly.c"
 #include "t-degree.c"
 #include "t-degrees_term_exp_fits_ui_si.c"
 #include "t-derivative.c"
@@ -65,6 +66,7 @@ test_struct tests[] =
     TEST_FUNCTION(fmpz_mod_mpoly_add_sub),
     TEST_FUNCTION(fmpz_mod_mpoly_add_sub_fmpz),
     TEST_FUNCTION(fmpz_mod_mpoly_add_sub_si),
+    TEST_FUNCTION(fmpz_mod_mpoly_compose),
     TEST_FUNCTION(fmpz_mod_mpoly_cmp),
     TEST_FUNCTION(fmpz_mod_mpoly_degree),
     TEST_FUNCTION(fmpz_mod_mpoly_degrees_term_exp_fits_ui_si),
diff --git a/src/fmpz_mod_mpoly/test/t-compose_fmpz_mod_mpoly.c b/src/fmpz_mod_mpoly/test/t-compose_fmpz_mod_mpoly.c
new file mode 100644
index 0000000000..ce8414f643
--- /dev/null
+++ b/src/fmpz_mod_mpoly/test/t-compose_fmpz_mod_mpoly.c
@@ -0,0 +1,385 @@
+/*
+    Copyright (C) 2018 Daniel Schultz
+
+    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 <https://www.gnu.org/licenses/>.
+*/
+
+#include "test_helpers.h"
+#include "fmpz_mod_mpoly.h"
+
+TEST_FUNCTION_START(fmpz_mod_mpoly_compose, state)
+{
+    slong i, j, v;
+
+    {
+        fmpz_t p;
+        slong nvarsAC = 2, nvarsB = 3;
+        fmpz_mod_mpoly_t A, A1, A2, B, B1;
+        fmpz_mod_mpoly_struct * Cp[3];
+        fmpz_mod_mpoly_struct C[3];
+        fmpz_mod_mpoly_ctx_t ctxAC, ctxB;
+        slong * c;
+
+        fmpz_init(p);
+        fmpz_randprime(p, state, n_randint(state, 15) + 2, 0);
+        fmpz_mod_mpoly_ctx_init(ctxB, nvarsB, ORD_LEX, p);
+        fmpz_mod_mpoly_ctx_init(ctxAC, nvarsAC, ORD_LEX, p);
+
+        fmpz_mod_mpoly_init(B, ctxB);
+        fmpz_mod_mpoly_init(A, ctxAC);
+        fmpz_mod_mpoly_init(A1, ctxAC);
+        fmpz_mod_mpoly_init(A2, ctxAC);
+        for (i = 0; i < 3; i++)
+        {
+            Cp[i] = C + i;
+            fmpz_mod_mpoly_init(C + i, ctxAC);
+        }
+
+        fmpz_mod_mpoly_set_str_pretty(B,
+                "1 + x1*x2^2 + x2^9999999999999999999999999*x3^9", NULL, ctxB);
+
+        fmpz_mod_mpoly_set_str_pretty(C + 0, "x1 + x2", NULL, ctxAC);
+        fmpz_mod_mpoly_set_str_pretty(C + 1, "x1 - x2", NULL, ctxAC);
+        fmpz_mod_mpoly_set_str_pretty(C + 2, "1", NULL, ctxAC);
+        if (fmpz_mod_mpoly_compose_fmpz_mod_mpoly(A, B, Cp, ctxB, ctxAC) ||
+            fmpz_mod_mpoly_compose_fmpz_mod_mpoly_horner(A1, B, Cp, ctxB, ctxAC) ||
+            fmpz_mod_mpoly_compose_fmpz_mod_mpoly_geobucket(A2, B, Cp, ctxB, ctxAC))
+            TEST_FUNCTION_FAIL("Check non-example 1\n");
+
+        fmpz_mod_mpoly_set_str_pretty(C + 0, "x1", NULL, ctxAC);
+        fmpz_mod_mpoly_set_str_pretty(C + 1, "2*x2", NULL, ctxAC);
+        fmpz_mod_mpoly_set_str_pretty(C + 2, "1", NULL, ctxAC);
+        if (!fmpz_mod_mpoly_compose_fmpz_mod_mpoly(A, B, Cp, ctxB, ctxAC) ||
+            !fmpz_mod_mpoly_compose_fmpz_mod_mpoly_horner(A1, B, Cp, ctxB, ctxAC) ||
+            !fmpz_mod_mpoly_compose_fmpz_mod_mpoly_geobucket(A2, B, Cp, ctxB, ctxAC))
+            TEST_FUNCTION_FAIL("Check non-example 2\n");
+
+        fmpz_mod_mpoly_set_str_pretty(C + 0, "2*x1", NULL, ctxAC);
+        fmpz_mod_mpoly_set_str_pretty(C + 1, "x2", NULL, ctxAC);
+        fmpz_mod_mpoly_set_str_pretty(C + 2, "1", NULL, ctxAC);
+        if (!fmpz_mod_mpoly_compose_fmpz_mod_mpoly(A, B, Cp, ctxB, ctxAC) ||
+            !fmpz_mod_mpoly_compose_fmpz_mod_mpoly_horner(A1, B, Cp, ctxB, ctxAC) ||
+            !fmpz_mod_mpoly_compose_fmpz_mod_mpoly_geobucket(A2, B, Cp, ctxB, ctxAC))
+            TEST_FUNCTION_FAIL("Check example 3\n");
+
+        /* Aliased generator composition */
+        c = (slong *) flint_malloc(nvarsB*sizeof(slong));
+        fmpz_mod_mpoly_init(B1, ctxB);
+        fmpz_mod_mpoly_set(B1, B, ctxB);
+        for (i = 0; i < nvarsB; i++)
+            c[i] = i;
+
+        fmpz_mod_mpoly_compose_fmpz_mod_mpoly_gen(B1, B1, c, ctxB, ctxB);
+        if (!fmpz_mod_mpoly_equal(B, B1, ctxB))
+            TEST_FUNCTION_FAIL("Check composition with aliased generators\n");
+
+        /* Reverse the generators, twice */
+        for (i = 0; i < nvarsB; i++)
+            c[i] = nvarsB - i - 1;
+
+        fmpz_mod_mpoly_compose_fmpz_mod_mpoly_gen(B1, B1, c, ctxB, ctxB);
+        if (fmpz_mod_mpoly_equal(B, B1, ctxB))
+            TEST_FUNCTION_FAIL("Check composition with reversed aliased generators\n");
+
+        fmpz_mod_mpoly_compose_fmpz_mod_mpoly_gen(B1, B1, c, ctxB, ctxB);
+        if (!fmpz_mod_mpoly_equal(B, B1, ctxB))
+            TEST_FUNCTION_FAIL("Check composition with un-reversed aliased generators\n");
+
+        /* Composition with zero polys */
+        fmpz_mod_mpoly_zero(B1, ctxB);
+
+        fmpz_mod_mpoly_compose_fmpz_mod_mpoly_gen(A, B1, c, ctxB, ctxAC);
+        if (!fmpz_mod_mpoly_is_zero(B1, ctxB))
+            TEST_FUNCTION_FAIL("Check composition with generators of zero poly\n");
+
+        fmpz_mod_mpoly_set_str_pretty(C + 0, "2*x1", NULL, ctxAC);
+        fmpz_mod_mpoly_set_str_pretty(C + 1, "x2", NULL, ctxAC);
+        fmpz_mod_mpoly_set_str_pretty(C + 2, "1", NULL, ctxAC);
+        if (!fmpz_mod_mpoly_compose_fmpz_mod_mpoly(A, B1, Cp, ctxB, ctxAC) ||
+            !fmpz_mod_mpoly_compose_fmpz_mod_mpoly_horner(A1, B1, Cp, ctxB, ctxAC) ||
+            !fmpz_mod_mpoly_compose_fmpz_mod_mpoly_geobucket(A2, B1, Cp, ctxB, ctxAC))
+            TEST_FUNCTION_FAIL("Check example 4\n");
+
+        if (!fmpz_mod_mpoly_is_zero(A, ctxAC) ||
+            !fmpz_mod_mpoly_is_zero(A1, ctxAC) ||
+            !fmpz_mod_mpoly_is_zero(A2, ctxAC))
+            TEST_FUNCTION_FAIL("Check composition with zero poly\n");
+
+        fmpz_mod_mpoly_clear(B, ctxB);
+        fmpz_mod_mpoly_clear(A, ctxAC);
+        fmpz_mod_mpoly_clear(A1, ctxAC);
+        fmpz_mod_mpoly_clear(A2, ctxAC);
+        for (i = 0; i < 3; i++)
+            fmpz_mod_mpoly_clear(C + i, ctxAC);
+
+        fmpz_mod_mpoly_ctx_clear(ctxB);
+        fmpz_mod_mpoly_ctx_clear(ctxAC);
+    }
+
+    /* check composition with generators */
+    for (i = 0; i < 20*flint_test_multiplier(); i++)
+    {
+        fmpz_t p;
+        fmpz_mod_mpoly_ctx_t ctxB, ctxAC;
+        fmpz_mod_mpoly_t B, A, A1, A2, A3;
+        fmpz_mod_mpoly_struct ** C;
+        slong * c;
+        slong nvarsB, nvarsAC;
+        slong len;
+        flint_bitcnt_t exp_bits;
+
+        fmpz_init(p);
+        fmpz_randprime(p, state, n_randint(state, 15) + 2, 0);
+        fmpz_mod_mpoly_ctx_init_rand(ctxB, state, 20, p);
+        fmpz_mod_mpoly_ctx_init_rand(ctxAC, state, 20, p);
+
+        nvarsB = ctxB->minfo->nvars;
+        nvarsAC = ctxAC->minfo->nvars;
+
+        fmpz_mod_mpoly_init(B, ctxB);
+        fmpz_mod_mpoly_init(A, ctxAC);
+        fmpz_mod_mpoly_init(A1, ctxAC);
+        fmpz_mod_mpoly_init(A2, ctxAC);
+        fmpz_mod_mpoly_init(A3, ctxAC);
+
+        c = (slong *) flint_malloc(nvarsB*sizeof(slong));
+        C = (fmpz_mod_mpoly_struct **) flint_malloc(nvarsB * sizeof(fmpz_mod_mpoly_struct *));
+        for (v = 0; v < nvarsB; v++)
+        {
+            C[v] = (fmpz_mod_mpoly_struct *) flint_malloc(sizeof(fmpz_mod_mpoly_struct));
+            fmpz_mod_mpoly_init(C[v], ctxAC);
+        }
+
+        for (j = 0; j < 4; j++)
+        {
+            len = n_randint(state, 200);
+            exp_bits = n_randint(state, 100) + 1;
+
+            for (v = 0; v < nvarsB; v++)
+            {
+                c[v] = n_randint(state, nvarsAC + 2) - 2;
+                if (c[v] >= 0)
+                    fmpz_mod_mpoly_gen(C[v], c[v], ctxAC);
+                else
+                    fmpz_mod_mpoly_zero(C[v], ctxAC);
+            }
+
+            fmpz_mod_mpoly_randtest_bits(B, state, len, exp_bits, ctxB);
+
+            fmpz_mod_mpoly_compose_fmpz_mod_mpoly_gen(A, B, c, ctxB, ctxAC);
+
+            if (!fmpz_mod_mpoly_compose_fmpz_mod_mpoly(A1, B, C, ctxB, ctxAC) ||
+                !fmpz_mod_mpoly_compose_fmpz_mod_mpoly_horner(A2, B, C, ctxB, ctxAC) ||
+                !fmpz_mod_mpoly_compose_fmpz_mod_mpoly_geobucket(A3, B, C, ctxB, ctxAC))
+                TEST_FUNCTION_FAIL("Check composition success with generators\n"
+                                   "i: %wd, j: %wd\n", i, j);
+
+            if (!fmpz_mod_mpoly_equal(A, A1, ctxAC) ||
+                !fmpz_mod_mpoly_equal(A, A2, ctxAC) ||
+                !fmpz_mod_mpoly_equal(A, A3, ctxAC))
+                TEST_FUNCTION_FAIL("Check composition with generators\n"
+                                   "i: %wd, j: %wd\n", i, j);
+
+            fmpz_mod_mpoly_assert_canonical(A, ctxAC);
+            fmpz_mod_mpoly_assert_canonical(A1, ctxAC);
+            fmpz_mod_mpoly_assert_canonical(A2, ctxAC);
+            fmpz_mod_mpoly_assert_canonical(A3, ctxAC);
+        }
+
+        for (v = 0; v < nvarsB; v++)
+        {
+            fmpz_mod_mpoly_clear(C[v], ctxAC);
+            flint_free(C[v]);
+        }
+        flint_free(C);
+        flint_free(c);
+
+        fmpz_mod_mpoly_clear(B, ctxB);
+        fmpz_mod_mpoly_clear(A, ctxAC);
+        fmpz_mod_mpoly_clear(A1, ctxAC);
+        fmpz_mod_mpoly_clear(A2, ctxAC);
+        fmpz_mod_mpoly_clear(A3, ctxAC);
+
+        fmpz_mod_mpoly_ctx_clear(ctxB);
+        fmpz_mod_mpoly_ctx_clear(ctxAC);
+    }
+
+    /* Check composition with identity */
+    for (i = 0; i < 20*flint_test_multiplier(); i++)
+    {
+        fmpz_t p;
+        slong len1;
+        flint_bitcnt_t exp_bits;
+        fmpz_mod_mpoly_struct ** vals1;
+        fmpz_mod_mpoly_t f, g, g1, g2;
+        fmpz_mod_mpoly_ctx_t ctx;
+
+        fmpz_init(p);
+        fmpz_randprime(p, state, n_randint(state, 15) + 2, 0);
+        fmpz_mod_mpoly_ctx_init_rand(ctx, state, 10, p);
+
+        vals1 = (fmpz_mod_mpoly_struct **) flint_malloc(ctx->minfo->nvars*
+                                                  sizeof(fmpz_mod_mpoly_struct *));
+        for (v = 0; v < ctx->minfo->nvars; v++)
+        {
+            vals1[v] = (fmpz_mod_mpoly_struct *) flint_malloc(sizeof(fmpz_mod_mpoly_struct));
+            fmpz_mod_mpoly_init(vals1[v], ctx);
+            fmpz_mod_mpoly_gen(vals1[v], v, ctx);
+        }
+
+        fmpz_mod_mpoly_init(f, ctx);
+        fmpz_mod_mpoly_init(g, ctx);
+        fmpz_mod_mpoly_init(g1, ctx);
+        fmpz_mod_mpoly_init(g2, ctx);
+
+        len1 = n_randint(state, 200);
+        exp_bits = n_randint(state, 100) + 1;
+        fmpz_mod_mpoly_randtest_bits(f, state, len1, exp_bits, ctx);
+
+        if (!fmpz_mod_mpoly_compose_fmpz_mod_mpoly(g, f, vals1, ctx, ctx) ||
+            !fmpz_mod_mpoly_compose_fmpz_mod_mpoly_horner(g1, f, vals1, ctx, ctx) ||
+            !fmpz_mod_mpoly_compose_fmpz_mod_mpoly_geobucket(g2, f, vals1, ctx, ctx) ||
+            !fmpz_mod_mpoly_equal(g, g1, ctx) ||
+            !fmpz_mod_mpoly_equal(g, g2, ctx))
+            TEST_FUNCTION_FAIL("Check composition success\ni: %wd\n", i);
+
+        fmpz_mod_mpoly_assert_canonical(g, ctx);
+        fmpz_mod_mpoly_assert_canonical(g1, ctx);
+        fmpz_mod_mpoly_assert_canonical(g2, ctx);
+
+        if (!fmpz_mod_mpoly_equal(f, g, ctx))
+            TEST_FUNCTION_FAIL("Check composition with identity\ni: %wd\n", i);
+
+        fmpz_mod_mpoly_clear(f, ctx);
+        fmpz_mod_mpoly_clear(g, ctx);
+        fmpz_mod_mpoly_clear(g1, ctx);
+        fmpz_mod_mpoly_clear(g2, ctx);
+
+        for (v = 0; v < ctx->minfo->nvars; v++)
+        {
+            fmpz_mod_mpoly_clear(vals1[v], ctx);
+            flint_free(vals1[v]);
+        }
+        flint_free(vals1);
+
+        fmpz_mod_mpoly_ctx_clear(ctx);
+    }
+
+    /* Check composition and evalall commute */
+    for (i = 0; i < 50 * flint_test_multiplier(); i++)
+    {
+        fmpz_t p;
+        fmpz_mod_mpoly_ctx_t ctx1, ctx2;
+        fmpz_mod_mpoly_t f, g, g1, g2;
+        fmpz_mod_mpoly_struct ** vals1;
+        fmpz_t fe, ge;
+        fmpz ** vals2, ** vals3;
+        slong nvars1, nvars2;
+        slong len1, len2;
+        slong exp_bound1, exp_bound2;
+
+        fmpz_init(p);
+        fmpz_randprime(p, state, n_randint(state, 15) + 2, 0);
+        fmpz_mod_mpoly_ctx_init_rand(ctx1, state, 6, p);
+        fmpz_mod_mpoly_ctx_init_rand(ctx2, state, 6, p);
+
+        nvars1 = ctx1->minfo->nvars;
+        nvars2 = ctx2->minfo->nvars;
+
+        fmpz_mod_mpoly_init(f, ctx1);
+        fmpz_mod_mpoly_init(g, ctx2);
+        fmpz_mod_mpoly_init(g1, ctx2);
+        fmpz_mod_mpoly_init(g2, ctx2);
+        fmpz_init(fe);
+        fmpz_init(ge);
+
+        len1 = n_randint(state, 50/FLINT_MAX(WORD(1), nvars1) + 1) + 1;
+        len2 = n_randint(state, 10/FLINT_MAX(WORD(1), nvars2) + 1) + 1;
+        exp_bound1 = n_randint(state, 15/FLINT_MAX(WORD(1), nvars1) + 2) + 2;
+        exp_bound2 = n_randint(state, 15/FLINT_MAX(WORD(1), nvars2) + 2) + 3;
+
+        vals1 = (fmpz_mod_mpoly_struct **) flint_malloc(nvars1
+                                                * sizeof(fmpz_mod_mpoly_struct *));
+        for (v = 0; v < nvars1; v++)
+        {
+            vals1[v] = (fmpz_mod_mpoly_struct *) flint_malloc(
+                                                    sizeof(fmpz_mod_mpoly_struct));
+            fmpz_mod_mpoly_init(vals1[v], ctx2);
+            fmpz_mod_mpoly_randtest_bound(vals1[v], state, len2, exp_bound2, ctx2);
+        }
+
+        vals2 = (fmpz **) flint_malloc(nvars2*sizeof(fmpz*));
+        for (v = 0; v < nvars2; v++)
+        {
+            vals2[v] = (fmpz *) flint_malloc(sizeof(fmpz));
+            fmpz_init(vals2[v]);
+            fmpz_randbits(vals2[v], state, 4);
+        }
+
+        vals3 = (fmpz **) flint_malloc(nvars1*sizeof(fmpz*));
+        for (v = 0; v < nvars1; v++)
+        {
+            vals3[v] = (fmpz *) flint_malloc(sizeof(fmpz));
+            fmpz_init(vals3[v]);
+            fmpz_mod_mpoly_evaluate_all_fmpz(vals3[v], vals1[v], vals2, ctx2);
+        }
+
+        fmpz_mod_mpoly_randtest_bound(f, state, len1, exp_bound1, ctx1);
+
+        if (!fmpz_mod_mpoly_compose_fmpz_mod_mpoly(g, f, vals1, ctx1, ctx2) ||
+            !fmpz_mod_mpoly_compose_fmpz_mod_mpoly_horner(g1, f, vals1, ctx1, ctx2) ||
+            !fmpz_mod_mpoly_compose_fmpz_mod_mpoly_geobucket(g2, f, vals1, ctx1, ctx2) ||
+            !fmpz_mod_mpoly_equal(g, g1, ctx2) ||
+            !fmpz_mod_mpoly_equal(g, g2, ctx2))
+            TEST_FUNCTION_FAIL("Check composition success\ni: %wd\n", i);
+
+        fmpz_mod_mpoly_assert_canonical(g, ctx2);
+        fmpz_mod_mpoly_assert_canonical(g1, ctx2);
+        fmpz_mod_mpoly_assert_canonical(g2, ctx2);
+
+        fmpz_mod_mpoly_evaluate_all_fmpz(fe, f, vals3, ctx1);
+        fmpz_mod_mpoly_evaluate_all_fmpz(ge, g, vals2, ctx2);
+
+        if (!fmpz_equal(fe, ge))
+            TEST_FUNCTION_FAIL("Check composition and evalall commute\ni: %wd\n", i);
+
+        for (v = 0; v < nvars1; v++)
+        {
+            fmpz_mod_mpoly_clear(vals1[v], ctx2);
+            flint_free(vals1[v]);
+        }
+        flint_free(vals1);
+
+        for (v = 0; v < nvars2; v++)
+        {
+            fmpz_clear(vals2[v]);
+            flint_free(vals2[v]);
+        }
+        flint_free(vals2);
+
+        for (v = 0; v < nvars1; v++)
+        {
+            fmpz_clear(vals3[v]);
+            flint_free(vals3[v]);
+        }
+        flint_free(vals3);
+
+        fmpz_mod_mpoly_clear(f, ctx1);
+        fmpz_mod_mpoly_clear(g, ctx2);
+        fmpz_mod_mpoly_clear(g1, ctx2);
+        fmpz_mod_mpoly_clear(g2, ctx2);
+
+        fmpz_clear(fe);
+        fmpz_clear(ge);
+
+        fmpz_mod_mpoly_ctx_clear(ctx1);
+        fmpz_mod_mpoly_ctx_clear(ctx2);
+    }
+
+    TEST_FUNCTION_END(state);
+}
diff --git a/src/fmpz_mpoly/test/t-compose_fmpz_mpoly.c b/src/fmpz_mpoly/test/t-compose_fmpz_mpoly.c
index af9d6978dc..957f0b18a1 100644
--- a/src/fmpz_mpoly/test/t-compose_fmpz_mpoly.c
+++ b/src/fmpz_mpoly/test/t-compose_fmpz_mpoly.c
@@ -17,13 +17,15 @@ TEST_FUNCTION_START(fmpz_mpoly_compose_fmpz_mpoly, state)
     slong i, j, v;
 
     {
-        fmpz_mpoly_t A, A1, A2, B;
+        slong nvarsAC = 2, nvarsB = 3;
+        fmpz_mpoly_t A, A1, A2, B, B1;
         fmpz_mpoly_struct * Cp[3];
         fmpz_mpoly_struct C[3];
         fmpz_mpoly_ctx_t ctxAC, ctxB;
+        slong * c;
 
-        fmpz_mpoly_ctx_init(ctxB, 3, ORD_LEX);
-        fmpz_mpoly_ctx_init(ctxAC, 2, ORD_LEX);
+        fmpz_mpoly_ctx_init(ctxB, nvarsB, ORD_LEX);
+        fmpz_mpoly_ctx_init(ctxAC, nvarsAC, ORD_LEX);
 
         fmpz_mpoly_init(B, ctxB);
         fmpz_mpoly_init(A, ctxAC);
@@ -44,12 +46,7 @@ TEST_FUNCTION_START(fmpz_mpoly_compose_fmpz_mpoly, state)
         if (fmpz_mpoly_compose_fmpz_mpoly(A, B, Cp, ctxB, ctxAC) ||
             fmpz_mpoly_compose_fmpz_mpoly_horner(A1, B, Cp, ctxB, ctxAC) ||
             fmpz_mpoly_compose_fmpz_mpoly_geobucket(A2, B, Cp, ctxB, ctxAC))
-        {
-            printf("FAIL\n");
-            flint_printf("Check non-example 1\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check non-example 1\n");
 
         fmpz_mpoly_set_str_pretty(C + 0, "x1", NULL, ctxAC);
         fmpz_mpoly_set_str_pretty(C + 1, "2*x2", NULL, ctxAC);
@@ -57,12 +54,7 @@ TEST_FUNCTION_START(fmpz_mpoly_compose_fmpz_mpoly, state)
         if (fmpz_mpoly_compose_fmpz_mpoly(A, B, Cp, ctxB, ctxAC) ||
             fmpz_mpoly_compose_fmpz_mpoly_horner(A1, B, Cp, ctxB, ctxAC) ||
             fmpz_mpoly_compose_fmpz_mpoly_geobucket(A2, B, Cp, ctxB, ctxAC))
-        {
-            printf("FAIL\n");
-            flint_printf("Check non-example 2\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check non-example 2\n");
 
         fmpz_mpoly_set_str_pretty(C + 0, "2*x1", NULL, ctxAC);
         fmpz_mpoly_set_str_pretty(C + 1, "x2", NULL, ctxAC);
@@ -70,12 +62,50 @@ TEST_FUNCTION_START(fmpz_mpoly_compose_fmpz_mpoly, state)
         if (!fmpz_mpoly_compose_fmpz_mpoly(A, B, Cp, ctxB, ctxAC) ||
             !fmpz_mpoly_compose_fmpz_mpoly_horner(A1, B, Cp, ctxB, ctxAC) ||
             !fmpz_mpoly_compose_fmpz_mpoly_geobucket(A2, B, Cp, ctxB, ctxAC))
-        {
-            printf("FAIL\n");
-            flint_printf("Check example 3\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check example 3\n");
+
+        /* Aliased generator composition */
+        c = (slong *) flint_malloc(nvarsB*sizeof(slong));
+        fmpz_mpoly_init(B1, ctxB);
+        fmpz_mpoly_set(B1, B, ctxB);
+        for (i = 0; i < nvarsB; i++)
+            c[i] = i;
+
+        fmpz_mpoly_compose_fmpz_mpoly_gen(B1, B1, c, ctxB, ctxB);
+        if (!fmpz_mpoly_equal(B, B1, ctxB))
+            TEST_FUNCTION_FAIL("Check composition with aliased generators\n");
+
+        /* Reverse the generators, twice */
+        for (i = 0; i < nvarsB; i++)
+            c[i] = nvarsB - i - 1;
+
+        fmpz_mpoly_compose_fmpz_mpoly_gen(B1, B1, c, ctxB, ctxB);
+        if (fmpz_mpoly_equal(B, B1, ctxB))
+            TEST_FUNCTION_FAIL("Check composition with reversed aliased generators\n");
+
+        fmpz_mpoly_compose_fmpz_mpoly_gen(B1, B1, c, ctxB, ctxB);
+        if (!fmpz_mpoly_equal(B, B1, ctxB))
+            TEST_FUNCTION_FAIL("Check composition with un-reversed aliased generators\n");
+
+        /* Composition with zero polys */
+        fmpz_mpoly_zero(B1, ctxB);
+
+        fmpz_mpoly_compose_fmpz_mpoly_gen(A, B1, c, ctxB, ctxAC);
+        if (!fmpz_mpoly_is_zero(B1, ctxB))
+            TEST_FUNCTION_FAIL("Check composition with generators of zero poly\n");
+
+        fmpz_mpoly_set_str_pretty(C + 0, "2*x1", NULL, ctxAC);
+        fmpz_mpoly_set_str_pretty(C + 1, "x2", NULL, ctxAC);
+        fmpz_mpoly_set_str_pretty(C + 2, "1", NULL, ctxAC);
+        if (!fmpz_mpoly_compose_fmpz_mpoly(A, B1, Cp, ctxB, ctxAC) ||
+            !fmpz_mpoly_compose_fmpz_mpoly_horner(A1, B1, Cp, ctxB, ctxAC) ||
+            !fmpz_mpoly_compose_fmpz_mpoly_geobucket(A2, B1, Cp, ctxB, ctxAC))
+            TEST_FUNCTION_FAIL("Check example 4\n");
+
+        if (!fmpz_mpoly_is_zero(A, ctxAC) ||
+            !fmpz_mpoly_is_zero(A1, ctxAC) ||
+            !fmpz_mpoly_is_zero(A2, ctxAC))
+            TEST_FUNCTION_FAIL("Check composition with zero poly\n");
 
         fmpz_mpoly_clear(B, ctxB);
         fmpz_mpoly_clear(A, ctxAC);
@@ -141,24 +171,14 @@ TEST_FUNCTION_START(fmpz_mpoly_compose_fmpz_mpoly, state)
             if (!fmpz_mpoly_compose_fmpz_mpoly(A1, B, C, ctxB, ctxAC) ||
                 !fmpz_mpoly_compose_fmpz_mpoly_horner(A2, B, C, ctxB, ctxAC) ||
                 !fmpz_mpoly_compose_fmpz_mpoly_geobucket(A3, B, C, ctxB, ctxAC))
-            {
-                printf("FAIL\n");
-                flint_printf("Check composition success with generators\n"
-                                                     "i: %wd, j: %wd\n", i, j);
-                fflush(stdout);
-                flint_abort();
-            }
+                TEST_FUNCTION_FAIL("Check composition success with generators\n"
+                                   "i: %wd, j: %wd\n", i, j);
 
             if (!fmpz_mpoly_equal(A, A1, ctxAC) ||
                 !fmpz_mpoly_equal(A, A2, ctxAC) ||
                 !fmpz_mpoly_equal(A, A3, ctxAC))
-            {
-                printf("FAIL\n");
-                flint_printf("Check composition with generators\n"
-                                                     "i: %wd, j: %wd\n", i, j);
-                fflush(stdout);
-                flint_abort();
-            }
+                TEST_FUNCTION_FAIL("Check composition with generators\n"
+                                   "i: %wd, j: %wd\n", i, j);
 
             fmpz_mpoly_assert_canonical(A, ctxAC);
             fmpz_mpoly_assert_canonical(A1, ctxAC);
@@ -219,24 +239,14 @@ TEST_FUNCTION_START(fmpz_mpoly_compose_fmpz_mpoly, state)
             !fmpz_mpoly_compose_fmpz_mpoly_geobucket(g2, f, vals1, ctx, ctx) ||
             !fmpz_mpoly_equal(g, g1, ctx) ||
             !fmpz_mpoly_equal(g, g2, ctx))
-        {
-            printf("FAIL\n");
-            flint_printf("Check composition success\ni: %wd\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check composition success\ni: %wd\n", i);
 
         fmpz_mpoly_assert_canonical(g, ctx);
         fmpz_mpoly_assert_canonical(g1, ctx);
         fmpz_mpoly_assert_canonical(g2, ctx);
 
         if (!fmpz_mpoly_equal(f, g, ctx))
-        {
-            printf("FAIL\n");
-            flint_printf("Check composition with identity\ni: %wd\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check composition with identity\ni: %wd\n", i);
 
         fmpz_mpoly_clear(f, ctx);
         fmpz_mpoly_clear(g, ctx);
@@ -311,12 +321,7 @@ TEST_FUNCTION_START(fmpz_mpoly_compose_fmpz_mpoly, state)
             vals3[v] = (fmpz *) flint_malloc(sizeof(fmpz));
             fmpz_init(vals3[v]);
             if (!fmpz_mpoly_evaluate_all_fmpz(vals3[v], vals1[v], vals2, ctx2))
-            {
-                printf("FAIL\n");
-                flint_printf("Check evaluation success\ni: %wd\n", i);
-                fflush(stdout);
-                flint_abort();
-            }
+                TEST_FUNCTION_FAIL("Check evaluation success\ni: %wd\n", i);
         }
 
         fmpz_mpoly_randtest_bound(f, state, len1, coeff_bits1, exp_bound1, ctx1);
@@ -326,12 +331,7 @@ TEST_FUNCTION_START(fmpz_mpoly_compose_fmpz_mpoly, state)
             !fmpz_mpoly_compose_fmpz_mpoly_geobucket(g2, f, vals1, ctx1, ctx2) ||
             !fmpz_mpoly_equal(g, g1, ctx2) ||
             !fmpz_mpoly_equal(g, g2, ctx2))
-        {
-            printf("FAIL\n");
-            flint_printf("Check composition success\ni: %wd\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check composition success\ni: %wd\n", i);
 
         fmpz_mpoly_assert_canonical(g, ctx2);
         fmpz_mpoly_assert_canonical(g1, ctx2);
@@ -339,20 +339,10 @@ TEST_FUNCTION_START(fmpz_mpoly_compose_fmpz_mpoly, state)
 
         if (!fmpz_mpoly_evaluate_all_fmpz(fe, f, vals3, ctx1) ||
             !fmpz_mpoly_evaluate_all_fmpz(ge, g, vals2, ctx2))
-        {
-            printf("FAIL\n");
-            flint_printf("Check evaluation success\ni: %wd\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check evaluation success\ni: %wd\n", i);
 
         if (!fmpz_equal(fe, ge))
-        {
-            printf("FAIL\n");
-            flint_printf("Check composition and evalall commute\ni: %wd\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check composition and evalall commute\ni: %wd\n", i);
 
         for (v = 0; v < nvars1; v++)
         {
diff --git a/src/nmod_mpoly/test/t-compose_nmod_mpoly.c b/src/nmod_mpoly/test/t-compose_nmod_mpoly.c
index 5127a5dc82..e56d175b30 100644
--- a/src/nmod_mpoly/test/t-compose_nmod_mpoly.c
+++ b/src/nmod_mpoly/test/t-compose_nmod_mpoly.c
@@ -18,13 +18,18 @@ TEST_FUNCTION_START(nmod_mpoly_compose_nmod_mpoly, state)
     slong i, j, v;
 
     {
-        nmod_mpoly_t A, A1, A2, B;
+        slong nvarsAC = 2, nvarsB = 3;
+        nmod_mpoly_t A, A1, A2, B, B1;
         nmod_mpoly_struct * Cp[3];
         nmod_mpoly_struct C[3];
         nmod_mpoly_ctx_t ctxAC, ctxB;
+        slong * c;
+        ulong modulus;
 
-        nmod_mpoly_ctx_init(ctxB, 3, ORD_LEX, 13);
-        nmod_mpoly_ctx_init(ctxAC, 2, ORD_LEX, 13);
+        modulus = n_randint(state, FLINT_BITS - 1) + 1;
+        modulus = n_randbits(state, modulus);
+        nmod_mpoly_ctx_init(ctxB, nvarsB, ORD_LEX, modulus);
+        nmod_mpoly_ctx_init(ctxAC, nvarsAC, ORD_LEX, modulus);
 
         nmod_mpoly_init(B, ctxB);
         nmod_mpoly_init(A, ctxAC);
@@ -45,12 +50,7 @@ TEST_FUNCTION_START(nmod_mpoly_compose_nmod_mpoly, state)
         if (nmod_mpoly_compose_nmod_mpoly(A, B, Cp, ctxB, ctxAC) ||
             nmod_mpoly_compose_nmod_mpoly_horner(A1, B, Cp, ctxB, ctxAC) ||
             nmod_mpoly_compose_nmod_mpoly_geobucket(A2, B, Cp, ctxB, ctxAC))
-        {
-            printf("FAIL\n");
-            flint_printf("Check non-example 1\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check non-example 1\n");
 
         nmod_mpoly_set_str_pretty(C + 0, "x1", NULL, ctxAC);
         nmod_mpoly_set_str_pretty(C + 1, "2*x2", NULL, ctxAC);
@@ -58,12 +58,7 @@ TEST_FUNCTION_START(nmod_mpoly_compose_nmod_mpoly, state)
         if (!nmod_mpoly_compose_nmod_mpoly(A, B, Cp, ctxB, ctxAC) ||
             !nmod_mpoly_compose_nmod_mpoly_horner(A1, B, Cp, ctxB, ctxAC) ||
             !nmod_mpoly_compose_nmod_mpoly_geobucket(A2, B, Cp, ctxB, ctxAC))
-        {
-            printf("FAIL\n");
-            flint_printf("Check non-example 2\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check non-example 2\n");
 
         nmod_mpoly_set_str_pretty(C + 0, "2*x1", NULL, ctxAC);
         nmod_mpoly_set_str_pretty(C + 1, "x2", NULL, ctxAC);
@@ -71,12 +66,50 @@ TEST_FUNCTION_START(nmod_mpoly_compose_nmod_mpoly, state)
         if (!nmod_mpoly_compose_nmod_mpoly(A, B, Cp, ctxB, ctxAC) ||
             !nmod_mpoly_compose_nmod_mpoly_horner(A1, B, Cp, ctxB, ctxAC) ||
             !nmod_mpoly_compose_nmod_mpoly_geobucket(A2, B, Cp, ctxB, ctxAC))
-        {
-            printf("FAIL\n");
-            flint_printf("Check example 3\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check example 3\n");
+
+        /* Aliased generator composition */
+        c = (slong *) flint_malloc(nvarsB*sizeof(slong));
+        nmod_mpoly_init(B1, ctxB);
+        nmod_mpoly_set(B1, B, ctxB);
+        for (i = 0; i < nvarsB; i++)
+            c[i] = i;
+
+        nmod_mpoly_compose_nmod_mpoly_gen(B1, B1, c, ctxB, ctxB);
+        if (!nmod_mpoly_equal(B, B1, ctxB))
+            TEST_FUNCTION_FAIL("Check composition with aliased generators\n");
+
+        /* Reverse the generators, twice */
+        for (i = 0; i < nvarsB; i++)
+            c[i] = nvarsB - i - 1;
+
+        nmod_mpoly_compose_nmod_mpoly_gen(B1, B1, c, ctxB, ctxB);
+        if (nmod_mpoly_equal(B, B1, ctxB))
+            TEST_FUNCTION_FAIL("Check composition with reversed aliased generators\n");
+
+        nmod_mpoly_compose_nmod_mpoly_gen(B1, B1, c, ctxB, ctxB);
+        if (!nmod_mpoly_equal(B, B1, ctxB))
+            TEST_FUNCTION_FAIL("Check composition with un-reversed aliased generators\n");
+
+        /* Composition with zero polys */
+        nmod_mpoly_zero(B1, ctxB);
+
+        nmod_mpoly_compose_nmod_mpoly_gen(A, B1, c, ctxB, ctxAC);
+        if (!nmod_mpoly_is_zero(B1, ctxB))
+            TEST_FUNCTION_FAIL("Check composition with generators of zero poly\n");
+
+        nmod_mpoly_set_str_pretty(C + 0, "2*x1", NULL, ctxAC);
+        nmod_mpoly_set_str_pretty(C + 1, "x2", NULL, ctxAC);
+        nmod_mpoly_set_str_pretty(C + 2, "1", NULL, ctxAC);
+        if (!nmod_mpoly_compose_nmod_mpoly(A, B1, Cp, ctxB, ctxAC) ||
+            !nmod_mpoly_compose_nmod_mpoly_horner(A1, B1, Cp, ctxB, ctxAC) ||
+            !nmod_mpoly_compose_nmod_mpoly_geobucket(A2, B1, Cp, ctxB, ctxAC))
+            TEST_FUNCTION_FAIL("Check example 4\n");
+
+        if (!nmod_mpoly_is_zero(A, ctxAC) ||
+            !nmod_mpoly_is_zero(A1, ctxAC) ||
+            !nmod_mpoly_is_zero(A2, ctxAC))
+            TEST_FUNCTION_FAIL("Check composition with zero poly\n");
 
         nmod_mpoly_clear(B, ctxB);
         nmod_mpoly_clear(A, ctxAC);
@@ -145,24 +178,14 @@ TEST_FUNCTION_START(nmod_mpoly_compose_nmod_mpoly, state)
             if (!nmod_mpoly_compose_nmod_mpoly(A1, B, C, ctxB, ctxAC) ||
                 !nmod_mpoly_compose_nmod_mpoly_horner(A2, B, C, ctxB, ctxAC) ||
                 !nmod_mpoly_compose_nmod_mpoly_geobucket(A3, B, C, ctxB, ctxAC))
-            {
-                printf("FAIL\n");
-                flint_printf("Check composition success with generators\n"
-                                                     "i: %wd, j: %wd\n", i, j);
-                fflush(stdout);
-                flint_abort();
-            }
+                TEST_FUNCTION_FAIL("Check composition success with generators\n"
+                                   "i: %wd, j: %wd\n", i, j);
 
             if (!nmod_mpoly_equal(A, A1, ctxAC) ||
                 !nmod_mpoly_equal(A, A2, ctxAC) ||
                 !nmod_mpoly_equal(A, A3, ctxAC))
-            {
-                printf("FAIL\n");
-                flint_printf("Check composition with generators\n"
-                                                     "i: %wd, j: %wd\n", i, j);
-                fflush(stdout);
-                flint_abort();
-            }
+                TEST_FUNCTION_FAIL("Check composition with generators\n"
+                                   "i: %wd, j: %wd\n", i, j);
 
             nmod_mpoly_assert_canonical(A, ctxAC);
             nmod_mpoly_assert_canonical(A1, ctxAC);
@@ -225,24 +248,14 @@ TEST_FUNCTION_START(nmod_mpoly_compose_nmod_mpoly, state)
             !nmod_mpoly_compose_nmod_mpoly_geobucket(g2, f, vals1, ctx, ctx) ||
             !nmod_mpoly_equal(g, g1, ctx) ||
             !nmod_mpoly_equal(g, g2, ctx))
-        {
-            printf("FAIL\n");
-            flint_printf("Check composition success\ni: %wd\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check composition success\ni: %wd\n", i);
 
         nmod_mpoly_assert_canonical(g, ctx);
         nmod_mpoly_assert_canonical(g1, ctx);
         nmod_mpoly_assert_canonical(g2, ctx);
 
         if (!nmod_mpoly_equal(f, g, ctx))
-        {
-            printf("FAIL\n");
-            flint_printf("Check composition with identity\ni: %wd\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check composition with identity\ni: %wd\n", i);
 
         nmod_mpoly_clear(f, ctx);
         nmod_mpoly_clear(g, ctx);
@@ -319,12 +332,7 @@ TEST_FUNCTION_START(nmod_mpoly_compose_nmod_mpoly, state)
             !nmod_mpoly_compose_nmod_mpoly_geobucket(g2, f, vals1, ctx1, ctx2) ||
             !nmod_mpoly_equal(g, g1, ctx2) ||
             !nmod_mpoly_equal(g, g2, ctx2))
-        {
-            printf("FAIL\n");
-            flint_printf("Check composition success\ni: %wd\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check composition success\ni: %wd\n", i);
 
         nmod_mpoly_assert_canonical(g, ctx2);
         nmod_mpoly_assert_canonical(g1, ctx2);
@@ -334,12 +342,7 @@ TEST_FUNCTION_START(nmod_mpoly_compose_nmod_mpoly, state)
         ge = nmod_mpoly_evaluate_all_ui(g, vals2, ctx2);
 
         if (fe != ge)
-        {
-            printf("FAIL\n");
-            flint_printf("Check composition and evalall commute\ni: %wd\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check composition and evalall commute\ni: %wd\n", i);
 
         for (v = 0; v < nvars1; v++)
         {
@@ -404,12 +407,7 @@ TEST_FUNCTION_START(nmod_mpoly_compose_nmod_mpoly, state)
         if (!nmod_mpoly_compose_nmod_mpoly(g, f, vals1, ctx1, ctx2) ||
             !nmod_mpoly_compose_nmod_mpoly_geobucket(g1, f, vals1, ctx1, ctx2) ||
             !nmod_mpoly_compose_nmod_mpoly_horner(g2, f, vals1, ctx1, ctx2))
-        {
-            printf("FAIL\n");
-            flint_printf("Check composition success\ni: %wd\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check composition success\ni: %wd\n", i);
 
         nmod_mpoly_assert_canonical(g, ctx2);
         nmod_mpoly_assert_canonical(g1, ctx2);
@@ -417,12 +415,7 @@ TEST_FUNCTION_START(nmod_mpoly_compose_nmod_mpoly, state)
 
         if (!nmod_mpoly_is_ui(g, ctx2) ||
             nmod_mpoly_get_ui(g, ctx2) != nmod_mpoly_evaluate_all_ui(f, vals2, ctx1))
-        {
-            printf("FAIL\n");
-            flint_printf("Check composition with constants matches evalall\ni: %wd\n", i);
-            fflush(stdout);
-            flint_abort();
-        }
+            TEST_FUNCTION_FAIL("Check composition with constants matches evalall\ni: %wd\n", i);
 
         for (v = 0; v < nvars1; v++)
         {