ECC Deep Dive, in depth: the arithmetic engine under elliptic-curve crypto, from field reduction and Jacobian formulas to batch inversion, wNAF and constant time

The cost model: M, S and I

Engine designers count field operations. M is a multiplication mod p, S is a squaring (about 0.8 M with dedicated code), and I is an inversion. With Fermat's little theorem, x⁻¹ = x^p−2, which is roughly 255 squarings and a few dozen multiplications. A constant-time binary extended Euclid variant is faster but still costs many M. Additions and subtractions are nearly free.

The textbook affine formula needs one inversion per point operation, and a scalar multiplication does hundreds of them. So engines keep a denominator in the point, never divide inside the loop, and pay one inversion at the end.

The cost model: M, S and I

Engine designers count field operations. M is a multiplication mod p, S is a squaring (about 0.8 M with dedicated code), and I is an inversion. With Fermat's little theorem, x⁻¹ = x^p−2, which is roughly 255 squarings and a few dozen multiplications. A constant-time binary extended Euclid variant is faster but still costs many M. Additions and subtractions are nearly free.

The textbook affine formula needs one inversion per point operation, and a scalar multiplication does hundreds of them. So engines keep a denominator in the point, never divide inside the loop, and pay one inversion at the end.

Field arithmetic and special primes

A 256-bit field element is stored as limbs: four 64-bit words, or five 51-bit words for Curve25519 so that partial products and their sums fit in 128-bit accumulators without immediate carries. Multiplying two elements gives a 512-bit product that must be reduced mod p. For a general p that needs Montgomery multiplication, which replaces division by shifts using a precomputed constant. Standard curves choose p so that reduction is cheap instead:

• Curve25519 uses p = 2²⁵⁵ − 19. Since 2²⁵⁵ ≡ 19 (mod p), the high half of a product folds back as 19 times itself, and a final conditional subtraction makes the result canonical.
• P-256 uses a Solinas prime, 2²⁵⁶ − 2²²⁴ + 2¹⁹² + 2⁹⁶ − 1. Its reduction is a fixed sum and difference of 32-bit-word rearrangements of the high half, with no multiplications.

def reduce25519(x):                  # x < 2**510, e.g. a product of two field elements
    q = 2**255 - 19
    while x >> 255:
        x = (x & (2**255 - 1)) + 19 * (x >> 255)   # 2**255 = 19 (mod q)
    return x - q if x >= q else x                    # canonical: 0 <= x < q

This fold matched x % q on 10,000 random products. In production the while loop becomes a fixed two-step fold, and the final subtraction becomes a masked select, so the work never depends on the data. Lazy reduction lets values exceed p between operations, which saves work. Every intermediate bound must be proved, because one overflowing carry gives wrong answers for rare inputs.

Jacobian coordinates

Jacobian coordinates represent the affine point (x, y) as (X, Y, Z) with x = X/Z² and y = Y/Z³. Each point has many representations, the point at infinity is any triple with Z = 0, and the formulas need no division. P-256 has a = −3, which lets doubling factor 3x² − 3z⁴ as 3(X−Z²)(X+Z²) and save operations. This speed-up is the usual explanation for that choice of a.

p = 2**256 - 2**224 + 2**192 + 2**96 - 1          # P-256, a = -3
b = 0x5AC635D8AA3A93E7B3EBBD55769886BC651D06B0CC53B0F63BCE3C3E27D2604B
n = 0xFFFFFFFF00000000FFFFFFFFFFFFFFFFBCE6FAADA7179E84F3B9CAC2FC632551
G = (0x6B17D1F2E12C4247F8BCE6E563A440F277037D812DEB33A0F4A13945D898C296,
     0x4FE342E2FE1A7F9B8EE7EB4A7C0F9E162BCE33576B315ECECBB6406837BF51F5)
INF = (1, 1, 0)                                    # Jacobian: x = X/Z^2, y = Y/Z^3

def dbl(P):                                        # 3M + 5S
    X1, Y1, Z1 = P
    if Y1 == 0 or Z1 == 0:
        return INF
    delta, gamma = Z1 * Z1 % p, Y1 * Y1 % p
    beta = X1 * gamma % p
    alpha = 3 * (X1 - delta) * (X1 + delta) % p   # uses a = -3
    X3 = (alpha * alpha - 8 * beta) % p
    Z3 = ((Y1 + Z1) ** 2 - gamma - delta) % p
    Y3 = (alpha * (4 * beta - X3) - 8 * gamma * gamma) % p
    return (X3, Y3, Z3)

def madd(P, Q):                                    # Jacobian + affine: 7M + 4S
    X1, Y1, Z1 = P
    x2, y2 = Q
    if Z1 == 0:
        return (x2, y2, 1)
    Z1Z1 = Z1 * Z1 % p
    U2, S2 = x2 * Z1Z1 % p, y2 * Z1 * Z1Z1 % p
    H, r = (U2 - X1) % p, 2 * (S2 - Y1) % p
    if H == 0:                                     # exceptional: P = Q or P = -Q
        return dbl(P) if r == 0 else INF
    HH = H * H % p
    I = 4 * HH
    J, V = H * I % p, X1 * I % p
    X3 = (r * r - J - 2 * V) % p
    Y3 = (r * (V - X3) - 2 * Y1 * J) % p
    Z3 = ((Z1 + H) ** 2 - Z1Z1 - HH) % p
    return (X3, Y3, Z3)

Instrumented, dbl costs 3M + 5S and madd costs 7M + 4S. Mixed addition takes one input in affine form (Z = 1), which removes several multiplications. That is the reason the precomputed table in the next sections is stored in affine form.

Note the H == 0 branch. The general addition formula is incomplete: for equal or opposite inputs it returns a meaningless triple. Normal scalar multiplications almost never hit this, so tests miss it, but attackers can choose inputs that do. The robust fix is complete formulas, which are correct for every input pair. Renes, Costello and Batina (2016) published complete formulas for prime-order short Weierstrass curves such as P-256, and Edwards curves such as Ed25519's are complete by construction.

Batch inversion

Sometimes you need many inversions at once, for example to convert a table of Jacobian points to affine. Montgomery's trick inverts k elements with one inversion. It computes the running products a₁, a₁a₂, ..., inverts the final product once, and then walks backwards, peeling off one factor per step.

def batch_inverse(xs):               # 1 inversion + about 3k multiplications
    prefix, acc = [], 1
    for x in xs:
        prefix.append(acc)
        acc = acc * x % p
    inv = pow(acc, p - 2, p)         # the only inversion
    out = [0] * len(xs)
    for i in range(len(xs) - 1, -1, -1):
        out[i] = inv * prefix[i] % p  # = 1 / xs[i]
        inv = inv * xs[i] % p         # drop xs[i] from the running inverse
    return out

def to_affine_many(pts):
    zinv = batch_inverse([Z for _, _, Z in pts])
    return [(X * zi * zi % p, Y * zi * zi * zi % p) for (X, Y, _), zi in zip(pts, zinv)]

On 50 elements it used 1 inversion and 150 multiplications instead of 50 inversions. One zero input corrupts every output, so check for zeros first.

Scalar multiplication with wNAF

Plain double-and-add uses one addition for each 1 bit in the scalar, about 128 of them. Two ideas cut that. First, subtraction costs the same as addition, because −(x, y) = (x, −y). Second, precomputed odd multiples P, 3P, 5P, ... let one addition handle several bits. The width-w NAF recoding uses both. It writes the scalar with digits in {0, ±1, ±3, ..., ±(2^w−1−1)}, and every nonzero digit is followed by at least w−1 zeros. For example, 93 = 1011101₂ has five 1 bits, while its width-3 NAF is 3, 0, 0, 0, 0, −3 (3·32 − 3), with just two nonzero digits.

def wnaf(k, w):                      # least significant digit first
    digits = []
    while k:
        d = 0
        if k & 1:
            d = k % (1 << w)
            if d >= 1 << (w - 1):
                d -= 1 << w          # choose the signed residue
            k -= d                   # now k is divisible by 2^w
        digits.append(d)
        k >>= 1
    return digits

def scalar_mult(k, P, w=5):
    jac = [(P[0], P[1], 1)]
    twoP = to_affine_many([dbl(jac[0])])[0]
    for _ in range((1 << (w - 2)) - 1):           # P, 3P, 5P, ..., 15P
        jac.append(madd(jac[-1], twoP))
    table = to_affine_many(jac)                     # one shared inversion
    acc = INF
    for d in reversed(wnaf(k % n, w)):
        acc = dbl(acc)
        if d > 0:
            acc = madd(acc, table[d >> 1])
        elif d < 0:
            x, y = table[(-d) >> 1]
            acc = madd(acc, (x, -y % p))            # negation is free
    if acc[2] == 0:
        raise ValueError("result is the point at infinity")
    return to_affine_many([acc])[0]

With w = 5, over 200 random 256-bit scalars, 16.9% of wNAF digits were nonzero, close to the theoretical 1/(w+1). Binary digits were 49.7% nonzero. One full multiplication measured 1,183 M + 1,483 S + 3 I. For five random private keys, scalar_mult(d, G) returned exactly the public key that OpenSSL derived.

Fixed-base multiplication by G (keygen, signing) uses large precomputed comb tables that remove most doublings. ECDH builds a small table per call, as above. ECDSA verification interleaves two wNAF expansions over one chain of doublings.

Constant time: what production engines change

The engine above leaks. The digit pattern decides when madd runs, which timing and power traces reveal, and the secret table index shows up in cache behaviour. That is fine for public scalars, as in verification. For secret scalars, production engines change four things:

1. Regular recoding. Signed fixed-window digits that are all nonzero, so every window does the same w doublings and one addition.
2. Constant-time table lookup. Read every entry and keep the right one with a mask computed from the index, so memory access never depends on it.
3. Branch-free field code. Fixed limb counts, masked conditional subtraction, and no early exits.
4. Complete formulas or a ladder. No exceptional-case branches. X25519 uses the Montgomery ladder with a constant-time conditional swap, as the Diffie–Hellman article shows.

Point decompression

A compressed point is x plus one bit giving the parity of y. Decompression solves y² = x³ − 3x + b. For P-256, p ≡ 3 (mod 4), so a square root is a single exponentiation, a^(p+1)/4. The step that matters is checking the result. About half of all x values are not on the curve, and for those the exponentiation returns a number whose square is not the input.

def decompress(x, y_is_odd):
    rhs = (x ** 3 - 3 * x + b) % p
    y = pow(rhs, (p + 1) // 4, p)
    if y * y % p != rhs:
        raise ValueError("x is not on the curve")    # e.g. x = 1
    return (x, y if y & 1 == y_is_odd else p - y)

On the five OpenSSL keys, decompression reproduced each y. For x = 1 it raised an error. Curve25519's prime is 5 mod 8, so it needs a different square-root routine. To map arbitrary data to a point, use the constant-time hash-to-curve methods of RFC 9380, not try-and-increment.

Testing an engine

Engine bugs fail on a tiny fraction of inputs, so random end-to-end tests miss them. Test in layers:

• Differential tests. Compare points with a trusted library, and every field routine with a big-integer reference on random and boundary values (0, 1, p−1, just above p).
• Edge scalars. 0 and n give infinity, 1 gives G, n−1 gives −G, and small k can be checked against repeated affine addition. Our engine passed k = 1 to 59 and n−1.
• Exceptional additions. Call addition on P + P and P + (−P) directly, since normal scalar multiplications rarely reach them.
• Adversarial vectors. Project Wycheproof publishes test vectors aimed at these bugs.
• Constant-time and verified code. Tools such as dudect and ctgrind detect secret-dependent timing or memory access. Formally verified field code, such as fiat-crypto's output, removes the carry-bug class.

Failure modes

• Converting infinity to affine. Z = 0 has no inverse. Before the guard was added, the toy engine silently returned (0, 0) for k = 0 and k = n. Check Z before converting and treat infinity as an error wherever a protocol forbids it.
• Incomplete addition. Doubling routed through the addition formula, or P + (−P), produces garbage. Use complete formulas or handle both cases explicitly.
• Missing on-curve check. Decompression without the squaring check, or accepting uncompressed points without checking the equation, enables invalid-curve attacks.
• Variable-time code on secrets. wNAF branches, table indexing and big-integer libraries are fine for verification and fatal for signing or ECDH.
• Non-canonical outputs. A value at or above p that never got a final reduction breaks comparisons and encodings.

What to do next

1. Assemble the code blocks into one file and check scalar_mult(d, G) against ec.generate_private_key(ec.SECP256R1()): compare private_numbers().private_value with public_key().public_numbers().
2. Wrap the field operations in counters and reproduce 3M + 5S, 7M + 4S and the totals. Then try w = 4 and w = 6 and find the cheapest width.
3. Replace madd's H == 0 branch with an assertion and find a scalar that triggers it.
4. Implement the 2²⁵⁵ − 19 fold with five 51-bit limbs and a fixed number of carry steps, and fuzz it against %.
5. Run Project Wycheproof's vectors against the production library you depend on.
6. Review the exponentiation pattern in modular exponentiation. Fermat inversion and point decompression are both exponentiations.