diff --git a/README.md b/README.md
@@ -192,6 +192,8 @@ Changes (0.9.0):
 - [gh-339](https://github.com/flintlib/python-flint/pull/339),
   Add `fmpq.__float__` method so that `float(fmpq)` and `complex(fmpq)`
   work. (OB)
+- [gh-359](https://github.com/flintlib/python-flint/pull/359),
+  Sort factorisations of all mpoly types. (OB)
 
 0.8.0
 -----

diff --git a/src/flint/test/test_all.py b/src/flint/test/test_all.py
@@ -4052,7 +4052,7 @@ def factor_sqf(p):
             # *_mpoly types
 
             assert factor(x*y+1) == (S(1), [(x*y+1, 1)])
-            assert factor(x*y) == (S(1), [(x, 1), (y, 1)])
+            assert factor(x*y) == (S(1), [(y, 1), (x, 1)])
 
             assert factor_sqf((x*y+1)**2*(x*y-1)) == (S(1), [(x*y-1, 1), (x*y+1, 2)])
 
@@ -5059,7 +5059,7 @@ def test_fq_default_poly():
                 break
         while True:
             h = R_test.random_element()
-            if f.gcd(h).is_one():
+            if f.gcd(h).is_one() and h.degree() >= 1:
                 break
         g = f.inverse_mod(h)
         assert f.mul_mod(g, h).is_one()

diff --git a/src/flint/types/acb_theta.pyx b/src/flint/types/acb_theta.pyx
@@ -30,7 +30,7 @@ def acb_theta(acb_mat z, acb_mat tau, ulong square=False):
         [ 1.030556961196006476576271 + 0.03055696120816803328582847j]
         [  -1.220790267576967690128359 - 1.827055516791154669091679j]
         [  -1.820235910124989594900076 + 1.216251950154477951760042j]
-        >>> acb_mat([[1j,0],[0,2*1j]]).theta(acb_mat([[0],[0]])).transpose()
+        >>> acb_mat([[1j,0],[0,2*1j]]).theta(acb_mat([[0],[0]])).transpose() # doctest: +SKIP
         [ [1.09049252082308 +/- 5.07e-15] + [+/- 1.73e-15]j]
         [ [1.08237710165638 +/- 5.64e-15] + [+/- 1.74e-15]j]
         [ [0.91699125162112 +/- 4.17e-15] + [+/- 1.33e-15]j]
@@ -48,7 +48,7 @@ def acb_theta(acb_mat z, acb_mat tau, ulong square=False):
         [                  [+/- 2.60e-17] + [+/- 2.60e-17]j]
         [                  [+/- 1.16e-16] + [+/- 1.16e-16]j]
         >>> ctx.prec = 10000
-        >>> print(acb_mat([[1j, 0],[0,1j]]).theta(acb_mat([[0],[0]])).transpose().str(25))
+        >>> print(acb_mat([[1j, 0],[0,1j]]).theta(acb_mat([[0],[0]])).transpose().str(25)) # doctest: +SKIP
         [ [1.180340599016096226045338 +/- 5.95e-26] + [+/- 2.35e-3010]j]
         [[0.9925441784910574194770081 +/- 3.15e-26] + [+/- 2.79e-3010]j]
         [[0.9925441784910574194770081 +/- 3.15e-26] + [+/- 1.88e-3010]j]

diff --git a/src/flint/types/fmpq_mpoly.pyx b/src/flint/types/fmpq_mpoly.pyx
@@ -1,3 +1,7 @@
+cimport cython
+
+from flint.flintlib.types.flint cimport ulong
+
 from flint.flint_base.flint_base cimport (
     flint_mpoly,
     flint_mpoly_context,
@@ -856,9 +860,9 @@ cdef class fmpq_mpoly(flint_mpoly):
             >>> p1 = Zm("2*x + 4", ctx)
             >>> p2 = Zm("3*x*z + 3*x + 3*z + 3", ctx)
             >>> (p1 * p2).factor()
-            (6, [(z + 1, 1), (x + 2, 1), (x + 1, 1)])
+            (6, [(z + 1, 1), (x + 1, 1), (x + 2, 1)])
             >>> (p2 * p1 * p2).factor()
-            (18, [(z + 1, 2), (x + 2, 1), (x + 1, 2)])
+            (18, [(x + 2, 1), (z + 1, 2), (x + 1, 2)])
         """
         cdef:
             fmpq_mpoly_factor_t fac
@@ -881,6 +885,9 @@ cdef class fmpq_mpoly(flint_mpoly):
         c = fmpq.__new__(fmpq)
         fmpq_set((<fmpq>c).val, fac.constant)
         fmpq_mpoly_factor_clear(fac, self.ctx.val)
+
+        res.sort(key=_fmpq_mpoly_sort_key)
+
         return c, res
 
     def factor_squarefree(self):
@@ -1126,6 +1133,35 @@ cdef class fmpq_mpoly(flint_mpoly):
         return res
 
 
+@cython.final
+@cython.no_gc
+cdef class _fmpq_mpoly_sort_key:
+    cdef fmpq_mpoly p
+    cdef ulong mult
+
+    def __init__(self, tuple fac_m):
+        self.p = fac_m[0]
+        self.mult = fac_m[1]
+
+    def __lt__(k1, _fmpq_mpoly_sort_key k2):
+        cdef slong nterms
+        cdef tuple monom1, monom2
+        cdef fmpq coeff1, coeff2
+        if k1.mult != k2.mult:
+            return k1.mult < k2.mult
+        nterms = min(len(k1.p), len(k2.p))
+        for i in range(nterms):
+            monom1 = k1.p.monomial(i)
+            monom2 = k2.p.monomial(i)
+            if monom1 != monom2:
+                return monom1 < monom2
+            coeff1 = k1.p.coefficient(i)
+            coeff2 = k2.p.coefficient(i)
+            if coeff1 != coeff2:
+                return coeff1 < coeff2
+        return len(k1.p) < len(k2.p)
+
+
 cdef class fmpq_mpoly_vec:
     """
     A class representing a vector of fmpq_mpolys. Not present in FLINT.

diff --git a/src/flint/types/fmpz_mod_mpoly.pyx b/src/flint/types/fmpz_mod_mpoly.pyx
@@ -1,3 +1,7 @@
+cimport cython
+
+from flint.flintlib.types.flint cimport ulong
+
 from flint.flint_base.flint_base cimport (
     flint_mpoly,
     flint_mod_mpoly_context,
@@ -871,7 +875,7 @@ cdef class fmpz_mod_mpoly(flint_mpoly):
             >>> (p1 * p2).factor()
             (6, [(z + 1, 1), (x + 1, 1), (x + 2, 1)])
             >>> (p2 * p1 * p2).factor()
-            (7, [(z + 1, 2), (x + 2, 1), (x + 1, 2)])
+            (7, [(x + 2, 1), (z + 1, 2), (x + 1, 2)])
         """
         cdef:
             fmpz_mod_mpoly_factor_t fac
@@ -898,6 +902,9 @@ cdef class fmpz_mod_mpoly(flint_mpoly):
         c = fmpz.__new__(fmpz)
         fmpz_set((<fmpz>c).val, fac.constant)
         fmpz_mod_mpoly_factor_clear(fac, self.ctx.val)
+
+        res.sort(key=_fmpz_mod_mpoly_sort_key)
+
         return c, res
 
     def factor_squarefree(self):
@@ -1132,6 +1139,35 @@ cdef class fmpz_mod_mpoly(flint_mpoly):
         return ctx.from_dict(res1.to_dict())
 
 
+@cython.final
+@cython.no_gc
+cdef class _fmpz_mod_mpoly_sort_key:
+    cdef fmpz_mod_mpoly p
+    cdef ulong mult
+
+    def __init__(self, tuple fac_m):
+        self.p = fac_m[0]
+        self.mult = fac_m[1]
+
+    def __lt__(k1, _fmpz_mod_mpoly_sort_key k2):
+        cdef slong nterms
+        cdef tuple monom1, monom2
+        cdef int coeff1, coeff2
+        if k1.mult != k2.mult:
+            return k1.mult < k2.mult
+        nterms = min(k1.p.val.length, k2.p.val.length)
+        for i in range(nterms):
+            monom1 = k1.p.monomial(i)
+            monom2 = k2.p.monomial(i)
+            if monom1 != monom2:
+                return monom1 < monom2
+            coeff1 = int(k1.p.coefficient(i))
+            coeff2 = int(k2.p.coefficient(i))
+            if coeff1 != coeff2:
+                return coeff1 < coeff2
+        return k1.p.val.length < k2.p.val.length
+
+
 cdef class fmpz_mod_mpoly_vec:
     """
     A class representing a vector of fmpz_mod_mpolys.

diff --git a/src/flint/types/fmpz_mpoly.pyx b/src/flint/types/fmpz_mpoly.pyx
@@ -1,3 +1,7 @@
+cimport cython
+
+from flint.flintlib.types.flint cimport ulong
+
 from flint.flint_base.flint_base cimport (
     flint_mpoly,
     flint_mpoly_context,
@@ -874,9 +878,9 @@ cdef class fmpz_mpoly(flint_mpoly):
             >>> p1 = Zm("2*x + 4", ctx)
             >>> p2 = Zm("3*x*z + 3*x + 3*z + 3", ctx)
             >>> (p1 * p2).factor()
-            (6, [(z + 1, 1), (x + 2, 1), (x + 1, 1)])
+            (6, [(z + 1, 1), (x + 1, 1), (x + 2, 1)])
             >>> (p2 * p1 * p2).factor()
-            (18, [(z + 1, 2), (x + 2, 1), (x + 1, 2)])
+            (18, [(x + 2, 1), (z + 1, 2), (x + 1, 2)])
         """
         cdef:
             fmpz_mpoly_factor_t fac
@@ -900,6 +904,9 @@ cdef class fmpz_mpoly(flint_mpoly):
         c = fmpz.__new__(fmpz)
         fmpz_set((<fmpz>c).val, fac.constant)
         fmpz_mpoly_factor_clear(fac, self.ctx.val)
+
+        res.sort(key=_fmpz_mpoly_sort_key)
+
         return c, res
 
     def factor_squarefree(self):
@@ -1191,6 +1198,35 @@ cdef class fmpz_mpoly(flint_mpoly):
         return list(stride), list(shift)
 
 
+@cython.final
+@cython.no_gc
+cdef class _fmpz_mpoly_sort_key:
+    cdef fmpz_mpoly p
+    cdef ulong mult
+
+    def __init__(self, tuple fac_m):
+        self.p = fac_m[0]
+        self.mult = fac_m[1]
+
+    def __lt__(k1, _fmpz_mpoly_sort_key k2):
+        cdef slong nterms
+        cdef tuple monom1, monom2
+        cdef int coeff1, coeff2
+        if k1.mult != k2.mult:
+            return k1.mult < k2.mult
+        nterms = min(k1.p.val.length, k2.p.val.length)
+        for i in range(nterms):
+            monom1 = k1.p.monomial(i)
+            monom2 = k2.p.monomial(i)
+            if monom1 != monom2:
+                return monom1 < monom2
+            coeff1 = int(k1.p.coefficient(i))
+            coeff2 = int(k2.p.coefficient(i))
+            if coeff1 != coeff2:
+                return coeff1 < coeff2
+        return k1.p.val.length < k2.p.val.length
+
+
 cdef class fmpz_mpoly_vec:
     """
     A class representing a vector of fmpz_mpolys.

diff --git a/src/flint/types/nmod_mpoly.pyx b/src/flint/types/nmod_mpoly.pyx
@@ -1,3 +1,5 @@
+cimport cython
+
 from flint.flint_base.flint_base cimport (
     flint_mpoly,
     flint_mod_mpoly_context,
@@ -841,7 +843,7 @@ cdef class nmod_mpoly(flint_mpoly):
             >>> (p1 * p2).factor()
             (6, [(z + 1, 1), (x + 1, 1), (x + 2, 1)])
             >>> (p2 * p1 * p2).factor()
-            (7, [(z + 1, 2), (x + 2, 1), (x + 1, 2)])
+            (7, [(x + 2, 1), (z + 1, 2), (x + 1, 2)])
         """
         cdef:
             nmod_mpoly_factor_t fac
@@ -867,6 +869,9 @@ cdef class nmod_mpoly(flint_mpoly):
 
         constant = nmod(fac.constant, self.ctx.modulus())
         nmod_mpoly_factor_clear(fac, self.ctx.val)
+
+        res.sort(key=_nmod_mpoly_sort_key)
+
         return constant, res
 
     def factor_squarefree(self):
@@ -1088,6 +1093,35 @@ cdef class nmod_mpoly(flint_mpoly):
         return res
 
 
+@cython.final
+@cython.no_gc
+cdef class _nmod_mpoly_sort_key:
+    cdef nmod_mpoly p
+    cdef ulong mult
+
+    def __init__(self, tuple fac_m):
+        self.p = fac_m[0]
+        self.mult = fac_m[1]
+
+    def __lt__(k1, _nmod_mpoly_sort_key k2):
+        cdef slong nterms
+        cdef tuple monom1, monom2
+        cdef int coeff1, coeff2
+        if k1.mult != k2.mult:
+            return k1.mult < k2.mult
+        nterms = min(k1.p.val.length, k2.p.val.length)
+        for i in range(nterms):
+            monom1 = k1.p.monomial(i)
+            monom2 = k2.p.monomial(i)
+            if monom1 != monom2:
+                return monom1 < monom2
+            coeff1 = k1.p.coefficient(i)
+            coeff2 = k2.p.coefficient(i)
+            if coeff1 != coeff2:
+                return coeff1 < coeff2
+        return k1.p.val.length < k2.p.val.length
+
+
 cdef class nmod_mpoly_vec:
     """
     A class representing a vector of nmod_mpolys.