# Copyright (c) The btclib developers
# Distributed under the MIT software license, see the accompanying
# LICENSE file or https://opensource.org/license/mit for the full text.
"""Elliptic CurveGroup class and functions.
CurveGroup is every point of the curve, a group that needs be neither
cyclic nor of prime order. Curve, in the ecc.curves.curve module, is
the cyclic subgroup of prime order generated by G, which is the group
everything built on curves multiplies in; the two are the whole of the
hierarchy, and the class docstring of Curve says what separates them.
"""
from __future__ import annotations
import functools
import heapq
import secrets
from collections.abc import Sequence
from math import ceil
from typing_extensions import override
from ecc._utils import assert_type, hex_string, int_from_integer, is_integer
from ecc.alias import INF, INFJ, Integer, JacPoint, Point
from ecc.exceptions import EccTypeError, EccValueError
from ecc.number_theory import (
_field_element,
_is_prime_var,
mod_inv_batch_var,
mod_inv_var,
mod_sqrt_var,
)
__all__ = [
"BOS_COSTER_THRESHOLD",
"HEX_THRESHOLD",
"MAX_W",
"CurveGroup",
"signed_odd_digits",
]
HEX_THRESHOLD = 0xFFFFFFFF
def _int_point(Q: Point) -> Point:
"""Return the affine point with int coordinates.
The group law computes in the type of the curve's `_modulus`, which may
be gmpy2's mpz, and a point leaves the package as int.
"""
return int(Q[0]), int(Q[1])
def _int_jac(Q: JacPoint) -> JacPoint:
"""Return the Jacobian point with int coordinates, as `_int_point`."""
return int(Q[0]), int(Q[1]), int(Q[2])
def _jac_from_aff(Q: Point) -> JacPoint:
"""Return the Jacobian representation of the affine point.
The input point is assumed to be on curve.
"""
return Q[0], Q[1], 1 if Q[1] else 0
def _blinded_jac(Q: JacPoint, ec: CurveGroup) -> JacPoint:
"""Return the same point with a Z nobody can predict.
A Jacobian point has a coordinate more than the curve needs: (X, Y, Z)
and (l^2*X, l^3*Y, l*Z) are the same affine point for every nonzero l.
Spending that freedom on a random l is what libsecp256k1's
`secp256k1_ecmult_gen_gej` does with `secp256k1_gej_rescale` on its
first table lookup, so that the values a multiplication forms are not
a function of the scalar alone: without it the same scalar produces
the same integers on every call and in every process, and a table of
them can be built offline.
It is not what makes a multiplication constant-time and does not
claim to be -- SECURITY.md publishes the Python path as
variable-time, and the recodings of curve_group_2 are what address
the duration. This addresses the values.
Three multiplications and a random field element, against the ~2500
multiplications of the multiplication it precedes. Infinity comes
back infinity, Z == 0 scaling to Z == 0.
"""
p = ec.p
# a nonzero l: zero would send the point to infinity, which is the one
# value of l that is not an isomorphism
lam = 1 + secrets.randbelow(p - 1)
lam2 = lam * lam % p
return Q[0] * lam2 % p, Q[1] * lam2 % p * lam % p, Q[2] * lam % p
[docs]
class CurveGroup:
"""Finite group of the points of an elliptic curve over Fp.
The elliptic curve is the set of points (x, y) that are solutions to
a Weierstrass equation y^2 = x^3 + a*x + b, with x, y, a, and b in
Fp (p being a prime), together with a point at infinity INF. The
constants a, b must satisfy the relationship 4 a^3 + 27 b^2 ā 0.
The group is defined by the point addition group law, INF being its
neutral element, and it is finite and abelian ā hence a product of at
most two cyclic groups, Z_n1 x Z_n2 with n2 dividing n1, so neither
cyclic nor of prime order in general.
How many points it holds is not among its data, and cannot be
computed from them: the order is p+1 minus the trace of Frobenius,
which Hasse's theorem bounds by twice the square root of p and
nothing here counts. Nor is any point of the group distinguished from
the others. The generator G, the prime order n of the cyclic subgroup
it generates, and the cofactor relating that n to the order of this
group are the parameters of Curve, in ecc.curves.curve, which is
what a caller of the library multiplies in; what this class is for is
the arithmetic underneath, shared by the two and defined by p, a and b
alone.
"""
# Set on an instance by `curve._catalogued_curve` alone, to skip the
# test that p (and n, for a Curve) is prime
_primes_are_trusted = False
def __init__(self, p: Integer, a: Integer, b: Integer) -> None:
# Parameters are checked according to SEC 1 v.2 3.1.1.2.1
p = int_from_integer(p)
a = int_from_integer(a)
b = int_from_integer(b)
# 1) check that p is a prime
if not self._primes_are_trusted and not _is_prime_var(p):
err_msg = "p is not prime: "
err_msg += f"'{hex_string(p)}'" if p > HEX_THRESHOLD else f"{p}"
raise EccValueError(err_msg)
plen = p.bit_length()
# byte-length
self.p_size = ceil(plen / 8)
# must be true to break symmetry using quadratic residue
self.p_is_3_mod_4 = p % 4 == 3
self.p = p
# what the Jacobian group law reduces by: p as gmpy2's mpz where the
# `gmpy2` extra is installed, so that its products are GMP's
self._modulus = _field_element(p)
# how many bits a scalar of this group has, which is what fixes
# the digit count of _mult_regular_window below -- the quantity a
# loop running once per digit of the scalar in hand puts on the
# clock. Hasse bounds the order of the group by p + 1 + 2*sqrt(p),
# so a scalar reduced mod that order is under 2^(plen+1); Curve
# narrows it to nlen, the order of its subgroup being the bound
# there. It is a count and not a limit: a scalar above it is
# multiplied all the same, in the digits it needs
self.scalar_len = plen + 1
# 2. check that a and b are integers in [0, pā1] # noqa: RUF003
if a < 0:
raise EccValueError(f"negative a: {a}")
if p <= a:
err_msg = "p <= a: " + (
f"'{hex_string(p)}' <= '{hex_string(a)}'"
if p > HEX_THRESHOLD
else f"{p} <= {a}"
)
raise EccValueError(err_msg)
if b < 0:
raise EccValueError(f"negative b: {b}")
if p <= b:
err_msg = "p <= b: " + (
f"'{hex_string(p)}' <= '{hex_string(b)}'"
if p > HEX_THRESHOLD
else f"{p} <= {b}"
)
raise EccValueError(err_msg)
# 3. Check that 4*a^3 + 27*b^2 ā 0 (mod p)
d = 4 * a * a * a + 27 * b * b
if d % p == 0:
raise EccValueError("zero discriminant")
self._a = a
self._b = b
# which spelling of the tangent law _double_jac_helper runs. The
# curve decides it and a curve does not change, so it is decided
# here rather than by a test inside the function every Jacobian
# doubling reaches -- double_jac's own, and add_jac's and
# add_jac_aff's when the two points they are given coincide -- the
# a*Z^4 term being zero for the k1 family, and the tangent
# numerator a difference of two squares for the a == p - 3 of most
# catalogued curves, costing one product instead of three
self._a_is_zero = a == 0
self._a_is_minus_3 = a == p - 3
# what add_jac feeds its formula in place of an operand at
# infinity: any value does, the arithmetic on it being thrown
# away, as long as it is the size of a real coordinate -- which is
# the whole of what it is for, so a fraction of p rather than a
# small constant. On a curve whose p is a handful of bits some of
# these are zero, and nothing there is worth timing anyway
self._stand_in_q = p // 3, p // 5, p // 7
self._stand_in_r = p // 11, p // 13, p // 17
# the points whose wNAF table is worth memoizing rather than
# rebuilding per call: the ones the curve itself names, which a
# verification hands the same on every call. A CurveGroup names
# none -- it has no generator -- and Curve fills this in
self._fixed_points: frozenset[JacPoint] = frozenset()
@override
def __str__(self) -> str:
# the class name, not the literal "Curve": the two classes hold
# different parameters, so a CurveGroup announcing itself as a
# Curve reads as a curve whose G and n went missing -- and the
# repr below, which is meant to eval back to an equal object,
# would name a constructor four arguments short
result = type(self).__name__
if self.p > HEX_THRESHOLD:
result += f"\n p = {hex_string(self.p)}"
else:
result += f"\n p = {self.p}"
if self._a > HEX_THRESHOLD or self._b > HEX_THRESHOLD:
result += f"\n a = {hex_string(self._a)}"
result += f"\n b = {hex_string(self._b)}"
else:
result += f"\n a = {self._a}"
result += f"\n b = {self._b}"
return result
@override
def __repr__(self) -> str:
result = f"{type(self).__name__}("
result += f"'{hex_string(self.p)}'" if self.p > HEX_THRESHOLD else f"{self.p}"
if self._a > HEX_THRESHOLD or self._b > HEX_THRESHOLD:
result += f", '{hex_string(self._a)}', '{hex_string(self._b)}'"
else:
result += f", {self._a}, {self._b}"
result += ")"
return result
def _eq_key(self) -> tuple[int, ...]:
"""Return the parameters that define the group, and nothing else.
Curve appends its own, the way __str__ and __repr__ chain: the
tuple is flat so that the subgroup cannot be equal to the group by
having a shorter key compare equal.
"""
return self.p, self._a, self._b
@override
def __eq__(self, other: object) -> bool:
# a curve is its parameters: two objects built from the same ones
# are the same curve, whatever the catalogue they came from. This
# is not cosmetic, as the libsecp256k1 dispatch tests ec against
# secp256k1, and the default identity comparison would silently
# send every other object holding the secp256k1 parameters down
# the slow Python path
if self is other:
return True
# an exact type test, not isinstance alone: a CurveGroup is not a
# Curve even when p, a and b are the same, one being every point
# of the curve and the other a subgroup of those points, and
# returning NotImplemented lets the reflected operand have its
# say before Python falls back to identity
if not isinstance(other, CurveGroup) or type(self) is not type(other):
return NotImplemented
return self._eq_key() == other._eq_key()
@override
def __hash__(self) -> int:
# __eq__ without __hash__ would set __hash__ to None, and curves
# are lru_cache keys in _cached_multiples below; equal curves
# share those cache entries rather than each filling its own
return hash(self._eq_key())
# methods using p: they could become functions
[docs]
def negate(self, Q: Point) -> Point:
"""Return the opposite point.
The input point is not checked to be on the curve.
"""
# % self.p is required to account for INF (i.e. Q[1]==0)
# so that negate(INF) = INF
if len(Q) == 2:
return Q[0], (self.p - Q[1]) % self.p
raise EccValueError("not a point")
[docs]
def negate_jac(self, Q: JacPoint) -> JacPoint:
"""Return the opposite Jacobian point.
The input point is not checked to be on the curve.
"""
# % self.p is required to account for INF (i.e. Q[1]==0)
# so that negate(INF) = INF
if len(Q) == 3:
return Q[0], (self.p - Q[1]) % self.p, Q[2]
raise EccValueError("not a Jacobian point")
[docs]
def aff_from_jac_var(self, Q: JacPoint) -> Point:
"""Return the affine point: one modular inversion.
The input point is assumed to be on the curve; infinity comes
back as INF, its affine spelling.
"""
if Q[2] == 0: # Infinity point in Jacobian coordinates
return INF
return _int_point(self._aff_from_z_inv(Q, mod_inv_var(Q[2], self._modulus)))
[docs]
def aff_from_jac_batch_var(self, Qs: Sequence[JacPoint]) -> list[Point]:
"""Return the affine points: one modular inversion for all of them.
`aff_from_jac_var` over a sequence, with `mod_inv_batch_var` in place of
the one inverse each: the conversion is two products a point once
the inverse is in hand, so a caller holding several Jacobian
points pays one extended Euclid instead of one per point.
The input points are assumed to be on the curve. Infinity is not
in the batch, having no Z to invert, and comes back as INF where
it stood.
"""
return [_int_point(Q) for Q in self._aff_from_jac_batch_var(Qs)]
def _aff_from_jac_batch_var(self, Qs: Sequence[JacPoint]) -> list[Point]:
# in the type of self._modulus: the tables of the multiplications
# are built here, and stay in the type their arithmetic runs on
inverses = iter(mod_inv_batch_var([Q[2] for Q in Qs if Q[2]], self._modulus))
return [
INF if Q[2] == 0 else self._aff_from_z_inv(Q, next(inverses)) for Q in Qs
]
def _aff_from_z_inv(self, Q: JacPoint, z_inv: int) -> Point:
# x is X/Z^2 and y is Y/Z^3, and both powers are built from one
# inverse: Z^-2 is a product away from Z^-1, and Z^-3 another.
# Inverting Z^2 and Z^3 separately is two extended Euclids where
# this is one and two multiplications, and an inverse costs what
# a hundred products cost -- which is also why the inverse itself
# is the caller's, one of them being what the batch above shares
p = self._modulus
z_inv2 = z_inv**2 % p
return Q[0] * z_inv2 % p, Q[1] * z_inv2 % p * z_inv % p
[docs]
def x_aff_from_jac_var(self, Q: JacPoint) -> int:
"""Return the affine x alone, without the products y costs.
One inversion, as `aff_from_jac_var`, of Z^2 rather than of Z: the
power x wants is the one inverted, so nothing is rebuilt from it.
The input point is assumed to be on the curve; infinity has no
x and is refused.
"""
if Q[2] == 0: # Infinity point in Jacobian coordinates
raise EccValueError("INF has no x-coordinate")
p = self._modulus
Z2 = Q[2] * Q[2]
return int(Q[0] * mod_inv_var(Z2, p) % p)
[docs]
def y_aff_from_jac_var(self, Q: JacPoint) -> int:
"""Return the affine y alone, without the product x costs.
One inversion, as `aff_from_jac_var`, of Z^3 rather than of Z, for
the same reason `x_aff_from_jac_var` inverts Z^2.
The input point is assumed to be on the curve; infinity has no
y and is refused.
"""
if Q[2] == 0: # Infinity point in Jacobian coordinates
raise EccValueError("INF has no y-coordinate")
p = self._modulus
Z2 = Q[2] * Q[2]
return int(Q[1] * mod_inv_var(Z2 * Q[2], p) % p)
[docs]
def is_jac_equal(self, QJ: JacPoint, PJ: JacPoint) -> bool:
"""Return True if Jacobian points are equal in affine coordinates.
The input points are assumed to be on curve.
"""
PJ2 = PJ[2] * PJ[2]
QJ2 = QJ[2] * QJ[2]
if QJ[0] * PJ2 % self.p != PJ[0] * QJ2 % self.p:
return False
PJ3 = PJ2 * PJ[2]
QJ3 = QJ2 * QJ[2]
return QJ[1] * PJ3 % self.p == PJ[1] * QJ3 % self.p
# methods using _a, _b, p
[docs]
def add_var(self, Q1: Point, Q2: Point) -> Point:
"""Return the sum of two points.
The input points must be on the curve.
"""
self.require_on_curve(Q1)
self.require_on_curve(Q2)
# affine arithmetic, and it is not the inversion that decides it:
# aff_from_jac_var makes one, add_aff_var makes one. The products
# decide, and reaching this sum through Jacobian coordinates and
# back makes more of them than add_aff_var does
return self.add_aff_var(Q1, Q2)
[docs]
def add_jac(self, Q: JacPoint, R: JacPoint) -> JacPoint:
"""Return the sum of two Jacobian points, branch-free.
The input points are assumed to be on the curve. One sequence
of operations whatever the operands -- infinity and doubling
included, the comment of `_add_jac` saying why that is
load-bearing.
"""
return _int_jac(self._add_jac(Q, R))
def _add_jac(self, Q: JacPoint, R: JacPoint) -> JacPoint:
# `add_jac` in the type of self._modulus, which is what the
# multiplications call: a point is converted to int once, where it
# leaves the package.
#
# the group law has four special cases, and they are of two kinds
# that want opposite treatment. An operand at infinity is
# bookkeeping rather than geometry: infinity is the identity, so
# every double-and-add starts its accumulator there and the
# windowed ones hold it in the table as the multiple 0*Q, which
# the digit 0 indexes. On secp256k1 that is not rare and not
# random -- over 2000 random scalars _mult reached it 3.98 times
# on average and only 23 of them not at all, and the count *is*
# the secret: exactly the number of zero base-16 digits of the
# scalar, or for _mult_jac_var the number of its low zero bits. An
# early return there puts that on the clock: skipping the
# addition costs a fraction of making it, so a scalar with no
# zero digit takes measurably longer through _mult than one with
# 63 of them.
#
# So no branch on infinity, and no arithmetic on it either: a
# Python integer costs what its size costs, and the zero
# coordinates of infinity send half the products to zero. A
# stand-in of full size takes its place, the table at the end
# answers for it, and an addition with no infinity in it pays two
# comparisons: 1.02x to 1.03x on every multiplication in the
# package, and 1.22x on _double_mult_var, whose Shamir-Strauss loop
# adds the infinity of a zero digit pair a quarter of the time --
# the additions an early return answers for free, at the price of
# putting their count on the clock. Dropping the early return
# without the stand-ins is not enough: P + INFJ still costs a
# third of a generic addition, and _mult_jac_var is still 25%
# faster on a scalar with 128 low zero bits
QS = (Q, self._stand_in_q)[Q[2] == 0]
RS = (R, self._stand_in_r)[R[2] == 0]
p = self._modulus
# every intermediate reduced as it is formed, which is what makes
# this the fast path rather than a transcription of the formula:
# left to grow, V3 and M*V2 reach p^9 and the products that close
# X and Y reach p^12 -- three thousand bits on secp256k1 -- and
# Python multiplies whatever it is handed. Worth between 2.0 and
# 3.0 times over every mult_* variant in the package, measured
# against the unreduced spelling on secp256k1 and secp256r1
RZ2 = RS[2] * RS[2] % p
RZ3 = RZ2 * RS[2] % p
QZ2 = QS[2] * QS[2] % p
QZ3 = QZ2 * QS[2] % p
# the two points in a common frame: V is the difference of their affine
# x-coordinates and W of their affine y ones, up to that frame. Which is
# what makes the two remaining cases free to detect, both being needed
# by the formula below -- the four multiplications issue
# btclib-org/btclib#171 counted were never their cost either, M, N, T
# and U all feeding them too
M = QS[0] * RZ2 % p
T = QS[1] * RZ3 % p
V = (RS[0] * QZ2 - M) % p
W = (RS[1] * QZ3 - T) % p
# the other kind, and a branch: V == 0 is one x-coordinate for
# both points, so they are the same point -- W == 0 too, the chord
# is the tangent, and the answer is the doubling -- or opposite
# ones, and it is infinity. That is geometry, and inside a
# multiplication it needs the accumulator to land on a table entry
# or its negation: 2^-250 on a curve with a real order, against
# the 3.98 times a random scalar reaches infinity above. What the
# branch tells an attacker is that two points coincided, which on
# secp256k1 they will not, and a caller spelling out P + P knows
# it did. It is the toy curves of the test suite that take it, and
# take it constantly, an accumulator wrapping around a group of
# eleven points; the exhaustive test over every pair of points of
# every one of them is what holds this branch correct. What would
# retire the case rather than branch on it is a complete formula:
# Renes-Costello-Batina, homogeneous coordinates, no exceptional
# case at all -- and measured on this very fixed window 34% slower
# for a == 0 and 59% slower for the a == p-3 of most catalogued
# curves, on top of JacPoint being public.
#
# Q[2] and R[2] because V says nothing when either is zero: Z == 0
# leaves X and Y unconstrained, and the x of INFJ is 7 while the x of
# _jac_from_aff(INF) is 5, picked for being invalid rather than for
# being zero, so on a curve of characteristic 7 or 5 they reduce to zero
# and V vanishes for an operand that is not a point at all. Without
# the two tests P + INFJ answers 2*P, which is what issue
# btclib-org/btclib#171 found making most of btclib-org/btclib's scalar
# multiplications wrong for most scalars on every p == 7 curve. The
# stand-ins do not make it safe to read: they are values, and V can be
# zero on a pair of them too, one time in p
if V == 0 and Q[2] and R[2]:
return self._double_jac_helper(Q, QZ2) if W == 0 else INFJ
V2 = V * V % p
V3 = V2 * V % p
MV2 = M * V2 % p
X = (W * W - V3 - 2 * MV2) % p
Y = (W * (MV2 - X) - T * V3) % p
Z = V * QS[2] % p * RS[2] % p
# and the infinity cases, answered rather than branched on. The
# arithmetic above was done on the stand-ins and is thrown away
# here, which is the whole point of it
return ((X, Y, Z), R, Q, INFJ)[(Q[2] == 0) + 2 * (R[2] == 0)]
[docs]
def add_jac_aff(self, Q: JacPoint, R: Point) -> JacPoint:
"""Return the sum of a Jacobian point and an affine one, branch-free.
`add_jac` with the second operand's Z known to be one, which is
what an affine point is: five of the sixteen products become
multiplications by one -- R's two powers of Z, the two that put Q
in R's frame, and one factor of the answer's Z. libsecp256k1 keeps
its tables in affine coordinates for exactly this and states
`secp256k1_gej_add_ge_var` as 8 mul and 3 sqr against the 12 and 4
of `secp256k1_gej_add_var`.
The affine operand is a Point and not a JacPoint whose Z happens
to be one, so the precondition is in the signature rather than in
a sentence nothing checks.
The input points are assumed to be on the curve. Every case
`add_jac` answers is answered here the same way and for the same
reasons, which its comments carry: infinity through a stand-in and
a selection at the end, spelled `R[1] == 0` because that is what
infinity is in affine coordinates, and the doubling and the sum
that is infinity through the one branch on V.
"""
return _int_jac(self._add_jac_aff(Q, R))
def _add_jac_aff(self, Q: JacPoint, R: Point) -> JacPoint:
# `add_jac_aff` in the type of self._modulus, as `_add_jac` is
QS = (Q, self._stand_in_q)[Q[2] == 0]
RS = (R, self._stand_in_r)[R[1] == 0]
p = self._modulus
QZ2 = QS[2] * QS[2] % p
QZ3 = QZ2 * QS[2] % p
# M and T are Q's coordinates as they stand: R's frame is the
# affine one, so putting Q in it is what QZ2 and QZ3 already did
M = QS[0]
T = QS[1]
V = (RS[0] * QZ2 - M) % p
W = (RS[1] * QZ3 - T) % p
if V == 0 and Q[2] and R[1]:
return self._double_jac_helper(Q, QZ2) if W == 0 else INFJ
V2 = V * V % p
V3 = V2 * V % p
MV2 = M * V2 % p
X = (W * W - V3 - 2 * MV2) % p
Y = (W * (MV2 - X) - T * V3) % p
Z = V * QS[2] % p
return ((X, Y, Z), (R[0], R[1], 1), Q, INFJ)[(Q[2] == 0) + 2 * (R[1] == 0)]
[docs]
def double_jac(self, Q: JacPoint) -> JacPoint:
"""Return twice the Jacobian point, assumed to be on the curve."""
return _int_jac(self._double_jac(Q))
def _double_jac(self, Q: JacPoint) -> JacPoint:
# `double_jac` in the type of self._modulus, as `_add_jac` is.
# Z^2 is what the a*Z^4 term is built from and the only thing that
# needs it, so a curve whose a is zero is handed a value the
# helper never reads rather than a squaring and its reduction
QZ2 = 0 if self._a_is_zero else Q[2] * Q[2] % self._modulus
return self._double_jac_helper(Q, QZ2)
def _double_jac_helper(self, Q: JacPoint, QZ2: int) -> JacPoint:
# QZ2 is Q[2]^2 reduced mod p, which add_jac has already formed;
# a curve whose a is zero never reads it, and double_jac above
# does not form it. The tangent law has no case to tell apart -- a
# point is never exceptional to itself -- so this is the whole of
# it: no test on Q, one sequence of operations for every point of
# the curve, and the same one for a secret point as for a public
# one. Which is also what lets add_jac call it on every addition
p = self._modulus
QY2 = Q[1] * Q[1] % p
# the a*Z^4 term, in the spelling this curve's a allows: the tests
# are on the curve and not on the point, so the operations a given
# curve makes are still the same for every point of it.
# libsecp256k1 carries the first spelling alone, secp256k1's a
# being zero, and states the doubling it makes as 3 mul and 4 sqr
# against the 12 and 4 of its addition
if self._a_is_zero:
W = 3 * Q[0] * Q[0] % p
elif self._a_is_minus_3:
# 3*(X - Z^2)*(X + Z^2) is 3*X^2 - 3*Z^4, one product where
# the general spelling below has two squarings and a product
W = 3 * (Q[0] - QZ2) * (Q[0] + QZ2) % p
else:
W = (3 * Q[0] * Q[0] + self._a * QZ2 * QZ2) % p
V = 4 * Q[0] * QY2 % p
X = (W * W - 2 * V) % p
Y = (W * (V - X) - 8 * QY2 * QY2) % p
# Q[1] == 0 is both the infinity point and the two-torsion, and
# both double to infinity, which is what Z == 0 says
Z = 2 * Q[1] * Q[2] % p
return X, Y, Z
[docs]
def add_aff_var(self, Q: Point, R: Point) -> Point:
"""Return the sum of two affine points, assumed on the curve.
One modular inversion; the special cases are branched on, the
comment below saying why affine coordinates leave no choice.
"""
# infinity is a special case here and cannot stop being one: an
# affine point is two field elements and the group has one point
# more than any pair of them can name, so INF is spelled y == 0
# (alias.py) and no formula reaches it. That spelling is also why
# a two-torsion point, the one real point with y == 0, has no
# affine form at all -- Jacobian coordinates do hold it, Z != 0
# telling it from infinity, and add_jac doubles it to infinity
# correctly; none of the test curves of cofactor 1 has one, and
# two catalogued curves of cofactor 4 and two test curves of
# cofactor 2 and 4 do (issue btclib-org/ecc#16). This method
# still assumes its operands on the curve and reads either one's
# y == 0 as infinity whatever it stands for, as add_jac's own
# stand-ins do; on a Curve, add_var -- the validated entry point --
# never reaches it with the real two-torsion point,
# Curve.is_on_curve refusing that one tuple ahead of this call for
# exactly this reason. A
# bare CurveGroup makes no such refusal, off āØGā© entirely.
# The order is load-bearing too, and it is this one rather than
# the doubling test first: INF is (5, 0), its x-coordinate
# arbitrary, so a doubling test reading it ahead of here answers
# "same x, different y, hence INF" for INF + P whenever P has
# x == 5 -- and ec23_19 of the test suite is generated by (5, 4),
# which is exactly that point (issue btclib-org/btclib#171)
if R[1] == 0: # Infinity point in affine coordinates
return Q
if Q[1] == 0: # Infinity point in affine coordinates
return R
# then doubling: two equal points make the chord a tangent and the
# slope below 0/0. Branches, where add_jac computes every case and
# selects, because there is nothing here for that to protect: the
# mod_inv_var this path exists to spend is an extended Euclid, whose
# iteration count follows the value it is inverting, so the affine
# law cannot be made uniform whatever is done to these three tests
if R[0] == Q[0]:
return self.double_aff_var(R) if R[1] == Q[1] else INF
p = self.p
# lam reduced before it is squared, as in add_jac, though here it
# is worth 6% of _mult_aff_var rather than a factor of two:
# mod_inv_var is what affine coordinates cost, and is why the
# ladders are not
lam = (R[1] - Q[1]) * mod_inv_var(R[0] - Q[0], p) % p
x = (lam * lam - Q[0] - R[0]) % p
y = (lam * (Q[0] - x) - Q[1]) % p
return x, y
[docs]
def double_aff_var(self, Q: Point) -> Point:
"""Return twice the affine point, assumed to be on the curve."""
if Q[1] == 0: # Infinity point in affine coordinates
return INF
p = self.p
lam = (3 * Q[0] * Q[0] + self._a) * mod_inv_var(2 * Q[1], p) % p
x = (lam * lam - Q[0] - Q[0]) % p
y = (lam * (Q[0] - x) - Q[1]) % p
return x, y
def _y2(self, x: int) -> int:
# skipping a crucial check here:
# if sqrt(y*y) does not exist, then x is not valid.
# This is a good reason to keep this method private
return ((x**2 + self._a) * x + self._b) % self.p
[docs]
def y_var(self, x: int) -> int:
"""Return the y coordinate from x, as in (x, y)."""
if not 0 <= x < self.p:
raise EccValueError("x-coordinate not in 0..p-1")
try:
return mod_sqrt_var(self._y2(x), self.p)
except EccValueError:
pass
# raised outside the handler, so that neither __cause__ nor __context__
# carries mod_sqrt_var's message, which names y**2 and so x
raise EccValueError("invalid x-coordinate")
[docs]
def require_on_curve(self, Q: Point) -> None:
"""Require the input curve Point to be on the curve.
An Error is raised if not.
"""
if not self.is_on_curve(Q):
raise EccValueError("point not on curve")
[docs]
def is_on_curve(self, Q: Point) -> bool:
"""Return True if the point is on the curve."""
# the type before the length: `len` of what is not sized is a
# TypeError about a builtin, where a Point is what this asks for
assert_type(Q, tuple, "point")
if len(Q) != 2:
raise EccValueError("point must be a tuple[int, int]")
# the type before the y==0 infinity check below: `False == 0`
# would read any x paired with a bool y as infinity, and
# `is_integer` is where this library says a bool is not a number
# (issue btclib-org/btclib#1249)
if not is_integer(Q[0]):
raise EccTypeError(f"non-integer x-coordinate: {type(Q[0]).__name__}")
if not is_integer(Q[1]):
raise EccTypeError(f"non-integer y-coordinate: {type(Q[1]).__name__}")
# before the y==0 check too: _y2 reduces x mod p, so an x outside
# 0..p-1 would otherwise be accepted -- x + p and x - p read back
# as a second, third encoding of the same point (issue
# btclib-org/ecc#7)
if not 0 <= Q[0] < self.p:
raise EccValueError("x-coordinate not in 0..p-1")
if Q[1] == 0: # Infinity point in affine coordinates
# true of every point this class enumerates on its own --
# find_all_points walks a bare CurveGroup with no G and no n,
# and a real order-2 point there is exactly as legitimate a
# member of the group as any other and is read as one. Curve
# overrides this method: āØGā© is what a cofactor makes
# ambiguous, a real two-torsion point outside it having no
# affine spelling distinct from INF's (issue
# btclib-org/ecc#16)
return True
if not 0 < Q[1] < self.p: # y cannot be zero
raise EccValueError("y-coordinate not in 1..p-1")
return self._y2(Q[0]) == (Q[1] * Q[1] % self.p)
# y-symmetry tiebreaker criteria: even/odd, low/high, or quadratic residue
[docs]
def y_even_var(self, x: int) -> int:
"""Return the odd/even affine y-coordinate associated to x."""
root = self.y_var(x)
# switch even/odd root as needed (XORing the conditions)
return self.p - root if root % 2 else root
[docs]
def y_low_var(self, x: int) -> int:
"""Return the low/high affine y-coordinate associated to x."""
root = self.y_var(x)
return root if root <= self.p // 2 else self.p - root
[docs]
def y_quadratic_residue_var(self, x: int) -> int:
"""Return the quadratic residue affine y-coordinate."""
if not self.p_is_3_mod_4:
err_msg = "field prime is not equal to 3 mod 4: "
err_msg += (
f"'{hex_string(self.p)}'" if self.p > HEX_THRESHOLD else f"{self.p}"
)
raise EccValueError(err_msg)
# for a p of this form self.y answers the residue root already:
# it is a ** ((p + 1) // 4), a power of the square a and so a
# square itself, which leaves p - root the non-residue, -1 being
# one modulo such a p. Asking legendre_symbol_var which of the two
# this is would spend a second exponentiation the size of the
# first to be told what the exponent settles
return self.y_var(x)
def _assert_valid_ec(ec: CurveGroup) -> None:
"""Refuse an ec that is not a group of curve points.
`hashes._assert_valid_hf` is the shape: one check where the parameter
is first read, rather than one per function that takes it. What it
catches is what an `ec` of no curve type otherwise becomes -- a field
read off it, `'NoneType' object has no attribute 'p'`, which is a
caller's own mistake reported from underneath the library.
Every public function whose `ec` is a `CurveGroup` asks here, and the
ones whose `ec` is a `Curve` ask `curve._assert_valid_ec` instead:
`n` and `G` are Curve's parameters and this class has neither, so
the group check would pass an ec that the multiplication then reads a
missing field off.
"""
assert_type(ec, CurveGroup, "ec")
def _mult_recursive_aff_var(m: int, Q: Point, ec: CurveGroup) -> Point:
"""Scalar multiplication of a curve point in affine coordinates.
This implementation uses a recursive version of 'double & add',
affine coordinates.
The input point is assumed to be on curve and the m coefficient is
assumed to have been reduced mod n if appropriate (e.g. cyclic
groups of order n).
"""
if m < 0:
raise EccValueError("negative m")
if m == 0:
return INF
if m % 2 == 1:
return ec.add_aff_var(Q, _mult_recursive_aff_var((m - 1), Q, ec))
return _mult_recursive_aff_var((m // 2), ec.double_aff_var(Q), ec)
def _mult_recursive_jac_var(m: int, Q: JacPoint, ec: CurveGroup) -> JacPoint:
"""Scalar multiplication of a curve point in affine coordinates.
This implementation uses a recursive version of 'double & add',
jacobian coordinates.
The input point is assumed to be on curve and the m coefficient is
assumed to have been reduced mod n if appropriate (e.g. cyclic
groups of order n).
"""
if m < 0:
raise EccValueError("negative m")
if m == 0:
return INFJ
if m % 2 == 1:
return ec._add_jac(Q, _mult_recursive_jac_var((m - 1), Q, ec))
return _mult_recursive_jac_var((m // 2), ec._double_jac(Q), ec)
def _mult_aff_var(m: int, Q: Point, ec: CurveGroup) -> Point:
"""Scalar multiplication of a curve point in affine coordinates.
This implementation uses 'double & add' algorithm, 'right-to-left'
binary decomposition of the m coefficient, affine coordinates.
The input point is assumed to be on curve and the m coefficient is
assumed to have been reduced mod n if appropriate (e.g. cyclic
groups of order n).
"""
if m < 0:
raise EccValueError("negative m")
# R[0] is the running result, R[1] = R[0] + Q is an ancillary variable
R = [INF, Q]
# if least significant bit of m is 1, then add Q to R[0]
R[0] = R[m & 1]
# remove the bit just accounted for
m >>= 1
while m > 0:
# the doubling part of 'double & add'
Q = ec.double_aff_var(Q)
# always perform the 'add', even if useless: one addition per bit
# whatever the bit is, as in _mult_jac_var. It buys less here,
# add_aff_var branching on its own cases and mod_inv_var taking the time
# its input asks for, which is why the affine variants are the
# readable ones rather than the ones to sign with
R[1] = ec.add_aff_var(R[0], Q)
# if least significant bit of m is 1, then add Q to R[0]
R[0] = R[m & 1]
m >>= 1
return R[0]
def _mult_jac_var(m: int, Q: JacPoint, ec: CurveGroup) -> JacPoint:
"""Scalar multiplication of a curve point in Jacobian coordinates.
This implementation uses 'double & add' algorithm, 'right-to-left'
binary decomposition of the m coefficient, Jacobian coordinates.
The input point is assumed to be on curve and the m coefficient is
assumed to have been reduced mod n if appropriate (e.g. cyclic
groups of order n).
"""
if m < 0:
raise EccValueError("negative m")
# R[0] is the running result, R[1] = R[0] + Q is an ancillary variable
R = [INFJ, Q]
# if least significant bit of m is 1, then add Q to R[0]
R[not m & 1] = Q
# remove the bit just accounted for
m >>= 1
while m > 0:
# the doubling part of 'double & add'
Q = ec._double_jac(Q)
# always perform the addition, even if useless, but use it as R[0]
# only if the least significant bit of m is 1: one doubling and
# one addition per bit, and with a branch-free add_jac under it
# the loop costs the same for every scalar of a given length
R[not m & 1] = ec._add_jac(R[0], Q)
m >>= 1
return R[0]
def _multiples(Q: JacPoint, size: int, ec: CurveGroup) -> list[JacPoint]:
"""Return {k_i * Q} for k_i in {0, ..., size-1}."""
if size < 2:
raise EccValueError(f"size too low: {size}")
k, odd = divmod(size, 2)
T = [INFJ, Q]
for i in range(3, k * 2, 2):
# not T.extend((double_jac(...), add_jac(T[-1], Q))): a tuple
# literal evaluates both elements before either is appended, so
# T[-1] there would still be the previous iteration's, not the
# double_jac just computed -- preview's FURB113 suggests it and
# unsafe-fixes applies it without noticing
T.append(ec._double_jac(T[(i - 1) // 2]))
T.append(ec._add_jac(T[-1], Q))
if odd:
T.append(ec._double_jac(T[(size - 1) // 2]))
return T
MAX_W = 5
@functools.lru_cache # results are cached to increase efficiency
def _cached_multiples(Q: JacPoint, ec: CurveGroup) -> list[JacPoint]:
"""Return the multiples 0..(2^MAX_W - 1) of Q, memoized on (Q, ec).
The table the fixed-window ladder indexes when asked to cache; the
memoization pays because the generator is the Q of most calls.
"""
T = [INFJ, Q]
for i in range(3, 2**MAX_W, 2):
# not extend(): _multiples() above has why
T.append(ec._double_jac(T[(i - 1) // 2]))
T.append(ec._add_jac(T[-1], Q))
return T
@functools.lru_cache # results are cached to increase efficiency
def _cached_multiples_fixwind(
Q: JacPoint, ec: CurveGroup, w: int
) -> list[list[JacPoint]]:
"""Made to precompute values for _mult_fixed_window_cached_var.
Do not use it for other functions. Made to be used for w=4, do not
use w.
"""
T = []
K = Q
for _ in range((ec.p_size * 8) // w + 1):
sublist = [INFJ, K]
for j in range(3, 2**w, 2):
# not extend(): _multiples() above has why
sublist.append(ec._double_jac(sublist[(j - 1) // 2]))
sublist.append(ec._add_jac(sublist[-1], K))
K = ec._double_jac(sublist[2 ** (w - 1)])
T.append(sublist)
return T
def _convert_number_to_base_var(i: int, base: int) -> list[int]:
"""Return the digits of an integer in the requested base."""
digits: list[int] = []
while i or not digits:
i, idx = divmod(i, base)
digits.append(idx)
return digits[::-1]
def _wNAF_of_m_var(m: int, w: int) -> list[int]:
"""WNAF (width-w Non-adjacent form) of number m.
Given an integer m, wNAF is a method of representation
with powers of 2, where the coefficients are odd or 0,
and where at most one of any w consecutive digits is nonzero.
It has the following properties:
- m has a unique width-w NAF.
-The length of wNAF(m) is at most one more than the length of the binary
representation of k.
-The average density of nonzero digits is approximately 1/(w + 1).
This recoding sits here next to _convert_number_to_base_var, the same kind
of integer-only helper, rather than with the wNAF multiplications of
curve_group_2: the interleaved _multi_mult_w_NAF_var below needs it, and
curve_group cannot import the module that imports it.
For complete reference see:
D. Hankerson, 'Guide to Elliptic Curve Cryptography' chapter 3
"""
w2 = 1 << w
M: list[int] = []
while m > 0:
# a run of zero bits goes in one shift, not one bit at a time
if zeros := (m & -m).bit_length() - 1:
M.extend([0] * zeros)
m >>= zeros
if w == 1:
# Computing binary NAF of m
d = 2 - (m & 3)
else:
# Computing wNAF of m
d = m & (w2 - 1)
if d >= w2 >> 1:
d -= w2
M.append(d)
m = (m - d) >> 1
return M
[docs]
def signed_odd_digits(m: int, w: int, size: int) -> list[int]:
"""Return the size signed odd base-2^w digits of an odd m.
Regular recoding, Joye-Tunstall: m = sum(digits[i] * 2^(w*i)) with
every digit odd and in {Ā±1, Ā±3, ..., Ā±(2^w - 1)}, least significant
first, as _wNAF_of_m_var returns its own. Two properties are what it is
for, and neither is the wNAF's:
- no digit is zero, so the multiplication that indexes them makes one
addition per digit whatever m is, where a wNAF adds on a nonzero
digit and so once per unit of the recoded weight
- the count is `size` and not the length of m, so a scalar's *size*
stops being visible too
The digits are signed, and that is what buys the first property: 2^w
odd values in -(2^w - 1)..2^w - 1 name 2^(w-1) points and their
opposites, and negating a Jacobian point is one modular subtraction.
m must be odd -- a sum of odd digits weighted by powers of 2 has the
parity of digits[0], so no even number has such a form at all -- and
must fit the size asked for, which is `m < 2^(w*size)`: the last digit
is what is left of m and is not reduced further.
"""
if not is_integer(m):
raise EccTypeError(f"non-integer m: {type(m).__name__}")
if not is_integer(w):
raise EccTypeError(f"non-integer w: {w}")
if not is_integer(size):
raise EccTypeError(f"non-integer size: {size}")
if m < 0:
raise EccValueError("negative m")
# a number cannot be written in basis 1 (ie w=0)
if w <= 0:
raise EccValueError(f"non positive w: {w}")
if m % 2 == 0:
raise EccValueError("even m")
if size < 1:
raise EccValueError(f"size too low: {size}")
if m >> (w * size):
raise EccValueError(f"m does not fit {size} digits")
w2 = 1 << w
digits: list[int] = []
for _ in range(size - 1):
# the low w+1 bits of m read as a signed value: odd, since m is,
# and in -(2^w - 1)..2^w - 1
digit = (m % (2 * w2)) - w2
digits.append(digit)
# m - digit is 2^w times an odd number -- m - (m mod 2^(w+1)) is a
# multiple of 2^(w+1) and 2^w is added back to it -- so the next
# digit is odd in its turn, and the recoding never needs a zero
m = (m - digit) >> w
# what is left, which the bound above makes an odd digit of its own:
# positive, so the accumulator can start at the entry it names
digits.append(m)
return digits
def _odd_multiples(Q: JacPoint, ec: CurveGroup, w: int) -> list[JacPoint]:
"""Return the odd multiples 1*Q, 3*Q, ..., (2^(w-1) - 1)*Q.
The table a width-w NAF indexes: its digits are odd and smaller than
2^(w-1) in absolute value, so a digit d picks entry (abs(d) - 1) // 2,
to be negated on the fly when d is negative. For w of 1 or 2 the
digits are only +-1 and the table is the point itself.
"""
T = [Q]
if w > 2:
Q2 = ec._double_jac(Q)
for _ in range(2 ** (w - 2) - 1):
T.append(ec._add_jac(T[-1], Q2))
return T
def _odd_multiples_aff(Q: JacPoint, ec: CurveGroup, w: int) -> list[Point]:
"""Return `_odd_multiples` in affine coordinates, one inversion for all.
What it buys is the loop that indexes it: `add_jac_aff` makes five
products fewer than `add_jac` on every addition, and a windowed
multiplication makes one addition per window. What it costs is one
extended Euclid, `aff_from_jac_batch_var` sharing it across the entries --
at w=4 one inverse and some forty products, against five saved on each
of the ~63 additions a 256-bit scalar makes.
The inversion is a function of the point and not of the scalar, so a
multiplication that is regular in its scalar stays regular in it.
"""
return ec._aff_from_jac_batch_var(_odd_multiples(Q, ec, w))
# the width a point of ec._fixed_points has its table built at, against
# the narrow one a point that arrives once is worth. The measurement is in
# _multi_mult_w_NAF_var's docstring
_FIXED_POINT_W = 10
@functools.lru_cache # a fixed point's table is the same on every call
def _cached_odd_multiples_aff(Q: JacPoint, ec: CurveGroup, w: int) -> list[Point]:
"""Return `_odd_multiples_aff` memoized, for a point known in advance.
What memoizing buys is not the building alone: a table that is not
rebuilt can be held at a width a per-call one could not pay for, and a
wNAF has one nonzero digit in w+1, so the wider window is fewer
additions in the loop that indexes it. libsecp256k1 splits the two the
same way, `WINDOW_A` of 5 for the variable point against a `WINDOW_G`
of 15 for the generator, whose table it generates at build time.
Only a point somebody has said will come back may be handed here:
`ec._fixed_points`, or one a caller prepared with `curve.PreparedPoint`. A
table memoized on whatever key arrives is a table built once, used once and
kept for ever, which is why the preparing is a caller's word and never
inferred -- issue btclib-org/btclib#287 is the same bound
`_cached_base58_decode` and `pedersen.second_generator` hold. What bounds it
when the word is given is this cache's own maxsize, 128 tables of 2^(w-2)
points.
"""
return _odd_multiples_aff(Q, ec, w)
def _signed_odd_multiples_aff(Q: JacPoint, ec: CurveGroup, w: int) -> list[Point]:
"""Return the 2^w multiples a signed odd width-w digit names, affine.
-(2^w - 1)*Q first and (2^w - 1)*Q last, so the digit d of
`signed_odd_digits` picks entry (d + 2^w - 1) // 2 and nothing has to
test its sign, where the wNAF loops negate an entry of
`_odd_multiples_aff` on the fly under an `if`. Cheaper, too, and not
only regular: 2^(w-1) modular subtractions once -- 8 of them at w=4 --
against the one per negative digit the loop would make, which is half
of the 64 digits a 256-bit scalar has at that width.
The table it is built on is `_odd_multiples_aff(Q, ec, w + 1)`: a
regular digit reaches 2^w - 1 where a width-w NAF digit stops at
2^(w-1) - 1, so the table is twice as long as that recoding's for the
same window. The negation is one modular subtraction in affine
coordinates as it is in Jacobian ones, `(x, p - y)`, so nothing is
paid for the entries being affine here rather than there.
"""
T = _odd_multiples_aff(Q, ec, w + 1)
return [ec.negate(P) for P in reversed(T)] + T
@functools.lru_cache # the generator is the Q of most calls, so it pays
def _cached_fixed_base_multiples(
Q: JacPoint, ec: CurveGroup, teeth: int, spacing: int
) -> list[list[Point]]:
"""Return the comb's tables: per block, 2^teeth affine points.
Tooth i of block j is the bit (j*teeth + i)*spacing of a scalar. Entry
e of block j is the sum over its teeth of +-2^((j*teeth + i)*spacing)
* Q, + where bit i of e is set and - where it is not. The blocks are
`ceil(ec.scalar_len / (teeth * spacing))`, enough teeth for every bit
of a scalar of the group.
Half of each table is the other half negated, which costs a modular
subtraction an entry. Keeping both halves is what lets the
multiplication index an entry without testing its sign, as
`_signed_odd_multiples_aff` does for the regular window.
Memoized on (Q, ec, teeth, spacing) because that is what makes it
worth building: the table is the point's and not the scalar's, so a
caller multiplying the generator pays for it once -- and so does a
caller multiplying its own point, once it has said with
`curve.PreparedPoint` that the point will come back.
"""
blocks = ceil(ec.scalar_len / (teeth * spacing))
# 2^p * Q for every tooth position p, and twice it, which is what
# turning that tooth from - to + adds: affine, one inversion for all,
# so that the entries below are built by the cheaper mixed addition
chain: list[JacPoint] = []
K = Q
for _ in range(blocks * teeth):
chain.append(K)
K = ec._double_jac(K)
chain.append(K)
for _ in range(spacing - 1):
K = ec._double_jac(K)
aff = ec._aff_from_jac_batch_var(chain)
powers, twice = aff[::2], aff[1::2]
jac: list[JacPoint] = []
for j in range(0, blocks * teeth, teeth):
# the first entry whose top tooth is +: every other tooth is -
entry = _jac_from_aff(powers[j + teeth - 1])
for P in powers[j : j + teeth - 1]:
entry = ec._add_jac_aff(entry, ec.negate(P))
block = [entry]
for i in range(j, j + teeth - 1):
block += [ec._add_jac_aff(P, twice[i]) for P in block]
jac += block
half = 1 << (teeth - 1)
aff = ec._aff_from_jac_batch_var(jac)
tables: list[list[Point]] = []
for at in range(0, len(aff), half):
upper = aff[at : at + half]
# entry e of the lower half is minus entry 2^teeth - 1 - e
tables.append([ec.negate(P) for P in reversed(upper)] + upper)
return tables
def _mult_fixed_base(
m: int, Q: JacPoint, ec: CurveGroup, teeth: int, spacing: int
) -> JacPoint:
"""Scalar multiplication by a signed-digit multi-comb.
Hamburg's comb, eprint 2012/309 section 3.3, which libsecp256k1's
`secp256k1_ecmult_gen_gej` runs for the generator. The teeth of a
block are `spacing` bits apart, and `_cached_fixed_base_multiples`
holds every sum they can name. Each comb offset adds one entry per
block, and a doubling separates one offset from the next. The first
entry starts the accumulator and the parity correction adds one
more, so a scalar costs blocks * spacing additions and spacing - 1
doublings, the same for every scalar of the curve. On secp256k1 at 8
teeth and spacing 8 that is 4 blocks, 32 additions and 7 doublings.
`_CountingGroup` in tests/curves/curve_group_test.py counts them.
Every bit of the recoded scalar is a sign and none is a zero. An odd
k below 2^L, L being the bits the blocks cover, is sum((2*d_i - 1) *
2^i) over the bits d_i of d = (k + 2^L - 1) / 2. That is integer
arithmetic and needs no group order, which `CurveGroup` does not
carry. libsecp256k1 computes d modulo n instead, over a table of G/2
and with a blinding term. So an even m is multiplied as m + 1 and Q
is subtracted at the end, as in `_mult_regular_window`.
The tables are affine, and an affine point of order 2, (x, 0), reads
as infinity. So Q must have no multiple of order 2. Every point of a
subgroup of odd prime order meets that, and `curve.PreparedPoint`
sends any other point to another arm.
8 teeth and spacing 8, by measurement on secp256k1; issue
btclib-org/ecc#194 holds the script and the table.
The reason for keeping a table is the point repeating and not the
point being the generator, so `curve.PreparedPoint` reaches here as
well, for a caller who has said its own point will come back. That
object is where the trade is written.
The accumulator is rescaled where `curves.mult` rescales the point on
its other arm: the table here is memoized and canonical, so it is the
running value that has to stop being a function of m alone.
`_blinded_jac` says what that buys and what it does not.
The input point is assumed to be on the curve, and m to be below 2^L,
which `ec.scalar_len` bounds and every entry point of the library
satisfies, each reducing mod n first.
"""
if m < 0:
raise EccValueError("negative m")
if teeth <= 0:
raise EccValueError(f"non positive teeth: {teeth}")
if spacing <= 0:
raise EccValueError(f"non positive spacing: {spacing}")
T = _cached_fixed_base_multiples(Q, ec, teeth, spacing)
size = teeth * spacing * len(T)
if m >> size:
raise EccValueError(f"m does not fit {size} bits")
d = ((m | 1) + (1 << size) - 1) >> 1
bits = format(d, f"0{size}b")
mask = (1 << teeth) - 1
R = INFJ
for off in range(spacing):
# bit q of the column is bit q*spacing + spacing-1-off of d: the
# tooth q % teeth of block q // teeth, at this offset
column = int(bits[off::spacing], 2)
if off:
R = ec._double_jac(R)
for Tj in T:
R = ec._add_jac_aff(R, Tj[column & mask])
column >>= teeth
else:
# the accumulator starts at a table entry, rescaled
R = _blinded_jac(_jac_from_aff(T[0][column & mask]), ec)
for Tj in T[1:]:
column >>= teeth
R = ec._add_jac_aff(R, Tj[column & mask])
# the parity correction of _mult_regular_window, made whatever the
# parity for the same reason, and what answers m == 0
return ec._add_jac(R, (INFJ, ec.negate_jac(Q))[not m & 1])
def _mult_mont_ladder_var(m: int, Q: JacPoint, ec: CurveGroup) -> JacPoint:
"""Scalar multiplication using 'Montgomery ladder' algorithm.
This implementation uses
'Montgomery ladder' algorithm,
'left-to-right' binary decomposition of the m coefficient,
Jacobian coordinates.
The loop body is one addition and one doubling for every bit of m,
and the bit only picks which of the two accumulators each is applied
to: no test on the scalar, which is what the FLUSH+RELOAD recovery of
OpenSSL's nonces read out of a double-and-add. With add_jac branching
on nothing either, this is the variant of the package that comes
closest to constant-time, and what is left is the loop itself: it
runs once per bit of m, so the length of a scalar is not hidden, only
its bits. See:
- https://eprint.iacr.org/2014/140
The input point is assumed to be on curve and
the m coefficient is assumed to have been reduced mod n
if appropriate (e.g. cyclic groups of order n).
"""
if m < 0:
raise EccValueError("negative m")
# R[0] is the running resultR[1] = R[0] + Q is an ancillary variable
R = [INFJ, Q]
for i in [int(i) for i in f"{m:b}"]:
R[not i] = ec._add_jac(R[i], R[not i])
R[i] = ec._double_jac(R[i])
return R[0]
def _mult_base_3_var(m: int, Q: JacPoint, ec: CurveGroup) -> JacPoint:
"""Scalar multiplication using ternary decomposition of the scalar.
This implementation uses 'triple & add' algorithm, 'left-to-right'
ternary decomposition of the m coefficient, Jacobian coordinates.
The input point is assumed to be on curve and the m coefficient is
assumed to have been reduced mod n if appropriate (e.g. cyclic
groups of order n).
"""
if m < 0:
raise EccValueError("negative m")
# at each step one of the points in T will be added
T = [INFJ, Q, ec._double_jac(Q)]
digits = _convert_number_to_base_var(m, 3)
R = T[digits[0]]
for i in digits[1:]:
# 'triple'
R2 = ec._double_jac(R)
R3 = ec._add_jac(R2, R)
# and 'add'
R = ec._add_jac(R3, T[i])
return R
def _mult_fixed_window_var(
m: int, Q: JacPoint, ec: CurveGroup, w: int, cached: bool
) -> JacPoint:
"""Scalar multiplication using "fixed window".
This implementation uses 'multiple-double & add' algorithm, 'left-to-right'
window decomposition of the m coefficient, Jacobian coordinates.
For 256-bit scalars it is suggested to choose w=4 or w=5.
The input point is assumed to be on curve and the m coefficient is
assumed to have been reduced mod n if appropriate (e.g. cyclic
groups of order n).
"""
if m < 0:
raise EccValueError("negative m")
# a number cannot be written in basis 1 (ie w=0)
if w <= 0:
raise EccValueError(f"non positive w: {w}")
# at each step one of the points in T will be added
T = _cached_multiples(Q, ec) if cached else _multiples(Q, 2**w, ec)
digits = _convert_number_to_base_var(m, 2**w)
R = T[digits[0]]
for i in digits[1:]:
# multiple 'double'
for _ in range(w):
R = ec._double_jac(R)
# and 'add'
R = ec._add_jac(R, T[i])
return R
def _mult_fixed_window_cached_var(
m: int, Q: JacPoint, ec: CurveGroup, w: int
) -> JacPoint:
"""Scalar multiplication using "fixed window" & cached values.
This implementation uses 'multiple-double & add' algorithm, 'left-to-right'
window decomposition of the m coefficient, Jacobian coordinates.
For 256-bit scalars it is suggested to choose w=4. Thanks to the
pre-computed values, it just needs addictions.
The input point is assumed to be on curve and the m coefficient is
assumed to have been reduced mod n if appropriate (e.g. cyclic
groups of order n).
"""
if m < 0:
raise EccValueError("negative m")
# a number cannot be written in basis 1 (ie w=0)
if w <= 0:
raise EccValueError(f"non positive w: {w}")
T = _cached_multiples_fixwind(Q, ec, w)
digits = _convert_number_to_base_var(m, 2**w)
k = len(digits) - 1
R = T[k][digits[0]]
for i in range(1, len(digits)):
k -= 1
# only 'add'
R = ec._add_jac(R, T[k][digits[i]])
return R
def _mult_regular_window(m: int, Q: JacPoint, ec: CurveGroup, w: int) -> JacPoint:
"""Scalar multiplication using a regular window: no digit is zero.
The fixed window above, with the quantity the scalar still decided taken
away from it (issue btclib-org/btclib#254). That one adds a table entry on
every digit, the infinity of T[0] included, so *which* digits a scalar has
does not change what it costs -- but how many it has does: the count is
ceil(m.bit_length() / w), and a scalar short of a full top window runs one
window less. Here the digits are the signed odd ones of signed_odd_digits,
none of them zero, and there are ceil(ec.scalar_len / w) of them whatever m
is: one addition and w doublings per window, for every scalar of every size.
On secp256k1 that is 71 additions and 253 doublings at w=4, the same for
every scalar and counted at the curve's own scalar_len of 256, where the
fixed window's count follows the size of the scalar it is given. What counts
them is `_CountingGroup` in tests/curves/curve_group_test.py, which carries
the command in its docstring and the reason for that length where it sets
it.
Uniformity for free, in other words -- indistinguishable in cost from
the fixed window over 30 random 256-bit scalars, best of five -- which
is what makes this the `_mult` the package reaches for rather than an
option beside it. The fixed window stays as what it is didactically,
the plain one, and the wNAF multiplications stay the fast irregular
ones for coefficients that are public: `_double_mult_w_NAF_var` is
what verification runs.
w=4 by measurement, as for the fixed window: w=5 and w=6 land within
a couple of percent either side of it, which buys nothing against a
table twice the size of the NAF's at the same width.
Two things follow for the group law underneath. The accumulator starts
at a table entry instead of at infinity and no digit names infinity, so
the identity is unreachable from the loop rather than merely uniform in
it; and the two special cases left in `add_jac`, a doubling and a sum
that is infinity, need the accumulator to land on a table entry or its
opposite, 2^-250 on a curve with a real order. On the toy curves of
the test suite they are reached constantly, an accumulator wrapping
around a group of eleven points, which is what holds them tested.
The input point is assumed to be on curve. m is *not* assumed to have
been reduced: a scalar above ec.scalar_len bits is multiplied in the
digits it needs, which is the one case where the count is the scalar's
own again -- and no entry point of the library forms one, each reducing
mod n first.
"""
if m < 0:
raise EccValueError("negative m")
# a number cannot be written in basis 1 (ie w=0)
if w <= 0:
raise EccValueError(f"non positive w: {w}")
size = ceil(max(ec.scalar_len, m.bit_length()) / w)
# the recoding needs an odd scalar, and m | 1 is m for an odd one and
# m + 1 for an even one. The final addition below takes the 1 back
# rather than n being added to make m odd: n is not on CurveGroup at
# all, and m + n would be the same point only for a point whose order
# it is -- which the group this multiplies in cannot say
digits = signed_odd_digits(m | 1, w, size)
T = _signed_odd_multiples_aff(Q, ec, w)
offset = (1 << w) - 1
# a table entry, the top digit being positive, and never infinity
R = _jac_from_aff(T[(digits[-1] + offset) // 2])
for digit in digits[-2::-1]:
# multiple 'double'
for _ in range(w):
R = ec._double_jac(R)
# and 'add', on every digit: there is no zero one to skip
R = ec._add_jac_aff(R, T[(digit + offset) // 2])
# the parity correction, made whatever the parity so that it cannot be
# read off the clock: one addition of infinity, which by the same
# property costs what the addition of -Q costs. It is also what
# answers m == 0, whose recoding is that of 1
return ec._add_jac(R, (INFJ, ec.negate_jac(Q))[not m & 1])
# the width _mult hands the regular window; the measurement behind it is
# in the docstring of the function it is passed to
_MULT_W = 4
def _mult(m: int, Q: JacPoint, ec: CurveGroup) -> JacPoint:
"""Return m*Q, with what the rest of the package multiplies with.
On every curve but secp256k1, where curves.mult has the GLV
endomorphism: the regular window, whose cost is the same for every
scalar of the curve, at the cost of nothing measurable against the
fixed window it recodes.
A function and not an alias of _mult_regular_window, the width being an
argument: a private signature carries no default, so the width this one
picks belongs here, where the choice is, rather than in the signature of
a variant every other caller passes its own width to.
"""
return _mult_regular_window(m, Q, ec, _MULT_W)
def _double_mult_var(
u: int, HJ: JacPoint, v: int, QJ: JacPoint, ec: CurveGroup
) -> JacPoint:
"""Double scalar multiplication (u*H + v*Q).
This implementation uses the Shamir-Strauss algorithm, 'left-to-right'
binary decomposition of the u and v coefficients, Jacobian coordinates.
Strauss algorithm consists of a single 'double & add' loop for the
parallel calculation of u*H and v*Q, efficiently using a single
'doubling' for both scalar multiplications (see
https://stackoverflow.com/questions/50993471/ec-scalar-multiplication-with-strauss-shamir-method).
The Shamir trick adds the precomputation of H+Q,
which is to be added in the loop when the binary digits
of u and v are both equal to 1 (on average 1/4 of the cases).
The input points are assumed to be on curve,
the u and v coefficients are assumed to have been reduced mod n
if appropriate (e.g. cyclic groups of order n).
"""
if u < 0:
raise EccValueError("negative first coefficient")
if v < 0:
raise EccValueError("negative second coefficient")
# at each step one of the following points will be added
T = [INFJ, HJ, QJ, ec._add_jac(HJ, QJ)]
# which one depends on binary digit for that step
ui = f"{u:b}"
vi = f"{v:b}".zfill(len(ui))
ui = ui.zfill(len(vi))
digits = [int(j) + 2 * int(k) for j, k in zip(ui, vi, strict=True)]
# R[0] is the running result, R[1] = R[0] + T[*] is an ancillary variable
R = T[digits[0]]
for i in digits[1:]:
# the doubling part of 'double & add'
R = ec._double_jac(R)
# always perform the 'add', even if the digit pair is 0 and the
# entry it names is infinity: one addition per step, whatever the
# coefficients. Which is where an add_jac that does not shortcut
# infinity is dearest -- a quarter of the pairs are 0 on average,
# each an addition a shortcut would answer for free -- and it is
# the one place in the package that measures 22% slower for it,
# against the 2% to 3% everywhere else
R = ec._add_jac(R, T[i])
return R
def _multi_mult_pairs(
scalars: Sequence[int], jac_points: Sequence[JacPoint]
) -> list[tuple[int, JacPoint]]:
"""Check the arguments of a multi multiplication, drop what adds nothing.
Shared by the two implementations below so that they answer the same
errors to the same arguments: whichever of them _multi_mult_var calls is
a performance decision, and a caller must not be able to tell which
it got. A zero scalar contributes nothing to the sum and is dropped
rather than carried, which Bos-Coster needs -- it would be a zero
divisor in its Euclidean step -- and interleaved wNAF is glad of,
a scalar of zero costing a table it would never index.
A point at infinity contributes nothing either, and libsecp256k1's
`secp256k1_ecmult_strauss_wnaf` skips it too. In the interleaved wNAF
a point kept costs a table of odd multiples like any other (issue
btclib-org/ecc#198 has the script that measures it).
"""
if len(scalars) != len(jac_points):
err_msg = "mismatch between number of scalars and points: "
err_msg += f"{len(scalars)} vs {len(jac_points)}"
raise EccValueError(err_msg)
if len(scalars) < 2:
raise EccValueError("not a multi_mult_var")
pairs: list[tuple[int, JacPoint]] = []
for n, PJ in zip(scalars, jac_points, strict=True):
if n < 0:
raise EccValueError("negative coefficient")
if n and PJ[2]:
pairs.append((n, PJ))
return pairs
# the width _multi_mult_var hands the interleaved wNAF; the measurement
# behind it is in the docstring of the function it is passed to
_MULTI_MULT_W = 5
def _additions_by_position_var(
nafs: Sequence[list[int]], tables: Sequence[list[Point]], ec: CurveGroup
) -> list[list[Point]]:
"""Return the points each bit position adds, indexed by position.
An entry per position, holding the table entry every wNAF with a nonzero
digit there names -- negated where the digit is, a signed digit naming a
point and its opposite. The alternative is every wNAF being asked for a
digit at every position: a wNAF has one nonzero digit in w+1, so a
secp256k1 double multiplication would ask a question at every position of
all four halves the GLV split makes -- 512 at the split's nominal 128 a
half -- and most of them would be answered by a `zip` step, a comparison
against a `len` recomputed per iteration, and an index that finds a zero,
rather than by an addition (issue btclib-org/btclib#906).
Read this way round the digits are walked once, where the alternative
reads each one at the position it belongs to and skips it at every other,
and the loop that consumes the answer is a doubling and its additions.
The table index and the negation are here rather than there for the same
reason: they are what a digit means, and the loop has no use for the
digit itself.
The doublings and the additions are the same operations in the same
order, this being where a digit is looked up and not what is done with
it -- `test_the_interleaved_loop_makes_the_operations_its_wnafs_name`
asserts the sequence of both, with the point each addition takes,
against the sequence the wNAFs name. A comparison of results would be
blind to a reordering that answers the same point, group addition
commuting.
The gain is CPython's: the loop control this avoids is what a tracing
JIT flattens on its own, and under PyPy the two spellings measure the
same. btclib-org/btclib's changelog/v2026.8.21.md has the CPython and
PyPy figures.
"""
at_position: list[list[Point]] = [[] for _ in range(max(len(naf) for naf in nafs))]
for naf, T in zip(nafs, tables, strict=True):
for j, d in enumerate(naf):
if d:
at_position[j].append(
T[(d - 1) // 2] if d > 0 else ec.negate(T[(-d - 1) // 2])
)
return at_position
def _multi_mult_w_NAF_var(
scalars: Sequence[int],
jac_points: Sequence[JacPoint],
ec: CurveGroup,
w: int,
fixed: frozenset[JacPoint],
) -> JacPoint:
"""Return the multi scalar multiplication u1*Q1 + ... + un*Qn.
Strauss' algorithm with interleaved wNAFs, the many-point form of
curve_group_2's _double_mult_w_NAF_var (HMV algorithm 3.51): each scalar
gets its own width-w NAF and its own table of odd multiples, and a
single left-to-right loop shares one doubling per bit position among
all of them, adding a (possibly negated) table entry wherever a NAF
has a nonzero digit.
Sharing the doublings is the whole of the idea, and what the dispatch
of _multi_mult_var weighs: the ~256 doublings a 256-bit scalar needs are
paid once for the batch instead of once per point, leaving some 50
point operations per scalar of its own. Bos-Coster has no fixed cost
to share and settles at 44.5 per scalar, so few scalars is where the
sharing pays and many is where the per-scalar figure is all there is.
Cost here is a function of the bit length of the scalars and of nothing
else, so no pair of magnitudes is pathological the way 10**6 next to 1 is
for a subtractive Bos-Coster step (issue btclib-org/btclib#175).
w=5 by measurement, on random 256-bit scalars: the table costs
2^(w-2) - 1 additions and a doubling per point while the NAF has one
nonzero digit in w+1, so the total falls to w=5 and rises again at
w=6 (per scalar at 128 scalars: 57.7 operations at w=4, 53.3 at w=5,
55.1 at w=6). In milliseconds over 10 scalars: 7.81 at w=4, 7.27 at
w=5, 7.38 at w=6, 8.68 at w=7.
No GLV splitting of its own, which would halve the doublings on
secp256k1: the split is curve-specific, this function serves every
curve, and with many scalars the doublings it would halve are already
shared. A caller that has only two of them has the opposite balance,
and curve_group_2's _double_mult_endomorphism_secp256k1_var is that
caller: it splits both coefficients itself and hands the four halves
here, so the split lives with the curve that has it and the
interleaving stays one implementation.
`fixed` is the points whose table is worth memoizing rather than
building here: `ec._fixed_points`, or that set with the endomorphism's
images added by the caller that forms them. A memoized table can be
held at a width a per-call one could not pay for, `_FIXED_POINT_W`
against the `w` every other point gets, and a wNAF has one nonzero
digit in w+1, so the wide window is fewer additions in the loop. That
is libsecp256k1's `WINDOW_G` of 15 beside its `WINDOW_A` of 5, on a
table it generates at build time.
Measured over 30 random 256-bit coefficient pairs, best of seven
alternating rounds, memoized against built per call: 1.24x on
`_double_mult_endomorphism_secp256k1_var`, whose four halves are two fixed
points and two of the caller's, and 1.14x on `_double_mult_w_NAF_var` over
secp256r1, which has one of each. With eight points of which one is
the generator it is 1.05x, the fixed point being an eighth of the
work. `_FIXED_POINT_W` is 10 by measurement on the first of those:
1.07x at 4, 1.12x at 6, 1.15x at 8, 1.18x at 10, 1.20x at 12, against
a table that doubles at every step -- and capped at `ec.scalar_len`,
a window wider than the scalars buying nothing and costing a table.
The input points are assumed to be on curve, the scalar coefficients
are assumed to have been reduced mod n if appropriate (e.g. cyclic
groups of order n).
"""
# a number cannot be written in basis 1 (ie w=0)
if w <= 0:
raise EccValueError(f"non positive w: {w}")
pairs = _multi_mult_pairs(scalars, jac_points)
if not pairs:
return INFJ
# a point the caller hands the same on every call gets the memoized
# table of _cached_odd_multiples_aff, at the wider window that only a
# memoized table can pay for; every other point gets its own, built
# here and dropped with the call
# a window wider than the scalars themselves buys nothing and costs a
# table, which is what the cap is for: the low-cardinality curves of
# the test suite have a scalar_len of a handful of bits
fixed_w = min(_FIXED_POINT_W, ec.scalar_len)
nafs = []
tables: list[list[Point]] = []
# the tables this call has to build, kept in Jacobian coordinates
# until they can be converted together
pending: list[tuple[int, list[JacPoint]]] = []
for i, (n, PJ) in enumerate(pairs):
width = fixed_w if PJ in fixed else w
nafs.append(_wNAF_of_m_var(n, width))
if PJ in fixed:
tables.append(_cached_odd_multiples_aff(PJ, ec, width))
else:
# a place in the list, filled by the conversion below
tables.append([])
pending.append((i, _odd_multiples(PJ, ec, width)))
# one extended Euclid for every table the call builds, where a table
# at a time is one apiece: the batch of `aff_from_jac_batch_var` over the
# concatenation rather than over one point's share of it. That is
# libsecp256k1's `secp256k1_ge_set_all_gej_var` over the whole batch;
# its own multiplication does not call it, its tables sharing one Z by
# `secp256k1_ge_table_set_globalz` instead. An empty concatenation is
# not a case to test for: every point being memoized leaves nothing to
# convert and nothing to do
aff = ec._aff_from_jac_batch_var([P for _, jac in pending for P in jac])
at = 0
for i, jac in pending:
tables[i] = aff[at : at + len(jac)]
at += len(jac)
# one doubling per position and the additions that position has:
# `_additions_by_position_var` groups them this way rather than asking
# every wNAF for a digit at every position, and is also where a
# digit becomes the point it names
R = INFJ
for additions in reversed(_additions_by_position_var(nafs, tables, ec)):
R = ec._double_jac(R)
for P in additions:
R = ec._add_jac_aff(R, P)
return R
def _multi_mult_bos_coster_var(
scalars: Sequence[int], jac_points: Sequence[JacPoint], ec: CurveGroup
) -> JacPoint:
"""Return the multi scalar multiplication u1*Q1 + ... + un*Qn.
Use Bos-Coster's algorithm for efficient computation.
A max-heap holds the (scalar, point) pairs still to be summed; each
step pops the two largest scalars and rewrites their contribution as
n1*P1 + n2*P2 = r*P1 + n2*(q*P1 + P2), with n1 = q*n2 + r
an identity, so the sum over the heap never changes: it is the
invariant, and the single pair left at the end is the answer. The
step also terminates, because r < n2 <= n1 makes the largest scalar
strictly decrease, and r == 0 costs the heap an entry outright.
Bos-Coster is usually written with the q == 1 case of that identity, n1*P1 +
n2*P2 = (n1-n2)*P1 + n2*(P1+P2), which is Euclid by repeated subtraction:
n1/n2 steps to do what one divmod does (issue btclib-org/btclib#175) --
multi_mult_var([10**6, 1], [G, H]) would take on the order of a million
steps of Euclid, and both multi_mult_var([n-1, 1], [G, H]) and
multi_mult_var([-1, 1], [G, H]) would never finish, on scalars a caller has
every right to pass. The whole of Euclid is what makes this a live path
rather than a reading exercise: _multi_mult_var calls it above its
threshold.
Where the interleaved wNAF above shares one doubling per bit among
all the points, this one has no fixed cost to amortize at all -- the
heap converges towards scalars of the same size, and the quotient is
1 about 40% of the time (Gauss-Kuzmin), so a step is usually one
addition. That is why it wins on many scalars and loses on few: 44.5
point operations per scalar over 128 random 256-bit scalars, against
267.5 over two.
The input points are assumed to be on curve, the scalar coefficients
are assumed to have been reduced mod n if appropriate (e.g. cyclic
groups of order n).
"""
# source: https://cr.yp.to/badbatch/boscoster2.py
pairs = _multi_mult_pairs(scalars, jac_points)
if not pairs:
return INFJ
# a max-heap, which Python spells as a min-heap of negated scalars
x: list[tuple[int, JacPoint]] = [(-n, PJ) for n, PJ in pairs]
heapq.heapify(x)
while len(x) > 1:
np1 = heapq.heappop(x)
np2 = heapq.heappop(x)
n_1, p_1 = -np1[0], np1[1]
n_2, p_2 = -np2[0], np2[1]
# _mult_jac_var, not the faster _mult: the quotient is 1 about 40% of
# the time (Gauss-Kuzmin), and _mult would rebuild its 16-entry
# table of multiples for each of those -- 10x slower over 128
# random scalars. _mult_jac_var(1, P) is P at no cost, so this is the
# subtractive step precisely when the subtractive step is right;
# measured, spelling that case out as a branch buys nothing, and
# subtracting for q up to 4 or up to 16 costs 3% and 11%.
q, n_1 = divmod(n_1, n_2)
p_2 = ec._add_jac(_mult_jac_var(q, p_1, ec), p_2)
if n_1 > 0:
heapq.heappush(x, (-n_1, p_1))
heapq.heappush(x, (-n_2, p_2))
np1 = heapq.heappop(x)
n_1, p_1 = -np1[0], np1[1]
# n_1 is left as it is rather than reduced mod n: _mult reduces, and
# the scalars reaching here are already below the order
return _mult(n_1, p_1, ec)
# the number of nonzero scalars from which Bos-Coster is the cheaper of
# the two, measured here rather than borrowed: libsecp256k1 dispatches
# at 88, but it dispatches between other algorithms, on machine words,
# in C, and a number that has not been measured on Python integers is
# decoration.
# Best of five interleaved reps, random 256-bit scalars on secp256k1,
# Python 3.14, macOS arm64, wNAF over Bos-Coster: 0.91 at 40 scalars,
# 0.94 at 48, 1.000 at 56, 1.001 at 64, 1.04 at 72, 1.09 at 96. The two
# meet at 56, so 56 it is; anywhere in 56..64 measures the same, and the
# cost of being wrong by a dozen scalars either way is under 1%.
BOS_COSTER_THRESHOLD = 56
def _multi_mult_var(
scalars: Sequence[int], jac_points: Sequence[JacPoint], ec: CurveGroup
) -> JacPoint:
"""Return the multi scalar multiplication u1*Q1 + ... + un*Qn.
Interleaved wNAF below BOS_COSTER_THRESHOLD kept pairs, each a nonzero
scalar on a finite point, and Bos-Coster from there on, the two crossing
where the shared doublings stop paying for themselves (issue
btclib-org/btclib#212). Both are kept, and not only to be
dispatched between: the library is didactic as much as it is a library, and
Bos-Coster is the one that can be read in twenty lines.
Pippenger, the third algorithm libsecp256k1 has, is *not* here, and
the measurement is the reason. Prototyped with sparse buckets and the
running-sum trick, best window per size, it lost at every size tried:
against Bos-Coster in ms, 162 vs 125 at 256 scalars, 284 vs 225 at
512, 500 vs 406 at 1024. Its bucket sums are additions of points that
cost the same as any other addition in Python, where in C they are
the cheap part, so the crossover it has there is beyond any batch
this package is asked to verify -- and a threshold nobody reaches is a
branch nobody takes.
The zeros and the points at infinity are counted out before the dispatch
rather than left to len(), because both are dropped downstream: 56
scalars of which 2 are nonzero is a batch of two, and sending it to
Bos-Coster on its length costs more than the wNAF does.
The input points are assumed to be on curve, the scalar coefficients
are assumed to have been reduced mod n if appropriate (e.g. cyclic
groups of order n).
"""
# a length mismatch falls through to _multi_mult_pairs, which raises it
kept = sum(bool(n and PJ[2]) for n, PJ in zip(scalars, jac_points, strict=False))
if kept < BOS_COSTER_THRESHOLD:
return _multi_mult_w_NAF_var(
scalars, jac_points, ec, _MULTI_MULT_W, ec._fixed_points
)
return _multi_mult_bos_coster_var(scalars, jac_points, ec)