diff --git a/library/bignum_core.c b/library/bignum_core.c
index 88582c2d38..8091694615 100644
--- a/library/bignum_core.c
+++ b/library/bignum_core.c
@@ -18,6 +18,7 @@
 #include "mbedtls/platform.h"
 
 #include "bignum_core.h"
+#include "bignum_core_invasive.h"
 #include "bn_mul.h"
 #include "constant_time_internal.h"
 
@@ -1019,4 +1020,221 @@ 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.
+ */
+MBEDTLS_STATIC_TESTABLE
+void mbedtls_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
+ *
+ * The paper gives two computationally equivalent algorithms: Alg 7 (readable)
+ * and Alg 8 (constant-time). We use a third version that's hopefully both:
+ *
+ *  u, v = A, N  # N is called p in the paper but doesn't have to be prime
+ *  q, r = 0, 1
+ *  repeat bits(A_limbs + N_limbs) times:
+ *      d = v - u  # t1 in Alg 7
+ *      t1 = (u and v both odd) ? u : d  # t1 in Alg 8
+ *      t2 = (u and v both odd) ? d : (u odd) ? v : u  # t2 in Alg 8
+ *      t2 >>= 1
+ *      swap = t1 > t2  # similar to s, z in Alg 8
+ *      u, v = (swap) ? t2, t1 : t1, t2
+ *
+ *      d = r - q mod N  # t2 in Alg 7
+ *      t1 = (u and v both odd) ? q : d  # t3 in Alg 8
+ *      t2 = (u and v both odd) ? d : (u odd) ? r : q  # t4 Alg 8
+ *      t2 /= 2 mod N  # see below (pre_com)
+ *      q, r = (swap) ? t2, t1 : t1, t2
+ *  return v, q  # v: GCD, see Alg 6; q: no mult by pre_com, see below
+ *
+ * The ternary operators in the above pseudo-code need to be realised in a
+ * constant-time fashion. We use conditional assign for t1, t2 and conditional
+ * swap for the final update. (Note: the similarity between branches of Alg 7
+ * are highlighted in tables 2 and 3 and the surrounding text.)
+ *
+ * Also, we re-order operations, grouping things related to the inverse, which
+ * facilitates making its computation optional, and requires fewer temporaries.
+ *
+ * The only actual change from the paper is dropping the trick with pre_com,
+ * which I think complicates things for no benefit.
+ * See the 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)
+{
+    /* 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
+     *
+     * We only write to G (aka v) after reading from inputs (A and N), which
+     * allows aliasing, except with N when I != NULL, as then we'll be operating
+     * mod N on q and r later - see the public documentation.
+     */
+    if (A_limbs > N_limbs) {
+        /* Violating this precondition should not result in memory errors. */
+        A_limbs = N_limbs;
+    }
+    memcpy(u, A, A_limbs * ciL);
+    memset((char *) u + A_limbs * ciL, 0, (N_limbs - A_limbs) * ciL);
+
+    /* Avoid possible UB with memcpy when src == dst. */
+    if (v != N) {
+        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 multiplying
+         * by pre_com, 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);
+
+            mbedtls_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..f044b33f93 100644
--- a/library/bignum_core.h
+++ b/library/bignum_core.h
@@ -822,4 +822,45 @@ 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, 0 <= A <= N and A_limbs <= N_limbs.
+ * When I != NULL, N (the modulus) must be greater than 1.
+ *
+ * A and N may not alias each other.
+ * When I == NULL (computing only the GCD), G may alias A or N.
+ * When I != NULL (computing the modular inverse), G or I may alias A
+ * but none of them may alias N (the modulus).
+ *
+ * If any of the above preconditions is not met, output values are unspecified.
+ *
+ * \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.
+ *                        This value must be odd.
+ *                        If I != NULL this value must be greater than 1.
+ * \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/library/bignum_core_invasive.h b/library/bignum_core_invasive.h
index 167099dc91..a9d447f792 100644
--- a/library/bignum_core_invasive.h
+++ b/library/bignum_core_invasive.h
@@ -13,11 +13,26 @@
 
 #include "bignum_core.h"
 
-#if defined(MBEDTLS_TEST_HOOKS) && !defined(MBEDTLS_THREADING_C)
+#if defined(MBEDTLS_TEST_HOOKS)
+
+#if !defined(MBEDTLS_THREADING_C)
 
 extern void (*mbedtls_safe_codepath_hook)(void);
 extern void (*mbedtls_unsafe_codepath_hook)(void);
 
-#endif /* MBEDTLS_TEST_HOOKS && !MBEDTLS_THREADING_C */
+#endif /* !MBEDTLS_THREADING_C */
+
+/** Divide X by 2 mod N in place, assuming N is odd.
+ *
+ * \param[in,out] X     The value to divide by 2 mod \p N.
+ * \param[in]     N     The modulus. Must be odd.
+ * \param[in]     limbs The number of limbs in \p X and \p N.
+ */
+MBEDTLS_STATIC_TESTABLE
+void mbedtls_mpi_core_div2_mod_odd(mbedtls_mpi_uint *X,
+                                   const mbedtls_mpi_uint *N,
+                                   size_t limbs);
+
+#endif /* MBEDTLS_TEST_HOOKS */
 
 #endif /* MBEDTLS_BIGNUM_CORE_INVASIVE_H */
diff --git a/tests/suites/test_suite_bignum_core.function b/tests/suites/test_suite_bignum_core.function
index c2b44bccdd..759801906f 100644
--- a/tests/suites/test_suite_bignum_core.function
+++ b/tests/suites/test_suite_bignum_core.function
@@ -2,6 +2,7 @@
 #include "mbedtls/bignum.h"
 #include "mbedtls/entropy.h"
 #include "bignum_core.h"
+#include "bignum_core_invasive.h"
 #include "constant_time_internal.h"
 #include "test/constant_flow.h"
 #include "test/bignum_codepath_check.h"
@@ -148,6 +149,29 @@ exit:
     return ret;
 }
 
