Add docs and tests to the ECTester target.

2 files changed, 124 insertions, 3 deletions

diff --git a/pyecsca/ec/curve.py b/pyecsca/ec/curve.py
index e8ae66c..a05ff3f 100644
--- a/pyecsca/ec/curve.py
+++ b/pyecsca/ec/curve.py
@@ -6,8 +6,8 @@ from public import public
 
 from .coordinates import CoordinateModel, AffineCoordinateModel
 from .mod import Mod
-from .model import CurveModel
-from .point import Point
+from .model import CurveModel, ShortWeierstrassModel
+from .point import Point, InfinityPoint
 
 
 @public
@@ -93,17 +93,62 @@ class EllipticCurve(object):
 
     @property
     def neutral_is_affine(self):
+        """Whether the neurtal point is an affine point."""
         return bool(self.model.base_neutral)
 
     def is_neutral(self, point: Point) -> bool:
+        """Check whether the point is the neutral point."""
         return self.neutral == point
 
     def is_on_curve(self, point: Point) -> bool:
+        """Check whether the point is on the curve."""
         if point.coordinate_model.curve_model != self.model:
             return False
+        if self.is_neutral(point):
+            return True
         loc = {**self.parameters, **point.to_affine().coords}
         return eval(compile(self.model.equation, "", mode="eval"), loc)
 
+    def to_affine(self) -> "EllipticCurve":
+        """Convert this curve into the affine coordinate model, if possible."""
+        coord_model = AffineCoordinateModel(self.model)
+        return EllipticCurve(self.model, coord_model, self.prime, self.neutral.to_affine(), self.parameters)
+
+    def decode_point(self, encoded: bytes) -> Point:
+        """Decode a point encoded as a sequence of bytes (ANSI X9.62)."""
+        if encoded[0] == 0x00 and len(encoded) == 1:
+            return InfinityPoint(self.coordinate_model)
+        coord_len = (self.prime.bit_length() + 7) // 8
+        if encoded[0] in (0x04, 0x06):
+            data = encoded[1:]
+            if len(data) != coord_len * len(self.coordinate_model.variables):
+                raise ValueError("Encoded point has bad length")
+            coords = {}
+            for var in sorted(self.coordinate_model.variables):
+                coords[var] = Mod(int.from_bytes(data[:coord_len], "big"), self.prime)
+                data = data[coord_len:]
+            return Point(self.coordinate_model, **coords)
+        elif encoded[0] in (0x02, 0x03):
+            if isinstance(self.coordinate_model, AffineCoordinateModel) and isinstance(self.model, ShortWeierstrassModel):
+                data = encoded[1:]
+                if len(data) != coord_len:
+                    raise ValueError("Encoded point has bad length")
+                x = Mod(int.from_bytes(data, "big"), self.prime)
+                rhs = x**3 + self.parameters["a"] * x + self.parameters["b"]
+                if not rhs.is_residue():
+                    raise ValueError("Point not on curve")
+                sqrt = rhs.sqrt()
+                yp = encoded[0] & 0x01
+                if int(sqrt) & 0x01 == yp:
+                    y = sqrt
+                else:
+                    y = -sqrt
+                return Point(self.coordinate_model, x=x, y=y)
+            else:
+                raise NotImplementedError
+        else:
+            raise ValueError(f"Wrong encoding type: {hex(encoded[0])}, should be one of 0x04, 0x06, 0x02, 0x03 or 0x00")
+
     def __eq__(self, other):
         if not isinstance(other, EllipticCurve):
             return False
diff --git a/pyecsca/ec/mod.py b/pyecsca/ec/mod.py
index ec5dfe6..421f521 100644
--- a/pyecsca/ec/mod.py
+++ b/pyecsca/ec/mod.py
@@ -1,5 +1,6 @@
+import random
 import secrets
-from functools import wraps
+from functools import wraps, lru_cache
 
 from public import public
 
@@ -39,6 +40,34 @@ def extgcd(a, b):
     return x2, y2, a
 
 
+@public
+@lru_cache
+def miller_rabin(n: int, rounds: int = 50) -> bool:
+    """Miller-Rabin probabilistic primality test."""
+    if n == 2 or n == 3:
+        return True
+
+    if n % 2 == 0:
+        return False
+
+    r, s = 0, n - 1
+    while s % 2 == 0:
+        r += 1
+        s //= 2
+    for _ in range(rounds):
+        a = random.randrange(2, n - 1)
+        x = pow(a, s, n)
+        if x == 1 or x == n - 1:
+            continue
+        for _ in range(r - 1):
+            x = pow(x, 2, n)
+            if x == n - 1:
+                break
+        else:
+            return False
+    return True
+
+
 def check(func):
     @wraps(func)
     def method(self, other):
@@ -99,6 +128,53 @@ class Mod(object):
     def __invert__(self):
         return self.inverse()
 
+    def is_residue(self):
+        """Whether this element is a quadratic residue (only implemented for prime modulus)."""
+        if not miller_rabin(self.n):
+            raise NotImplementedError
+        if self.x == 0:
+            return True
+        if self.n == 2:
+            return self.x in (0, 1)
+        legendre = self ** ((self.n - 1) // 2)
+        return legendre == 1
+
+    def sqrt(self):
+        """
+        The modular square root of this element (only implemented for prime modulus).
+
+        Uses the `Tonelli-Shanks <https://en.wikipedia.org/wiki/Tonelli–Shanks_algorithm>`_ algorithm.
+        """
+        if not miller_rabin(self.n):
+            raise NotImplementedError
+        q = self.n - 1
+        s = 0
+        while q % 2 == 0:
+            q //= 2
+            s += 1
+
+        z = 2
+        while Mod(z, self.n).is_residue():
+            z += 1
+
+        m = s
+        c = Mod(z, self.n) ** q
+        t = self ** q
+        r_exp = (q + 1) // 2
+        r = self ** r_exp
+
+        while t != 1:
+            i = 1
+            while not (t ** (2**i)) == 1:
+                i += 1
+            two_exp = m - (i + 1)
+            b = c ** int(Mod(2, self.n)**two_exp)
+            m = int(Mod(i, self.n))
+            c = b ** 2
+            t *= c
+            r *= b
+        return r
+
     @check
     def __mul__(self, other):
         return Mod((self.x * other.x) % self.n, self.n)