diff --git a/library/bignum_core.c b/library/bignum_core.c
index 88582c2d38..1e80a5bf5d 100644
--- a/library/bignum_core.c
+++ b/library/bignum_core.c
@@ -1019,4 +1019,189 @@ void mbedtls_mpi_core_from_mont_rep(mbedtls_mpi_uint *X,
     mbedtls_mpi_core_montmul(X, A, &Rinv, 1, N, AN_limbs, mm, T);
 }
 
+/*
+ * Compute X = A - B mod N.
+ * Both A and B must be in [0, N) and so will the output.
+ */
+static void mpi_core_sub_mod(mbedtls_mpi_uint *X,
+                             const mbedtls_mpi_uint *A,
+                             const mbedtls_mpi_uint *B,
+                             const mbedtls_mpi_uint *N,
+                             size_t limbs)
+{
+    mbedtls_mpi_uint c = mbedtls_mpi_core_sub(X, A, B, limbs);
+    (void) mbedtls_mpi_core_add_if(X, N, limbs, (unsigned) c);
+}
+
+/*
+ * Divide X by 2 mod N in place, assuming N is odd.
+ * The input must be in [0, N) and so will the output.
+ */
+static void mpi_core_div2_mod_odd(mbedtls_mpi_uint *X,
+                                  const mbedtls_mpi_uint *N,
+                                  size_t limbs)
+{
+    /* If X is odd, add N to make it even before shifting. */
+    unsigned odd = (unsigned) X[0] & 1;
+    mbedtls_mpi_uint c = mbedtls_mpi_core_add_if(X, N, limbs, odd);
+    mbedtls_mpi_core_shift_r(X, limbs, 1);
+    X[limbs - 1] |= c << (biL - 1);
+}
+
+/*
+ * Constant-time GCD and modular inversion - odd modulus.
+ *
+ * Pre-conditions: see public documentation.
+ *
+ * See https://www.jstage.jst.go.jp/article/transinf/E106.D/9/E106.D_2022ICP0009/_pdf
+ * This is an adaptation of Alg 7 / Alg 8:
+ * - Alg 7 is readable but not constant-time, Alg 8 is constant-time but not
+ *   readable (and uses signed arithmetic). We mostly follow Alg 7 and make it
+ *   constant-time by using our usual primitives (conditional assign,
+ *   conditional swap) rather than re-inventing them. We only take a few
+ *   notations from Alg 8 for temporaries.
+ * - Compared to both, we skip the trick with pre_comm: I think this trick
+ *   complicates things for no benefit (see comment on the big I != NULL block
+ *   below for details).
+ */
+void mbedtls_mpi_core_gcd_modinv_odd(mbedtls_mpi_uint *G,
+                                     mbedtls_mpi_uint *I,
+                                     const mbedtls_mpi_uint *A,
+                                     size_t A_limbs,
+                                     const mbedtls_mpi_uint *N,
+                                     size_t N_limbs,
+                                     mbedtls_mpi_uint *T)
+{
+    /* Note: N is called p in the paper, but doesn't need to be prime, only odd.
+     */
+
+    /* GCD and modinv, names common to Alg 7 and Alg 8 */
+    mbedtls_mpi_uint *u = T + 0 * N_limbs;
+    mbedtls_mpi_uint *v = G;
+
+    /* GCD and modinv, my name (t1, t2 from Alg 7) */
+    mbedtls_mpi_uint *d = T + 1 * N_limbs;
+
+    /* GCD and modinv, names from Alg 8 (note: t1, t2 from Alg 7 are d above) */
+    mbedtls_mpi_uint *t1 = T + 2 * N_limbs;
+    mbedtls_mpi_uint *t2 = T + 3 * N_limbs;
+
+    /* modinv only, names common to Alg 7 and Alg 8 */
+    mbedtls_mpi_uint *q = I;
+    mbedtls_mpi_uint *r = I != NULL ? T + 4 * N_limbs : NULL;
+
+    /*
+     * Initial values:
+     * u, v = A, N
+     * q, r = 0, 1
+     */
+    memcpy(u, A, A_limbs * ciL);
+    memset((char *) u + A_limbs * ciL, 0, (N_limbs - A_limbs) * ciL);
+
+    memcpy(v, N, N_limbs * ciL);
+
+    if (I != NULL) {
+        memset(q, 0, N_limbs * ciL);
+
+        memset(r, 0, N_limbs * ciL);
+        r[0] = 1;
+    }
+
+    /*
+     * At each step, out of u, v, v - u we keep one, shift another, and discard
+     * the third, then update (u, v) with the ordered result.
+     * Then we mirror those actions with q, r, r - q mod N.
+     *
+     * Loop invariants:
+     *  u <= v                  (on entry: A <= N)
+     *  GCD(u, v) == GCD(A, N)  (on entry: trivial)
+     *  v = A * q mod N         (on entry: N = A * 0 mod N)
+     *  u = A * r mod N         (on entry: A = A * 1 mod N)
+     *  q, r in [0, N)          (on entry: 0, 1)
+     *
+     * On exit:
+     *  u = 0
+     *  v = GCD(A, N) = A * q mod N
+     *  if v == 1 then 1 = A * q mod N ie q is A's inverse mod N
+     *  r = 0
+     *
+     * The exit state is a fixed point of the loop's body.
+     * Alg 7 and Alg 8 use 2 * bitlen(N) iterations but Theorem 2 (above in the
+     * paper) says bitlen(A) + bitlen(N) is actually enough.
+     */
+    for (size_t i = 0; i < (A_limbs + N_limbs) * biL; i++) {
+        /* s, z in Alg 8 - use meaningful names instead */
+        mbedtls_ct_condition_t u_odd = mbedtls_ct_bool(u[0] & 1);
+        mbedtls_ct_condition_t v_odd = mbedtls_ct_bool(v[0] & 1);
+
+        /* Other conditions that will be useful below */
+        mbedtls_ct_condition_t u_odd_v_odd = mbedtls_ct_bool_and(u_odd, v_odd);
+        mbedtls_ct_condition_t v_even = mbedtls_ct_bool_not(v_odd);
+        mbedtls_ct_condition_t u_odd_v_even = mbedtls_ct_bool_and(u_odd, v_even);
+
+        /* This is called t1 in Alg 7 (no name in Alg 8).
+         * We know that u <= v so there is no carry */
+        (void) mbedtls_mpi_core_sub(d, v, u, N_limbs);
+
+        /* t1 (the thing that's kept) can be d (default) or u (if t2 is d) */
+        memcpy(t1, d, N_limbs * ciL);
+        mbedtls_mpi_core_cond_assign(t1, u, N_limbs, u_odd_v_odd);
+
+        /* t2 (the thing that's shifted) can be u (if even), or v (if even),
+         * or d (which is even if both u and v were odd) */
+        memcpy(t2, u, N_limbs * ciL);
+        mbedtls_mpi_core_cond_assign(t2, v, N_limbs, u_odd_v_even);
+        mbedtls_mpi_core_cond_assign(t2, d, N_limbs, u_odd_v_odd);
+
+        mbedtls_mpi_core_shift_r(t2, N_limbs, 1); // t2 is even
+
+        /* Update u, v and re-order them if needed */
+        memcpy(u, t1, N_limbs * ciL);
+        memcpy(v, t2, N_limbs * ciL);
+        mbedtls_ct_condition_t swap = mbedtls_mpi_core_lt_ct(v, u, N_limbs);
+        mbedtls_mpi_core_cond_swap(u, v, N_limbs, swap);
+
+        /* Now, if modinv was requested, do the same with q, r, but:
+         * - decisions still based on u and v (their initial values);
+         * - operations are now mod N;
+         * - we re-use t1, t2 for what the paper calls t3, t4 in Alg 8.
+         *
+         * Here we slightly diverge from the paper and instead do the obvious
+         * thing that preserves the invariants involving q and r: mirror
+         * operations on u and v, ie also divide by 2 here (mod N).
+         *
+         * The paper uses a trick where it replaces division by 2 with
+         * multiplication by 2 here, and compensates in the end by doing a
+         * final multiplication, which is probably intended as an optimisation.
+         *
+         * However I believe it's not actually an optimisation, since
+         * constant-time modular multiplication by 2 (left-shift + conditional
+         * subtract) is just as costly as constant-time modular division by 2
+         * (conditional add + right-shift). So, skip it and keep things simple.
+         */
+        if (I != NULL) {
+            /* This is called t2 in Alg 7 (no name in Alg 8). */
+            mpi_core_sub_mod(d, q, r, N, N_limbs);
+
+            /* t3 (the thing that's kept) */
+            memcpy(t1, d, N_limbs * ciL);
+            mbedtls_mpi_core_cond_assign(t1, r, N_limbs, u_odd_v_odd);
+
+            /* t4 (the thing that's shifted) */
+            memcpy(t2, r, N_limbs * ciL);
+            mbedtls_mpi_core_cond_assign(t2, q, N_limbs, u_odd_v_even);
+            mbedtls_mpi_core_cond_assign(t2, d, N_limbs, u_odd_v_odd);
+
+            mpi_core_div2_mod_odd(t2, N, N_limbs);
+
+            /* Update and possibly swap */
+            memcpy(r, t1, N_limbs * ciL);
+            memcpy(q, t2, N_limbs * ciL);
+            mbedtls_mpi_core_cond_swap(r, q, N_limbs, swap);
+        }
+    }
+
+    /* G and I already hold the correct values by virtue of being aliased */
+}
+
 #endif /* MBEDTLS_BIGNUM_C */
