Notebook

Sage is not required

In [3]:

from tqdm import tqdm

In [4]:

def legendre(a, p):
    return pow(a, (p - 1) // 2, p)
 
def tonelli(n, p):
    assert legendre(n, p) == 1, "not a square (mod p)"
    q = p - 1
    s = 0
    while q % 2 == 0:
        q //= 2
        s += 1
    if s == 1:
        return pow(n, (p + 1) // 4, p)
    for z in range(2, p):
        if p - 1 == legendre(z, p):
            break
    c = pow(z, q, p)
    r = pow(n, (q + 1) // 2, p)
    t = pow(n, q, p)
    m = s
    t2 = 0
    while (t - 1) % p != 0:
        t2 = (t * t) % p
        for i in range(1, m):
            if (t2 - 1) % p == 0:
                break
            t2 = (t2 * t2) % p
        b = pow(c, 1 << (m - i - 1), p)
        r = (r * b) % p
        c = (b * b) % p
        t = (t * c) % p
        m = i
    return r

Prerequisites¶

Elliptic curves

Reminder
Weierstrass curve $$y^2 = x^3 + ax + b$$ with $4a^3 + 27b^2 \neq 0$

Theory¶

Costello and Smith paper - must read, I followed from here

Definition

A Montgomery curve over a field $K$ is defined by the equation $$\boxed{M_{A,B}:By^{2}=x^{3}+Ax^{2}+x}$$ for 2 parameters $A, B \in K$ and with $B(A^2 − 4) \neq 0$.

Or in projective coordinates $$BY^2Z = X^3 + AX^2Z + XZ^2 \subseteq \mathbb P^2$$ With $\mathcal O = (0: 1: 0)$ being the point at infinity

J-invariant $$J(M_{A, B}) = \dfrac {256(A^2 - 3)^3} {A^2 - 4}$$

Notice that the $\overline K$ isomorphism class is determined only by $A$

In [26]:

def is_on(P, C):
    """
    Checks if a point P is on the montgomery curve C
    """
    A, B, p = C
    x, y = P
    return (B* y**2) % p == (x**3 + A * x**2 + x) % p

In [14]:

def lift(x, C):
    """
    Lifts a coordinate `x` to a point on the curve `C` if possible
    """
    A, B, p = C
    y = tonelli((x**3 + A * x**2 + x) % p, p)
    return (x, y)

The group structure¶

Identity: $\mathcal{O}$ acts as identity
Negation: $(x:y:1) \mapsto (x:-y:1)$

Addition
Let $P = (x_P, y_P), \ Q = (x_Q, y_Q)$ be a two points on our curve. Then for $R = P + Q = (x_R, y_R)$ we have

$x_R = B\lambda ^2 - (x_P + x_Q) - A$
$y_R = (2x_P + x_Q + A)\lambda - B\lambda ^3 - y_P

= \lambda(x_P - x_R) - yP$

For $\lambda = \begin{cases} \cfrac {y_Q - y_P} {x_Q - x_P} \text{ if } P \neq Q, \ P \neq -Q \\ \cfrac {3x^2_P + 2Ax_P + 1)} {2By_P} \text{ if } P = Q \end{cases} $

In [23]:

from Crypto.Util.number import inverse

def montgomery_addition(P1, P2, C):
    """
    Adds 2 points P1, P2 on the curve C
    """
    
    assert(is_on(P1, C)), f"The point {P1} is not on curve {C}"
    assert(is_on(P2, C)), f"The point {P2} is not on curve {C}"
    
    A, B, p = C
    x1, y1 = P1
    x2, y2 = P2
    
    #assert(x1 != x2 and y1 != y2), "P must not be equal to Q"
    lam = ((y2 - y1) * inverse(x2 - x1, p)) % p
    x3 = (B * lam ** 2 - A - x1 - x2) % p
    y3 = (lam * (x1 - x3) - y1) % p
    
    return (x3, y3)

def montgomery_doubling(P, C):
    """
    Double a point P on the curve C
    """
    
    assert(is_on(P, C)), f"The point {P} is not on curve {C}"
    
    A, B, p = C
    x, y = P
    
    lam = ((3*x**2 + 2*A*x + 1) * inverse(2 * B * y, p)) % p
    x3 = (B * lam ** 2 - A - 2 * x) % p
    y3 = (lam * (x - x3) - y) % p
    
    return (x3, y3)

Fast arithmetic¶

The paper presents pseudo-addition and pseudo-doubling which perform x-addition and x-doubling using $x(P), x(Q), x(P-Q)$ with the xADD and xDBL algorithms

$+$ efficient x-line arithmetic

Montgomery ladder¶

This is algorithm 3 from the paper and it is used to compute $P \mapsto [k]P$ on a curve

In [24]:

def montgomery_binary(P, k, C):
    """
    Return k*P on C given P, k, C
    """
    assert(is_on(P, C)), f"The point {P} is not on curve {C}"

    R0 = P
    R1 = montgomery_doubling(P, C)
    l = k.bit_length()
    assert((k >> (l-1)) & 1 == 1), "the l-bit k must have the l-1 bit equal to 1"
    for i in range(l-2, -1, -1):
        if (k>>i) & 1 == 0:
            R0, R1= montgomery_doubling(R0, C), montgomery_addition(R0, R1, C)
        else:
            R0, R1= montgomery_addition(R0, R1, C), montgomery_doubling(R1, C)
    
    return R0

In [36]:

p = (1<<255) - 19
B = 1
A = 486662
curve = (A, B, p)
x = 4
P = lift(x, curve)
k = 0xfeed

Q = montgomery_binary(P, k, curve)

print(Q)

(1681274011185213290829845720121217409632760205352840302183329577742232042267, 11000771751458358810062840331211779070151159325452818737845305294782515669468)

Correspondence with other models¶

Short Weierstrass models
From Montgomery to Weierstrass map $ M_{A,B}\rightarrow E_{a,b}$

$$(x,y)\mapsto (t,v)=\left({\frac {x}{B}}+{\frac {A}{3B}},{\frac {y}{B}}\right),a={\frac {3-A^{2}}{3B^{2}}},b={\frac {2A^{3}-9A}{27B^{3}}}$$

The reverse map does not always exist. In fact, in order to transform an Weierstrass form curve into an Montgomery curve we have to satisfy the following conditions:

$E_{a, b}$ has order divisible by 4
$z^{3}+az+b=0$ has at least one root $ \alpha \in \mathbb {F}$
$3\alpha ^{2}+a$ is a quadratic residue in $ \mathbb {F}$

In [54]:

def is_on_w(P, C):
    """
    Checks if a point P is on the Weierstrass curve C
    """
    a, b, p = C
    x, y = P
    return y**2 % p == (x**3 + a * x + b) % p

In [52]:

def montgomery_to_weierstrass_curve(C):
    """
    Returns the weierstrass form curve from a given Montgomery curve
    """
    A, B, p = C
    
    a = ((3 - A**2) * inverse(3 * B ** 2, p)) % p
    b = ((2 * A**3 - 9*A) * inverse(27 * B**3, p)) % p
    
    return (a, b, p)

def montgomery_to_weierstrass_point(P, CM, CW):
    """
    Maps a point on a montgomery curve CM onto one on the weierstrass curve CW
    """
    
    assert(CW == montgomery_to_weierstrass_curve(CM)), f"The Weierstrass curve {CW} is not mapped from the Montgomery curve {CM}"
    A, B, p = CM
    
    x, y = P
    t = (x * inverse(B, p) + A * inverse(3 * B, p)) % p
    v = (y * inverse(B, p)) % p
    
    Q = (t, v)
    assert(is_on_w(Q, CW))
    
    return Q

In [35]:

p = (1<<255) - 19
B = 1
A = 486662
curve = (A, B, p)
x = 4
P = lift(x, curve)
k = 0xfeed

Q = montgomery_binary(P, k, curve)

In [82]:

curve_w = montgomery_to_weierstrass_curve(curve)
print(curve_w)

T = montgomery_to_weierstrass_point(P, curve, curve_w)
print(T)

S = montgomery_to_weierstrass_point(Q, curve, curve_w)
print(S)



# Can't be bothered to implement double and add in Weierstrass form
from ecdsa.ellipticcurve import CurveFp, Point

a, b, p = curve_w
CW = CurveFp(p, a, b)
T_ = Point(CW, *T)
S_ = Point(CW, *S)

print(S_ == k * T_)

(19298681539552699237261830834781317975544997444273427339909597334573241639236, 55751746669818908907645289078257140818241103727901012315294400837956729358436, 57896044618658097711785492504343953926634992332820282019728792003956564819949)
(19298681539552699237261830834781317975544997444273427339909597334652188435541, 47499954730490638715091885595963621956489259355261559690379252041373947974816)
(20979955550737912528091676554902535385177757649626267642092926912394420477804, 11000771751458358810062840331211779070151159325452818737845305294782515669468)
True

Resources¶

In [ ]: