diff --git a/tests/api/test_wolfmath.c b/tests/api/test_wolfmath.c
index e5d1eebf934..bfe30df229a 100644
--- a/tests/api/test_wolfmath.c
+++ b/tests/api/test_wolfmath.c
@@ -217,6 +217,69 @@ int test_mp_exptmod_even_2k(void)
     return EXPECT_RESULT();
 }
 
+/*
+ * Testing TFM exponentiation with a modulus of FP_SIZE / 2 digits - the largest
+ * accepted. Constant time exponentiation keeps working values at the length of
+ * the modulus so products use all FP_SIZE digits.
+ */
+int test_mp_exptmod_fp_size_boundary(void)
+{
+    EXPECT_DECLS;
+#if defined(USE_FAST_MATH) && defined(WOLFSSL_PUBLIC_MP) && \
+    defined(TFM_TIMING_RESISTANT) && !defined(WC_NO_MODEXP_TIMING_RESISTANT)
+    mp_int g, x, p, y, e;
+    int k = (FP_SIZE / 2) * DIGIT_BIT;
+
+    XMEMSET(&g, 0, sizeof(g));
+    XMEMSET(&x, 0, sizeof(x));
+    XMEMSET(&p, 0, sizeof(p));
+    XMEMSET(&y, 0, sizeof(y));
+    XMEMSET(&e, 0, sizeof(e));
+    ExpectIntEQ(mp_init_multi(&g, &x, &p, &y, &e, NULL), MP_OKAY);
+
+    /* p = 2^k - 1: odd and all FP_SIZE / 2 digits set. */
+    ExpectIntEQ(mp_set(&p, 1), MP_OKAY);
+    ExpectIntEQ(mp_mul_2d(&p, k, &p), MP_OKAY);
+    ExpectIntEQ(mp_sub_d(&p, 1, &p), MP_OKAY);
+    ExpectIntEQ(p.used, FP_SIZE / 2);
+
+    /* Base 2: 2^(k.2^64 + 7) mod (2^k - 1) = 2^7 as 2^k = 1 mod p. */
+    ExpectIntEQ(mp_set(&g, 2), MP_OKAY);
+    ExpectIntEQ(mp_set_int(&x, (unsigned long)k), MP_OKAY);
+    ExpectIntEQ(mp_mul_2d(&x, 64, &x), MP_OKAY);
+    ExpectIntEQ(mp_add_d(&x, 7, &x), MP_OKAY);
+    ExpectIntEQ(mp_exptmod(&g, &x, &p, &y), MP_OKAY);
+    ExpectIntEQ(mp_cmp_d(&y, 128), MP_EQ);
+
+    /* General base: (p - 1)^x mod p = p - 1 when x odd and 1 when x even.
+     * Full width base and exponent. */
+    ExpectIntEQ(mp_sub_d(&p, 1, &g), MP_OKAY);
+    ExpectIntEQ(mp_sub_d(&p, 2, &x), MP_OKAY);
+    ExpectIntEQ(mp_exptmod(&g, &x, &p, &y), MP_OKAY);
+    ExpectIntEQ(mp_cmp(&y, &g), MP_EQ);
+    ExpectIntEQ(mp_sub_d(&p, 3, &x), MP_OKAY);
+    ExpectIntEQ(mp_exptmod(&g, &x, &p, &y), MP_OKAY);
+    ExpectIntEQ(mp_cmp_d(&y, 1), MP_EQ);
+
+    /* Exponent processed as a fixed number of digits. */
+    ExpectIntEQ(mp_sub_d(&p, 2, &x), MP_OKAY);
+    ExpectIntEQ(mp_exptmod_ex(&g, &x, x.used, &p, &y), MP_OKAY);
+    ExpectIntEQ(mp_cmp(&y, &g), MP_EQ);
+
+    /* One digit more than the boundary is rejected. */
+    ExpectIntEQ(mp_mul_2d(&p, DIGIT_BIT, &e), MP_OKAY);
+    ExpectIntEQ(mp_add_d(&e, 1, &e), MP_OKAY);
+    ExpectIntEQ(mp_exptmod(&g, &x, &e, &y), WC_NO_ERR_TRACE(MP_VAL));
+
+    mp_clear(&g);
+    mp_clear(&x);
+    mp_clear(&p);
+    mp_clear(&y);
+    mp_clear(&e);
+#endif
+    return EXPECT_RESULT();
+}
+
 /*
  * Testing wc_export_int
  */
diff --git a/tests/api/test_wolfmath.h b/tests/api/test_wolfmath.h
index 39de04f040d..1cd580d072d 100644
--- a/tests/api/test_wolfmath.h
+++ b/tests/api/test_wolfmath.h
@@ -30,6 +30,7 @@ int test_mp_get_rand_digit(void);
 int test_mp_cond_copy(void);
 int test_mp_rand(void);
 int test_mp_exptmod_even_2k(void);
+int test_mp_exptmod_fp_size_boundary(void);
 int test_wc_export_int(void);
 int test_wc_SpIntSizeDecisionCoverage(void);
 int test_wc_SpIntShiftDecisionCoverage(void);
@@ -45,6 +46,7 @@ int test_wc_SpIntExptGcdDecisionCoverage(void);
     TEST_DECL_GROUP("wolfmath", test_mp_cond_copy),     \
     TEST_DECL_GROUP("wolfmath", test_mp_rand),          \
     TEST_DECL_GROUP("wolfmath", test_mp_exptmod_even_2k), \
+    TEST_DECL_GROUP("wolfmath", test_mp_exptmod_fp_size_boundary), \
     TEST_DECL_GROUP("wolfmath", test_wc_export_int),    \
     TEST_DECL_GROUP("wolfmath", test_wc_SpIntSizeDecisionCoverage),      \
     TEST_DECL_GROUP("wolfmath", test_wc_SpIntShiftDecisionCoverage),     \
diff --git a/tests/unit-mcdc/test_ecc_fault_whitebox.c b/tests/unit-mcdc/test_ecc_fault_whitebox.c
index 44f7f1676e5..91a148f053c 100644
--- a/tests/unit-mcdc/test_ecc_fault_whitebox.c
+++ b/tests/unit-mcdc/test_ecc_fault_whitebox.c
@@ -794,38 +794,38 @@ static void wb_degenerate_points(void)
         wb_fail = 1;
         goto out;
     }
-    (void)_ecc_projective_add_point(P, Q, R, &a, &prime, mp);
+    (void)_ecc_projective_add_point(P, Q, R, &a, &prime, mp, 0);
 
     /* Q = -P: x and z equal, y differs -> the general add gives Z == 0 */
     (void)mp_sub(&prime, P->y, Q->y);
-    (void)_ecc_projective_add_point(P, Q, R, &a, &prime, mp);
+    (void)_ecc_projective_add_point(P, Q, R, &a, &prime, mp, 0);
 
     /* Q->z == 0: the digit-count operand FALSE with x still equal */
     (void)wc_ecc_copy_point(P, Q);
     (void)mp_set(Q->z, 0);
-    (void)_ecc_projective_add_point(P, Q, R, &a, &prime, mp);
+    (void)_ecc_projective_add_point(P, Q, R, &a, &prime, mp, 0);
 
     /* Q->z differs (non-zero): the z comparison FALSE */
     (void)mp_set(Q->z, 2);
-    (void)_ecc_projective_add_point(P, Q, R, &a, &prime, mp);
+    (void)_ecc_projective_add_point(P, Q, R, &a, &prime, mp, 0);
 
     /* Q->x differs: the first operand FALSE */
     (void)wc_ecc_copy_point(P, Q);
     (void)mp_set(Q->x, 3);
-    (void)_ecc_projective_add_point(P, Q, R, &a, &prime, mp);
+    (void)_ecc_projective_add_point(P, Q, R, &a, &prime, mp, 0);
 
     /* y differs and Q is not -P either: the last operand FALSE */
     (void)wc_ecc_copy_point(P, Q);
     (void)mp_set(Q->y, 5);
-    (void)_ecc_projective_add_point(P, Q, R, &a, &prime, mp);
+    (void)_ecc_projective_add_point(P, Q, R, &a, &prime, mp, 0);
 
     /* curve coefficient zero -> "modulus - a is zero" TRUE in the doubling
      * formula's coefficient selection; the real a gives the FALSE row. */
-    (void)_ecc_projective_dbl_point(P, R, &zero, &prime, mp);
-    (void)_ecc_projective_dbl_point(P, R, &a, &prime, mp);
+    (void)_ecc_projective_dbl_point(P, R, &zero, &prime, mp, 0);
+    (void)_ecc_projective_dbl_point(P, R, &a, &prime, mp, 0);
     for (n = 1; n <= 10; n++) {
         mcdc_fm_arm(n);
-        (void)_ecc_projective_dbl_point(P, R, &zero, &prime, mp);
+        (void)_ecc_projective_dbl_point(P, R, &zero, &prime, mp, 0);
         mcdc_fm_disarm();
     }
 
diff --git a/wolfcrypt/src/ecc.c b/wolfcrypt/src/ecc.c
index 0c060d34517..c1e853b45cf 100644
--- a/wolfcrypt/src/ecc.c
+++ b/wolfcrypt/src/ecc.c
@@ -2084,15 +2084,117 @@ static void alt_fp_init(mp_int* a)
 }
 #endif /* ALT_ECC_SIZE */
 
+/* Constant time Montgomery multiplication available from math library. */
+#if defined(USE_FAST_MATH) || defined(WOLFSSL_SP_MATH_ALL) || \
+    defined(USE_INTEGER_HEAP_MATH)
+    #define ECC_HAVE_MONT_MUL_CT
+#endif
+
+#ifdef USE_INTEGER_HEAP_MATH
+/* Heap math has no constant time Montgomery multiplication. */
+WC_MAYBE_UNUSED static int mp_mont_mul_ct(mp_int* a, mp_int* b, mp_int* m,
+    mp_digit mp, mp_int* r)
+{
+    int err = mp_mul(a, b, r);
+    if (err == MP_OKAY)
+        err = mp_montgomery_reduce(r, m, mp);
+    return err;
+}
+#endif
+
 
 #if !defined(WOLFSSL_ATECC508A) && !defined(WOLFSSL_ATECC608A) && \
     !defined(WOLFSSL_NO_ECC_SW) && \
     (!defined(WOLF_CRYPTO_CB_ONLY_ECC) || defined(WOLFSSL_QNX_CAAM) || \
       defined(WOLFSSL_IMXRT1170_CAAM))
 