diff --git a/library/bignum_core.h b/library/bignum_core.h
index 264ee63550..29e05cd3c9 100644
--- a/library/bignum_core.h
+++ b/library/bignum_core.h
@@ -822,4 +822,38 @@ void mbedtls_mpi_core_from_mont_rep(mbedtls_mpi_uint *X,
                                     mbedtls_mpi_uint mm,
                                     mbedtls_mpi_uint *T);
 
+/** Compute GCD(A, N) and optionally the inverse of A mod N if it exists.
+ *
+ * Requires N to be odd, and 0 <= A <= N.
+ *
+ * None of the parameters may alias or overlap another.
+ *
+ * \param[out]    G       The GCD of \p A and \p N.
+ *                        Must have the same number of limbs as \p N.
+ * \param[out]    I       The inverse of \p A modulo \p N if it exists (that is,
+ *                        if \p G above is 1 on exit); indeterminate otherwise.
+ *                        This must either be NULL (to only compute the GCD),
+ *                        or have the same number of limbs as \p N.
+ * \param[in]     A       The 1st operand of GCD and number to invert.
+ *                        This value must be less than or equal to \p N.
+ * \param         A_limbs The number of limbs of \p A.
+ *                        Must be less than or equal to \p N_limbs.
+ * \param[in]     N       The 2nd operand of GCD and modulus for inversion.
+ *                        Must be odd or the results are indeterminate.
+ * \param         N_limbs The number of limbs of \p N.
+ * \param[in,out] T       Temporary storage of size at least 5 * N_limbs limbs,
+ *                        or 4 * N_limbs if \p I is NULL (GCD only).
+ *                        Its initial content is unused and
+ *                        its final content is indeterminate.
+ *                        It must not alias or otherwise overlap any of the
+ *                        other parameters.
+ */
+void mbedtls_mpi_core_gcd_modinv_odd(mbedtls_mpi_uint *G,
+                                     mbedtls_mpi_uint *I,
+                                     const mbedtls_mpi_uint *A,
+                                     size_t A_limbs,
+                                     const mbedtls_mpi_uint *N,
+                                     size_t N_limbs,
+                                     mbedtls_mpi_uint *T);
+
 #endif /* MBEDTLS_BIGNUM_CORE_H */
diff --git a/tests/suites/test_suite_bignum_core.function b/tests/suites/test_suite_bignum_core.function
index c2b44bccdd..81c2242997 100644
--- a/tests/suites/test_suite_bignum_core.function
+++ b/tests/suites/test_suite_bignum_core.function
@@ -1356,3 +1356,94 @@ exit:
     mbedtls_free(X);
 }
 /* END_CASE */
+
+/* BEGIN_CASE */
+void mpi_core_gcd_modinv_odd(char *input_A, char *input_N,
+                             char *input_exp_G, char *input_exp_I)
+{
+    mbedtls_mpi_uint *A = NULL;
+    size_t A_limbs = 0;
+    mbedtls_mpi_uint *N = NULL;
+    size_t N_limbs = 0;
+    mbedtls_mpi_uint *exp_G = NULL;
+    size_t exp_G_limbs = 0;
+    mbedtls_mpi_uint *exp_I = NULL;
+    size_t exp_I_limbs = 0;
+    mbedtls_mpi_uint *G = NULL, *I = NULL, *T = NULL;
+
+    /* Read test parameters into MPI structures */
+    TEST_EQUAL(0, mbedtls_test_read_mpi_core(&A, &A_limbs, input_A));
+    TEST_EQUAL(0, mbedtls_test_read_mpi_core(&N, &N_limbs, input_N));
+    TEST_EQUAL(0, mbedtls_test_read_mpi_core(&exp_G, &exp_G_limbs, input_exp_G));
+    const unsigned char got_I = strlen(input_exp_I) != 0;
+    if (got_I) {
+        TEST_EQUAL(0, mbedtls_test_read_mpi_core(&exp_I, &exp_I_limbs, input_exp_I));
+    }
+
+    /* The function under test wants this */
+    TEST_EQUAL(N[0] & 1, 1);
+    TEST_LE_U(A_limbs, N_limbs);
+    if (A_limbs == N_limbs) {
+        TEST_EQUAL(mbedtls_mpi_core_lt_ct(N, A, N_limbs), MBEDTLS_CT_FALSE);
+    }
+    /* Other things we want from test data, for our convenience */
+    TEST_EQUAL(exp_G_limbs, N_limbs);
+    if (got_I) {
+        TEST_EQUAL(exp_I_limbs, N_limbs);
+    }
+
+    const size_t limbs = N_limbs;
+    const size_t bytes = limbs * sizeof(mbedtls_mpi_uint);
+
+    TEST_CF_SECRET(A, A_limbs * sizeof(mbedtls_mpi_uint));
+    TEST_CF_SECRET(N, N_limbs * sizeof(mbedtls_mpi_uint));
+
+    /*
+     * Test GCD only (I == NULL)
+     */
+    TEST_CALLOC(G, N_limbs);
+    memset(G, 'G', bytes);
+
+    TEST_CALLOC(T, 4 * N_limbs);
+    memset(T, 'T', bytes);
+
+    mbedtls_mpi_core_gcd_modinv_odd(G, NULL, A, A_limbs, N, N_limbs, T);
+    TEST_CF_PUBLIC(G, bytes);
+    TEST_MEMORY_COMPARE(G, bytes, exp_G, bytes);
+
+    mbedtls_free(G);
+    G = NULL;
+    mbedtls_free(T);
+    T = NULL;
+
+    /*
+     * Test GCD + modinv
+     */
+    TEST_CALLOC(G, N_limbs);
+    memset(G, 'G', bytes);
+
+    TEST_CALLOC(I, N_limbs);
+    memset(I, 'I', bytes);
+
+    TEST_CALLOC(T, 5 * N_limbs);
+    memset(T, 'T', bytes);
+
+    mbedtls_mpi_core_gcd_modinv_odd(G, I, A, A_limbs, N, N_limbs, T);
+
+    TEST_CF_PUBLIC(G, bytes);
+    TEST_MEMORY_COMPARE(G, bytes, exp_G, bytes);
+    if (got_I) {
+        TEST_CF_PUBLIC(I, bytes);
+        TEST_MEMORY_COMPARE(I, bytes, exp_I, bytes);
+    }
+
+exit:
+    mbedtls_free(A);
+    mbedtls_free(N);
+    mbedtls_free(exp_G);
+    mbedtls_free(exp_I);
+    mbedtls_free(G);
+    mbedtls_free(I);
+    mbedtls_free(T);
+}
+/* END_CASE */
diff --git a/tests/suites/test_suite_bignum_core.misc.data b/tests/suites/test_suite_bignum_core.misc.data
index ba86029977..118ca07e65 100644
--- a/tests/suites/test_suite_bignum_core.misc.data
+++ b/tests/suites/test_suite_bignum_core.misc.data
@@ -523,3 +523,9 @@ mpi_core_clz:64:0
 
 CLZ: 100000 0: skip overly long input
 mpi_core_clz:100000:0
+
+GCD-modinv random 80-bit, non-trivial GCD -> no inverse
+mpi_core_gcd_modinv_odd:"e4518a1900fce698fa3":"1a84113636607520200d":"00000000000000000003":""
+
+GCD-modinv random 80-bit, trivial GCD -> inverse
+mpi_core_gcd_modinv_odd:"7f2405d6de7db80a7bc":"1a84113636607520200d":"00000000000000000001":"15f158844a59cd7a3ed2"