+/*
+ * Return -1 if A < B, +1 if A > B and 0 if A == B
+ */
+static int mpi_core_cmp(const mbedtls_mpi_uint *A, size_t A_limbs,
+                        const mbedtls_mpi_uint *B, size_t B_limbs)
+{
+    TEST_CF_PUBLIC(A, A_limbs * sizeof(mbedtls_mpi_uint));
+    TEST_CF_PUBLIC(B, B_limbs * sizeof(mbedtls_mpi_uint));
+
+    const mbedtls_mpi AA = {
+        .p = (mbedtls_mpi_uint *) A,
+        .s = 1,
+        .n = (unsigned short) A_limbs,
+    };
+    const mbedtls_mpi BB = {
+        .p = (mbedtls_mpi_uint *) B,
+        .s = 1,
+        .n = (unsigned short) B_limbs,
+    };
+
+    return mbedtls_mpi_cmp_mpi(&AA, &BB);
+}
+
 /* END_HEADER */
 
 /* BEGIN_DEPENDENCIES
@@ -1356,3 +1380,285 @@ 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);
+    TEST_ASSERT(mpi_core_cmp(A, A_limbs, N, N_limbs) <= 0);
+
+    const size_t N_bytes = N_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));
+
+#define FREE_G_I_T      \
+    mbedtls_free(G);    \
+    G = NULL;           \
+    mbedtls_free(I);    \
+    I = NULL;           \
+    mbedtls_free(T);    \
+    T = NULL
+
+    /*
+     * Test GCD only (I == NULL)
+     */
+    TEST_CALLOC(G, N_limbs);
+    memset(G, 'G', N_bytes);
+
+    TEST_CALLOC(T, 4 * N_limbs);
+    memset(T, 'T', 4 * N_bytes);
+
+    mbedtls_mpi_core_gcd_modinv_odd(G, NULL, A, A_limbs, N, N_limbs, T);
+    TEST_EQUAL(mpi_core_cmp(G, N_limbs, exp_G, exp_G_limbs), 0);
+
+    FREE_G_I_T;
+
+    /* GCD only, G aliased to N */
+    TEST_CALLOC(G, N_limbs);
+    memcpy(G, N, N_bytes);
+
+    TEST_CALLOC(T, 4 * N_limbs);
+    memset(T, 'T', 4 * N_bytes);
+
+    mbedtls_mpi_core_gcd_modinv_odd(G, NULL, A, A_limbs, /* N */ G, N_limbs, T);
+    TEST_EQUAL(mpi_core_cmp(G, N_limbs, exp_G, exp_G_limbs), 0);
+
+    FREE_G_I_T;
+
+    /* GCD only, G aliased to A (hence A_limbs = N_limbs) */
+    TEST_CALLOC(G, N_limbs);
+    memcpy(G, A, A_limbs * sizeof(mbedtls_mpi_uint));
+
+    TEST_CALLOC(T, 4 * N_limbs);
+    memset(T, 'T', 4 * N_bytes);
+
+    mbedtls_mpi_core_gcd_modinv_odd(G, NULL, /* A */ G, N_limbs, N, N_limbs, T);
+    TEST_EQUAL(mpi_core_cmp(G, N_limbs, exp_G, exp_G_limbs), 0);
+
+    FREE_G_I_T;
+
+    /*
+     * Test GCD + modinv
+     */
+    TEST_CALLOC(G, N_limbs);
+    memset(G, 'G', N_bytes);
+
+    TEST_CALLOC(I, N_limbs);
+    memset(I, 'I', N_bytes);
+
+    TEST_CALLOC(T, 5 * N_limbs);
+    memset(T, 'T', 5 * N_bytes);
+
+    mbedtls_mpi_core_gcd_modinv_odd(G, I, A, A_limbs, N, N_limbs, T);
+
+    TEST_EQUAL(mpi_core_cmp(G, N_limbs, exp_G, exp_G_limbs), 0);
+    if (got_I) {
+        TEST_EQUAL(mpi_core_cmp(I, N_limbs, exp_I, exp_I_limbs), 0);
+    }
+
+    FREE_G_I_T;
+
+    /* GCD + modinv, G aliased to A */
+    TEST_CALLOC(G, N_limbs);
+    memcpy(G, A, A_limbs * sizeof(mbedtls_mpi_uint));
+
+    TEST_CALLOC(I, N_limbs);
+    memset(I, 'I', N_bytes);
+
+    TEST_CALLOC(T, 5 * N_limbs);
+    memset(T, 'T', 5 * N_bytes);
+
+    mbedtls_mpi_core_gcd_modinv_odd(G, I, /* A */ G, N_limbs, N, N_limbs, T);
+
+    TEST_EQUAL(mpi_core_cmp(G, N_limbs, exp_G, exp_G_limbs), 0);
+    if (got_I) {
+        TEST_EQUAL(mpi_core_cmp(I, N_limbs, exp_I, exp_I_limbs), 0);
+    }
+
+    FREE_G_I_T;
+
+    /* GCD + modinv, I aliased to A */
+    TEST_CALLOC(G, N_limbs);
+    memset(G, 'G', N_bytes);
+
+    TEST_CALLOC(I, N_limbs);
+    memcpy(I, A, A_limbs * sizeof(mbedtls_mpi_uint));
+
+    TEST_CALLOC(T, 5 * N_limbs);
+    memset(T, 'T', 5 * N_bytes);
+
+    mbedtls_mpi_core_gcd_modinv_odd(G, I, /* A */ I, N_limbs, N, N_limbs, T);
+
+    TEST_EQUAL(mpi_core_cmp(G, N_limbs, exp_G, exp_G_limbs), 0);
+    if (got_I) {
+        TEST_EQUAL(mpi_core_cmp(I, N_limbs, exp_I, exp_I_limbs), 0);
+    }
+
+    FREE_G_I_T;
+
+#undef FREE_G_I_T
+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 */
+
+/* BEGIN_CASE */
+void mpi_core_gcd_modinv_odd_preconditions()
+{
+    /*
+     * The purpose of this test function is to ensure that the function doesn't
+     * crash (but just outputs garbage) when preconditions are not met.
+     */
+
+    mbedtls_mpi_uint A[1];
+    mbedtls_mpi_uint N[2];
+    mbedtls_mpi_uint AAA[3];
+    mbedtls_mpi_uint *G = NULL, *I = NULL, *TG = NULL, *TI = NULL;
+
+    /* We'll always use a two-limbs N */
+    TEST_CALLOC(G, 2);
+    TEST_CALLOC(I, 2);
+    TEST_CALLOC(TG, 4 * 2);  // For I == NULL
+    TEST_CALLOC(TI, 5 * 2);  // For I != NULL
+
+    /*
+     * Input values
+     */
+
+    /* N is not odd */
+    N[0] = 2;       // N = 2^n + 2
+    N[1] = 1;
+    A[0] = 42;      // A = 42
+    mbedtls_mpi_core_gcd_modinv_odd(G, NULL, A, 1, N, 2, TG);
+    mbedtls_mpi_core_gcd_modinv_odd(G, I, A, 1, N, 2, TI);
+
+    /* A > N */
+    N[0] = 3;       // N = 3
+    N[1] = 0;
+    A[0] = 42;      // A = 42
+    mbedtls_mpi_core_gcd_modinv_odd(G, NULL, A, 1, N, 2, TG);
+    mbedtls_mpi_core_gcd_modinv_odd(G, I, A, 1, N, 2, TI);
+
+    /* A_limbs > N_limbs (but A <= N) */
+    N[0] = 3;       // N = 3
+    N[1] = 0;
+    AAA[0] = 1;     // A = 1
+    AAA[1] = 0;
+    AAA[2] = 0;
+    mbedtls_mpi_core_gcd_modinv_odd(G, NULL, AAA, 3, N, 2, TG);
+    mbedtls_mpi_core_gcd_modinv_odd(G, I, AAA, 3, N, 2, TI);
+
+    /* A_limbs > N_limbs (and A > N) */
+    N[0] = 3;       // N = 3
+    N[1] = 0;
+    AAA[0] = 0;     // A = 2^2n
+    AAA[1] = 0;
+    AAA[2] = 1;
+    mbedtls_mpi_core_gcd_modinv_odd(G, NULL, AAA, 3, N, 2, TG);
+    mbedtls_mpi_core_gcd_modinv_odd(G, I, AAA, 3, N, 2, TI);
+
+    /* I != NULL but N is 1 */
+    N[0] = 1;       // N = 1
+    N[1] = 0;
+    A[0] = 1;       // A = 1
+    mbedtls_mpi_core_gcd_modinv_odd(G, I, N, 2, A, 1, TI);
+
+    /*
+     * Aliasing
+     */
+
+    /* Now use valid values for remaining tests */
+    N[0] = 1;       // N = 2^n + 1
+    N[1] = 1;
+    A[0] = 42;      // A = 42
+
+    /* A aliased to N */
+    mbedtls_mpi_core_gcd_modinv_odd(G, NULL, N, 2, N, 2, TG);
+    mbedtls_mpi_core_gcd_modinv_odd(G, I, N, 2, N, 2, TI);
+
+    /* G aliased to A and N */
+    memcpy(G, N, 2 * sizeof(mbedtls_mpi_uint));
+    mbedtls_mpi_core_gcd_modinv_odd(G, NULL, G, 2, G, 2, TG);
+    mbedtls_mpi_core_gcd_modinv_odd(G, I, G, 2, G, 2, TI);
+
+    /* I != NULL, G aliased to N */
+    memcpy(G, N, 2 * sizeof(mbedtls_mpi_uint));
+    mbedtls_mpi_core_gcd_modinv_odd(G, I, A, 1, G, 2, TI);
+
+    /* I != NULL, I aliased to N */
+    memcpy(I, N, 2 * sizeof(mbedtls_mpi_uint));
+    mbedtls_mpi_core_gcd_modinv_odd(G, I, A, 1, I, 2, TI);
+
+    /* I aliased to A and N */
+    memcpy(I, N, 2 * sizeof(mbedtls_mpi_uint));
+    mbedtls_mpi_core_gcd_modinv_odd(G, I, I, 2, I, 2, TI);
+
+exit:
+    mbedtls_free(G);
+    mbedtls_free(I);
+    mbedtls_free(TG);
+    mbedtls_free(TI);
+}
+/* END_CASE */
+
+/* BEGIN_CASE depends_on:MBEDTLS_TEST_HOOKS */
+void mpi_core_div2_mod_odd(char *input_X, char *input_N, char *input_exp_X)
+{
+    mbedtls_mpi_uint *X = NULL;
+    size_t X_limbs = 0;
+    mbedtls_mpi_uint *N = NULL;
+    size_t N_limbs = 0;
+    mbedtls_mpi_uint *exp_X = NULL;
+    size_t exp_X_limbs = 0;
+
+    /* Read test parameters into MPI structures */
+    TEST_EQUAL(0, mbedtls_test_read_mpi_core(&X, &X_limbs, input_X));
+    TEST_EQUAL(0, mbedtls_test_read_mpi_core(&N, &N_limbs, input_N));
+    TEST_EQUAL(0, mbedtls_test_read_mpi_core(&exp_X, &exp_X_limbs, input_exp_X));
+
+    /* The function under test requires this */
+    TEST_EQUAL(X_limbs, N_limbs);
+
+    TEST_CF_SECRET(X, X_limbs * sizeof(mbedtls_mpi_uint));
+    TEST_CF_SECRET(N, N_limbs * sizeof(mbedtls_mpi_uint));
+
+    mbedtls_mpi_core_div2_mod_odd(X, N, N_limbs);
+
+    TEST_CF_PUBLIC(X, X_limbs * sizeof(mbedtls_mpi_uint));
+    TEST_EQUAL(0, mpi_core_cmp(X, X_limbs, exp_X, exp_X_limbs));
+
+exit:
+    mbedtls_free(X);
+    mbedtls_free(N);
+    mbedtls_free(exp_X);
+}
+/* END_CASE */
diff --git a/tests/suites/test_suite_bignum_core.misc.data b/tests/suites/test_suite_bignum_core.misc.data
index ba86029977..e2ba5a61f7 100644
--- a/tests/suites/test_suite_bignum_core.misc.data
+++ b/tests/suites/test_suite_bignum_core.misc.data
@@ -523,3 +523,29 @@ 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":"3":""
+
+GCD-modinv random 80-bit, trivial GCD -> inverse
+mpi_core_gcd_modinv_odd:"7f2405d6de7db80a7bc":"1a84113636607520200d":"1":"15f158844a59cd7a3ed2"
+
+# This data results in the gcd-modinv loop converging to its final state
+# only in the last iteration. See python script in commit message.
+GCD-modinv (almost) max iterations
+mpi_core_gcd_modinv_odd:"8000000000000000":"b26eb5721a2cb24c36acb4550b176671":"1":"77e1dd63583a6b3c8deefe7737862c89"
+
+GCD-modinv preconditions not met
+mpi_core_gcd_modinv_odd_preconditions:
+
+Div2 mod odd: even value
+mpi_core_div2_mod_odd:"4":"7":"2"
+
+Div2 mod odd: odd value, no carry
+mpi_core_div2_mod_odd:"5":"7":"6"
+
+Div2 mod odd: odd value with carry
+mpi_core_div2_mod_odd:"8000000000000001":"8000000000000003":"8000000000000002"
+
+Div2 mod odd: even value with top bit set
+mpi_core_div2_mod_odd:"8000000000000002":"8000000000000003":"4000000000000001"