diff --git a/doc/source/fmpz_factor.rst b/doc/source/fmpz_factor.rst
index b8bca58aa2..7085a213dc 100644
--- a/doc/source/fmpz_factor.rst
+++ b/doc/source/fmpz_factor.rst
@@ -176,6 +176,35 @@ A separate ``int`` field holds the sign, which may be `-1`, `0` or `1`.
     https://maths-people.anu.edu.au/~brent/pd/rpb051i.pdf
 
 
+Utility
+--------------------------------------------------------------------------------
+
+.. function:: void fmpz_factor_sort(fmpz_factor_t factor)
+
+    Sorts ``factor`` in-place such that the factors are in ascending order.
+
+Arithmetic
+--------------------------------------------------------------------------------
+
+.. function:: void _fmpz_factor_mul(fmpz_factor_t factor, const fmpz_factor_t a, const fmpz_factor_t b)
+
+    Sets ``factor`` to the factored representation of `a \times b` assuming ``a`` and ``b``
+    are sorted.
+
+.. function:: void fmpz_factor_mul(fmpz_factor_t factor, const fmpz_factor_t a, const fmpz_factor_t b)
+
+    Sets ``factor`` to the factored representation of `a \times b` by first sorting ``a`` and ``b``.
+
+.. function:: void _fmpz_factor_gcd(fmpz_factor_t factor, const fmpz_factor_t a, const fmpz_factor_t b)
+
+    Sets ``factor`` to the factored representation of `\gcd(a, b)` assuming ``a`` and ``b``
+    are sorted.
+
+.. function:: void fmpz_factor_gcd(fmpz_factor_t factor, const fmpz_factor_t a, const fmpz_factor_t b)
+
+    Sets ``factor`` to the factored representation of `\gcd(a, b)` by first sorting ``a`` and ``b``.
+
+
 Input and output
 --------------------------------------------------------------------------------
 
diff --git a/src/fmpz_factor.h b/src/fmpz_factor.h
index bfafd94fa1..35cfa56563 100644
--- a/src/fmpz_factor.h
+++ b/src/fmpz_factor.h
@@ -38,6 +38,16 @@ void _fmpz_factor_append_ui(fmpz_factor_t factor, ulong p, ulong exp);
 
 void _fmpz_factor_concat(fmpz_factor_t factor1, fmpz_factor_t factor2, ulong exp);
 
+void fmpz_factor_sort(fmpz_factor_t factor);
+
+/* Arithmetic ****************************************************************/
+
+void _fmpz_factor_mul(fmpz_factor_t factor, const fmpz_factor_t a, const fmpz_factor_t b);
+void fmpz_factor_gcd(fmpz_factor_t factor, fmpz_factor_t a, fmpz_factor_t b);
+
+void _fmpz_factor_gcd(fmpz_factor_t factor, const fmpz_factor_t a, const fmpz_factor_t b);
+void fmpz_factor_mul(fmpz_factor_t factor, fmpz_factor_t a, fmpz_factor_t b);
+
 /* I/O ***********************************************************************/
 
 #ifdef FLINT_HAVE_FILE