+/* Montgomery reduction used by the point arithmetic.
+ * With timing resistance, the ordinates being reduced may be derived from a
+ * private scalar (or the blinding value), so use the constant-time
+ * implementation when ct is set. The variable-time version clamps the result
+ * and only conditionally subtracts the modulus.
+ * ct is only set to 0 for operations on public values, such as verification.
+ */
+#ifdef ECC_TIMING_RESISTANT
+    #define ECC_MONT_RED(a, m, mp, ct)                         \
+        ((ct) ? mp_montgomery_reduce_ct(a, m, mp) :            \
+                mp_montgomery_reduce(a, m, mp))
+#else
+    #define ECC_MONT_RED(a, m, mp, ct)  mp_montgomery_reduce(a, m, mp)
+#endif
+
+#if !defined(WOLFSSL_SP_MATH) || defined(WOLFSSL_PUBLIC_ECC_ADD_DBL)
+/* Montgomery multiplication used by the point arithmetic: r = a * b / R mod m
+ *
+ * When ct is set, the operands are padded to the length of the modulus and
+ * the product is not clamped so that the operations performed don't depend on
+ * the number of leading zero digits. a and b must be less than m.
+ *
+ * @param [in]  a   First operand.
+ * @param [in]  b   Second operand.
+ * @param [in]  m   Modulus.
+ * @param [in]  mp  Montgomery multiplier for m.
+ * @param [out] r   Result. May be the same as a or b.
+ * @param [in]  ct  Use constant time operations - values may be secret.
+ * @return  MP_OKAY on success.
+ */
+WC_MAYBE_UNUSED static int ecc_mont_mul(mp_int* a, mp_int* b, mp_int* m,
+    mp_digit mp, mp_int* r, int ct)
+{
+    int err;
+
+    (void)ct;
+#if defined(ECC_TIMING_RESISTANT) && defined(ECC_HAVE_MONT_MUL_CT)
+    if (ct) {
+        err = mp_mont_mul_ct(a, b, m, mp, r);
+    }
+    else
+#endif
+    {
+        err = mp_mul(a, b, r);
+        if (err == MP_OKAY)
+            err = ECC_MONT_RED(r, m, mp, ct);
+    }
+
+    return err;
+}
+
+/* Montgomery squaring used by the point arithmetic: r = a * a / R mod m
+ *
+ * See ecc_mont_mul().
+ *
+ * @param [in]  a   Operand.
+ * @param [in]  m   Modulus.
+ * @param [in]  mp  Montgomery multiplier for m.
+ * @param [out] r   Result. May be the same as a.
+ * @param [in]  ct  Use constant time operations - values may be secret.
+ * @return  MP_OKAY on success.
+ */
+WC_MAYBE_UNUSED static int ecc_mont_sqr(mp_int* a, mp_int* m, mp_digit mp,
+    mp_int* r, int ct)
+{
+    int err;
+
+    (void)ct;
+#if defined(ECC_TIMING_RESISTANT) && defined(ECC_HAVE_MONT_MUL_CT)
+    if (ct) {
+        err = mp_mont_mul_ct(a, a, m, mp, r);
+    }
+    else
+#endif
+    {
+        err = mp_sqr(a, r);
+        if (err == MP_OKAY)
+            err = ECC_MONT_RED(r, m, mp, ct);
+    }
+
+    return err;
+}
+#endif /* !WOLFSSL_SP_MATH || WOLFSSL_PUBLIC_ECC_ADD_DBL */
+
 #if !defined(WOLFSSL_SP_MATH) || defined(WOLFSSL_PUBLIC_ECC_ADD_DBL)
 static int _ecc_projective_dbl_point(ecc_point *P, ecc_point *R, mp_int* a,
-                                     mp_int* modulus, mp_digit mp);
+                                     mp_int* modulus, mp_digit mp, int ct);
 
 /**
    Add two ECC points
@@ -2102,10 +2204,12 @@ static int _ecc_projective_dbl_point(ecc_point *P, ecc_point *R, mp_int* a,
    a        ECC curve parameter a
    modulus  The modulus of the field the ECC curve is in
    mp       The "b" value from montgomery_setup()
+   ct       Use constant time operations - values may be secret
    return   MP_OKAY on success
 */
 static int _ecc_projective_add_point(ecc_point* P, ecc_point* Q, ecc_point* R,
-                                     mp_int* a, mp_int* modulus, mp_digit mp)
+                                     mp_int* a, mp_int* modulus, mp_digit mp,
+                                     int ct)
 {
 #if !defined(WOLFSSL_SP_MATH)
    DECL_MP_INT_SIZE_DYN(t1, mp_bitsused(modulus), MAX_ECC_BITS_USE);
@@ -2118,6 +2222,8 @@ static int _ecc_projective_add_point(ecc_point* P, ecc_point* Q, ecc_point* R,
    mp_int  *x, *y, *z;
    int     err;
 
+   (void)ct;
+
    if (MP_BITS_OVER_MAX(mp_bitsused(modulus), MAX_ECC_BITS_USE)) {
        return WC_KEY_SIZE_E;
    }
@@ -2194,9 +2300,25 @@ static int _ecc_projective_add_point(ecc_point* P, ecc_point* Q, ecc_point* R,
 
    /* should we dbl instead? */
    if (err == MP_OKAY) {
-      if ( (mp_cmp(P->x, Q->x) == MP_EQ) &&
-           (mp_get_digit_count(Q->z) && mp_cmp(P->z, Q->z) == MP_EQ) &&
-           (mp_cmp(P->y, Q->y) == MP_EQ)) {
+      int eq;
+
+   #ifdef ECC_TIMING_RESISTANT
+      if (ct) {
+         /* Ordinates may be secret - compare all digits of every ordinate
+          * without early exit or short-circuit evaluation. */
+         eq  = (mp_cmp_ct(P->x, Q->x, modulus->used) == MP_EQ);
+         eq &= (mp_iszero(Q->z) == MP_NO);
+         eq &= (mp_cmp_ct(P->z, Q->z, modulus->used) == MP_EQ);
+         eq &= (mp_cmp_ct(P->y, Q->y, modulus->used) == MP_EQ);
+      }
+      else
+   #endif
+      {
+         eq = (mp_cmp(P->x, Q->x) == MP_EQ) &&
+              (mp_get_digit_count(Q->z) && mp_cmp(P->z, Q->z) == MP_EQ) &&
+              (mp_cmp(P->y, Q->y) == MP_EQ);
+      }
+      if (eq) {
           mp_clear(t1);
           mp_clear(t2);
    #ifdef WOLFSSL_SMALL_STACK
@@ -2213,7 +2335,7 @@ static int _ecc_projective_add_point(ecc_point* P, ecc_point* Q, ecc_point* R,
             FREE_MP_INT_SIZE(t2, NULL, DYNAMIC_TYPE_ECC);
             FREE_MP_INT_SIZE(t1, NULL, DYNAMIC_TYPE_ECC);
          }
-         return _ecc_projective_dbl_point(P, R, a, modulus, mp);
+         return _ecc_projective_dbl_point(P, R, a, modulus, mp, ct);
       }
    }
 
@@ -2257,53 +2379,37 @@ static int _ecc_projective_add_point(ecc_point* P, ecc_point* Q, ecc_point* R,
    if (err == MP_OKAY) {
        if (!mp_iszero(Q->z)) {
            /* T1 = Z' * Z' */
-           err = mp_sqr(Q->z, t1);
-           if (err == MP_OKAY)
-               err = mp_montgomery_reduce(t1, modulus, mp);
+           err = ecc_mont_sqr(Q->z, modulus, mp, t1, ct);
 
            /* X = X * T1 */
            if (err == MP_OKAY)
-               err = mp_mul(t1, x, x);
-           if (err == MP_OKAY)
-               err = mp_montgomery_reduce(x, modulus, mp);
+               err = ecc_mont_mul(t1, x, modulus, mp, x, ct);
 
            /* T1 = Z' * T1 */
            if (err == MP_OKAY)
-               err = mp_mul(Q->z, t1, t1);
-           if (err == MP_OKAY)
-               err = mp_montgomery_reduce(t1, modulus, mp);
+               err = ecc_mont_mul(Q->z, t1, modulus, mp, t1, ct);
 
            /* Y = Y * T1 */
            if (err == MP_OKAY)
-               err = mp_mul(t1, y, y);
-           if (err == MP_OKAY)
-               err = mp_montgomery_reduce(y, modulus, mp);
+               err = ecc_mont_mul(t1, y, modulus, mp, y, ct);
        }
    }
 
    /* T1 = Z*Z */
    if (err == MP_OKAY)
-       err = mp_sqr(z, t1);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(t1, modulus, mp);
+       err = ecc_mont_sqr(z, modulus, mp, t1, ct);
 
    /* T2 = X' * T1 */
    if (err == MP_OKAY)
-       err = mp_mul(Q->x, t1, t2);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(t2, modulus, mp);
+       err = ecc_mont_mul(Q->x, t1, modulus, mp, t2, ct);
 
    /* T1 = Z * T1 */
    if (err == MP_OKAY)
-       err = mp_mul(z, t1, t1);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(t1, modulus, mp);
+       err = ecc_mont_mul(z, t1, modulus, mp, t1, ct);
 
    /* T1 = Y' * T1 */
    if (err == MP_OKAY)
-       err = mp_mul(Q->y, t1, t1);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(t1, modulus, mp);
+       err = ecc_mont_mul(Q->y, t1, modulus, mp, t1, ct);
 
    /* Y = Y - T1 */
    if (err == MP_OKAY)
@@ -2327,47 +2433,33 @@ static int _ecc_projective_add_point(ecc_point* P, ecc_point* Q, ecc_point* R,
    if (err == MP_OKAY) {
        if (!mp_iszero(Q->z)) {
            /* Z = Z * Z' */
-           err = mp_mul(z, Q->z, z);
-           if (err == MP_OKAY)
-               err = mp_montgomery_reduce(z, modulus, mp);
+           err = ecc_mont_mul(z, Q->z, modulus, mp, z, ct);
        }
    }
 
    /* Z = Z * X */
    if (err == MP_OKAY)
-       err = mp_mul(z, x, z);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(z, modulus, mp);
+       err = ecc_mont_mul(z, x, modulus, mp, z, ct);
 
    /* T1 = T1 * X  */
    if (err == MP_OKAY)
-       err = mp_mul(t1, x, t1);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(t1, modulus, mp);
+       err = ecc_mont_mul(t1, x, modulus, mp, t1, ct);
 
    /* X = X * X */
    if (err == MP_OKAY)
-       err = mp_sqr(x, x);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(x, modulus, mp);
+       err = ecc_mont_sqr(x, modulus, mp, x, ct);
 
    /* T2 = T2 * x */
    if (err == MP_OKAY)
-       err = mp_mul(t2, x, t2);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(t2, modulus, mp);
+       err = ecc_mont_mul(t2, x, modulus, mp, t2, ct);
 
    /* T1 = T1 * X  */
    if (err == MP_OKAY)
-       err = mp_mul(t1, x, t1);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(t1, modulus, mp);
+       err = ecc_mont_mul(t1, x, modulus, mp, t1, ct);
 
    /* X = Y*Y */
    if (err == MP_OKAY)
-       err = mp_sqr(y, x);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(x, modulus, mp);
+       err = ecc_mont_sqr(y, modulus, mp, x, ct);
 
    /* X = X - T2 */
    if (err == MP_OKAY)
@@ -2380,9 +2472,7 @@ static int _ecc_projective_add_point(ecc_point* P, ecc_point* Q, ecc_point* R,
        err = mp_submod_ct(t2, x, modulus, t2);
    /* T2 = T2 * Y */
    if (err == MP_OKAY)
-       err = mp_mul(t2, y, t2);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(t2, modulus, mp);
+       err = ecc_mont_mul(t2, y, modulus, mp, t2, ct);
 
    /* Y = T2 - T1 */
    if (err == MP_OKAY)
@@ -2426,6 +2516,7 @@ static int _ecc_projective_add_point(ecc_point* P, ecc_point* Q, ecc_point* R,
 
     (void)a;
     (void)mp;
+    (void)ct;
 
 #if defined(WOLFSSL_SM2) && defined(WOLFSSL_SP_SM2)
     if ((modBits == 256) && (!mp_is_bit_set(modulus, 224))) {
@@ -2471,7 +2562,7 @@ int ecc_projective_add_point(ecc_point* P, ecc_point* Q, ecc_point* R,
         return ECC_OUT_OF_RANGE_E;
     }
 
-    return _ecc_projective_add_point(P, Q, R, a, modulus, mp);
+    return _ecc_projective_add_point(P, Q, R, a, modulus, mp, 1);
 }
 
 /* ### Point doubling in Jacobian coordinate system ###
@@ -2499,10 +2590,11 @@ int ecc_projective_add_point(ecc_point* P, ecc_point* Q, ecc_point* R,
    a   ECC curve parameter a
    modulus  The modulus of the field the ECC curve is in
    mp       The "b" value from montgomery_setup()
+   ct       Use constant time operations - values may be secret
    return   MP_OKAY on success
 */
 static int _ecc_projective_dbl_point(ecc_point *P, ecc_point *R, mp_int* a,
-                                     mp_int* modulus, mp_digit mp)
+                                     mp_int* modulus, mp_digit mp, int ct)
 {
 #if !defined(WOLFSSL_SP_MATH)
    DECL_MP_INT_SIZE_DYN(t1, mp_bitsused(modulus), MAX_ECC_BITS_USE);
@@ -2515,6 +2607,8 @@ static int _ecc_projective_dbl_point(ecc_point *P, ecc_point *R, mp_int* a,
    mp_int *x, *y, *z;
    int    err;
 
+   (void)ct;
+
    if (MP_BITS_OVER_MAX(mp_bitsused(modulus), MAX_ECC_BITS_USE)) {
        return WC_KEY_SIZE_E;
    }
@@ -2630,15 +2724,11 @@ static int _ecc_projective_dbl_point(ecc_point *P, ecc_point *R, mp_int* a,
 
    /* T1 = Z * Z */
    if (err == MP_OKAY)
-       err = mp_sqr(z, t1);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(t1, modulus, mp);
+       err = ecc_mont_sqr(z, modulus, mp, t1, ct);
 
    /* Z = Y * Z */
    if (err == MP_OKAY)
-       err = mp_mul(z, y, z);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(z, modulus, mp);
+       err = ecc_mont_mul(z, y, modulus, mp, z, ct);
 
    /* Z = 2Z */
    if (err == MP_OKAY)
@@ -2652,9 +2742,7 @@ static int _ecc_projective_dbl_point(ecc_point *P, ecc_point *R, mp_int* a,
    }
    if (err == MP_OKAY && mp_iszero((MP_INT_SIZE*)t2)) {
       /* T2 = X * X */
-      err = mp_sqr(x, t2);
-      if (err == MP_OKAY)
-          err = mp_montgomery_reduce(t2, modulus, mp);
+      err = ecc_mont_sqr(x, modulus, mp, t2, ct);
       /* T1 = T2 + T1 */
       if (err == MP_OKAY)
           err = mp_addmod_ct(t2, t2, modulus, t1);
@@ -2666,17 +2754,29 @@ static int _ecc_projective_dbl_point(ecc_point *P, ecc_point *R, mp_int* a,
       /* use "a" in calc */
 
       /* T2 = T1 * T1 */
-      err = mp_sqr(t1, t2);
-      if (err == MP_OKAY)
-          err = mp_montgomery_reduce(t2, modulus, mp);
+      err = ecc_mont_sqr(t1, modulus, mp, t2, ct);
       /* T1 = T2 * a */
-      if (err == MP_OKAY)
-          err = mp_mulmod(t2, a, modulus, t1);
+   #ifdef ECC_TIMING_RESISTANT
+      if (ct) {
+          /* T2 may be secret - mp_mulmod() is not constant time.
+           * T1 = a in Montgomery form - a is public so variable time ok. */
+          if (err == MP_OKAY)
+              err = mp_montgomery_calc_normalization(t1, modulus);
+          if (err == MP_OKAY)
+              err = mp_mulmod(t1, a, modulus, t1);
+          /* T1 = T2 * T1 - constant time Montgomery multiplication */
+          if (err == MP_OKAY)
+              err = ecc_mont_mul(t2, t1, modulus, mp, t1, ct);
+      }
+      else
+   #endif
+      {
+          if (err == MP_OKAY)
+              err = mp_mulmod(t2, a, modulus, t1);
+      }
       /* T2 = X * X */
       if (err == MP_OKAY)
-          err = mp_sqr(x, t2);
-      if (err == MP_OKAY)
-          err = mp_montgomery_reduce(t2, modulus, mp);
+          err = ecc_mont_sqr(x, modulus, mp, t2, ct);
       /* T1 = T2 + T1 */
       if (err == MP_OKAY)
           err = mp_addmod_ct(t1, t2, modulus, t1);
@@ -2701,9 +2801,7 @@ static int _ecc_projective_dbl_point(ecc_point *P, ecc_point *R, mp_int* a,
           err = mp_addmod_ct(t1, x, modulus, t1);
       /* T2 = T1 * T2 */
       if (err == MP_OKAY)
-          err = mp_mul(t1, t2, t2);
-      if (err == MP_OKAY)
-          err = mp_montgomery_reduce(t2, modulus, mp);
+          err = ecc_mont_mul(t1, t2, modulus, mp, t2, ct);
 
       /* T1 = 2T2 */
       if (err == MP_OKAY)
@@ -2718,15 +2816,11 @@ static int _ecc_projective_dbl_point(ecc_point *P, ecc_point *R, mp_int* a,
        err = mp_addmod_ct(y, y, modulus, y);
    /* Y = Y * Y */
    if (err == MP_OKAY)
-       err = mp_sqr(y, y);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(y, modulus, mp);
+       err = ecc_mont_sqr(y, modulus, mp, y, ct);
 
    /* T2 = Y * Y */
    if (err == MP_OKAY)
-       err = mp_sqr(y, t2);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(t2, modulus, mp);
+       err = ecc_mont_sqr(y, modulus, mp, t2, ct);
 
    /* T2 = T2/2 */
    if (err == MP_OKAY)
@@ -2734,15 +2828,11 @@ static int _ecc_projective_dbl_point(ecc_point *P, ecc_point *R, mp_int* a,
 
    /* Y = Y * X */
    if (err == MP_OKAY)
-       err = mp_mul(y, x, y);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(y, modulus, mp);
+       err = ecc_mont_mul(y, x, modulus, mp, y, ct);
 
    /* X = T1 * T1 */
    if (err == MP_OKAY)
-       err = mp_sqr(t1, x);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(x, modulus, mp);
+       err = ecc_mont_sqr(t1, modulus, mp, x, ct);
 
    /* X = X - Y */
    if (err == MP_OKAY)
@@ -2756,9 +2846,7 @@ static int _ecc_projective_dbl_point(ecc_point *P, ecc_point *R, mp_int* a,
        err = mp_submod_ct(y, x, modulus, y);
    /* Y = Y * T1 */
    if (err == MP_OKAY)
-       err = mp_mul(y, t1, y);
-   if (err == MP_OKAY)
-       err = mp_montgomery_reduce(y, modulus, mp);
+       err = ecc_mont_mul(y, t1, modulus, mp, y, ct);
 
    /* Y = Y - T2 */
    if (err == MP_OKAY)
@@ -2798,6 +2886,7 @@ static int _ecc_projective_dbl_point(ecc_point *P, ecc_point *R, mp_int* a,
 
     (void)a;
     (void)mp;
+    (void)ct;
 
 #if defined(WOLFSSL_SM2) && defined(WOLFSSL_SP_SM2)
     if ((modBits == 256) && (!mp_is_bit_set(modulus, 224))) {
@@ -2835,7 +2924,7 @@ int ecc_projective_dbl_point(ecc_point *P, ecc_point *R, mp_int* a,
         return ECC_OUT_OF_RANGE_E;
     }
 
-    return _ecc_projective_dbl_point(P, R, a, modulus, mp);
+    return _ecc_projective_dbl_point(P, R, a, modulus, mp, 1);
 }
 
 #if !defined(FREESCALE_LTC_ECC) && \
@@ -2980,44 +3069,54 @@ int ecc_map_ex(ecc_point* P, mp_int* modulus, mp_digit mp, int ct)
         z = P->z;
         #endif
 
-        /* get 1/z */
+        /* get 1/z, 1/z^2 and 1/z^3 */
         if (err == MP_OKAY) {
         #if defined(ECC_TIMING_RESISTANT) && (defined(USE_FAST_MATH) || \
                         defined(WOLFSSL_SP_MATH) || defined(WOLFSSL_SP_MATH_ALL))
             if (ct) {
+                /* Values derived from the private scalar - keep all
+                 * reductions constant time. Stay in Montgomery form so that
+                 * Montgomery reduction can be used instead of mp_mod().
+                 */
+                /* t1 = 1/z (Montgomery form) */
                 err = mp_invmod_mont_ct(z, modulus, t1, mp);
+                /* t2 = 1/z^2 (Montgomery form) */
+                if (err == MP_OKAY)
+                    err = ecc_mont_sqr(t1, modulus, mp, t2, ct);
+                /* t1 = 1/z^3 (Montgomery form) */
                 if (err == MP_OKAY)
-                    err = mp_montgomery_reduce(t1, modulus, mp);
+                    err = ecc_mont_mul(t1, t2, modulus, mp, t1, ct);
+                /* Convert 1/z^2 and 1/z^3 to normal form. */
+                if (err == MP_OKAY)
+                    err = ECC_MONT_RED(t2, modulus, mp, ct);
+                if (err == MP_OKAY)
+                    err = ECC_MONT_RED(t1, modulus, mp, ct);
             }
             else
         #endif
             {
                 /* first map z back to normal */
-                err = mp_montgomery_reduce(z, modulus, mp);
+                err = ECC_MONT_RED(z, modulus, mp, ct);
                 if (err == MP_OKAY)
                     err = mp_invmod(z, modulus, t1);
+
+                /* get 1/z^2 and 1/z^3 */
+                if (err == MP_OKAY)
+                    err = mp_sqr(t1, t2);
+                if (err == MP_OKAY)
+                    err = mp_mod(t2, modulus, t2);
+                if (err == MP_OKAY)
+                    err = mp_mul(t1, t2, t1);
+                if (err == MP_OKAY)
+                    err = mp_mod(t1, modulus, t1);
             }
         }
 
-        /* get 1/z^2 and 1/z^3 */
-        if (err == MP_OKAY)
-            err = mp_sqr(t1, t2);
-        if (err == MP_OKAY)
-            err = mp_mod(t2, modulus, t2);
-        if (err == MP_OKAY)
-            err = mp_mul(t1, t2, t1);
-        if (err == MP_OKAY)
-            err = mp_mod(t1, modulus, t1);
-
         /* multiply against x/y */
         if (err == MP_OKAY)
-            err = mp_mul(x, t2, x);
-        if (err == MP_OKAY)
-            err = mp_montgomery_reduce(x, modulus, mp);
+            err = ecc_mont_mul(x, t2, modulus, mp, x, ct);
         if (err == MP_OKAY)
-            err = mp_mul(y, t1, y);
-        if (err == MP_OKAY)
-            err = mp_montgomery_reduce(y, modulus, mp);
+            err = ecc_mont_mul(y, t1, modulus, mp, y, ct);
 
         if (err == MP_OKAY)
             err = mp_set(z, 1);
@@ -3270,27 +3369,27 @@ static int wc_ecc_gen_z(WC_RNG* rng, int size, ecc_point* p, mp_int* modulus,
     if (err == MP_OKAY)
         err = mp_mul(p->z, ty, p->z);
     if (err == MP_OKAY)
-        err = mp_montgomery_reduce(p->z, modulus, mp);
+        err = ECC_MONT_RED(p->z, modulus, mp, 1);
     /* Square random value for X (X' = X / Z^2). */
     if (err == MP_OKAY)
         err = mp_sqr(ty, tx);
     if (err == MP_OKAY)
-        err = mp_montgomery_reduce(tx, modulus, mp);
+        err = ECC_MONT_RED(tx, modulus, mp, 1);
     /* Multiply square of random by random value for Y. */
     if (err == MP_OKAY)
         err = mp_mul(ty, tx, ty);
     if (err == MP_OKAY)
-        err = mp_montgomery_reduce(ty, modulus, mp);
+        err = ECC_MONT_RED(ty, modulus, mp, 1);
     /* Multiply square into X. */
     if (err == MP_OKAY)
         err = mp_mul(p->x, tx, p->x);
     if (err == MP_OKAY)
-        err = mp_montgomery_reduce(p->x, modulus, mp);
+        err = ECC_MONT_RED(p->x, modulus, mp, 1);
     /* Multiply cube into Y (Y' = Y / Z^3). */
     if (err == MP_OKAY)
         err = mp_mul(p->y, ty, p->y);
     if (err == MP_OKAY)
-        err = mp_montgomery_reduce(p->y, modulus, mp);
+        err = ECC_MONT_RED(p->y, modulus, mp, 1);
 
     return err;
 }
@@ -5300,15 +5399,15 @@ int wc_ecc_shared_secret_gen_sync(ecc_key* private_key, ecc_point* point,
         }
         if (err == MP_OKAY) {
             x = mp_unsigned_bin_size(curve->prime);
-            if (*outlen < (word32)x || x < mp_unsigned_bin_size(result->x)) {
+            if (*outlen < (word32)x) {
                 err = BUFFER_E;
             }
         }
 
         if (err == MP_OKAY) {
-            XMEMSET(out, 0, (size_t)x);
-            err = mp_to_unsigned_bin(result->x, out +
-                                     (x - mp_unsigned_bin_size(result->x)));
+            /* Shared secret is secret - don't make the encoding depend on its
+             * length. x ordinate is less than the prime after mapping. */
+            err = mp_to_unsigned_bin_len_ct(result->x, out, x);
         }
         *outlen = (word32)x;
 
@@ -7503,28 +7602,84 @@ int wc_ecc_sign_hash_ex(const byte* in, word32 inlen, WC_RNG* rng,
       !defined(WOLFSSL_KCAPI_ECC) && \
       !defined(WOLFSSL_NO_ECC_SW)
 #ifndef WOLFSSL_SP_MATH
+/* Multiply two numbers modulo the order of the curve: r = a * b mod n
+ *
+ * Constant time - values are derived from the private key and nonce.
+ * Uses constant time Montgomery multiplication twice instead of a division
+ * based reduction. Operands are padded to the length of the order so that the
+ * operations don't depend on the number of leading zero digits:
+ *   r = (a * b) / R mod n
+ *   r = (r * R^2) / R mod n = a * b mod n
+ *
+ * a and b must be less than n.
+ *
+ * @param [in]  a      First operand.
+ * @param [in]  b      Second operand.
+ * @param [in]  order  Order of the curve - n.
+ * @param [in]  mp     Montgomery multiplier for order.
+ * @param [in]  norm2  R^2 mod n.
+ * @param [out] r      Result. May be the same as a or b.
+ * @return  MP_OKAY on success.
+ */
+static int ecc_mulmod_order_ct(mp_int* a, mp_int* b, mp_int* order,
+    mp_digit mp, mp_int* norm2, mp_int* r)
+{
+    /* r = a.b / R mod n */
+    int err = mp_mont_mul_ct(a, b, order, mp, r);
+    /* r = (a.b / R).R^2 / R mod n = a.b mod n */
+    if (err == MP_OKAY)
+        err = mp_mont_mul_ct(r, norm2, order, mp, r);
+
+    return err;
+}
+
 static int ecc_sign_hash_sw(ecc_key* key, ecc_key* pubkey, WC_RNG* rng,
                             ecc_curve_spec* curve, mp_int* e, mp_int* r,
                             mp_int* s)
 {
     int err = MP_OKAY;
     int loop_check = 0;
+    mp_digit mp = 0;
 
     DECL_MP_INT_SIZE_DYN(b, ECC_KEY_MAX_BITS_NONULLCHECK(key), MAX_ECC_BITS_USE);
+    DECL_MP_INT_SIZE_DYN(norm2, ECC_KEY_MAX_BITS_NONULLCHECK(key),
+        MAX_ECC_BITS_USE);
 
     if (MP_BITS_OVER_MAX(ECC_KEY_MAX_BITS_NONULLCHECK(key), MAX_ECC_BITS_USE)) {
         return WC_KEY_SIZE_E;
     }
 
     NEW_MP_INT_SIZE(b, ECC_KEY_MAX_BITS_NONULLCHECK(key), key->heap, DYNAMIC_TYPE_ECC);
+    NEW_MP_INT_SIZE(norm2, ECC_KEY_MAX_BITS_NONULLCHECK(key), key->heap,
+        DYNAMIC_TYPE_ECC);
 #ifdef MP_INT_SIZE_CHECK_NULL
-    if (b == NULL)
+    if ((b == NULL) || (norm2 == NULL))
         err = MEMORY_E;
 #endif
 
     if (err == MP_OKAY) {
         err = INIT_MP_INT_SIZE(b, ECC_KEY_MAX_BITS_NONULLCHECK(key));
     }
+    if (err == MP_OKAY) {
+        err = INIT_MP_INT_SIZE(norm2, ECC_KEY_MAX_BITS_NONULLCHECK(key));
+    }
+    if (err != MP_OKAY) {
+        /* Temporaries not all initialized - nothing secret stored yet. */
+        FREE_MP_INT_SIZE(norm2, key->heap, DYNAMIC_TYPE_ECC);
+        FREE_MP_INT_SIZE(b, key->heap, DYNAMIC_TYPE_ECC);
+        return err;
+    }
+    /* Setup for constant time Montgomery multiplication modulo the order.
+     * norm2 = R^2 mod n - order is public. */
+    if (err == MP_OKAY) {
+        err = mp_montgomery_setup(curve->order, &mp);
+    }
+    if (err == MP_OKAY) {
+        err = mp_montgomery_calc_normalization(norm2, curve->order);
+    }
+    if (err == MP_OKAY) {
+        err = mp_sqrmod(norm2, curve->order, norm2);
+    }
 
 #ifdef WOLFSSL_CUSTOM_CURVES
     /* if custom curve, apply params to pubkey */
@@ -7632,7 +7787,7 @@ static int ecc_sign_hash_sw(ecc_key* key, ecc_key* pubkey, WC_RNG* rng,
                       = b.(e/k.b + x.r/k.b) */
 
             /* k' = k.b */
-            err = mp_mulmod(kp, b, curve->order, kp);
+            err = ecc_mulmod_order_ct(kp, b, curve->order, mp, norm2, kp);
             if (err != MP_OKAY) break;
 
             /* k' = 1/k.b
@@ -7641,17 +7796,18 @@ static int ecc_sign_hash_sw(ecc_key* key, ecc_key* pubkey, WC_RNG* rng,
             if (err != MP_OKAY) break;
 
             /* s = x.r */
-            err = mp_mulmod(x, r, curve->order, s);
+            err = ecc_mulmod_order_ct(x, r, curve->order, mp, norm2, s);
             if (err != MP_OKAY) break;
 
             /* s = x.r/k.b
                  = k'.s */
-            err = mp_mulmod(kp, s, curve->order, s);
+            err = ecc_mulmod_order_ct(kp, s, curve->order, mp, norm2, s);
             if (err != MP_OKAY) break;
 
             /* e' = e/k.b
                   = e.k' */
-            err = mp_mulmod(kp, e, curve->order, ep);
+            err = ecc_mulmod_order_ct(kp, e, curve->order, mp, norm2,
+                ep);
             if (err != MP_OKAY) break;
 
             /* s = e/k.b + x.r/k.b = (e + x.r)/k.b
@@ -7661,7 +7817,7 @@ static int ecc_sign_hash_sw(ecc_key* key, ecc_key* pubkey, WC_RNG* rng,
 
             /* s = b.(e + x.r)/k.b = (e + x.r)/k
                  = b.s */
-            err = mp_mulmod(s, b, curve->order, s);
+            err = ecc_mulmod_order_ct(s, b, curve->order, mp, norm2, s);
             if (err != MP_OKAY) break;
 
             if (mp_iszero(s) == MP_NO) {
@@ -7682,7 +7838,9 @@ static int ecc_sign_hash_sw(ecc_key* key, ecc_key* pubkey, WC_RNG* rng,
          WC_RELAX_LONG_LOOP();
     }
     ecc_put_k(key);
+    mp_clear(norm2);
     mp_forcezero(b);
+    FREE_MP_INT_SIZE(norm2, key->heap, DYNAMIC_TYPE_ECC);
     FREE_MP_INT_SIZE(b, key->heap, DYNAMIC_TYPE_ECC);
 #if !defined(WOLFSSL_SMALL_STACK) && defined(WOLFSSL_CHECK_MEM_ZERO)
     mp_memzero_check(b);
@@ -8986,9 +9144,13 @@ static int ecc_cond_copy_point(ecc_point* S, int copy, ecc_point* D)
  * The infinity out-flag is set only for the last case - a finite pair whose
  * sum is infinity - matching the original branching version, which never
  * touched the flag when copying an operand out.
+ *
+ * ct indicates the values may be secret and constant time operations are
+ * required.
  */
-int ecc_projective_add_point_safe(ecc_point* A, ecc_point* B, ecc_point* R,
-    mp_int* a, mp_int* modulus, mp_digit mp, int* infinity)
+static int _ecc_projective_add_point_safe(ecc_point* A, ecc_point* B,
+    ecc_point* R, mp_int* a, mp_int* modulus, mp_digit mp, int* infinity,
+    int ct)
 {
     int err;
     int aInf;
@@ -9027,7 +9189,7 @@ int ecc_projective_add_point_safe(ecc_point* A, ecc_point* B, ecc_point* R,
     T->key = R->key;
 #endif
 
-    err = _ecc_projective_add_point(A, B, T, a, modulus, mp);
+    err = _ecc_projective_add_point(A, B, T, a, modulus, mp, ct);
     if ((err == MP_OKAY) && mp_iszero(T->z)) {
         /* When all zero then should have done a double */
         if (mp_iszero(T->x) && mp_iszero(T->y)) {
@@ -9037,11 +9199,12 @@ int ecc_projective_add_point_safe(ecc_point* A, ecc_point* B, ecc_point* R,
                     err = mp_montgomery_calc_normalization(T->z, modulus);
                 }
                 if (err == MP_OKAY) {
-                    err = _ecc_projective_dbl_point(T, T, a, modulus, mp);
+                    err = _ecc_projective_dbl_point(T, T, a, modulus, mp,
+                        ct);
                 }
             }
             else {
-                err = _ecc_projective_dbl_point(B, T, a, modulus, mp);
+                err = _ecc_projective_dbl_point(B, T, a, modulus, mp, ct);
             }
         }
         /* When only Z zero then result is infinity */
@@ -9081,13 +9244,27 @@ int ecc_projective_add_point_safe(ecc_point* A, ecc_point* B, ecc_point* R,
     return err;
 }
 
+/* Add A and B, handling the exceptional cases. Constant time.
+ *
+ * See _ecc_projective_add_point_safe().
+ */
+int ecc_projective_add_point_safe(ecc_point* A, ecc_point* B, ecc_point* R,
+    mp_int* a, mp_int* modulus, mp_digit mp, int* infinity)
+{
+    return _ecc_projective_add_point_safe(A, B, R, a, modulus, mp, infinity,
+        1);
+}
+
 /* Handles when P is the infinity point.
  *
  * Double infinity -> infinity.
  * Otherwise do normal double - which can't lead to infinity as odd order.
+ *
+ * ct indicates the values may be secret and constant time operations are
+ * required.
  */
-int ecc_projective_dbl_point_safe(ecc_point *P, ecc_point *R, mp_int* a,
-                                  mp_int* modulus, mp_digit mp)
+static int _ecc_projective_dbl_point_safe(ecc_point *P, ecc_point *R,
+    mp_int* a, mp_int* modulus, mp_digit mp, int ct)
 {
     int err;
     int inf;
@@ -9098,7 +9275,7 @@ int ecc_projective_dbl_point_safe(ecc_point *P, ecc_point *R, mp_int* a,
     inf  = (mp_iszero(P->x) == MP_YES);
     inf &= (mp_iszero(P->y) == MP_YES);
 
-    err = _ecc_projective_dbl_point(P, R, a, modulus, mp);
+    err = _ecc_projective_dbl_point(P, R, a, modulus, mp, ct);
     if (err == MP_OKAY) {
         inf |= (mp_iszero(R->z) == MP_YES);
         if (inf) {
@@ -9108,6 +9285,16 @@ int ecc_projective_dbl_point_safe(ecc_point *P, ecc_point *R, mp_int* a,
 
     return err;
 }
+
+/* Double P, handling when P is the infinity point. Constant time.
+ *
+ * See _ecc_projective_dbl_point_safe().
+ */
+int ecc_projective_dbl_point_safe(ecc_point *P, ecc_point *R, mp_int* a,
+                                  mp_int* modulus, mp_digit mp)
+{
+    return _ecc_projective_dbl_point_safe(P, R, a, modulus, mp, 1);
+}
 #endif /* !(WOLFSSL_ATECC508A) && !(WOLFSSL_ATECC608A) && \
           !(WOLFSSL_SP_MATH) && \
           (!(WOLF_CRYPTO_CB_ONLY_ECC) || (WOLFSSL_QNX_CAAM) || \
@@ -9343,20 +9530,22 @@ int ecc_mul2add(ecc_point* A, mp_int* kA,
 
   if (err == MP_OKAY) {
     /* precomp [i,0](A + B) table */
-    err = ecc_projective_dbl_point_safe(precomp[1], precomp[2], a, modulus, mp);
+    err = _ecc_projective_dbl_point_safe(precomp[1], precomp[2], a, modulus, mp,
+                                         0);
   }
   if (err == MP_OKAY) {
-    err = ecc_projective_add_point_safe(precomp[1], precomp[2], precomp[3],
-                                                          a, modulus, mp, NULL);
+    err = _ecc_projective_add_point_safe(precomp[1], precomp[2], precomp[3],
+                                                       a, modulus, mp, NULL, 0);
   }
 
   if (err == MP_OKAY) {
     /* precomp [0,i](A + B) table */
-    err = ecc_projective_dbl_point_safe(precomp[4], precomp[8], a, modulus, mp);
+    err = _ecc_projective_dbl_point_safe(precomp[4], precomp[8], a, modulus, mp,
+                                         0);
   }
   if (err == MP_OKAY) {
-    err = ecc_projective_add_point_safe(precomp[4], precomp[8], precomp[12], a,
-                                                             modulus, mp, NULL);
+    err = _ecc_projective_add_point_safe(precomp[4], precomp[8], precomp[12], a,
+                                                          modulus, mp, NULL, 0);
   }
 
   if (err == MP_OKAY) {
@@ -9364,9 +9553,9 @@ int ecc_mul2add(ecc_point* A, mp_int* kA,
     for (x = 1; x < 4; x++) {
       for (y = 1; y < 4; y++) {
         if (err == MP_OKAY) {
-          err = ecc_projective_add_point_safe(precomp[x], precomp[(y<<2)],
+          err = _ecc_projective_add_point_safe(precomp[x], precomp[(y<<2)],
                                                   precomp[x+(y<<2)], a, modulus,
-                                                  mp, NULL);
+                                                  mp, NULL, 0);
         }
       }
     }
@@ -9404,9 +9593,9 @@ int ecc_mul2add(ecc_point* A, mp_int* kA,
         if (first == 0) {
             /* double twice */
             if (err == MP_OKAY)
-                err = ecc_projective_dbl_point_safe(C, C, a, modulus, mp);
+                err = _ecc_projective_dbl_point_safe(C, C, a, modulus, mp, 0);
             if (err == MP_OKAY)
-                err = ecc_projective_dbl_point_safe(C, C, a, modulus, mp);
+                err = _ecc_projective_dbl_point_safe(C, C, a, modulus, mp, 0);
             else
                 break;
         }
@@ -9430,9 +9619,9 @@ int ecc_mul2add(ecc_point* A, mp_int* kA,
             } else {
                 /* if not first, add from table */
                 if (err == MP_OKAY)
-                    err = ecc_projective_add_point_safe(C, precomp[i],
+                    err = _ecc_projective_add_point_safe(C, precomp[i],
                                                         C, a, modulus, mp,
-                                                        &first);
+                                                        &first, 0);
                 if (err != MP_OKAY)
                     break;
             }
@@ -10181,8 +10370,8 @@ static int ecc_verify_hash(mp_int *r, mp_int *s, const byte* hash,
              * recovery for a missed-double doubles B into R, which requires
              * B and R to be different points. */
             if (err == MP_OKAY)
-                err = ecc_projective_add_point_safe(mG, mQ, mG, curve->Af,
-                                                        curve->prime, mp, NULL);
+                err = _ecc_projective_add_point_safe(mG, mQ, mG, curve->Af,
+                                                     curve->prime, mp, NULL, 0);
         #ifdef WOLFSSL_CHECK_VER_FAULTS
             if (err == MP_OKAY && wc_ecc_cmp_point(mG, mG1) == MP_EQ) {
                 err = BAD_STATE_E;
diff --git a/wolfcrypt/src/sp_int.c b/wolfcrypt/src/sp_int.c
index bdc67b67221..3a4f79a7367 100644
--- a/wolfcrypt/src/sp_int.c
+++ b/wolfcrypt/src/sp_int.c
@@ -5199,7 +5199,8 @@ WOLFSSL_LOCAL int sp_ModExp_4096(const sp_int* base, const sp_int* exp,
 #if defined(WOLFSSL_SP_MATH_ALL) || defined(WOLFSSL_HAVE_SP_DH) || \
     defined(OPENSSL_ALL)
 static int _sp_mont_red(sp_int* a, const sp_int* m, sp_int_digit mp, int ct);
-static void _sp_mont_red_sub_ct(sp_int* a, const sp_int* m, unsigned int n);
+WC_MAYBE_UNUSED static void _sp_mont_red_sub_ct(sp_int* a, const sp_int* m,
+    unsigned int n);
 #endif
 #if defined(WOLFSSL_SP_MATH_ALL) || defined(WOLFSSL_HAVE_SP_DH) || \
     defined(WOLFCRYPT_HAVE_ECCSI) || defined(WOLFCRYPT_HAVE_SAKKE) || \
@@ -5897,8 +5898,8 @@ int sp_cmp(const sp_int* a, const sp_int* b)
 }
 #endif
 
-#if defined(HAVE_ECC) && !defined(WC_NO_RNG) && \
-    defined(WOLFSSL_ECC_GEN_REJECT_SAMPLING)
+#if defined(HAVE_ECC) && ((!defined(WC_NO_RNG) && \
+    defined(WOLFSSL_ECC_GEN_REJECT_SAMPLING)) || defined(ECC_TIMING_RESISTANT))
 /* Compare two multi-precision numbers in constant time.
  *
  * Assumes a and b are not NULL.
@@ -5918,7 +5919,7 @@ static int _sp_cmp_ct(const sp_int* a, const sp_int* b, unsigned int n)
     int i;
     volatile int mask = -1;
 
-    for (i = n - 1; i >= 0; i--) {
+    for (i = (int)n - 1; i >= 0; i--) {
         sp_int_digit ad = a->dp[i] & ((sp_int_digit)0 - (i < (int)a->used));
         sp_int_digit bd = b->dp[i] & ((sp_int_digit)0 - (i < (int)b->used));
 
@@ -5969,7 +5970,8 @@ int sp_cmp_ct(const sp_int* a, const sp_int* b, unsigned int n)
 
     return ret;
 }
-#endif /* HAVE_ECC && !WC_NO_RNG && WOLFSSL_ECC_GEN_REJECT_SAMPLING */
+#endif /* HAVE_ECC && ((!WC_NO_RNG && WOLFSSL_ECC_GEN_REJECT_SAMPLING) ||
+        * ECC_TIMING_RESISTANT) */
 
 /* Constant time clamping.
  *
@@ -13021,14 +13023,14 @@ static int _sp_invmod_mont_ct(const sp_int* a, const sp_int* m, sp_int* r,
             _sp_init_size(pre[i], (sp_size_t)(m->used * 2 + 1));
             err = sp_sqr(pre[i-1], pre[i]);
             if (err == MP_OKAY) {
-                err = _sp_mont_red(pre[i], m, mp, 0);
+                err = _sp_mont_red(pre[i], m, mp, 1);
             }
             /* ..10 -> ..11 */
             if (err == MP_OKAY) {
                 err = sp_mul(pre[i], a, pre[i]);
             }
             if (err == MP_OKAY) {
-                err = _sp_mont_red(pre[i], m, mp, 0);
+                err = _sp_mont_red(pre[i], m, mp, 1);
             }
         }
     }
@@ -13082,7 +13084,7 @@ static int _sp_invmod_mont_ct(const sp_int* a, const sp_int* m, sp_int* r,
                     /* 6.4.2.1. t = (t ^ 2) mod m */
                     err = sp_sqr(t, t);
                     if (err == MP_OKAY) {
-                        err = _sp_mont_red(t, m, mp, 0);
+                        err = _sp_mont_red(t, m, mp, 1);
                     }
                 }
                 /* 6.4.3. s = 1 - bit */
@@ -13093,7 +13095,7 @@ static int _sp_invmod_mont_ct(const sp_int* a, const sp_int* m, sp_int* r,
                     err = sp_mul(t, pre[j-1], t);
                 }
                 if (err == MP_OKAY) {
-                    err = _sp_mont_red(t, m, mp, 0);
+                    err = _sp_mont_red(t, m, mp, 1);
                 }
                 /* 6.4.5. j = 0
                  *        Reset number of 1 bits seen.
@@ -13109,21 +13111,29 @@ static int _sp_invmod_mont_ct(const sp_int* a, const sp_int* m, sp_int* r,
             /* 7.1. t = (t ^ 2) mod m */
             err = sp_sqr(t, t);
             if (err == MP_OKAY) {
-                err = _sp_mont_red(t, m, mp, 0);
+                err = _sp_mont_red(t, m, mp, 1);
             }
         }
     }
     if (err == MP_OKAY) {
         /* 8. If j > 0 then r = (t * pre[j-1]) mod m */
         if (j > 0) {
-            err = sp_mul(t, pre[j-1], r);
+            /* Reduce in t - reduction needs 2 * m->used + 1 digits. */
+            err = sp_mul(t, pre[j-1], t);
             if (err == MP_OKAY) {
-                err = _sp_mont_red(r, m, mp, 0);
+                err = _sp_mont_red(t, m, mp, 1);
             }
         }
         /* 9. Else r = t */
-        else {
-            _sp_copy(t, r);
+        if (err == MP_OKAY) {
+            /* Copy fixed number of digits - result is less than m. */
+            for (i = 0; i < (int)m->used; i++) {
+                r->dp[i] = t->dp[i];
+            }
+            r->used = t->used;
+        #ifdef WOLFSSL_SP_INT_NEGATIVE
+            r->sign = MP_ZPOS;
+        #endif
         }
     }
 
@@ -13183,6 +13193,119 @@ int sp_invmod_mont_ct(const sp_int* a, const sp_int* m, sp_int* r,
 
 #endif /* WOLFSSL_SP_INVMOD_MONT_CT */
 
+#if defined(WOLFSSL_SP_MATH_ALL) && defined(HAVE_ECC)
+static int _sp_sqr_nc(const sp_int* a, sp_int* r);
+
+/* Montgomery multiplication in constant time: r = a * b / R mod m
+ *
+ * Operands are padded to the length of the modulus, the product is not clamped
+ * and the result is copied with a fixed number of digits. The operations
+ * performed depend only on the sizes of a, b, m and r.
+ *
+ * a * b must be less than m * R - true when a is less than m and b is less
+ * than R.
+ *
+ * @param [in]  a   SP integer to multiply. Has no more digits than m.
+ * @param [in]  b   SP integer to multiply. Has no more digits than m.
+ * @param [in]  m   SP integer that is the modulus. Must be odd.
+ * @param [in]  mp  SP integer digit that is the bottom digit of inv(-m).
+ * @param [out] r   SP integer result. May be the same as a or b.
+ *
+ * @return  MP_OKAY on success.
+ * @return  MP_VAL when a, b, m or r is NULL, an operand has more digits than
+ *          m, m is zero or too big, or r is too small.
+ * @return  MP_MEM when dynamic memory allocation fails.
+ */
+int sp_mont_mul_ct(const sp_int* a, const sp_int* b, const sp_int* m,
+    sp_int_digit mp, sp_int* r)
+{
+    int err = MP_OKAY;
+    unsigned int i;
+    unsigned int n = 0;
+    sp_int* ta = NULL;
+    sp_int* tb = NULL;
+    sp_int* t = NULL;
+    /* Called for each field multiplication in ECC point arithmetic - only
+     * allocate when small stack. */
+    DECL_SP_INT_ARRAY(pre, m == NULL ? 1 : m->used * 2 + 1, 3);
+
+    /* Validate parameters. Lengths are public. */
+    if ((a == NULL) || (b == NULL) || (m == NULL) || (r == NULL)) {
+        err = MP_VAL;
+    }
+    else if (sp_iszero(m) || (m->used * 2 >= SP_INT_DIGITS) ||
+             (a->used > m->used) || (b->used > m->used) ||
+             (m->used > r->size)) {
+        err = MP_VAL;
+    }
+#ifdef WOLFSSL_SP_INT_NEGATIVE
+    else if ((a->sign == MP_NEG) || (b->sign == MP_NEG) ||
+             (m->sign == MP_NEG)) {
+        err = MP_VAL;
+    }
+#endif
+
+    if (err == MP_OKAY) {
+        n = m->used;
+        ALLOC_SP_INT_ARRAY(pre, n * 2U + 1U, 3, err, NULL);
+    }
+    if (err == MP_OKAY) {
+        ta = pre[0];
+        tb = pre[1];
+        t  = pre[2];
+        _sp_init_size(ta, (sp_size_t)(n * 2 + 1));
+        _sp_init_size(tb, (sp_size_t)(n * 2 + 1));
+        _sp_init_size(t, (sp_size_t)(n * 2 + 1));
+
+        /* Pad operands to the length of the modulus. Digits above used are
+         * masked out. Digits read depend only on the sizes. */
+        for (i = 0; i < n; i++) {
+            sp_int_digit d;
+
+            d = (i < (unsigned int)a->size) ? a->dp[i] : 0;
+            ta->dp[i] = d & ((sp_int_digit)0 -
+                             (sp_int_digit)(i < (unsigned int)a->used));
+            d = (i < (unsigned int)b->size) ? b->dp[i] : 0;
+            tb->dp[i] = d & ((sp_int_digit)0 -
+                             (sp_int_digit)(i < (unsigned int)b->used));
+        }
+        ta->used = (sp_size_t)n;
+        tb->used = (sp_size_t)n;
+
+        /* t = a * b, not clamped, and reduced in constant time.
+         * Same pointer is public - square when operands are the same. */
+        if (a == b) {
+            err = _sp_sqr_nc(ta, t);
+        }
+        else {
+            err = _sp_mul_nc(ta, tb, t);
+        }
+    }
+    if (err == MP_OKAY) {
+        err = _sp_mont_red(t, m, mp, 1);
+    }
+    if (err == MP_OKAY) {
+        /* Copy a fixed number of digits - result is less than m. */
+        for (i = 0; i < n; i++) {
+            r->dp[i] = t->dp[i];
+        }
+        r->used = t->used;
+    #ifdef WOLFSSL_SP_INT_NEGATIVE
+        r->sign = MP_ZPOS;
+    #endif
+    }
+
+    if (ta != NULL) {
+        /* Operands use n digits and the product 2 * n + 1 digits. */
+        ForceZero(ta->dp, n * (word32)SP_WORD_SIZEOF);
+        ForceZero(tb->dp, n * (word32)SP_WORD_SIZEOF);
+        ForceZero(t->dp, (n * 2 + 1) * (word32)SP_WORD_SIZEOF);
+    }
+    FREE_SP_INT_ARRAY(pre, NULL);
+    return err;
+}
+#endif /* WOLFSSL_SP_MATH_ALL && HAVE_ECC */
+
 
 /**************************
  * Exponentiation functions
@@ -13533,7 +13656,10 @@ static void _sp_pad_to_mod(sp_int* a, const sp_int* m)
 #if defined(WOLFSSL_SP_MATH) && defined(WOLFSSL_SP_SMALL)
 #define SP_EXPTMOD_SQR(a, r)            _sp_mul_nc(a, a, r)
 #else
+#if !defined(WOLFSSL_SP_MATH_ALL) || !defined(HAVE_ECC)
+/* Otherwise declared before sp_mont_mul_ct(). */
 static int _sp_sqr_nc(const sp_int* a, sp_int* r);
+#endif
 #define SP_EXPTMOD_SQR(a, r)            _sp_sqr_nc(a, r)
 #endif
 #else
@@ -17941,6 +18067,114 @@ int sp_sqrmod(const sp_int* a, const sp_int* m, sp_int* r)
 #if defined(WOLFSSL_SP_MATH_ALL) || defined(WOLFSSL_HAVE_SP_DH) || \
     defined(WOLFCRYPT_HAVE_ECCSI) || defined(WOLFCRYPT_HAVE_SAKKE) || \
     defined(OPENSSL_ALL)
+#ifndef WOLFSSL_NO_CT_OPS
+/* Constant time zeroing of the digits of a from used up to n2.
+ *
+ * Inline so that the loop is unrolled when n2 is a constant.
+ *
+ * @param [in, out] a   SP integer to zero top digits of.
+ * @param [in]      n2  Number of digits that must be valid.
+ */
+static WC_INLINE void _sp_mont_red_zero_ct(sp_int* a, unsigned int n2)
+{
+    unsigned int i;
+    int top = (int)a->used - 1;
+
+    for (i = 0; i < n2; i++) {
+        a->dp[i] &= (sp_int_digit)(sp_int_sdigit)ctMaskIntGTE(top, (int)i);
+    }
+}
+#endif /* !WOLFSSL_NO_CT_OPS */
+
+/* Constant time final subtraction of Montgomery reduction.
+ *
+ * a = a - m when a >= m, otherwise a is unchanged.
+ * Always calculates a - m and then selects, in constant time, the result.
+ *
+ * Assumes a is less than 2 * m, digits 0..n of a are set, a has space for at
+ * least 2 * n + 1 digits and n is m->used.
+ * Digits n+1..2n of a are used as scratch - they are not part of the result.
+ * Inline so that the loops are unrolled when n is a constant.
+ *
+ * @param [in, out] a  SP integer to reduce.
+ * @param [in]      m  SP integer that is the modulus.
+ * @param [in]      n  Number of digits in modulus.
+ */
+static WC_INLINE void _sp_mont_red_sub_ct_n(sp_int* a, const sp_int* m,
+    unsigned int n)
+{
+#ifndef SQR_MUL_ASM
+    sp_int_sword w;
+#else
+    sp_int_digit l;
+    sp_int_digit h;
+    sp_int_digit b;
+#endif
+    sp_int_digit t;
+    sp_int_digit* d = a->dp + n + 1;
+    volatile sp_int_digit mask;
+    sp_int_digit mk;
+    sp_int_digit z;
+    sp_size_t used = (sp_size_t)n;
+    sp_size_t um = (sp_size_t)-1;
+    unsigned int i;
+
+    /* d = a - m, in constant time. */
+#ifndef SQR_MUL_ASM
+    w = 0;
+    for (i = 0; i < n; i++) {
+        w    += a->dp[i];
+        w    -= m->dp[i];
+        d[i]  = (sp_int_digit)w;
+        w   >>= DIGIT_BIT;
+    }
+    /* Top digit of a with borrow - 0 when a >= m and -1 otherwise. */
+    w += a->dp[n];
+    mask = (sp_int_digit)0 - (sp_int_digit)(w < 0);
+#else
+    /* b is the borrow: 0 or 1. */
+    b = 0;
+    for (i = 0; i < n; i++) {
+        l = a->dp[i];
+        h = 0;
+        t = m->dp[i];
+        SP_ASM_SUBB_REG(l, h, t);
+        /* Can't borrow twice: when l - t borrowed then l is at least 1. */
+        SP_ASM_SUBB_REG(l, h, b);
+        d[i] = l;
+        /* h is 0 or -1. */
+        b = (sp_int_digit)0 - h;
+    }
+    /* Top digit of a with borrow - 0 when a >= m and -1 otherwise. */
+    l = a->dp[n];
+    h = 0;
+    SP_ASM_SUBB_REG(l, h, b);
+    mask = h;
+#endif
+
+    /* Select a when a < m (mask set) and a - m otherwise.
+     * Work down from the top digit, counting leading zero digits in constant
+     * time to determine used.
+     */
+    mk = mask;
+    for (i = n; i > 0; i--) {
+        t = d[i - 1] ^ ((a->dp[i - 1] ^ d[i - 1]) & mk);
+        a->dp[i - 1] = t;
+        /* z is all ones when digit is zero. */
+        z = (sp_int_digit)((sp_int_sdigit)(~t & (t - 1)) >>
+                           (SP_WORD_SIZE - 1));
+        /* um stays all ones while only zero digits seen. */
+        /* Sign extend - sp_size_t may be wider than sp_int_digit. */
+        um &= (sp_size_t)(sp_int_sdigit)z;
+        used = (sp_size_t)(used + um);
+    }
+    /* Result is less than m. */
+    a->used = used;
+#ifdef WOLFSSL_SP_INT_NEGATIVE
+    a->sign = MP_ZPOS;
+#endif
+}
+
 /* Final subtraction of Montgomery reduction in constant time: a = a mod m
  *
  * a is in n digits with the carry in digit n and is less than 2 * m.
@@ -17951,7 +18185,8 @@ int sp_sqrmod(const sp_int* a, const sp_int* m, sp_int* r)
  * @param [in]      m  SP integer that is the modulus.
  * @param [in]      n  Number of digits in modulus.
  */
-static void _sp_mont_red_sub_ct(sp_int* a, const sp_int* m, unsigned int n)
+WC_MAYBE_UNUSED static void _sp_mont_red_sub_ct(sp_int* a, const sp_int* m,
+    unsigned int n)
 {
     unsigned int i;
     volatile sp_int_digit vmask;
@@ -18056,10 +18291,7 @@ static int _sp_mont_red(sp_int* a, const sp_int* m, sp_int_digit mp, int ct)
     /* Adding numbers into m->used * 2 digits - zero out unused digits. */
 #ifndef WOLFSSL_NO_CT_OPS
     if (ct == 1) {
-        for (i = 0; i < (unsigned int)m->used * 2; i++) {
-            a->dp[i] &= (sp_int_digit)
-                (sp_int_sdigit)ctMaskIntGTE((int)a->used - 1, (int)i);
-        }
+        _sp_mont_red_zero_ct(a, (unsigned int)m->used * 2U);
     }
     else
 #endif /* !WOLFSSL_NO_CT_OPS */
@@ -18151,9 +18383,7 @@ static int _sp_mont_red(sp_int* a, const sp_int* m, sp_int_digit mp, int ct)
         /* 3. a >>= NumBits(m) */
         (void)sp_rshb(a, bits, a);
         /* 4. a = a mod m */
-        _sp_mont_red_sub_ct(a, m, m->used);
-        /* Constant time clamping. */
-        sp_clamp_ct(a);
+        _sp_mont_red_sub_ct_n(a, m, m->used);
     }
 
 
@@ -18180,9 +18410,18 @@ static int _sp_mont_red(sp_int* a, const sp_int* m, sp_int_digit mp, int ct)
 
 #ifndef WOLFSSL_NO_CT_OPS
     if (ct == 1) {
-        for (i = 0; i < (unsigned int)m->used * 2; i++) {
-            a->dp[i] &= (sp_int_digit)
-                (sp_int_sdigit)ctMaskIntGTE((int)a->used - 1, (int)i);
+        /* Constant digit counts for common sizes so loop is unrolled. */
+    #if !defined(WOLFSSL_SP_MATH) && defined(HAVE_ECC) && (SP_WORD_SIZE == 64)
+        if (m->used == 4) {
+            _sp_mont_red_zero_ct(a, 8U);
+        }
+        else if (m->used == 6) {
+            _sp_mont_red_zero_ct(a, 12U);
+        }
+        else
+    #endif
+        {
+            _sp_mont_red_zero_ct(a, (unsigned int)m->used * 2U);
         }
     }
     else
@@ -18261,9 +18500,7 @@ static int _sp_mont_red(sp_int* a, const sp_int* m, sp_int_digit mp, int ct)
 
         if (ct) {
             /* a = a mod m - constant time */
-            _sp_mont_red_sub_ct(a, m, 4);
-            /* Remove leading zeros. */
-            sp_clamp_ct(a);
+            _sp_mont_red_sub_ct_n(a, m, 4U);
         }
         else {
             /* Remove leading zeros. */
@@ -18333,9 +18570,7 @@ static int _sp_mont_red(sp_int* a, const sp_int* m, sp_int_digit mp, int ct)
 
         if (ct) {
             /* a = a mod m - constant time */
-            _sp_mont_red_sub_ct(a, m, 6);
-            /* Remove leading zeros. */
-            sp_clamp_ct(a);
+            _sp_mont_red_sub_ct_n(a, m, 6U);
         }
         else {
             /* Remove leading zeros. */
@@ -18404,9 +18639,7 @@ static int _sp_mont_red(sp_int* a, const sp_int* m, sp_int_digit mp, int ct)
 
         if (ct) {
             /* a = a mod m - constant time */
-            _sp_mont_red_sub_ct(a, m, m->used);
-            /* Remove leading zeros. */
-            sp_clamp_ct(a);
+            _sp_mont_red_sub_ct_n(a, m, m->used);
         }
         else {
             /* Remove leading zeros. */
@@ -18496,9 +18729,7 @@ static int _sp_mont_red(sp_int* a, const sp_int* m, sp_int_digit mp, int ct)
     else {
         (void)sp_rshb(a, bits, a);
         /* a = a mod m */
-        _sp_mont_red_sub_ct(a, m, m->used);
-        /* Constant time clamping. */
-        sp_clamp_ct(a);
+        _sp_mont_red_sub_ct_n(a, m, m->used);
     }
 
 #if 0
diff --git a/wolfcrypt/src/tfm.c b/wolfcrypt/src/tfm.c
index fd47a60bce2..6b5ca077faf 100644
--- a/wolfcrypt/src/tfm.c
+++ b/wolfcrypt/src/tfm.c
@@ -46,6 +46,14 @@
 #include <wolfssl/wolfcrypt/wolfmath.h> /* common functions */
 #include <wolfssl/wolfcrypt/logging.h>
 
+/* Constant time Montgomery reduction and multiplication - for exponentiation
+ * with timing resistance, when requested and for ECC with timing resistance as
+ * point arithmetic operates on secret values. */
+#if defined(TFM_TIMING_RESISTANT) || defined(WOLFSSL_MONT_RED_CT) || \
+    (defined(HAVE_ECC) && defined(ECC_TIMING_RESISTANT))
+    #define TFM_MONT_RED_CT_AVAILABLE
+#endif
+
 #ifdef WOLFSSL_ESPIDF
     #include <esp_log.h>
     #include <wolfssl/wolfcrypt/port/Espressif/esp32-crypt.h>
@@ -571,7 +579,8 @@ int fp_mul_2d(fp_int *a, int b, fp_int *c)
 /* generic PxQ multiplier */
 #if defined(HAVE_INTEL_MULX)
 
-WC_INLINE static int fp_mul_comba_mulx(fp_int *A, fp_int *B, fp_int *C)
+WC_INLINE static int fp_mul_comba_mulx(fp_int *A, fp_int *B, fp_int *C,
+    int clamp)
 
 {
    int       ix, iy, iz, pa;
@@ -585,9 +594,9 @@ WC_INLINE static int fp_mul_comba_mulx(fp_int *A, fp_int *B, fp_int *C)
    WC_ALLOC_VAR_EX(tmp, fp_int, 1, NULL, DYNAMIC_TYPE_BIGINT,
        return FP_MEM);
 
-   /* get size of output and trim */
+   /* get size of output and trim - unclamped result may use all digits */
    pa = A->used + B->used;
-   if (pa >= FP_SIZE) {
+   if (pa >= FP_SIZE + (clamp == 0)) {
       pa = FP_SIZE-1;
    }
 
@@ -602,7 +611,9 @@ WC_INLINE static int fp_mul_comba_mulx(fp_int *A, fp_int *B, fp_int *C)
 
   dst->used = pa;
   dst->sign = A->sign ^ B->sign;
-  fp_clamp(dst);
+  if (clamp) {
+    fp_clamp(dst);
+  }
   fp_copy(dst, C);
 
   WC_FREE_VAR_EX(tmp, NULL, DYNAMIC_TYPE_BIGINT);
@@ -611,8 +622,12 @@ WC_INLINE static int fp_mul_comba_mulx(fp_int *A, fp_int *B, fp_int *C)
 }
 #endif
 
-/*  C = (A * B)   */
-int fp_mul_comba(fp_int *A, fp_int *B, fp_int *C)
+/*  C = (A * B)
+ *
+ * When clamp is zero, the result has used set to A->used + B->used and is not
+ * clamped. The number of operations then doesn't depend on the value.
+ */
+static WC_INLINE int _fp_mul_comba(fp_int *A, fp_int *B, fp_int *C, int clamp)
 {
    int       ret = 0;
    int       ix, iy, iz, tx, ty, pa;
@@ -620,9 +635,14 @@ int fp_mul_comba(fp_int *A, fp_int *B, fp_int *C)
    fp_int    *dst;
    WC_DECLARE_VAR(tmp, fp_int, 1, 0);
 
-   if (A->used + B->used >= FP_SIZE) return FP_VAL;
+   /* Unclamped result, used in exponentiation, may use all digits. */
+   if (A->used + B->used >= FP_SIZE + (clamp == 0)) return FP_VAL;
 
-   IF_HAVE_INTEL_MULX(ret = fp_mul_comba_mulx(A, B, C), return ret) ;
+   /* MULX processes B in groups of four digits and writes a carry after each
+    * row - only use when the rounded up product fits. Lengths are public. */
+   if (A->used + ((B->used + 3) & ~3) <= FP_SIZE) {
+      IF_HAVE_INTEL_MULX(ret = fp_mul_comba_mulx(A, B, C, clamp), return ret) ;
+   }
 
    WC_ALLOC_VAR_EX(tmp, fp_int, 1, NULL, DYNAMIC_TYPE_BIGINT,
        return FP_MEM);
@@ -632,7 +652,7 @@ int fp_mul_comba(fp_int *A, fp_int *B, fp_int *C)
 
    /* get size of output and trim */
    pa = A->used + B->used;
-   if (pa >= FP_SIZE) {
+   if (pa >= FP_SIZE + (clamp == 0)) {
       pa = FP_SIZE-1;
    }
 
@@ -674,7 +694,9 @@ int fp_mul_comba(fp_int *A, fp_int *B, fp_int *C)
 
   /* warning: WOLFSSL_SP_INT_NEGATIVE may disable negative numbers */
   dst->sign = A->sign ^ B->sign;
-  fp_clamp(dst);
+  if (clamp) {
+    fp_clamp(dst);
+  }
   fp_copy(dst, C);
 
   /* Variables used but not seen by cppcheck. */
@@ -684,6 +706,120 @@ int fp_mul_comba(fp_int *A, fp_int *B, fp_int *C)
   return ret;
 }
 
+#if defined(HAVE_ECC) && defined(TFM_MONT_RED_CT_AVAILABLE)
+/* r = a * b where a and b have n digits - r has 2 * n digits.
+ *
+ * Number of operations depends only on n. r must not overlap a or b.
+ *
+ * @param [in]  a  First operand of n digits.
+ * @param [in]  b  Second operand of n digits.
+ * @param [in]  n  Number of digits in operands. Must be at least 1.
+ * @param [out] r  Product of 2 * n digits.
+ */
+static void fp_mul_comba_n(const fp_digit* a, const fp_digit* b, int n,
+                           fp_digit* r)
+{
+  int       ix, iy, iz, tx, ty;
+  fp_digit  c0, c1, c2;
+  const fp_digit *tmpx, *tmpy;
+
+  COMBA_START;
+  COMBA_CLEAR;
+
+  for (ix = 0; ix < 2 * n; ix++) {
+    /* get offsets into the two bignums */
+    ty = MIN(ix, n - 1);
+    tx = ix - ty;
+
+    /* setup temp aliases */
+    tmpx = a + tx;
+    tmpy = b + ty;
+
+    /* number of products in this column */
+    iy = MIN(n - tx, ty + 1);
+
+    /* execute loop */
+    COMBA_FORWARD;
+    for (iz = 0; iz < iy; ++iz) {
+      fp_digit _tmpx = *tmpx++;
+      fp_digit _tmpy = *tmpy--;
+      MULADD(_tmpx, _tmpy);
+    }
+
+    /* store term */
+    COMBA_STORE(r[ix]);
+  }
+  COMBA_FINI;
+
+  /* Variables used but not seen by cppcheck. */
+  (void)c0; (void)c1; (void)c2;
+}
+
+/* r = a * a where a has n digits - r has 2 * n digits.
+ *
+ * Number of operations depends only on n. r must not overlap a.
+ *
+ * @param [in]  a  Operand of n digits.
+ * @param [in]  n  Number of digits in operand. Must be at least 1.
+ * @param [out] r  Square of 2 * n digits.
+ */
+static void fp_sqr_comba_n(const fp_digit* a, int n, fp_digit* r)
+{
+  int       ix, iy, iz, tx, ty;
+  fp_digit  c0, c1, c2;
+#ifdef TFM_ISO
+  fp_word   tt = 0;
+#endif
+  const fp_digit *tmpx, *tmpy;
+
+  COMBA_START;
+  COMBA_CLEAR;
+
+  for (ix = 0; ix < 2 * n; ix++) {
+    /* get offsets into the two bignums */
+    ty = MIN(n - 1, ix);
+    tx = ix - ty;
+
+    /* setup temp aliases */
+    tmpx = a + tx;
+    tmpy = a + ty;
+
+    /* number of products in this column, each doubled */
+    iy = MIN(n - tx, ty + 1);
+    iy = MIN(iy, (ty - tx + 1) >> 1);
+
+    /* forward carries */
+    COMBA_FORWARD;
+
+    /* execute loop */
+    for (iz = 0; iz < iy; iz++) {
+      SQRADD2(*tmpx++, *tmpy--);
+    }
+
+    /* even columns have the square term in them */
+    if ((ix & 1) == 0) {
+      SQRADD(a[ix >> 1], a[ix >> 1]);
+    }
+
+    /* store term */
+    COMBA_STORE(r[ix]);
+  }
+  COMBA_FINI;
+
+  /* Variables used but not seen by cppcheck. */
+  (void)c0; (void)c1; (void)c2;
+#ifdef TFM_ISO
+  (void)tt;
+#endif
+}
+#endif
+
+/*  C = (A * B)   */
+int fp_mul_comba(fp_int *A, fp_int *B, fp_int *C)
+{
+  return _fp_mul_comba(A, B, C, 1);
+}
+
 /* a/b => cb + d == a */
 int fp_div(fp_int *a, fp_int *b, fp_int *c, fp_int *d)
 {
@@ -1456,9 +1592,9 @@ int fp_invmod_mont_ct(fp_int *a, fp_int *b, fp_int *c, fp_digit mp)
   for (i = 1; i < CT_INV_MOD_PRE_CNT; i++) {
     fp_init(&pre[i]);
     err |= fp_sqr(&pre[i-1], &pre[i]);
-    err |= fp_montgomery_reduce(&pre[i], b, mp);
+    err |= fp_montgomery_reduce_ex(&pre[i], b, mp, 1);
     err |= fp_mul(&pre[i], a, &pre[i]);
-    err |= fp_montgomery_reduce(&pre[i], b, mp);
+    err |= fp_montgomery_reduce_ex(&pre[i], b, mp, 1);
   }
 
   err |= fp_sub_d(b, 2, e);
@@ -1476,16 +1612,16 @@ int fp_invmod_mont_ct(fp_int *a, fp_int *b, fp_int *c, fp_digit mp)
 
     if ((j == CT_INV_MOD_PRE_CNT) || (!set && j > 0)) {
       err |= fp_mul(t, &pre[j-1], t);
-      err |= fp_montgomery_reduce(t, b, mp);
+      err |= fp_montgomery_reduce_ex(t, b, mp, 1);
       j = 0;
     }
     err |= fp_sqr(t, t);
-    err |= fp_montgomery_reduce(t, b, mp);
+    err |= fp_montgomery_reduce_ex(t, b, mp, 1);
     j += set;
   }
   if (j > 0) {
     err |= fp_mul(t, &pre[j-1], c);
-    err |= fp_montgomery_reduce(c, b, mp);
+    err |= fp_montgomery_reduce_ex(c, b, mp, 1);
   }
   else
     fp_copy(t, c);
@@ -1696,6 +1832,44 @@ static void fp_cond_copy(const fp_int* a, int copy, fp_int* r)
 }
 #endif
 
+/* Harden the exponentiation ladders against attackers that can observe fine
+ * grained timing or caches (e.g. a co-resident process). Montgomery reductions
+ * are constant time and working values are kept at the length of the modulus
+ * so that the multiply and square operations don't depend on secret values.
+ * Define WC_NO_MODEXP_TIMING_RESISTANT to disable, as for SP math.
+ */
+#ifndef WC_NO_MODEXP_TIMING_RESISTANT
+static int fp_montgomery_reduce_ct(fp_int *a, fp_int *m, fp_digit mp,
+                                   int fixed);
+/* Montgomery reduction during exponentiation - constant time, fixed length. */
+#define FP_EXPTMOD_RED(a, m, mp)        fp_montgomery_reduce_ct(a, m, mp, 1)
+/* Final Montgomery reduction - constant time and result clamped. */
+#define FP_EXPTMOD_RED_LAST(a, m, mp)   fp_montgomery_reduce_ct(a, m, mp, 0)
+static WC_INLINE int _fp_sqr_comba(fp_int *A, fp_int *B, int clamp);
+/* Square without removing leading zeros - number of leading zero digits is
+ * value dependent. Result is reduced with FP_EXPTMOD_RED. */
+static int fp_sqr_nc(fp_int *a, fp_int *r)
+{
+  /* Modulus may be FP_SIZE/2 digits - product uses all digits. */
+  if (a->used + a->used > FP_SIZE) {
+    return FP_VAL;
+  }
+  return _fp_sqr_comba(a, r, 0);
+}
+/* Multiply and square without removing leading zeros. */
+#define FP_EXPTMOD_MUL(a, b, r)         _fp_mul_comba(a, b, r, 0)
+#define FP_EXPTMOD_SQR(a, r)            fp_sqr_nc(a, r)
+/* Pad non-secret starting value, less than modulus, to length of modulus.
+ * Digits at and above used are zero. */
+#define FP_EXPTMOD_PAD(a, m)            ((a)->used = (m)->used)
+#else
+#define FP_EXPTMOD_RED(a, m, mp)        fp_montgomery_reduce(a, m, mp)
+#define FP_EXPTMOD_RED_LAST(a, m, mp)   fp_montgomery_reduce(a, m, mp)
+#define FP_EXPTMOD_PAD(a, m)            WC_DO_NOTHING
+#define FP_EXPTMOD_MUL(a, b, r)         fp_mul(a, b, r)
+#define FP_EXPTMOD_SQR(a, r)            fp_sqr(a, r)
+#endif
+
 #ifdef WC_RSA_NONBLOCK
 
 #ifdef WC_RSA_NONBLOCK_TIME
@@ -1848,6 +2022,11 @@ int fp_exptmod_nb(exptModNb_t* nb, fp_int* G, fp_int* X, fp_int* P, fp_int* Y)
        }
     }
 
+    /* Values are not secret - pad to modulus length so that all working
+     * values have the same length. */
+    FP_EXPTMOD_PAD(&nb->R[0], P);
+    FP_EXPTMOD_PAD(&nb->R[1], P);
+
     /* set initial mode and bit cnt */
     nb->bitcnt = 1;
     nb->buf    = 0;
@@ -1876,12 +2055,16 @@ int fp_exptmod_nb(exptModNb_t* nb, fp_int* G, fp_int* X, fp_int* P, fp_int* Y)
     FALL_THROUGH;
 
   case TFM_EXPTMOD_NB_MUL:
-    fp_mul(&nb->R[0], &nb->R[1], &nb->R[nb->y^1]);
+    err = FP_EXPTMOD_MUL(&nb->R[0], &nb->R[1], &nb->R[nb->y^1]);
+    if (err != FP_OKAY) {
+      nb->state = TFM_EXPTMOD_NB_INIT;
+      return err;
+    }
     nb->state = TFM_EXPTMOD_NB_MUL_RED;
     break;
 
   case TFM_EXPTMOD_NB_MUL_RED:
-    err = fp_montgomery_reduce(&nb->R[nb->y^1], P, nb->mp);
+    err = FP_EXPTMOD_RED(&nb->R[nb->y^1], P, nb->mp);
     if (err != FP_OKAY) {
       nb->state = TFM_EXPTMOD_NB_INIT;
       return err;
@@ -1891,16 +2074,16 @@ int fp_exptmod_nb(exptModNb_t* nb, fp_int* G, fp_int* X, fp_int* P, fp_int* Y)
 
   case TFM_EXPTMOD_NB_SQR:
   #ifdef WC_NO_CACHE_RESISTANT
-    err = fp_sqr(&nb->R[nb->y], &nb->R[nb->y]);
+    err = FP_EXPTMOD_SQR(&nb->R[nb->y], &nb->R[nb->y]);
   #elif defined(WC_NO_PTR_INT_CAST)
     fp_cond_copy(&nb->R[0], nb->y^1, &nb->R[2]);
     fp_cond_copy(&nb->R[1], nb->y,   &nb->R[2]);
-    err = fp_sqr(&nb->R[2], &nb->R[2]);
+    err = FP_EXPTMOD_SQR(&nb->R[2], &nb->R[2]);
   #else
     fp_copy((fp_int*) ( ((wc_ptr_t)&nb->R[0] & wc_off_on_addr[nb->y^1]) +
                         ((wc_ptr_t)&nb->R[1] & wc_off_on_addr[nb->y]) ),
             &nb->R[2]);
-    err = fp_sqr(&nb->R[2], &nb->R[2]);
+    err = FP_EXPTMOD_SQR(&nb->R[2], &nb->R[2]);
   #endif /* WC_NO_CACHE_RESISTANT */
     if (err != FP_OKAY) {
       nb->state = TFM_EXPTMOD_NB_INIT;
@@ -1912,13 +2095,13 @@ int fp_exptmod_nb(exptModNb_t* nb, fp_int* G, fp_int* X, fp_int* P, fp_int* Y)
 
   case TFM_EXPTMOD_NB_SQR_RED:
   #ifdef WC_NO_CACHE_RESISTANT
-    err = fp_montgomery_reduce(&nb->R[nb->y], P, nb->mp);
+    err = FP_EXPTMOD_RED(&nb->R[nb->y], P, nb->mp);
   #elif defined(WC_NO_PTR_INT_CAST)
-    err = fp_montgomery_reduce(&nb->R[2], P, nb->mp);
+    err = FP_EXPTMOD_RED(&nb->R[2], P, nb->mp);
     fp_cond_copy(&nb->R[2], nb->y^1, &nb->R[0]);
     fp_cond_copy(&nb->R[2], nb->y,   &nb->R[1]);
   #else
-    err = fp_montgomery_reduce(&nb->R[2], P, nb->mp);
+    err = FP_EXPTMOD_RED(&nb->R[2], P, nb->mp);
     fp_copy(&nb->R[2],
             (fp_int*) ( ((wc_ptr_t)&nb->R[0] & wc_off_on_addr[nb->y^1]) +
                         ((wc_ptr_t)&nb->R[1] & wc_off_on_addr[nb->y]) ) );
@@ -1933,7 +2116,7 @@ int fp_exptmod_nb(exptModNb_t* nb, fp_int* G, fp_int* X, fp_int* P, fp_int* Y)
 
   case TFM_EXPTMOD_NB_RED:
     /* final reduce */
-    err = fp_montgomery_reduce(&nb->R[0], P, nb->mp);
+    err = FP_EXPTMOD_RED_LAST(&nb->R[0], P, nb->mp);
     if (err != FP_OKAY) {
       nb->state = TFM_EXPTMOD_NB_INIT;
       return err;
@@ -2021,6 +2204,10 @@ static int _fp_exptmod_ct(fp_int * G, fp_int * X, int digits, fp_int * P,
       WC_FREE_VAR_EX(R, NULL, DYNAMIC_TYPE_BIGINT);
       return err;
   }
+  /* Values are not secret - pad to modulus length so that all working values
+   * have the same length. */
+  FP_EXPTMOD_PAD(&R[0], P);
+  FP_EXPTMOD_PAD(&R[1], P);
 
   /* for j = t-1 downto 0 do
         r_!k = R0*R1; r_k = r_k^2
@@ -2049,35 +2236,35 @@ static int _fp_exptmod_ct(fp_int * G, fp_int * X, int digits, fp_int * P,
 
 #ifdef WC_NO_CACHE_RESISTANT
     /* do ops */
-    err = fp_mul(&R[0], &R[1], &R[y^1]);
+    err = FP_EXPTMOD_MUL(&R[0], &R[1], &R[y^1]);
     if (err != FP_OKAY) {
       WC_FREE_VAR_EX(R, NULL, DYNAMIC_TYPE_BIGINT);
       return err;
     }
-    err = fp_montgomery_reduce(&R[y^1], P, mp);
+    err = FP_EXPTMOD_RED(&R[y^1], P, mp);
     if (err != FP_OKAY) {
       WC_FREE_VAR_EX(R, NULL, DYNAMIC_TYPE_BIGINT);
       return err;
     }
 
-    err = fp_sqr(&R[y], &R[y]);
+    err = FP_EXPTMOD_SQR(&R[y], &R[y]);
     if (err != FP_OKAY) {
       WC_FREE_VAR_EX(R, NULL, DYNAMIC_TYPE_BIGINT);
       return err;
     }
-    err = fp_montgomery_reduce(&R[y], P, mp);
+    err = FP_EXPTMOD_RED(&R[y], P, mp);
     if (err != FP_OKAY) {
       WC_FREE_VAR_EX(R, NULL, DYNAMIC_TYPE_BIGINT);
       return err;
     }
 #else
     /* do ops */
-    err = fp_mul(&R[0], &R[1], &R[2]);
+    err = FP_EXPTMOD_MUL(&R[0], &R[1], &R[2]);
     if (err != FP_OKAY) {
       WC_FREE_VAR_EX(R, NULL, DYNAMIC_TYPE_BIGINT);
       return err;
     }
-    err = fp_montgomery_reduce(&R[2], P, mp);
+    err = FP_EXPTMOD_RED(&R[2], P, mp);
     if (err != FP_OKAY) {
       WC_FREE_VAR_EX(R, NULL, DYNAMIC_TYPE_BIGINT);
       return err;
@@ -2105,12 +2292,12 @@ static int _fp_exptmod_ct(fp_int * G, fp_int * X, int digits, fp_int * P,
                         ((wc_ptr_t)&R[1] & wc_off_on_addr[y]) ),
             &R[2]);
 #endif
-    err = fp_sqr(&R[2], &R[2]);
+    err = FP_EXPTMOD_SQR(&R[2], &R[2]);
     if (err != FP_OKAY) {
       WC_FREE_VAR_EX(R, NULL, DYNAMIC_TYPE_BIGINT);
       return err;
     }
-    err = fp_montgomery_reduce(&R[2], P, mp);
+    err = FP_EXPTMOD_RED(&R[2], P, mp);
     if (err != FP_OKAY) {
       WC_FREE_VAR_EX(R, NULL, DYNAMIC_TYPE_BIGINT);
       return err;
@@ -2126,7 +2313,7 @@ static int _fp_exptmod_ct(fp_int * G, fp_int * X, int digits, fp_int * P,
 #endif /* WC_NO_CACHE_RESISTANT */
   }
 
-   err = fp_montgomery_reduce(&R[0], P, mp);
+   err = FP_EXPTMOD_RED_LAST(&R[0], P, mp);
    fp_copy(&R[0], Y);
    WC_FREE_VAR_EX(R, NULL, DYNAMIC_TYPE_BIGINT);
 
@@ -2213,6 +2400,10 @@ static int _fp_exptmod_ct(fp_int * G, fp_int * X, int digits, fp_int * P,
       WC_FREE_VAR_EX(R, NULL, DYNAMIC_TYPE_BIGINT);
       return err;
   }
+  /* Values are not secret - pad to modulus length so that all working values
+   * have the same length. */
+  FP_EXPTMOD_PAD(&R[0], P);
+  FP_EXPTMOD_PAD(&R[1], P);
 
   /* for j = t-1 downto 0 do
         r_!k = R0*R1; r_k = r_k^2
@@ -2240,12 +2431,12 @@ static int _fp_exptmod_ct(fp_int * G, fp_int * X, int digits, fp_int * P,
     buf <<= (fp_digit)1;
 
     /* do ops */
-    err = fp_mul(&R[0], &R[1], &R[2]);
+    err = FP_EXPTMOD_MUL(&R[0], &R[1], &R[2]);
     if (err != FP_OKAY) {
       WC_FREE_VAR_EX(R, NULL, DYNAMIC_TYPE_BIGINT);
       return err;
     }
-    err = fp_montgomery_reduce(&R[2], P, mp);
+    err = FP_EXPTMOD_RED(&R[2], P, mp);
     if (err != FP_OKAY) {
       WC_FREE_VAR_EX(R, NULL, DYNAMIC_TYPE_BIGINT);
       return err;
@@ -2262,12 +2453,12 @@ static int _fp_exptmod_ct(fp_int * G, fp_int * X, int digits, fp_int * P,
                         ((wc_ptr_t)&R[1] & wc_off_on_addr[y]) ),
             &R[3]);
 #endif
-    err = fp_sqr(&R[3], &R[3]);
+    err = FP_EXPTMOD_SQR(&R[3], &R[3]);
     if (err != FP_OKAY) {
       WC_FREE_VAR_EX(R, NULL, DYNAMIC_TYPE_BIGINT);
       return err;
     }
-    err = fp_montgomery_reduce(&R[3], P, mp);
+    err = FP_EXPTMOD_RED(&R[3], P, mp);
     if (err != FP_OKAY) {
       WC_FREE_VAR_EX(R, NULL, DYNAMIC_TYPE_BIGINT);
       return err;
@@ -2275,7 +2466,7 @@ static int _fp_exptmod_ct(fp_int * G, fp_int * X, int digits, fp_int * P,
     fp_copy_2_ct(&R[2], &R[3], &R[0], &R[1], y, P->used);
   }
 
-  err = fp_montgomery_reduce(&R[0], P, mp);
+  err = FP_EXPTMOD_RED_LAST(&R[0], P, mp);
   fp_copy(&R[0], Y);
   WC_FREE_VAR_EX(R, NULL, DYNAMIC_TYPE_BIGINT);
   return err;
@@ -2569,7 +2760,108 @@ static int _fp_exptmod_nct(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
 }
 
 
-#ifdef TFM_TIMING_RESISTANT
+#if defined(TFM_TIMING_RESISTANT) && !defined(WC_NO_MODEXP_TIMING_RESISTANT)
+/* r = (2 * a) mod m in constant time.
+ *
+ * a must be less than m and have used equal to m->used.
+ * r has used set to m->used and is NOT clamped.
+ */
+static void fp_dbl_mod_ct(const fp_int* a, const fp_int* m, fp_int* r)
+{
+  int      i;
+  fp_digit c = 0;
+  fp_digit mask;
+  fp_sword w = 0;
+
+  /* r = (a << 1) - m, with top bit of a in c. */
+  for (i = 0; i < m->used; i++) {
+    fp_digit v = a->dp[i];
+    w        += (fp_digit)((v << 1) | c);
+    w        -= m->dp[i];
+    r->dp[i]  = (fp_digit)w;
+    w       >>= DIGIT_BIT;
+    c         = v >> (DIGIT_BIT - 1);
+  }
+  /* 2 * a < 2 * m - negative only when 2 * a was less than m. */
+  w += c;
+  mask = (fp_digit)0 - (fp_digit)(w < 0);
+  /* Add modulus back when negative. */
+  w = 0;
+  for (i = 0; i < m->used; i++) {
+    w        += r->dp[i];
+    w        += m->dp[i] & mask;
+    r->dp[i]  = (fp_digit)w;
+    w       >>= DIGIT_BIT;
+  }
+  r->used = m->used;
+  r->sign = FP_ZPOS;
+}
+
+/* y = 2**x (mod b) in constant time.
+ *
+ * Processes one exponent bit at a time with result in Montgomery form.
+ * Doubling is always calculated and selected in constant time when the
+ * exponent bit is set. Working values are kept at the length of the modulus.
+ */
+static int _fp_exptmod_base_2(fp_int * X, int digits, fp_int * P,
+                              fp_int * Y)
+{
+  fp_digit mp;
+  int      err, i, j;
+  fp_int*  tr;
+  fp_int*  t;
+  WC_DECLARE_VAR(res, fp_int, 2, 0);
+
+  /* now setup montgomery  */
+  if ((err = fp_montgomery_setup(P, &mp)) != FP_OKAY) {
+    return err;
+  }
+
+  WC_ALLOC_VAR_EX(res, fp_int, 2, NULL, DYNAMIC_TYPE_TMP_BUFFER,
+      return FP_MEM);
+  tr = &res[0];
+  t  = &res[1];
+  fp_init(tr);
+  fp_init(t);
+
+  /* tr = R mod m = 1 in Montgomery form - not secret. */
+  err = fp_montgomery_calc_normalization(tr, P);
+  if (err == FP_OKAY) {
+    FP_EXPTMOD_PAD(tr, P);
+  }
+
+  for (i = digits * (int)DIGIT_BIT - 1; (err == FP_OKAY) && (i >= 0); i--) {
+    volatile fp_digit mask;
+
+    err = FP_EXPTMOD_SQR(tr, tr);
+    if (err == FP_OKAY) {
+      err = FP_EXPTMOD_RED(tr, P, mp);
+    }
+    if (err == FP_OKAY) {
+      fp_dbl_mod_ct(tr, P, t);
+      /* tr = t when exponent bit is set. */
+      mask = (fp_digit)0 -
+             ((X->dp[i / (int)DIGIT_BIT] >> (i % (int)DIGIT_BIT)) & 1);
+      for (j = 0; j < P->used; j++) {
+        tr->dp[j] ^= (tr->dp[j] ^ t->dp[j]) & mask;
+      }
+    }
+  }
+
+  /* convert out of Montgomery form */
+  if (err == FP_OKAY) {
+    err = FP_EXPTMOD_RED_LAST(tr, P, mp);
+  }
+  if (err == FP_OKAY) {
+    fp_copy(tr, Y);
+  }
+
+  ForceZero(res, sizeof(fp_int) * 2);
+  WC_FREE_VAR_EX(res, NULL, DYNAMIC_TYPE_TMP_BUFFER);
+  return err;
+}
+
+#elif defined(TFM_TIMING_RESISTANT)
 #if DIGIT_BIT <= 16
     #define WINSIZE    2
     #define WINMASK    0x3
@@ -3358,8 +3650,12 @@ int fp_sqr(fp_int *A, fp_int *B)
   return err;
 }
 
-/* generic comba squarer */
-int fp_sqr_comba(fp_int *A, fp_int *B)
+/* generic comba squarer
+ *
+ * When clamp is zero, the result has used set to A->used * 2 and is not
+ * clamped. The number of operations then doesn't depend on the value.
+ */
+static WC_INLINE int _fp_sqr_comba(fp_int *A, fp_int *B, int clamp)
 {
   int       pa, ix, iz;
   fp_digit  c0, c1, c2;
@@ -3372,9 +3668,9 @@ int fp_sqr_comba(fp_int *A, fp_int *B)
    WC_ALLOC_VAR_EX(tmp, fp_int, 1, NULL, DYNAMIC_TYPE_BIGINT,
        return FP_MEM);
 
-  /* get size of output and trim */
+  /* get size of output and trim - unclamped result may use all digits */
   pa = A->used + A->used;
-  if (pa >= FP_SIZE) {
+  if (pa >= FP_SIZE + (clamp == 0)) {
      pa = FP_SIZE-1;
   }
 
@@ -3436,7 +3732,9 @@ int fp_sqr_comba(fp_int *A, fp_int *B)
 
   /* setup dest */
   dst->used = pa;
-  fp_clamp (dst);
+  if (clamp) {
+    fp_clamp (dst);
+  }
   if (dst != B) {
      fp_copy(dst, B);
   }
@@ -3451,6 +3749,12 @@ int fp_sqr_comba(fp_int *A, fp_int *B)
   return FP_OKAY;
 }
 
+/* generic comba squarer */
+int fp_sqr_comba(fp_int *A, fp_int *B)
+{
+  return _fp_sqr_comba(A, B, 1);
+}
+
 int fp_cmp(fp_int *a, fp_int *b)
 {
    if (a->sign == FP_NEG && b->sign == FP_ZPOS) {
@@ -3517,6 +3821,36 @@ int fp_cmp_mag(fp_int *a, fp_int *b)
 }
 
 
+/* Unsigned comparison of the first n digits in constant time.
+ *
+ * Digits at and above used are treated as zero.
+ *
+ * @param [in] a  First number to compare.
+ * @param [in] b  Second number to compare.
+ * @param [in] n  Number of digits to compare. Must not exceed FP_SIZE.
+ * @return  FP_GT when a is greater than b.
+ * @return  FP_LT when a is less than b.
+ * @return  FP_EQ when a is equal to b.
+ */
+int fp_cmp_ct(const fp_int *a, const fp_int *b, int n)
+{
+   int ret = FP_EQ;
+   int i;
+   volatile int mask = -1;
+
+   for (i = n - 1; i >= 0; i--) {
+      fp_digit ad = a->dp[i] & ((fp_digit)0 - (fp_digit)(i < a->used));
+      fp_digit bd = b->dp[i] & ((fp_digit)0 - (fp_digit)(i < b->used));
+
+      ret |= mask & ((0 - (ad < bd)) & FP_LT);
+      mask &= 0 - (ret == FP_EQ);
+      ret |= mask & ((0 - (ad > bd)) & FP_GT);
+      mask &= 0 - (ret == FP_EQ);
+   }
+
+   return ret;
+}
+
 /* sets up the montgomery reduction */
 int fp_montgomery_setup(fp_int *a, fp_digit *rho)
 {
@@ -3698,6 +4032,267 @@ static int fp_montgomery_reduce_mulx(fp_int *a, fp_int *m, fp_digit mp, int ct)
 }
 #endif
 
+/* Number of digits the integer can hold - an alt_fp_int under ALT_ECC_SIZE
+ * has fewer digits than FP_SIZE. */
+#if defined(ALT_ECC_SIZE) || defined(HAVE_WOLF_BIGINT)
+    #define FP_DIGITS_CAP(a)    ((a)->size)
+#else
+    #define FP_DIGITS_CAP(a)    FP_SIZE
+#endif
+
+#ifdef TFM_MONT_RED_CT_AVAILABLE
+/* Remove leading zero digits in constant time.
+ *
+ * All a->used digits are always checked - no early exit.
+ *
+ * @param [in, out] a  Number to clamp.
+ */
+static void fp_clamp_ct(fp_int* a)
+{
+  int      i;
+  int      used = 0;
+
+  for (i = 0; i < a->used; i++) {
+    fp_digit d = a->dp[i];
+    /* Top bit set only when digit is non-zero - calculated in fp_digit width
+     * as small digit types are promoted to int. */
+    fp_digit nzd = (fp_digit)(d | (fp_digit)((fp_digit)0 - d)) >>
+                   (DIGIT_BIT - 1);
+    /* All bits set when digit is non-zero. */
+    int nz = (int)0 - (int)nzd;
+    /* When digit is non-zero, used is this digit's index + 1. */
+    used ^= (used ^ (i + 1)) & nz;
+  }
+  a->used = used;
+  /* Zero has positive sign - sign of a is always positive here. */
+  a->sign = FP_ZPOS;
+}
+
+/* computes x/R == x (mod N) via Montgomery Reduction in constant time.
+ *
+ * Separate from fp_montgomery_reduce_ex() so that the variable time version is
+ * not slowed down. The operations performed depend only on the modulus:
+ *  - all m->used * 2 digits of a are used (digits above used are zero),
+ *  - the carry of each round is added to the next digit only, with the
+ *    overflow kept and added in the next round,
+ *  - the final subtraction is always performed and undone when negative.
+ *
+ * When fixed is set, the result is not clamped and has used set to m->used.
+ * This keeps the length of working values in an exponentiation independent of
+ * their value. Such a result must only be used as an intermediate value.
+ */
+static int fp_montgomery_reduce_ct(fp_int *a, fp_int *m, fp_digit mp,
+                                   int fixed)
+{
+   WC_DECLARE_VAR(c, fp_digit, FP_SIZE+1, 0);
+   fp_digit *_c, *tmpm, mu = 0;
+   fp_digit h = 0;
+   int      oldused, x, y, pa;
+
+   /* bail if too large */
+   if (m->used > (FP_SIZE/2)) {
+      (void)mu;                     /* shut up compiler */
+      return FP_VAL;
+   }
+
+   WC_ALLOC_VAR_EX(c, fp_digit, (FP_SIZE+1), NULL, DYNAMIC_TYPE_BIGINT,
+       return FP_MEM);
+
+   pa = m->used;
+   /* a must hold the m->used + 1 digit intermediate result. Capacity is
+    * public - an alt_fp_int under ALT_ECC_SIZE may hold fewer than
+    * 2 * m->used digits. */
+   if (FP_DIGITS_CAP(a) < pa + 1) {
+      WC_FREE_VAR_EX(c, NULL, DYNAMIC_TYPE_BIGINT);
+      return FP_VAL;
+   }
+
+   /* now zero the buff - only 2 * pa + 1 digits are used */
+   XMEMSET(c, 0, sizeof(fp_digit) * (size_t)(2 * pa + 1));
+
+   /* copy the input - fixed length bounded by the capacity of a. Missing high
+    * digits are zero in c. */
+   oldused = MIN(m->used * 2, FP_DIGITS_CAP(a));
+   for (x = 0; x < oldused; x++) {
+       c[x] = a->dp[x];
+   }
+   MONT_START;
+
+   for (x = 0; x < pa; x++) {
+       fp_digit cy = 0;
+       /* get Mu for this round */
+       LOOP_START;
+       _c   = c + x;
+       tmpm = m->dp;
+       y = 0;
+#if defined(INNERMUL8)
+        for (; y < (pa & ~7); y += 8) {
+              INNERMUL8 ;
+              _c   += 8;
+              tmpm += 8;
+           }
+#endif
+       for (; y < pa; y++) {
+          INNERMUL;
+          ++_c;
+       }
+       LOOP_END;
+       /* Add carry and the overflow of the previous round into the next
+        * digit. Overflow is at most one and is added in the next round - the
+        * same number of operations every round. */
+       {
+           fp_word t = (fp_word)c[x + pa] + cy + h;
+           c[x + pa] = (fp_digit)t;
+           h = (fp_digit)(t >> DIGIT_BIT);
+       }
+  }
+  c[2 * pa] = h;
+
+  /* now copy out */
+  _c   = c + pa;
+  tmpm = a->dp;
+  for (x = 0; x < pa+1; x++) {
+     *tmpm++ = *_c++;
+  }
+
+  /* zero any excess digits on the destination that we didn't write to */
+  for (; x < oldused; x++) {
+     *tmpm++ = 0;
+  }
+
+  MONT_FINI;
+
+  /* a is less than 2 * m with all m->used + 1 digits set.
+   * Always subtract m and add it back when the result is negative. */
+  {
+    fp_sword w = 0;
+    fp_digit mask;
+
+    for (x = 0; x < pa; x++) {
+      w        += a->dp[x];
+      w        -= m->dp[x];
+      a->dp[x]  = (fp_digit)w;
+      w       >>= DIGIT_BIT;
+    }
+    w += a->dp[pa];
+    /* Negative when a was less than m - add modulus back. */
+    mask = (fp_digit)0 - (fp_digit)(w < 0);
+    w = 0;
+    for (x = 0; x < pa; x++) {
+      w        += a->dp[x];
+      w        += m->dp[x] & mask;
+      a->dp[x]  = (fp_digit)w;
+      w       >>= DIGIT_BIT;
+    }
+  }
+  /* Result fits in m->used digits - clear the carry digit as digits above used
+   * must be zero. */
+  a->dp[pa] = 0;
+  a->used = pa;
+  a->sign = FP_ZPOS;
+  if (!fixed) {
+    fp_clamp_ct(a);
+  }
+
+  WC_FREE_VAR_EX(c, NULL, DYNAMIC_TYPE_BIGINT);
+  return FP_OKAY;
+}
+#endif /* TFM_MONT_RED_CT_AVAILABLE */
+
+#ifdef HAVE_ECC
+/* r = a * b / R mod m - Montgomery multiplication.
+ *
+ * Constant time when TFM_MONT_RED_CT_AVAILABLE is defined:
+ * operands are padded to the length of the modulus, the product is not
+ * clamped and the result is copied with a fixed number of digits.
+ * Called for each field multiplication in ECC point arithmetic - only the
+ * digits needed are written and zeroized.
+ *
+ * a * b must be less than m * R - true when a is less than m and b is less
+ * than R.
+ *
+ * @param [in]  a   First operand. Has no more digits than m.
+ * @param [in]  b   Second operand. Has no more digits than m.
+ * @param [in]  m   Modulus. Must be odd.
+ * @param [in]  mp  Montgomery multiplier for m.
+ * @param [out] r   Result. May be the same as a or b.
+ * @return  FP_OKAY on success.
+ * @return  FP_VAL when an operand is longer than m or m is too long.
+ * @return  FP_MEM when dynamic memory allocation fails.
+ */
+int fp_mont_mul_ct(fp_int* a, fp_int* b, fp_int* m, fp_digit mp, fp_int* r)
+{
+  int err = FP_OKAY;
+#ifdef TFM_MONT_RED_CT_AVAILABLE
+  int i;
+  int n = m->used;
+  fp_digit* ta;
+  fp_digit* tb;
+  WC_DECLARE_VAR(t, fp_int, 2, 0);
+
+  /* Lengths and capacities are public - operands are less than the
+   * modulus. */
+  if ((n == 0) || (a->used > n) || (b->used > n) || (n + n > FP_SIZE) ||
+      (FP_DIGITS_CAP(r) < n)) {
+    return FP_VAL;
+  }
+
+  WC_ALLOC_VAR_EX(t, fp_int, 2, NULL, DYNAMIC_TYPE_BIGINT, return FP_MEM);
+
+  /* Pad operands to the length of the modulus in the digits of t[0]. Digits
+   * read depend only on the capacities - digits at and above used are masked
+   * out. */
+  ta = t[0].dp;
+  tb = t[0].dp + n;
+  for (i = 0; i < n; i++) {
+    fp_digit d;
+
+    d = (i < FP_DIGITS_CAP(a)) ? a->dp[i] : 0;
+    ta[i] = d & ((fp_digit)0 - (fp_digit)(i < a->used));
+    d = (i < FP_DIGITS_CAP(b)) ? b->dp[i] : 0;
+    tb[i] = d & ((fp_digit)0 - (fp_digit)(i < b->used));
+  }
+
+  /* t[1] = a * b, not clamped, and reduced in constant time.
+   * Same pointer is public - square when operands are the same. */
+  if (a == b) {
+    fp_sqr_comba_n(ta, n, t[1].dp);
+  }
+  else {
+    fp_mul_comba_n(ta, tb, n, t[1].dp);
+  }
+  t[1].used = n + n;
+  t[1].sign = FP_ZPOS;
+#if defined(ALT_ECC_SIZE) || defined(HAVE_WOLF_BIGINT)
+  t[1].size = FP_SIZE;
+#endif
+  err = fp_montgomery_reduce_ct(&t[1], m, mp, 0);
+  if (err == FP_OKAY) {
+    /* Copy a fixed number of digits and clear the rest of r's digits. */
+    for (i = 0; i < n; i++) {
+      r->dp[i] = t[1].dp[i];
+    }
+    XMEMSET(r->dp + n, 0,
+            sizeof(fp_digit) * (size_t)(FP_DIGITS_CAP(r) - n));
+    r->used = t[1].used;
+    r->sign = FP_ZPOS;
+  }
+
+  /* Only the first 2 * n digits of each temporary were used. */
+  ForceZero(t[0].dp, sizeof(fp_digit) * (size_t)(n + n));
+  ForceZero(t[1].dp, sizeof(fp_digit) * (size_t)(n + n));
+  WC_FREE_VAR_EX(t, NULL, DYNAMIC_TYPE_BIGINT);
+#else
+  err = fp_mul(a, b, r);
+  if (err == FP_OKAY) {
+    err = fp_montgomery_reduce(r, m, mp);
+  }
+#endif
+  return err;
+}
+#endif /* HAVE_ECC */
+
+
 /* computes x/R == x (mod N) via Montgomery Reduction */
 int fp_montgomery_reduce_ex(fp_int *a, fp_int *m, fp_digit mp, int ct)
 {
@@ -3705,6 +4300,11 @@ int fp_montgomery_reduce_ex(fp_int *a, fp_int *m, fp_digit mp, int ct)
    fp_digit *_c, *tmpm, mu = 0;
    int      oldused, x, y, pa, err = 0;
 
+#ifdef TFM_MONT_RED_CT_AVAILABLE
+   if (ct) {
+      return fp_montgomery_reduce_ct(a, m, mp, 0);
+   }
+#endif
    IF_HAVE_INTEL_MULX(err=fp_montgomery_reduce_mulx(a, m, mp, ct), return err) ;
    (void)err;
 
@@ -3807,7 +4407,12 @@ int fp_montgomery_reduce_ex(fp_int *a, fp_int *m, fp_digit mp, int ct)
 
 int fp_montgomery_reduce(fp_int *a, fp_int *m, fp_digit mp)
 {
+#ifdef WOLFSSL_MONT_RED_CT
   return fp_montgomery_reduce_ex(a, m, mp, 1);
+#else
+  /* Constant time only when requested with fp_montgomery_reduce_ex(). */
+  return fp_montgomery_reduce_ex(a, m, mp, 0);
+#endif
 }
 
 int fp_read_unsigned_bin(fp_int *a, const unsigned char *b, int c)
@@ -4627,6 +5232,25 @@ int mp_exptmod_nct (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
 
 
 /* compare two ints (signed)*/
+/* Compare first n digits of a and b in constant time. */
+int mp_cmp_ct(const mp_int *a, const mp_int *b, int n)
+{
+  /* Same as SP math: both NULL equal, NULL is less than non-NULL. */
+  if (a == b) {
+    return FP_EQ;
+  }
+  if (a == NULL) {
+    return FP_LT;
+  }
+  if (b == NULL) {
+    return FP_GT;
+  }
+  if ((n < 0) || (n > FP_SIZE)) {
+    n = FP_SIZE;
+  }
+  return fp_cmp_ct(a, b, n);
+}
+
 int mp_cmp (const mp_int * a, const mp_int * b)
 {
   return fp_cmp((mp_int *)a, (mp_int *)b);
diff --git a/wolfssl/wolfcrypt/sp_int.h b/wolfssl/wolfcrypt/sp_int.h
index 4149f1c0e62..10e8cc922bb 100644
--- a/wolfssl/wolfcrypt/sp_int.h
+++ b/wolfssl/wolfcrypt/sp_int.h
@@ -1144,6 +1144,7 @@ typedef sp_int_digit mp_digit;
     #define sp_invmod wc_sp_invmod
     #if defined(WOLFSSL_SP_MATH_ALL) && defined(HAVE_ECC)
         #define sp_invmod_mont_ct wc_sp_invmod_mont_ct
+        #define sp_mont_mul_ct wc_sp_mont_mul_ct
     #endif
 
     #define sp_exptmod_ex wc_sp_exptmod_ex
@@ -1300,6 +1301,8 @@ MP_API int sp_invmod(const sp_int* a, const sp_int* m, sp_int* r);
 #if defined(WOLFSSL_SP_MATH_ALL) && defined(HAVE_ECC)
 MP_API int sp_invmod_mont_ct(const sp_int* a, const sp_int* m, sp_int* r,
     sp_int_digit mp);
+WOLFSSL_LOCAL int sp_mont_mul_ct(const sp_int* a, const sp_int* b,
+    const sp_int* m, sp_int_digit mp, sp_int* r);
 #endif
 
 MP_API int sp_exptmod_ex(const sp_int* b, const sp_int* e, int digits,
@@ -1432,6 +1435,7 @@ WOLFSSL_LOCAL void sp_memzero_check(sp_int* sp);
 #define mp_mul                              sp_mul
 #define mp_invmod                           sp_invmod
 #define mp_invmod_mont_ct                   sp_invmod_mont_ct
+#define mp_mont_mul_ct                      sp_mont_mul_ct
 #define mp_exptmod_ex                       sp_exptmod_ex
 #define mp_exptmod_nct                      sp_exptmod_nct
 #define mp_div_2d                           sp_div_2d
diff --git a/wolfssl/wolfcrypt/tfm.h b/wolfssl/wolfcrypt/tfm.h
index a17a13dcfc6..fa310d67df9 100644
--- a/wolfssl/wolfcrypt/tfm.h
+++ b/wolfssl/wolfcrypt/tfm.h
@@ -555,6 +555,9 @@ int fp_cmp(fp_int *a, fp_int *b);
 /* unsigned comparison */
 int fp_cmp_mag(fp_int *a, fp_int *b);
 
+/* unsigned comparison of n digits in constant time */
+int fp_cmp_ct(const fp_int *a, const fp_int *b, int n);
+
 /* power of 2 operations */
 void fp_div_2d(fp_int *a, int b, fp_int *c, fp_int *d);
 void fp_mod_2d(fp_int *a, int b, fp_int *c);
@@ -651,6 +654,11 @@ int fp_montgomery_calc_normalization(fp_int *a, fp_int *b);
 /* computes x/R == x (mod N) via Montgomery Reduction */
 int fp_montgomery_reduce(fp_int *a, fp_int *m, fp_digit mp);
 int fp_montgomery_reduce_ex(fp_int *a, fp_int *m, fp_digit mp, int ct);
+#ifdef HAVE_ECC
+/* r = a * b / R mod m - Montgomery multiplication in constant time */
+WOLFSSL_LOCAL int fp_mont_mul_ct(fp_int* a, fp_int* b, fp_int* m, fp_digit mp,
+    fp_int* r);
+#endif
 
 /* d = a**b (mod c) */
 int fp_exptmod(fp_int *G, fp_int *X, fp_int *P, fp_int *Y);
@@ -853,7 +861,7 @@ MP_API int  mp_2expt(mp_int* a, int b);
 MP_API int  mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d);
 
 MP_API int  mp_cmp(const mp_int *a, const mp_int *b);
-#define mp_cmp_ct(a, b, n) mp_cmp(a, b)
+MP_API int  mp_cmp_ct(const mp_int *a, const mp_int *b, int n);
 MP_API int  mp_cmp_d(mp_int *a, mp_digit b);
 
 MP_API int  mp_unsigned_bin_size(const mp_int * a);
@@ -889,6 +897,9 @@ MP_API int mp_radix_size (mp_int * a, int radix, int *size);
 
 #define mp_montgomery_reduce_ct(a, m, mp) \
     mp_montgomery_reduce_ex(a, m, mp, 1)
+#ifdef HAVE_ECC
+#define mp_mont_mul_ct(a, b, m, mp, r)  fp_mont_mul_ct(a, b, m, mp, r)
+#endif
 MP_API int mp_montgomery_reduce(fp_int *a, fp_int *m, fp_digit mp);
 MP_API int mp_montgomery_reduce_ex(fp_int *a, fp_int *m, fp_digit mp, int ct);
 MP_API int mp_montgomery_setup(fp_int *a, fp_digit *rho);