diff --git a/src/fmpz_factor/gcd.c b/src/fmpz_factor/gcd.c
new file mode 100644
index 0000000000..7338f9593c
--- /dev/null
+++ b/src/fmpz_factor/gcd.c
@@ -0,0 +1,59 @@
+/*
+    Copyright (C) 2025 Mael Hostettler
+
+    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_factor.h"
+
+void
+_fmpz_factor_gcd(fmpz_factor_t factor, const fmpz_factor_t a, const fmpz_factor_t b)
+{
+    slong i, j;
+    ulong num_len;
+    int cmp;
+
+    num_len = FLINT_MIN(a->num, b->num);
+    _fmpz_factor_fit_length(factor, num_len);
+    factor->sign = 1;
+    factor->num = 0;
+
+    i = 0;
+    j = 0;
+    factor->num = 0;
+
+    while ((i < a->num) && (j < b->num))
+    {
+        cmp = fmpz_cmp(a->p + i, b->p + j);
+        if (cmp < 0) // a[i] < b[i]
+        {
+            i++;
+        }
+        else if (cmp > 0) // b[i] < a[i]
+        {
+            j++;
+        }
+        else // b[i] = a[i]
+        {
+            fmpz_set(factor->p + factor->num, a->p + i);
+            factor->exp[factor->num] = FLINT_MIN(a->exp[i], b->exp[j]);
+            factor->num++; 
+            i++;
+            j++;
+        }
+    }
+}
+
+void fmpz_factor_gcd(fmpz_factor_t factor, fmpz_factor_t a, fmpz_factor_t b)
+{
+    fmpz_factor_sort(a);
+    fmpz_factor_sort(b);
+
+    _fmpz_factor_gcd(factor, a, b);
+}
diff --git a/src/fmpz_factor/mul.c b/src/fmpz_factor/mul.c
new file mode 100644
index 0000000000..376e660092
--- /dev/null
+++ b/src/fmpz_factor/mul.c
@@ -0,0 +1,78 @@
+/*
+    Copyright (C) 2025 Mael Hostettler
+
+    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_factor.h"
+
+void
+_fmpz_factor_mul(fmpz_factor_t factor, const fmpz_factor_t a, const fmpz_factor_t b)
+{
+    slong i, j;
+    int cmp;
+
+    _fmpz_factor_fit_length(factor, a->num + b->num);
+    factor->sign = a->sign * b->sign;
+    
+    i = 0;
+    j = 0;
+    factor->num = 0;
+
+    while ((i < a->num) && (j < b->num))
+    {
+        cmp = fmpz_cmp(a->p + i, b->p + j);
+        if (cmp < 0) // a[i] < b[i]
+        {
+            fmpz_set(factor->p + factor->num, a->p + i);
+            factor->exp[factor->num] = a->exp[i];
+            factor->num++; 
+            i++;
+        }
+        else if (cmp > 0) // b[i] < a[i]
+        {
+            fmpz_set(factor->p + factor->num, b->p + j);
+            factor->exp[factor->num] = b->exp[j];
+            factor->num++; 
+            j++;
+        }
+        else // b[i] = a[i]
+        {
+            fmpz_set(factor->p + factor->num, a->p + i);
+            factor->exp[factor->num] = a->exp[i] + b->exp[j];
+            factor->num++; 
+            i++;
+            j++;
+        }
+    }
+
+    while (i < a->num)
+    {
+        fmpz_set(factor->p + factor->num, a->p + i);
+        factor->exp[factor->num] = a->exp[i];
+        factor->num++; 
+        i++;
+    }
+
+    while (j < b->num)
+    {
+        fmpz_set(factor->p + factor->num, b->p + j);
+        factor->exp[factor->num] = b->exp[j];
+        factor->num++; 
+        j++;
+    }
+}
+
+void fmpz_factor_mul(fmpz_factor_t factor, fmpz_factor_t a, fmpz_factor_t b)
+{
+    fmpz_factor_sort(a);
+    fmpz_factor_sort(b);
+
+    _fmpz_factor_mul(factor, a, b);
+}
diff --git a/src/fmpz_factor/sort.c b/src/fmpz_factor/sort.c
new file mode 100644
index 0000000000..a93b446b92
--- /dev/null
+++ b/src/fmpz_factor/sort.c
@@ -0,0 +1,43 @@
+/*
+    Copyright (C) 2025 Mael Hostettler
+
+    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_factor.h"
+
+// Tiny implem of quick sort, qsort is not usable because we also want to permute exp
+// and qsort_r / qsort_s have portability issues
+static void _fmpz_factor_sort_part(fmpz* p, ulong* exp, slong n)
+{
+    slong i, j;
+    if (n <= 1) return;
+
+    fmpz_swap(p + n / 2, p + n - 1);
+    FLINT_SWAP(ulong, exp[n / 2], exp[n - 1]);
+    for (i = 0, j = 0; i < n - 1; i++)
+    {
+        if (fmpz_cmp(p + i, p + n - 1) < 0)
+        {
+            fmpz_swap(p + i, p + j);
+            FLINT_SWAP(ulong, exp[i], exp[j]);
+            j++;
+        }
+    }
+    fmpz_swap(p + j, p + n - 1);
+    FLINT_SWAP(ulong, exp[j], exp[n - 1]);
+
+    _fmpz_factor_sort_part(p, exp, j);
+    _fmpz_factor_sort_part(p + j + 1, exp + j + 1, n - j - 1);
+}
+
+void
+fmpz_factor_sort(fmpz_factor_t factor)
+{
+    _fmpz_factor_sort_part(factor->p, factor->exp, factor->num);
+}
diff --git a/src/fmpz_factor/test/main.c b/src/fmpz_factor/test/main.c
index 8a09a8519f..b6b5fe2bac 100644
--- a/src/fmpz_factor/test/main.c
+++ b/src/fmpz_factor/test/main.c
@@ -22,7 +22,9 @@
 #include "t-pollard_brent.c"
 #include "t-pollard_brent_single.c"
 #include "t-refine.c"
-
+#include "t-sort.c"
+#include "t-mul.c"
+#include "t-gcd.c"
 /* Array of test functions ***************************************************/
 
 test_struct tests[] =
@@ -35,7 +37,10 @@ test_struct tests[] =
     TEST_FUNCTION(fmpz_factor_trial),
     TEST_FUNCTION(fmpz_factor_pollard_brent),
     TEST_FUNCTION(fmpz_factor_pollard_brent_single),
-    TEST_FUNCTION(fmpz_factor_refine)
+    TEST_FUNCTION(fmpz_factor_refine),
+    TEST_FUNCTION(fmpz_factor_sort),
+    TEST_FUNCTION(fmpz_factor_mul),
+    TEST_FUNCTION(fmpz_factor_gcd)
 };
 
 /* main function *************************************************************/
diff --git a/src/fmpz_factor/test/t-gcd.c b/src/fmpz_factor/test/t-gcd.c
new file mode 100644
index 0000000000..b7a52108ad
--- /dev/null
+++ b/src/fmpz_factor/test/t-gcd.c
@@ -0,0 +1,65 @@
+/*
+    Copyright (C) 2025 Mael Hostettler
+
+    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.h"
+#include "fmpz_factor.h"
+
+TEST_FUNCTION_START(fmpz_factor_gcd, state)
+{
+    slong i, b;
+    fmpz_t x, y;
+    fmpz_factor_t fx, fy, fz;
+
+    fmpz_init(x);
+    fmpz_init(y);
+    fmpz_factor_init(fx);
+    fmpz_factor_init(fy);
+    fmpz_factor_init(fz);
+
+    for (i = 0; i < 1000 * flint_test_multiplier(); i++)
+    {
+        b = n_randint(state, 1000000);
+        if (n_randint(state, 2))
+            b = -b;
+        fmpz_set_si(x, b);
+
+        b = n_randint(state, 1000000);
+        if (n_randint(state, 2))
+            b = -b;
+        fmpz_set_si(y, b);
+
+        fmpz_factor(fx, x);
+        fmpz_factor(fy, y);
+
+        fmpz_factor_gcd(fz, fx, fy);
+
+        fmpz_gcd(x, x, y);
+        fmpz_factor_expand(y, fz);
+
+        if (fmpz_cmp(x, y) != 0)
+        {
+            flint_printf("FAIL (factor gcd)\n");
+            fmpz_print(x); flint_printf(" ");
+            fmpz_print(y); flint_printf("\n");
+            fflush(stdout);
+            flint_abort();
+        }
+    }
+
+    fmpz_clear(x);
+    fmpz_clear(y);
+    fmpz_factor_clear(fx);
+    fmpz_factor_clear(fy);
+    fmpz_factor_clear(fz);
+
+    TEST_FUNCTION_END(state);
+}
diff --git a/src/fmpz_factor/test/t-mul.c b/src/fmpz_factor/test/t-mul.c
new file mode 100644
index 0000000000..2d23b85cb2
--- /dev/null
+++ b/src/fmpz_factor/test/t-mul.c
@@ -0,0 +1,65 @@
+/*
+    Copyright (C) 2025 Mael Hostettler
+
+    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.h"
+#include "fmpz_factor.h"
+
+TEST_FUNCTION_START(fmpz_factor_mul, state)
+{
+    slong i, b;
+    fmpz_t x, y;
+    fmpz_factor_t fx, fy, fz;
+
+    fmpz_init(x);
+    fmpz_init(y);
+    fmpz_factor_init(fx);
+    fmpz_factor_init(fy);
+    fmpz_factor_init(fz);
+
+    for (i = 0; i < 1000 * flint_test_multiplier(); i++)
+    {
+        b = n_randint(state, 1000000);
+        if (n_randint(state, 2))
+            b = -b;
+        fmpz_set_si(x, b);
+
+        b = n_randint(state, 1000000);
+        if (n_randint(state, 2))
+            b = -b;
+        fmpz_set_si(y, b);
+
+        fmpz_factor(fx, x);
+        fmpz_factor(fy, y);
+
+        fmpz_factor_mul(fz, fx, fy);
+
+        fmpz_mul(x, x, y);
+        fmpz_factor_expand(y, fz);
+
+        if (fmpz_cmp(x, y) != 0)
+        {
+            flint_printf("FAIL (factor mul)\n");
+            fmpz_print(x); flint_printf(" ");
+            fmpz_print(y); flint_printf("\n");
+            fflush(stdout);
+            flint_abort();
+        }
+    }
+
+    fmpz_clear(x);
+    fmpz_clear(y);
+    fmpz_factor_clear(fx);
+    fmpz_factor_clear(fy);
+    fmpz_factor_clear(fz);
+
+    TEST_FUNCTION_END(state);
+}
diff --git a/src/fmpz_factor/test/t-sort.c b/src/fmpz_factor/test/t-sort.c
new file mode 100644
index 0000000000..4af2953a63
--- /dev/null
+++ b/src/fmpz_factor/test/t-sort.c
@@ -0,0 +1,64 @@
+/*
+    Copyright (C) 2025 Mael Hostettler
+
+    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.h"
+#include "fmpz_factor.h"
+
+TEST_FUNCTION_START(fmpz_factor_sort, state)
+{
+    slong i, j;
+    fmpz_t x, y, z;
+    fmpz_factor_t factor;
+
+    fmpz_init(x);
+    fmpz_init(y);
+    fmpz_init(z);
+
+    for (i = 0; i < 1000 * flint_test_multiplier(); i++)
+    {
+        fmpz_factor_init(factor);
+
+        for(j = 0; j < 10 * flint_test_multiplier(); j++)
+        {
+            fmpz_randprime(x, state, 2 + j, 1);
+            _fmpz_factor_append(factor, x, 1 + n_randint(state, 5));
+        }
+        fmpz_factor_expand(y, factor);
+        fmpz_factor_sort(factor);
+        fmpz_factor_expand(z, factor);
+
+        // check order
+        int ord = 1;
+        for (j = 0; j < factor->num - 1; j++)
+        {
+            ord &= (fmpz_cmp(factor->p + j, factor->p + j + 1) < 0);
+        }
+
+        // check equality
+        if (!fmpz_equal(y, z) || !ord)
+        {
+            flint_printf("FAIL (factor sort)\n");
+            fmpz_factor_print(factor); flint_printf("\n");
+            fmpz_print(y); flint_printf(" ");
+            fmpz_print(z); flint_printf("\n");
+            fflush(stdout);
+            flint_abort();
+        }
+        fmpz_factor_clear(factor);
+    }
+
+    fmpz_clear(x);
+    fmpz_clear(y);
+    fmpz_clear(z);
+
+    TEST_FUNCTION_END(state);
+}