From e4d02c940ac99719503f8dfd678e2a2cc2ae555f Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Sun, 6 Oct 2019 16:58:07 +0200
Subject: [PATCH 01/30] add implementation of EC using Jacobi representation

---
 src/ecdsa/ellipticcurve.py | 125 ++++++++++++++++++++++++++++++++++++-
 src/ecdsa/test_jacobi.py   |  59 +++++++++++++++++
 2 files changed, 183 insertions(+), 1 deletion(-)
 create mode 100644 src/ecdsa/test_jacobi.py

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index e78314cc..3d23db6a 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -57,11 +57,130 @@ def b(self):
 
   def contains_point(self, x, y):
     """Is the point (x,y) on this curve?"""
-    return (y * y - (x * x * x + self.__a * x + self.__b)) % self.__p == 0
+    return ((y**2) - (x**3 + self.__a * x + self.__b)) % self.__p == 0
 
   def __str__(self):
     return "CurveFp(p=%d, a=%d, b=%d)" % (self.__p, self.__a, self.__b)
 
+
+class PointJacobi(object):
+  """Point on an elliptic curve. Uses Jacobi representation"""
+  def __init__(self, curve, x, y, z, order=None):
+      self.__curve = curve
+      self.__x = x
+      self.__y = y
+      self.__z = z
+      self.__order = order
+
+  def __eq__(self, other):
+      """Compare two points with each-other."""
+      if self.__y == 0 and other == INFINITY:
+          return True
+      if isinstance(other, Point):
+          return self.to_weierstrass() == other
+      return self.x() == other.x() and self.y() == other.y()
+
+  def order(self):
+      return self.__order
+
+  def curve(self):
+      return self.__curve
+
+  def x(self):
+      """Return x coordinate in Weierstrass form"""
+      p = self.__curve.p()
+      z = numbertheory.inverse_mod(self.__z, p)
+      return self.__x * z**2 % p
+
+  def y(self):
+      """Return y coordinate in Weierstrass form"""
+      p = self.__curve.p()
+      z = numbertheory.inverse_mod(self.__z, p)
+      return self.__y * z**3 % p
+
+  def to_weierstrass(self):
+      """Return point in Weierstrass form"""
+      p = self.__curve.p()
+      z = numbertheory.inverse_mod(self.__z, p)
+      return Point(self.__curve, self.__x * z**2 % p,
+                   self.__y * z**3 % p, self.__order)
+
+  @staticmethod
+  def from_weierstrass(point):
+      """Create from a Weierstrass form point"""
+      return PointJacobi(point.curve(), point.x(), point.y(), 1,
+                         point.order())
+
+  def double(self):
+      """Add a point to itself."""
+      if not self.__y:
+          return PointJacobi(self.__curve, 0, 0, 0)
+
+      p = self.__curve.p()
+      a = self.__curve.a()
+      x, y, z = self.__x, self.__y, self.__z
+      ysq = y**2 % p
+      s = 4 * x * ysq % p
+      m = (3 * x**2 + a * pow(z, 4, p)) % p
+      nx = (m**2 - 2 * s) % p
+      ny = (m * (s - nx) - 8 * ysq**2) % p
+      nz = 2 * y * z % p
+      return PointJacobi(self.__curve, nx, ny, nz, self.__order)
+
+  def __add__(self, other):
+      """Add two points on elliptic curve."""
+      if self == INFINITY:
+          return other
+      if other == INFINITY:
+          return self
+      if isinstance(other, Point):
+          other = PointJacobi.from_weierstrass(other)
+      if self.__curve != other.__curve:
+          raise ValueError("The other point is on different curve")
+
+      p = self.__curve.p()
+      s_x, s_y, s_z = self.__x, self.__y, self.__z
+      o_x, o_y, o_z = other.__x, other.__y, other.__z
+      o_z2 = o_z**2 % p
+      s_z2 = s_z**2 % p
+      u1 = s_x * o_z2 % p
+      u2 = o_x * s_z2 % p
+      s1 = s_y * o_z2 * o_z % p
+      s2 = o_y * s_z2 * s_z % p
+
+      if u1 == u2:
+          if s1 != s2:
+              return PointJacobi(self.__curve, 0, 0, 1)
+          return self.double()
+
+      h = u2 - u1
+      r = s2 - s1
+      h2 = h**2 % p
+      h3 = h * h2 % p
+      u1h2 = u1 * h2 % p
+      nx = (r ** 2 - h3 - 2 * u1h2) % p
+      ny = (r * (u1h2 - nx) - s1 * h3) % p
+      nz = h * s_z * o_z % p
+      return PointJacobi(self.__curve, nx, ny, nz)
+
+  def __rmul__(self, other):
+      """Multiply point by an integer."""
+      return self * other
+
+  def __mul__(self, other):
+      """Multiply point by an integer."""
+      if self.__y == 0 or other == 0:
+          return PointJacobi(self.__curve, 0, 0, 1)
+      if other == 1:
+          return self
+      if self.__order:
+          other = other % self.__order
+
+      if other % 2:
+          return (self * (other // 2)).double() + self
+      return (self * (other // 2)).double()
+
+
 class Point(object):
   """A point on an elliptic curve. Altering x and y is forbidding,
      but they can be read by the x() and y() methods."""
@@ -79,6 +198,8 @@ def __init__(self, curve, x, y, order=None):
 
   def __eq__(self, other):
     """Return True if the points are identical, False otherwise."""
+    if isinstance(other, PointJacobi):
+        return self == other.to_weierstrass()
     if self.__curve == other.__curve \
        and self.__x == other.__x \
        and self.__y == other.__y:
@@ -98,6 +219,8 @@ def __add__(self, other):
       return self
     if self == INFINITY:
       return other
+    if not isinstance(other, Point):
+        return NotImplemented
     assert self.__curve == other.__curve
     if self.__x == other.__x:
       if (self.__y + other.__y) % self.__curve.p() == 0:
diff --git a/src/ecdsa/test_jacobi.py b/src/ecdsa/test_jacobi.py
new file mode 100644
index 00000000..a54ccbc4
--- /dev/null
+++ b/src/ecdsa/test_jacobi.py
@@ -0,0 +1,59 @@
+
+try:
+    import unittest2 as unittest
+except ImportError:
+    import unittest
+
+from .ellipticcurve import Point, PointJacobi
+from .ecdsa import generator_256, curve_256
+
+class TestJacobi(unittest.TestCase):
+    def test_conversion(self):
+        pj = PointJacobi.from_weierstrass(generator_256)
+        pw = pj.to_weierstrass()
+
+        self.assertEqual(generator_256, pw)
+
+    def test_single_double(self):
+        pj = PointJacobi.from_weierstrass(generator_256)
+        pw = generator_256.double()
+
+        pj = pj.double()
+
+        self.assertEqual(pj.x(), pw.x())
+        self.assertEqual(pj.y(), pj.y())
+
+    def test_multiply_by_one(self):
+        pj = PointJacobi.from_weierstrass(generator_256)
+        pw = generator_256 * 1
+
+        pj = pj * 1
+
+        self.assertEqual(pj.x(), pw.x())
+        self.assertEqual(pj.y(), pj.y())
+
+    def test_multiply_by_two(self):
+        pj = PointJacobi.from_weierstrass(generator_256)
+        pw = generator_256 * 2
+
+        pj = pj * 2
+
+        self.assertEqual(pj.x(), pw.x())
+        self.assertEqual(pj.y(), pj.y())
+
+    def test_compare_double_with_multiply(self):
+        pj = PointJacobi.from_weierstrass(generator_256)
+        dbl = pj.double()
+        mlpl = pj * 2
+
+        self.assertEqual(dbl, mlpl)
+
+    def test_simple_multiplications(self):
+        for i in range(1, 256):
+            pj = PointJacobi.from_weierstrass(generator_256)
+            pw = generator_256 * i
+
+            pj = pj * i
+
+            self.assertEqual((pj.x(), pj.y()), (pw.x(), pw.y()))
+

From a336e9bd3e34ab14b7a9b1e2df009bcac48d1148 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Sun, 6 Oct 2019 19:22:58 +0200
Subject: [PATCH 02/30] recursion-less scalar multiplication

---
 src/ecdsa/ellipticcurve.py | 26 ++++++++++++++++++++------
 1 file changed, 20 insertions(+), 6 deletions(-)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index 3d23db6a..8cc1d671 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -173,12 +173,26 @@ def __mul__(self, other):
           return PointJacobi(self.__curve, 0, 0, 1)
       if other == 1:
           return self
-      if self.__order:
-          other = other % self.__order
-
-      if other % 2:
-          return (self * (other // 2)).double() + self
-      return (self * (other // 2)).double()
+      # makes vulnerable to Minerva
+      #if self.__order:
+      #    other = other % self.__order
+
+      def leftmost_bit(x):
+        assert x > 0
+        result = 1
+        while result <= x:
+          result = 2 * result
+        return result // 2
+
+      e = other
+      i = leftmost_bit(e)
+      result = self
+      while i > 1:
+          result = result.double()
+          if e & (i // 2) != 0:
+              result = result + self
+          i = i // 2
+      return result
 
 
 class Point(object):

From dc6d2d6fa89d6e2f46cfd2552b248f82236c93d7 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Sun, 6 Oct 2019 19:23:21 +0200
Subject: [PATCH 03/30] faster point addition

---
 src/ecdsa/ellipticcurve.py | 104 +++++++++++++++++++++++++++++--------
 1 file changed, 81 insertions(+), 23 deletions(-)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index 8cc1d671..17532358 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -127,6 +127,55 @@ def double(self):
       nz = 2 * y * z % p
       return PointJacobi(self.__curve, nx, ny, nz, self.__order)
 
+  def _add_with_z_1(self, X1, Y1, X2, Y2):
+      """add points when both Z1 and Z2 equal 1"""
+      p = self.__curve.p()
+      H = X2-X1
+      HH = H**2
+      I = 4*HH % p
+      J = H*I
+      r = 2*(Y2-Y1)
+      V = X1*I
+      X3 = (r**2-J-2*V)%p
+      Y3 = (r*(V-X3)-2*Y1*J)%p
+      Z3 = 2*H%p
+      if Y3 == 0 or Z3 == 0:
+          return INFINITY
+      return PointJacobi(self.__curve, X3, Y3, Z3, self.__order)
+
+  def _add_with_z_eq(self, X1, Y1, Z1, X2, Y2):
+      """add points when Z1 == Z2"""
+      p = self.__curve.p()
+      A = (X2-X1)**2%p
+      B = X1*A%p
+      C = X2*A
+      D = (Y2-Y1)**2 %p
+      X3 = (D-B-C)%p
+      Y3 = ((Y2-Y1)*(B-X3)-Y1*(C-B)) %p
+      Z3 = Z1*(X2-X1) %p
+      if Y3 == 0 or Z3 == 0:
+          return INFINITY
+      return PointJacobi(self.__curve, X3, Y3, Z3, self.__order)
+
+  def _add_with_z2_1(self, X1, Y1, Z1, X2, Y2):
+      """add points when Z2 == 1"""
+      p = self.__curve.p()
+      Z1Z1 = Z1**2 % p
+      U2 = X2*Z1Z1
+      S2 = Y2*Z1*Z1Z1
+      H = U2-X1
+      HH = H**2%p
+      I = 4*HH%p
+      J = H*I
+      r = 2*(S2-Y1)
+      V = X1*I
+      X3 = (r**2-J-2*V)%p
+      Y3 = (r*(V-X3)-2*Y1*J)%p
+      Z3 = ((Z1+H)**2-Z1Z1-HH)%p
+      if Y3 == 0 or Z3 == 0:
+          return INFINITY
+      return PointJacobi(self.__curve, X3, Y3, Z3, self.__order)
+
   def __add__(self, other):
       """Add two points on elliptic curve."""
       if self == INFINITY:
@@ -139,29 +188,38 @@ def __add__(self, other):
           raise ValueError("The other point is on different curve")
 
       p = self.__curve.p()
-      s_x, s_y, s_z = self.__x, self.__y, self.__z
-      o_x, o_y, o_z = other.__x, other.__y, other.__z
-      o_z2 = o_z**2 % p
-      s_z2 = s_z**2 % p
-      u1 = s_x * o_z2 % p
-      u2 = o_x * s_z2 % p
-      s1 = s_y * o_z2 * o_z % p
-      s2 = o_y * s_z2 * s_z % p
-
-      if u1 == u2:
-          if s1 != s2:
-              return PointJacobi(self.__curve, 0, 0, 1)
-          return self.double()
-
-      h = u2 - u1
-      r = s2 - s1
-      h2 = h**2 % p
-      h3 = h * h2 % p
-      u1h2 = u1 * h2 % p
-      nx = (r ** 2 - h3 - 2 * u1h2) % p
-      ny = (r * (u1h2 - nx) - s1 * h3) % p
-      nz = h * s_z * o_z % p
-      return PointJacobi(self.__curve, nx, ny, nz)
+      X1, Y1, Z1 = self.__x, self.__y, self.__z
+      X2, Y2, Z2 = other.__x, other.__y, other.__z
+      if Z1 == Z2:
+          if Z1 == 1:
+              return self._add_with_z_1(X1, Y1, X2, Y2)
+          return self._add_with_z_eq(X1, Y1, Z1, X2, Y2)
+      if Z1 == 1 or Z2 == 1:
+          if Z1 == 1:
+              return self._add_with_z2_1(X2, Y2, Z2, X1, Y1)
+          return self._add_with_z2_1(X1, Y1, Z1, X2, Y2)
+
+      # after:
+      # http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-add-2007-bl
+      Z1Z1 = Z1**2 % p
+      Z2Z2 = Z2**2 % p
+      U1 = X1*Z2Z2 % p
+      U2 = X2*Z1Z1 % p
+      S1 = Y1*Z2*Z2Z2 %p
+      S2 = Y2*Z1*Z1Z1 % p
+      H = U2-U1
+      I = (2*H)**2 % p
+      J = H*I % p
+      r = 2*(S2-S1) % p
+      V = U1*I
+      X3 = (r**2-J-2*V) % p
+      Y3 = (r*(V-X3)-2*S1*J) % p
+      Z3 = ((Z1+Z2)**2-Z1Z1-Z2Z2)*H % p
+
+      if Y3 == 0 or Z3 == 0:
+          return INFINITY
+
+      return PointJacobi(self.__curve, X3, Y3, Z3)
 
   def __rmul__(self, other):
       """Multiply point by an integer."""

From d660eb41be5165f3e2072e3e87f7e1052b543fcf Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Sun, 6 Oct 2019 19:53:42 +0200
Subject: [PATCH 04/30] precompute for generators

---
 src/ecdsa/curves.py        | 13 +++++++------
 src/ecdsa/ellipticcurve.py | 36 +++++++++++++++++++++++++++++++++---
 src/ecdsa/test_jacobi.py   | 11 +++++++++++
 3 files changed, 51 insertions(+), 9 deletions(-)

diff --git a/src/ecdsa/curves.py b/src/ecdsa/curves.py
index b61f4b45..2d6877f1 100644
--- a/src/ecdsa/curves.py
+++ b/src/ecdsa/curves.py
@@ -2,6 +2,7 @@
 
 from . import der, ecdsa
 from .util import orderlen
+from .ellipticcurve import PointJacobi
 
 
 # orderlen was defined in this module previously, so keep it in __all__,
@@ -30,22 +31,22 @@ def __init__(self, name, curve, generator, oid, openssl_name=None):
         self.encoded_oid = der.encode_oid(*oid)
 
 NIST192p = Curve("NIST192p", ecdsa.curve_192,
-                 ecdsa.generator_192,
+                 PointJacobi.from_weierstrass(ecdsa.generator_192, True),
                  (1, 2, 840, 10045, 3, 1, 1), "prime192v1")
 NIST224p = Curve("NIST224p", ecdsa.curve_224,
-                 ecdsa.generator_224,
+                 PointJacobi.from_weierstrass(ecdsa.generator_224, True),
                  (1, 3, 132, 0, 33), "secp224r1")
 NIST256p = Curve("NIST256p", ecdsa.curve_256,
-                 ecdsa.generator_256,
+                 PointJacobi.from_weierstrass(ecdsa.generator_256, True),
                  (1, 2, 840, 10045, 3, 1, 7), "prime256v1")
 NIST384p = Curve("NIST384p", ecdsa.curve_384,
-                 ecdsa.generator_384,
+                 PointJacobi.from_weierstrass(ecdsa.generator_384, True),
                  (1, 3, 132, 0, 34), "secp384r1")
 NIST521p = Curve("NIST521p", ecdsa.curve_521,
-                 ecdsa.generator_521,
+                 PointJacobi.from_weierstrass(ecdsa.generator_521, True),
                  (1, 3, 132, 0, 35), "secp521r1")
 SECP256k1 = Curve("SECP256k1", ecdsa.curve_secp256k1,
-                  ecdsa.generator_secp256k1,
+                  PointJacobi.from_weierstrass(ecdsa.generator_secp256k1, True),
                   (1, 3, 132, 0, 10), "secp256k1")
 
 curves = [NIST192p, NIST224p, NIST256p, NIST384p, NIST521p, SECP256k1]
diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index 17532358..a69f9338 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -65,12 +65,24 @@ def __str__(self):
 
 class PointJacobi(object):
   """Point on an elliptic curve. Uses Jacobi representation"""
-  def __init__(self, curve, x, y, z, order=None):
+  def __init__(self, curve, x, y, z, order=None, generator=False):
       self.__curve = curve
       self.__x = x
       self.__y = y
       self.__z = z
       self.__order = order
+      self.__precompute={}
+      if generator:
+          assert order
+          i = 1
+          doubler = PointJacobi(curve, x, y, z, order)
+          self.__precompute[i] = doubler
+          doubler = doubler.double().normalise()
+
+          while i < order:
+              i *= 2
+              self.__precompute[i] = doubler
+              doubler = doubler.double().normalise()
 
   def __eq__(self, other):
       """Compare two points with each-other."""
@@ -98,6 +110,10 @@ def y(self):
       z = numbertheory.inverse_mod(self.__z, p)
       return self.__y * z**3 % p
 
+  def normalise(self):
+      """Return point with z == 1."""
+      return self.from_weierstrass(self.to_weierstrass())
+
   def to_weierstrass(self):
       """Return point in Weierstrass form"""
       p = self.__curve.p()
@@ -106,10 +122,10 @@ def to_weierstrass(self):
                    self.__y * z**3 % p, self.__order)
 
   @staticmethod
-  def from_weierstrass(point):
+  def from_weierstrass(point, generator=False):
       """Create from a Weierstrass form point"""
       return PointJacobi(point.curve(), point.x(), point.y(), 1,
-                         point.order())
+                         point.order(), generator)
 
   def double(self):
       """Add a point to itself."""
@@ -225,12 +241,26 @@ def __rmul__(self, other):
       """Multiply point by an integer."""
       return self * other
 
+  def _mul_precompute(self, other):
+      i = 1
+      result = None
+      while i <= other:
+          if i & other:
+              if result:
+                  result = result + self.__precompute[i]
+              else:
+                  result = self.__precompute[i]
+          i *= 2
+      return result
+
   def __mul__(self, other):
       """Multiply point by an integer."""
       if self.__y == 0 or other == 0:
           return PointJacobi(self.__curve, 0, 0, 1)
       if other == 1:
           return self
+      if self.__precompute:
+          return self._mul_precompute(other)
       # makes vulnerable to Minerva
       #if self.__order:
       #    other = other % self.__order
diff --git a/src/ecdsa/test_jacobi.py b/src/ecdsa/test_jacobi.py
index a54ccbc4..94a483c6 100644
--- a/src/ecdsa/test_jacobi.py
+++ b/src/ecdsa/test_jacobi.py
@@ -56,4 +56,15 @@ def test_simple_multiplications(self):
             pj = pj * i
 
             self.assertEqual((pj.x(), pj.y()), (pw.x(), pw.y()))
+            self.assertEqual(pj, PointJacobi.from_weierstrass(pw))
 
+    def test_precompute(self):
+        precomp = PointJacobi.from_weierstrass(generator_256, True)
+        pj = PointJacobi.from_weierstrass(generator_256)
+
+
+        for i in range(255):
+            print(i)
+            a = precomp * (2 << i)
+            b = pj * (2 << i)
+            self.assertEqual(precomp * (2 << i), pj * (2 << i))

From f86b0b23d8cdf89a24e8d8b4b66859729f8a2adc Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Sun, 6 Oct 2019 20:26:52 +0200
Subject: [PATCH 05/30] use sliding window for precomputation

---
 src/ecdsa/ellipticcurve.py | 43 +++++++++++++++++++++++---------------
 1 file changed, 26 insertions(+), 17 deletions(-)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index a69f9338..9b3658b2 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -36,6 +36,20 @@
 
 from six import python_2_unicode_compatible
 from . import numbertheory
+import binascii
+
+
+def orderlen(order):
+    return (1+len("%x" % order))//2  # bytes
+
+
+def number_to_string(num, order):
+    l = orderlen(order)
+    fmt_str = "%0" + str(2 * l) + "x"
+    string = binascii.unhexlify((fmt_str % num).encode())
+    assert len(string) == l, (len(string), l)
+    return string
+
 
 @python_2_unicode_compatible
 class CurveFp(object):
@@ -73,16 +87,9 @@ def __init__(self, curve, x, y, z, order=None, generator=False):
       self.__order = order
       self.__precompute={}
       if generator:
-          assert order
-          i = 1
-          doubler = PointJacobi(curve, x, y, z, order)
-          self.__precompute[i] = doubler
-          doubler = doubler.double().normalise()
-
-          while i < order:
-              i *= 2
-              self.__precompute[i] = doubler
-              doubler = doubler.double().normalise()
+          point = PointJacobi(curve, x, y, z, order)
+          for i in range(1, 256):
+              self.__precompute[i] = (point * i).normalise()
 
   def __eq__(self, other):
       """Compare two points with each-other."""
@@ -244,13 +251,15 @@ def __rmul__(self, other):
   def _mul_precompute(self, other):
       i = 1
       result = None
-      while i <= other:
-          if i & other:
-              if result:
-                  result = result + self.__precompute[i]
-              else:
-                  result = self.__precompute[i]
-          i *= 2
+      bit_mask = number_to_string(other, other)
+      bit_mask = bytearray(bit_mask)
+      result = self.__precompute[bit_mask[0]]
+      for b in bit_mask[1:]:
+          for _ in range(8):
+              result = result.double()
+          if b:
+              result = result + self.__precompute[b]
+
       return result
 
   def __mul__(self, other):

From 991f1d1e0433e17c6ba5c025320619aaf71083a2 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Mon, 7 Oct 2019 00:09:31 +0200
Subject: [PATCH 06/30] more efficient doubling algorithm

---
 src/ecdsa/ellipticcurve.py | 33 +++++++++++++++++++++++----------
 1 file changed, 23 insertions(+), 10 deletions(-)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index 9b3658b2..60747810 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -119,6 +119,8 @@ def y(self):
 
   def normalise(self):
       """Return point with z == 1."""
+      if self.__z == 1:
+          return self
       return self.from_weierstrass(self.to_weierstrass())
 
   def to_weierstrass(self):
@@ -141,14 +143,24 @@ def double(self):
 
       p = self.__curve.p()
       a = self.__curve.a()
-      x, y, z = self.__x, self.__y, self.__z
-      ysq = y**2 % p
-      s = 4 * x * ysq % p
-      m = (3 * x**2 + a * pow(z, 4, p)) % p
-      nx = (m**2 - 2 * s) % p
-      ny = (m * (s - nx) - 8 * ysq**2) % p
-      nz = 2 * y * z % p
-      return PointJacobi(self.__curve, nx, ny, nz, self.__order)
+
+      X1, Y1, Z1 = self.__x, self.__y, self.__z
+
+      # after:
+      # http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#doubling-dbl-2007-bl
+
+      XX = X1**2%p
+      YY = Y1**2%p
+      YYYY = YY**2%p
+      ZZ = Z1**2%p
+      S = 2*((X1+YY)**2-XX-YYYY)%p
+      M = (3*XX+a*ZZ**2)%p
+      T = (M**2-2*S)%p
+      X3 = T
+      Y3 = (M*(S-T)-8*YYYY)%p
+      Z3 = ((Y1+Z1)**2-YY-ZZ)%p
+
+      return PointJacobi(self.__curve, X3, Y3, Z3, self.__order)
 
   def _add_with_z_1(self, X1, Y1, X2, Y2):
       """add points when both Z1 and Z2 equal 1"""
@@ -283,12 +295,13 @@ def leftmost_bit(x):
 
       e = other
       i = leftmost_bit(e)
+      self = self.normalise()
       result = self
       while i > 1:
           result = result.double()
-          if e & (i // 2) != 0:
-              result = result + self
           i = i // 2
+          if e & i != 0:
+              result = result + self
       return result
 
 

From c00bd47392a77495a61a93cb978da1dfce69e5c7 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Mon, 7 Oct 2019 00:09:42 +0200
Subject: [PATCH 07/30] references for algorithms

---
 src/ecdsa/ellipticcurve.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index 60747810..8ff70381 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -164,6 +164,8 @@ def double(self):
 
   def _add_with_z_1(self, X1, Y1, X2, Y2):
       """add points when both Z1 and Z2 equal 1"""
+      # after:
+      # http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-mmadd-2007-bl
       p = self.__curve.p()
       H = X2-X1
       HH = H**2
@@ -180,6 +182,8 @@ def _add_with_z_1(self, X1, Y1, X2, Y2):
 
   def _add_with_z_eq(self, X1, Y1, Z1, X2, Y2):
       """add points when Z1 == Z2"""
+      # after:
+      # http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-zadd-2007-m
       p = self.__curve.p()
       A = (X2-X1)**2%p
       B = X1*A%p
@@ -194,6 +198,8 @@ def _add_with_z_eq(self, X1, Y1, Z1, X2, Y2):
 
   def _add_with_z2_1(self, X1, Y1, Z1, X2, Y2):
       """add points when Z2 == 1"""
+      # after:
+      # http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-madd-2007-bl
       p = self.__curve.p()
       Z1Z1 = Z1**2 % p
       U2 = X2*Z1Z1

From 77ef6002ba3b3d0c5ab8ac504c7962000f273732 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Mon, 7 Oct 2019 00:12:03 +0200
Subject: [PATCH 08/30] faster inverse modulo

turns out the divmod() call is very expensive
---
 src/ecdsa/numbertheory.py | 30 +++++++++---------------------
 1 file changed, 9 insertions(+), 21 deletions(-)

diff --git a/src/ecdsa/numbertheory.py b/src/ecdsa/numbertheory.py
index 27a9b0f1..1a34c8b2 100644
--- a/src/ecdsa/numbertheory.py
+++ b/src/ecdsa/numbertheory.py
@@ -207,27 +207,15 @@ def square_root_mod_prime(a, p):
 
 
 def inverse_mod(a, m):
-  """Inverse of a mod m."""
-
-  if a < 0 or m <= a:
-    a = a % m
-
-  # From Ferguson and Schneier, roughly:
-
-  c, d = a, m
-  uc, vc, ud, vd = 1, 0, 0, 1
-  while c != 0:
-    q, c, d = divmod(d, c) + (c,)
-    uc, vc, ud, vd = ud - q * uc, vd - q * vc, uc, vc
-
-  # At this point, d is the GCD, and ud*a+vd*m = d.
-  # If d == 1, this means that ud is a inverse.
-
-  assert d == 1
-  if ud > 0:
-    return ud
-  else:
-    return ud + m
+    """Inverse of a mod m."""
+    if a == 0:
+        return 0
+    lm, hm = 1, 0
+    low, high = a % m, m
+    while low > 1:
+        r = high // low
+        lm, low, hm, high = hm - lm * r, high-low *r, lm, low
+    return lm % m
 
 
 def gcd2(a, b):

From 73a5f9f78a0a670ac0edd10aaef3444a9fdc7f4f Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Mon, 7 Oct 2019 00:20:03 +0200
Subject: [PATCH 09/30] ensure that jacobi implementation is used for user
 methods

---
 src/ecdsa/keys.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/ecdsa/keys.py b/src/ecdsa/keys.py
index 4f926429..c380952b 100644
--- a/src/ecdsa/keys.py
+++ b/src/ecdsa/keys.py
@@ -146,6 +146,8 @@ def from_public_point(cls, point, curve=NIST192p, hashfunc=sha1):
         :rtype: VerifyingKey
         """
         self = cls(_error__please_use_generate=True)
+        if not isinstance(point, ellipticcurve.PointJacobi):
+            point = ellipticcurve.PointJacobi.from_weierstrass(point)
         self.curve = curve
         self.default_hashfunc = hashfunc
         self.pubkey = ecdsa.Public_key(curve.generator, point)
@@ -682,6 +684,7 @@ def from_secret_exponent(cls, secexp, curve=NIST192p, hashfunc=sha1):
                 "Invalid value for secexp, expected integer between 1 and {0}"
                 .format(n))
         pubkey_point = curve.generator * secexp
+        pubkey_point = pubkey_point.normalise()
         pubkey = ecdsa.Public_key(curve.generator, pubkey_point)
         pubkey.order = n
         self.verifying_key = VerifyingKey.from_public_point(pubkey_point, curve,

From bd6d745e5142bb55849442e9385f79735933b410 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Mon, 7 Oct 2019 00:21:13 +0200
Subject: [PATCH 10/30] use jacobi implementation also for pubkey recovery

---
 src/ecdsa/ecdsa.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/ecdsa/ecdsa.py b/src/ecdsa/ecdsa.py
index 213c67d8..90493f48 100644
--- a/src/ecdsa/ecdsa.py
+++ b/src/ecdsa/ecdsa.py
@@ -88,12 +88,12 @@ def recover_public_keys(self, hash, generator):
     y = beta if beta % 2 == 0 else curve.p() - beta
 
     # Compute the public key
-    R1 = ellipticcurve.Point(curve, x, y, n)
+    R1 = ellipticcurve.PointJacobi(curve, x, y, 1, n)
     Q1 = numbertheory.inverse_mod(r, n) * (s * R1 + (-e % n) * generator)
     Pk1 = Public_key(generator, Q1)
 
     # And the second solution
-    R2 = ellipticcurve.Point(curve, x, -y, n)
+    R2 = ellipticcurve.PointJacobi(curve, x, -y, 1, n)
     Q2 = numbertheory.inverse_mod(r, n) * (s * R2 + (-e % n) * generator)
     Pk2 = Public_key(generator, Q2)
 

From 1c3b33ff02b2af85dbe4a5b581baa1b8c6cd6374 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Mon, 7 Oct 2019 00:22:40 +0200
Subject: [PATCH 11/30] we're so fast we need wider output :)

---
 speed.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/speed.py b/speed.py
index 3f05309b..686ebf52 100644
--- a/speed.py
+++ b/speed.py
@@ -31,5 +31,7 @@ def do(setup_statements, statement):
     import ecdsa
     c = getattr(ecdsa, curve)
     sig = ecdsa.SigningKey.generate(c).sign(six.b("msg"))
-    print("%9s: siglen=%3d, keygen=%.3fs, sign=%.3fs, verify=%.3fs" \
-          % (curve, len(sig), keygen, sign, verf))
+    print("%9s: siglen=%3d, keygen=%.4fs, %5.1f/s, "
+          "sign=%.4fs, %5.1f/s, verify=%.4fs, %5.1f/s" \
+          % (curve, len(sig), keygen, 1.0/keygen, sign, 1.0/sign,
+             verf, 1.0/verf))

From c4452e4a6925bcb1e149d4cc0770dd036a24bca5 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Mon, 7 Oct 2019 01:02:09 +0200
Subject: [PATCH 12/30] more documentation

---
 src/ecdsa/ellipticcurve.py | 34 ++++++++++++++++++++++++++++++----
 1 file changed, 30 insertions(+), 4 deletions(-)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index 8ff70381..0bc54114 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -1,4 +1,5 @@
 #! /usr/bin/env python
+# -*- coding: utf-8 -*-
 #
 # Implementation of elliptic curves, for cryptographic applications.
 #
@@ -78,8 +79,21 @@ def __str__(self):
 
 
 class PointJacobi(object):
-  """Point on an elliptic curve. Uses Jacobi representation"""
+  """
+  Point on an elliptic curve. Uses Jacobi representation
+
+  In Jacobian coordinates, there are three parameters, X, Y and Z.
+  They correspond to Weierstrass representation points 'x' and 'y' like so:
+
+  x = X / Z²
+  y = Y / Z³
+  """
   def __init__(self, curve, x, y, z, order=None, generator=False):
+      """
+      :param bool generator: the point provided is a curve generator, as
+        such, it will be commonly used with scalar multiplication. This will
+        cause to precompute multiplication table for it
+      """
       self.__curve = curve
       self.__x = x
       self.__y = y
@@ -106,19 +120,31 @@ def curve(self):
       return self.__curve
 
   def x(self):
-      """Return x coordinate in Weierstrass form"""
+      """
+      Return x coordinate in Weierstrass form.
+
+      This method should be used only when the 'y' coordinate is not needed.
+      It's computationally more efficient to use `to_weierstrass()` and then
+      call x() and y() on the returned instance.
+      """
       p = self.__curve.p()
       z = numbertheory.inverse_mod(self.__z, p)
       return self.__x * z**2 % p
 
   def y(self):
-      """Return y coordinate in Weierstrass form"""
+      """
+      Return y coordinate in Weierstrass form.
+
+      This method should be used only when the 'x' coordinate is not needed.
+      It's computationally more efficient to use `to_weierstrass()` and then
+      call x() and y() on the returned instance.
+      """
       p = self.__curve.p()
       z = numbertheory.inverse_mod(self.__z, p)
       return self.__y * z**3 % p
 
   def normalise(self):
-      """Return point with z == 1."""
+      """Return point scaled so that z == 1."""
       if self.__z == 1:
           return self
       return self.from_weierstrass(self.to_weierstrass())

From 6d5b2671b1f3e8b13b9ac4631c3793e0900a1a26 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Wed, 9 Oct 2019 22:55:24 +0200
Subject: [PATCH 13/30] unify INFINITY

---
 src/ecdsa/ellipticcurve.py | 21 ++++++++++++---------
 1 file changed, 12 insertions(+), 9 deletions(-)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index 0bc54114..b7653750 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -107,7 +107,7 @@ def __init__(self, curve, x, y, z, order=None, generator=False):
 
   def __eq__(self, other):
       """Compare two points with each-other."""
-      if self.__y == 0 and other == INFINITY:
+      if (not self.__y or not self.__z) and other is INFINITY:
           return True
       if isinstance(other, Point):
           return self.to_weierstrass() == other
@@ -165,7 +165,7 @@ def from_weierstrass(point, generator=False):
   def double(self):
       """Add a point to itself."""
       if not self.__y:
-          return PointJacobi(self.__curve, 0, 0, 0)
+          return INFINITY
 
       p = self.__curve.p()
       a = self.__curve.a()
@@ -186,6 +186,9 @@ def double(self):
       Y3 = (M*(S-T)-8*YYYY)%p
       Z3 = ((Y1+Z1)**2-YY-ZZ)%p
 
+      if not Y3:
+          return INFINITY
+
       return PointJacobi(self.__curve, X3, Y3, Z3, self.__order)
 
   def _add_with_z_1(self, X1, Y1, X2, Y2):
@@ -202,7 +205,7 @@ def _add_with_z_1(self, X1, Y1, X2, Y2):
       X3 = (r**2-J-2*V)%p
       Y3 = (r*(V-X3)-2*Y1*J)%p
       Z3 = 2*H%p
-      if Y3 == 0 or Z3 == 0:
+      if not Y3:
           return INFINITY
       return PointJacobi(self.__curve, X3, Y3, Z3, self.__order)
 
@@ -218,7 +221,7 @@ def _add_with_z_eq(self, X1, Y1, Z1, X2, Y2):
       X3 = (D-B-C)%p
       Y3 = ((Y2-Y1)*(B-X3)-Y1*(C-B)) %p
       Z3 = Z1*(X2-X1) %p
-      if Y3 == 0 or Z3 == 0:
+      if not Y3:
           return INFINITY
       return PointJacobi(self.__curve, X3, Y3, Z3, self.__order)
 
@@ -239,15 +242,15 @@ def _add_with_z2_1(self, X1, Y1, Z1, X2, Y2):
       X3 = (r**2-J-2*V)%p
       Y3 = (r*(V-X3)-2*Y1*J)%p
       Z3 = ((Z1+H)**2-Z1Z1-HH)%p
-      if Y3 == 0 or Z3 == 0:
+      if not Y3:
           return INFINITY
       return PointJacobi(self.__curve, X3, Y3, Z3, self.__order)
 
   def __add__(self, other):
       """Add two points on elliptic curve."""
-      if self == INFINITY:
+      if self.__y == 0:
           return other
-      if other == INFINITY:
+      if other is INFINITY:
           return self
       if isinstance(other, Point):
           other = PointJacobi.from_weierstrass(other)
@@ -283,7 +286,7 @@ def __add__(self, other):
       Y3 = (r*(V-X3)-2*S1*J) % p
       Z3 = ((Z1+Z2)**2-Z1Z1-Z2Z2)*H % p
 
-      if Y3 == 0 or Z3 == 0:
+      if not Y3:
           return INFINITY
 
       return PointJacobi(self.__curve, X3, Y3, Z3)
@@ -309,7 +312,7 @@ def _mul_precompute(self, other):
   def __mul__(self, other):
       """Multiply point by an integer."""
       if self.__y == 0 or other == 0:
-          return PointJacobi(self.__curve, 0, 0, 1)
+          return INFINITY
       if other == 1:
           return self
       if self.__precompute:

From 36b14bc957ab2be84d4031dbc100c6dc02a58f1d Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Wed, 9 Oct 2019 22:59:17 +0200
Subject: [PATCH 14/30] speed up scaling, name it in standard way

in curve arithmetic it's called scaling, not normalisation

also change one cubing into multiplication
---
 src/ecdsa/ellipticcurve.py | 29 ++++++++++++++++++++---------
 src/ecdsa/keys.py          |  2 +-
 2 files changed, 21 insertions(+), 10 deletions(-)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index b7653750..4d9d4f20 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -103,7 +103,7 @@ def __init__(self, curve, x, y, z, order=None, generator=False):
       if generator:
           point = PointJacobi(curve, x, y, z, order)
           for i in range(1, 256):
-              self.__precompute[i] = (point * i).normalise()
+              self.__precompute[i] = (point * i).scale()
 
   def __eq__(self, other):
       """Compare two points with each-other."""
@@ -143,18 +143,29 @@ def y(self):
       z = numbertheory.inverse_mod(self.__z, p)
       return self.__y * z**3 % p
 
-  def normalise(self):
-      """Return point scaled so that z == 1."""
+  def scale(self):
+      """
+      Return point scaled so that z == 1.
+
+      Modifies point in place, returns self.
+      """
       if self.__z == 1:
           return self
-      return self.from_weierstrass(self.to_weierstrass())
+      p = self.__curve.p()
+      z_inv = numbertheory.inverse_mod(self.__z, p)
+      zz_inv = z_inv ** 2 % p
+      self.__x = self.__x * zz_inv % p
+      self.__y = self.__y * zz_inv * z_inv % p
+      self.__z = 1
+      return self
 
   def to_weierstrass(self):
       """Return point in Weierstrass form"""
-      p = self.__curve.p()
-      z = numbertheory.inverse_mod(self.__z, p)
-      return Point(self.__curve, self.__x * z**2 % p,
-                   self.__y * z**3 % p, self.__order)
+      if not self.__y or not self.__z:
+          return INFINITY
+      self.scale()
+      return Point(self.__curve, self.__x,
+                   self.__y, self.__order)
 
   @staticmethod
   def from_weierstrass(point, generator=False):
@@ -330,7 +341,7 @@ def leftmost_bit(x):
 
       e = other
       i = leftmost_bit(e)
-      self = self.normalise()
+      self = self.scale()
       result = self
       while i > 1:
           result = result.double()
diff --git a/src/ecdsa/keys.py b/src/ecdsa/keys.py
index c380952b..d0c24d2e 100644
--- a/src/ecdsa/keys.py
+++ b/src/ecdsa/keys.py
@@ -684,7 +684,7 @@ def from_secret_exponent(cls, secexp, curve=NIST192p, hashfunc=sha1):
                 "Invalid value for secexp, expected integer between 1 and {0}"
                 .format(n))
         pubkey_point = curve.generator * secexp
-        pubkey_point = pubkey_point.normalise()
+        pubkey_point = pubkey_point.scale()
         pubkey = ecdsa.Public_key(curve.generator, pubkey_point)
         pubkey.order = n
         self.verifying_key = VerifyingKey.from_public_point(pubkey_point, curve,

From 38bb4dabc8ff4bb98e370ad380062694079d562d Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Thu, 10 Oct 2019 00:03:20 +0200
Subject: [PATCH 15/30] implement equality without inverse_mod

since inverse_mod is very computationally expensive (around 100
multiplications) it's cheaper to just bring the fractions to the
same denominator
---
 src/ecdsa/ellipticcurve.py | 17 +++++++++++++++--
 1 file changed, 15 insertions(+), 2 deletions(-)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index 4d9d4f20..e1ad8e0b 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -110,8 +110,21 @@ def __eq__(self, other):
       if (not self.__y or not self.__z) and other is INFINITY:
           return True
       if isinstance(other, Point):
-          return self.to_weierstrass() == other
-      return self.x() == other.x() and self.y() == other.y()
+          x2, y2, z2 = other.x(), other.y(), 1
+      elif isinstance(other, PointJacobi):
+          x2, y2, z2 = other.__x, other.__y, other.__z
+      else:
+          return NotImplemented
+      x1, y1, z1 = self.__x, self.__y, self.__z
+      p = self.__curve.p()
+
+      zz1 = z1 ** 2 % p
+      zz2 = z2 ** 2 % p
+
+      # compare the fractions by bringing them to same denominator
+      # depend on short-circuit to save 4 multiplications in case of inequality
+      return (x1 * zz2 - x2 * zz1) % p == 0 and \
+              (y1 * zz2 * z2 - y2 * zz1 * z1) % p == 0
 
   def order(self):
       return self.__order

From 672bdbccae3281ab807906b928f4d91d86f7679c Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Thu, 10 Oct 2019 00:04:21 +0200
Subject: [PATCH 16/30] don't re-scale whole point when just one coordinate is
 needed

---
 src/ecdsa/ellipticcurve.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index e1ad8e0b..6ca3fef8 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -140,6 +140,8 @@ def x(self):
       It's computationally more efficient to use `to_weierstrass()` and then
       call x() and y() on the returned instance.
       """
+      if self.__z == 1:
+          return self.__x
       p = self.__curve.p()
       z = numbertheory.inverse_mod(self.__z, p)
       return self.__x * z**2 % p
@@ -152,6 +154,8 @@ def y(self):
       It's computationally more efficient to use `to_weierstrass()` and then
       call x() and y() on the returned instance.
       """
+      if self.__z == 1:
+          return self.__y
       p = self.__curve.p()
       z = numbertheory.inverse_mod(self.__z, p)
       return self.__y * z**3 % p

From cdecb334613f00c8110df063dab2e768c679095b Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Thu, 10 Oct 2019 00:04:37 +0200
Subject: [PATCH 17/30] microoptimisation for double()

---
 src/ecdsa/ellipticcurve.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index 6ca3fef8..a66bb285 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -210,14 +210,14 @@ def double(self):
       S = 2*((X1+YY)**2-XX-YYYY)%p
       M = (3*XX+a*ZZ**2)%p
       T = (M**2-2*S)%p
-      X3 = T
+      #X3 = T
       Y3 = (M*(S-T)-8*YYYY)%p
       Z3 = ((Y1+Z1)**2-YY-ZZ)%p
 
       if not Y3:
           return INFINITY
 
-      return PointJacobi(self.__curve, X3, Y3, Z3, self.__order)
+      return PointJacobi(self.__curve, T, Y3, Z3, self.__order)
 
   def _add_with_z_1(self, X1, Y1, X2, Y2):
       """add points when both Z1 and Z2 equal 1"""

From 68342024c7b92c41be137b9b0ae84d067a07f9a8 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Thu, 10 Oct 2019 00:05:47 +0200
Subject: [PATCH 18/30] don't create Public_key object twice

since Public_key verifies correctness of point, it's expensive to
create, do it just once by reusing one created for VerifyingKey
---
 src/ecdsa/keys.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/ecdsa/keys.py b/src/ecdsa/keys.py
index d0c24d2e..3d6e5f0e 100644
--- a/src/ecdsa/keys.py
+++ b/src/ecdsa/keys.py
@@ -685,10 +685,9 @@ def from_secret_exponent(cls, secexp, curve=NIST192p, hashfunc=sha1):
                 .format(n))
         pubkey_point = curve.generator * secexp
         pubkey_point = pubkey_point.scale()
-        pubkey = ecdsa.Public_key(curve.generator, pubkey_point)
-        pubkey.order = n
         self.verifying_key = VerifyingKey.from_public_point(pubkey_point, curve,
                                                             hashfunc)
+        pubkey = self.verifying_key.pubkey
         self.privkey = ecdsa.Private_key(pubkey, secexp)
         self.privkey.order = n
         return self

From 7f4e03057f395ce25a0369a266d0a3f99c0624d0 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Thu, 10 Oct 2019 09:02:01 +0200
Subject: [PATCH 19/30] Revert "use sliding window for precomputation"

This reverts commit 92c7422d851e0b654a31bdc0fd3dc415a3b8f45d.
---
 src/ecdsa/ellipticcurve.py | 43 +++++++++++++++-----------------------
 1 file changed, 17 insertions(+), 26 deletions(-)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index a66bb285..3ea87b01 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -37,20 +37,6 @@
 
 from six import python_2_unicode_compatible
 from . import numbertheory
-import binascii
-
-
-def orderlen(order):
-    return (1+len("%x" % order))//2  # bytes
-
-
-def number_to_string(num, order):
-    l = orderlen(order)
-    fmt_str = "%0" + str(2 * l) + "x"
-    string = binascii.unhexlify((fmt_str % num).encode())
-    assert len(string) == l, (len(string), l)
-    return string
-
 
 @python_2_unicode_compatible
 class CurveFp(object):
@@ -101,9 +87,16 @@ def __init__(self, curve, x, y, z, order=None, generator=False):
       self.__order = order
       self.__precompute={}
       if generator:
-          point = PointJacobi(curve, x, y, z, order)
-          for i in range(1, 256):
-              self.__precompute[i] = (point * i).scale()
+          assert order
+          i = 1
+          doubler = PointJacobi(curve, x, y, z, order)
+          self.__precompute[i] = doubler
+          doubler = doubler.double().scale()
+
+          while i < order:
+              i *= 2
+              self.__precompute[i] = doubler
+              doubler = doubler.double().scale()
 
   def __eq__(self, other):
       """Compare two points with each-other."""
@@ -326,15 +319,13 @@ def __rmul__(self, other):
   def _mul_precompute(self, other):
       i = 1
       result = None
-      bit_mask = number_to_string(other, other)
-      bit_mask = bytearray(bit_mask)
-      result = self.__precompute[bit_mask[0]]
-      for b in bit_mask[1:]:
-          for _ in range(8):
-              result = result.double()
-          if b:
-              result = result + self.__precompute[b]
-
+      while i <= other:
+          if i & other:
+              if result:
+                  result = result + self.__precompute[i]
+              else:
+                  result = self.__precompute[i]
+          i *= 2
       return result
 
   def __mul__(self, other):

From 1f6110fd621b58da2229353e681041f0e6351887 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Thu, 10 Oct 2019 09:32:57 +0200
Subject: [PATCH 20/30] add ability to precompute for verification

---
 speed.py          | 5 +++--
 src/ecdsa/keys.py | 4 ++++
 2 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/speed.py b/speed.py
index 686ebf52..a74fcd7f 100644
--- a/speed.py
+++ b/speed.py
@@ -20,14 +20,15 @@ def do(setup_statements, statement):
     S3 = "msg = six.b('msg')"
     S4 = "sig = sk.sign(msg)"
     S5 = "vk = sk.get_verifying_key()"
-    S6 = "vk.verify(sig, msg)"
+    S6 = "vk.precompute()"
+    S7 = "vk.verify(sig, msg)"
     # We happen to know that .generate() also calculates the
     # verifying key, which is the time-consuming part. If the code
     # were changed to lazily calculate vk, we'd need to change this
     # benchmark to loop over S5 instead of S2
     keygen = do([S1], S2)
     sign = do([S1,S2,S3], S4)
-    verf = do([S1,S2,S3,S4,S5], S6)
+    verf = do([S1,S2,S3,S4,S5,S6], S7)
     import ecdsa
     c = getattr(ecdsa, curve)
     sig = ecdsa.SigningKey.generate(c).sign(six.b("msg"))
diff --git a/src/ecdsa/keys.py b/src/ecdsa/keys.py
index 3d6e5f0e..65d9f819 100644
--- a/src/ecdsa/keys.py
+++ b/src/ecdsa/keys.py
@@ -154,6 +154,10 @@ def from_public_point(cls, point, curve=NIST192p, hashfunc=sha1):
         self.pubkey.order = curve.order
         return self
 
+    def precompute(self):
+        self.pubkey.point = ellipticcurve.PointJacobi.from_weierstrass(
+            self.pubkey.point, True)
+
     @staticmethod
     def _from_raw_encoding(string, curve, validate_point):
         """

From 5c65e2380e1c7a40cd64994597b5ba4726ed29a3 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Fri, 11 Oct 2019 02:12:19 +0200
Subject: [PATCH 21/30] verify in one place, skip for self-created points

when loading public keys, perform the point verification just once
when loading private keys, do not verify the derived public point
---
 src/ecdsa/ecdsa.py | 30 ++++++++++++++++++++++--------
 src/ecdsa/keys.py  | 44 +++++++++++++++++++++++++++-----------------
 2 files changed, 49 insertions(+), 25 deletions(-)

diff --git a/src/ecdsa/ecdsa.py b/src/ecdsa/ecdsa.py
index 90493f48..3cfc5d18 100644
--- a/src/ecdsa/ecdsa.py
+++ b/src/ecdsa/ecdsa.py
@@ -63,6 +63,10 @@ class RSZeroError(RuntimeError):
     pass
 
 
+class InvalidPointError(RuntimeError):
+    pass
+
+
 class Signature(object):
   """ECDSA signature.
   """
@@ -99,25 +103,35 @@ def recover_public_keys(self, hash, generator):
 
     return [Pk1, Pk2]
 
+
 class Public_key(object):
   """Public key for ECDSA.
   """
 
-  def __init__(self, generator, point):
-    """generator is the Point that generates the group,
-    point is the Point that defines the public key.
+  def __init__(self, generator, point, verify=True):
+    """
+    Low level ECDSA public key object.
+
+    :param generator: the Point that generates the group
+    :param point: the Point that defines the public key
+    :param bool verify: if True check if point is valid point on curve
+
+    :raises InvalidPointError: if the point parameters are invalid or
+        point does not lie on the curve
     """
 
     self.curve = generator.curve()
     self.generator = generator
     self.point = point
     n = generator.order()
-    if not n:
-      raise RuntimeError("Generator point must have order.")
-    if not n * point == ellipticcurve.INFINITY:
-      raise RuntimeError("Generator point order is bad.")
     if point.x() < 0 or n <= point.x() or point.y() < 0 or n <= point.y():
-      raise RuntimeError("Generator point has x or y out of range.")
+      raise InvalidPointError("Generator point has x or y out of range.")
+    if verify and not self.curve.contains_point(point.x(), point.y()):
+        raise InvalidPointError("Point does not lie on the curve")
+    if not n:
+      raise InvalidPointError("Generator point must have order.")
+    if verify and not n * point == ellipticcurve.INFINITY:
+      raise InvalidPointError("Generator point order is bad.")
 
   def verifies(self, hash, signature):
     """Verify that signature is a valid signature of hash.
diff --git a/src/ecdsa/keys.py b/src/ecdsa/keys.py
index 65d9f819..132b6654 100644
--- a/src/ecdsa/keys.py
+++ b/src/ecdsa/keys.py
@@ -126,7 +126,8 @@ def __init__(self, _error__please_use_generate=None):
         self.pubkey = None
 
     @classmethod
-    def from_public_point(cls, point, curve=NIST192p, hashfunc=sha1):
+    def from_public_point(cls, point, curve=NIST192p, hashfunc=sha1,
+                          validate_point=True):
         """
         Initialise the object from a Point object.
 
@@ -141,6 +142,12 @@ def from_public_point(cls, point, curve=NIST192p, hashfunc=sha1):
             verification, needs to implement the same interface
             as hashlib.sha1
         :type hashfunc: callable
+        :type bool validate_point: whether to check if the point lies on curve
+            should always be used if the public point is not a result
+            of our own calculation
+
+        :raises MalformedPointError: if the public point does not lie on the
+            curve
 
         :return: Initialised VerifyingKey object
         :rtype: VerifyingKey
@@ -150,7 +157,11 @@ def from_public_point(cls, point, curve=NIST192p, hashfunc=sha1):
             point = ellipticcurve.PointJacobi.from_weierstrass(point)
         self.curve = curve
         self.default_hashfunc = hashfunc
-        self.pubkey = ecdsa.Public_key(curve.generator, point)
+        try:
+            self.pubkey = ecdsa.Public_key(curve.generator, point,
+                                           validate_point)
+        except ecdsa.InvalidPointError:
+            raise MalformedPointError("Point does not lie on the curve")
         self.pubkey.order = curve.order
         return self
 
@@ -159,7 +170,7 @@ def precompute(self):
             self.pubkey.point, True)
 
     @staticmethod
-    def _from_raw_encoding(string, curve, validate_point):
+    def _from_raw_encoding(string, curve):
         """
         Decode public point from :term:`raw encoding`.
 
@@ -177,13 +188,11 @@ def _from_raw_encoding(string, curve, validate_point):
             raise MalformedPointError("Unexpected length of encoded y")
         x = string_to_number(xs)
         y = string_to_number(ys)
-        if validate_point and not ecdsa.point_is_valid(curve.generator, x, y):
-            raise MalformedPointError("Point does not lie on the curve")
 
-        return ellipticcurve.Point(curve.curve, x, y, order)
+        return ellipticcurve.PointJacobi(curve.curve, x, y, 1, order)
 
     @staticmethod
-    def _from_compressed(string, curve, validate_point):
+    def _from_compressed(string, curve):
         """Decode public point from compressed encoding."""
         if string[:1] not in (b('\x02'), b('\x03')):
             raise MalformedPointError("Malformed compressed point encoding")
@@ -202,9 +211,7 @@ def _from_compressed(string, curve, validate_point):
             y = p - beta
         else:
             y = beta
-        if validate_point and not ecdsa.point_is_valid(curve.generator, x, y):
-            raise MalformedPointError("Point does not lie on curve")
-        return ellipticcurve.Point(curve.curve, x, y, order)
+        return ellipticcurve.PointJacobi(curve.curve, x, y, 1, order)
 
     @classmethod
     def _from_hybrid(cls, string, curve, validate_point):
@@ -213,7 +220,7 @@ def _from_hybrid(cls, string, curve, validate_point):
         assert string[:1] in (b('\x06'), b('\x07'))
 
         # primarily use the uncompressed as it's easiest to handle
-        point = cls._from_raw_encoding(string[1:], curve, validate_point)
+        point = cls._from_raw_encoding(string[1:], curve)
 
         # but validate if it's self-consistent if we're asked to do that
         if validate_point and \
@@ -248,29 +255,32 @@ def from_string(cls, string, curve=NIST192p, hashfunc=sha1,
             provided curve or not, defaults to True
         :type validate_point: bool
 
+        :raises MalformedPointError: if the public point does not lie on the
+            curve or the encoding is invalid
+
         :return: Initialised VerifyingKey object
         :rtype: VerifyingKey
         """
         sig_len = len(string)
         if sig_len == curve.verifying_key_length:
-            point = cls._from_raw_encoding(string, curve, validate_point)
+            point = cls._from_raw_encoding(string, curve)
         elif sig_len == curve.verifying_key_length + 1:
             if string[:1] in (b('\x06'), b('\x07')):
                 point = cls._from_hybrid(string, curve, validate_point)
             elif string[:1] == b('\x04'):
-                point = cls._from_raw_encoding(string[1:], curve,
-                                               validate_point)
+                point = cls._from_raw_encoding(string[1:], curve)
             else:
                 raise MalformedPointError(
                     "Invalid X9.62 encoding of the public point")
         elif sig_len == curve.baselen + 1:
-            point = cls._from_compressed(string, curve, validate_point)
+            point = cls._from_compressed(string, curve)
         else:
             raise MalformedPointError(
                 "Length of string does not match lengths of "
                 "any of the supported encodings of {0} "
                 "curve.".format(curve.name))
-        return cls.from_public_point(point, curve, hashfunc)
+        return cls.from_public_point(point, curve, hashfunc,
+                                     validate_point)
 
     @classmethod
     def from_pem(cls, string):
@@ -690,7 +700,7 @@ def from_secret_exponent(cls, secexp, curve=NIST192p, hashfunc=sha1):
         pubkey_point = curve.generator * secexp
         pubkey_point = pubkey_point.scale()
         self.verifying_key = VerifyingKey.from_public_point(pubkey_point, curve,
-                                                            hashfunc)
+                                                            hashfunc, False)
         pubkey = self.verifying_key.pubkey
         self.privkey = ecdsa.Private_key(pubkey, secexp)
         self.privkey.order = n

From ea4787e92c6ae2c7bef7e9555dac76e54cdecb64 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Sat, 12 Oct 2019 12:56:39 +0200
Subject: [PATCH 22/30] don't waste point double in precomputation

---
 src/ecdsa/ellipticcurve.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index 3ea87b01..db25380c 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -91,12 +91,11 @@ def __init__(self, curve, x, y, z, order=None, generator=False):
           i = 1
           doubler = PointJacobi(curve, x, y, z, order)
           self.__precompute[i] = doubler
-          doubler = doubler.double().scale()
 
           while i < order:
               i *= 2
-              self.__precompute[i] = doubler
               doubler = doubler.double().scale()
+              self.__precompute[i] = doubler
 
   def __eq__(self, other):
       """Compare two points with each-other."""

From e14ae61027168e6567b0d79acceb48f801c6f64c Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Sat, 12 Oct 2019 12:57:56 +0200
Subject: [PATCH 23/30] handle addition of the same point

---
 src/ecdsa/ellipticcurve.py | 10 ++++++++++
 src/ecdsa/test_jacobi.py   |  3 +++
 2 files changed, 13 insertions(+)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index db25380c..9550fa45 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -221,6 +221,8 @@ def _add_with_z_1(self, X1, Y1, X2, Y2):
       I = 4*HH % p
       J = H*I
       r = 2*(Y2-Y1)
+      if not H and not r:
+          return self.double()
       V = X1*I
       X3 = (r**2-J-2*V)%p
       Y3 = (r*(V-X3)-2*Y1*J)%p
@@ -238,6 +240,8 @@ def _add_with_z_eq(self, X1, Y1, Z1, X2, Y2):
       B = X1*A%p
       C = X2*A
       D = (Y2-Y1)**2 %p
+      if not A and not D:
+          return self.double()
       X3 = (D-B-C)%p
       Y3 = ((Y2-Y1)*(B-X3)-Y1*(C-B)) %p
       Z3 = Z1*(X2-X1) %p
@@ -258,6 +262,8 @@ def _add_with_z2_1(self, X1, Y1, Z1, X2, Y2):
       I = 4*HH%p
       J = H*I
       r = 2*(S2-Y1)
+      if not r and not H:
+          return self.double()
       V = X1*I
       X3 = (r**2-J-2*V)%p
       Y3 = (r*(V-X3)-2*Y1*J)%p
@@ -280,6 +286,8 @@ def __add__(self, other):
       p = self.__curve.p()
       X1, Y1, Z1 = self.__x, self.__y, self.__z
       X2, Y2, Z2 = other.__x, other.__y, other.__z
+      if (X1, Y1, Z1) == (X2, Y2, Z2):
+          return self.double()
       if Z1 == Z2:
           if Z1 == 1:
               return self._add_with_z_1(X1, Y1, X2, Y2)
@@ -301,6 +309,8 @@ def __add__(self, other):
       I = (2*H)**2 % p
       J = H*I % p
       r = 2*(S2-S1) % p
+      if not H and not r:
+          return self.double()
       V = U1*I
       X3 = (r**2-J-2*V) % p
       Y3 = (r*(V-X3)-2*S1*J) % p
diff --git a/src/ecdsa/test_jacobi.py b/src/ecdsa/test_jacobi.py
index 94a483c6..e6357e86 100644
--- a/src/ecdsa/test_jacobi.py
+++ b/src/ecdsa/test_jacobi.py
@@ -68,3 +68,6 @@ def test_precompute(self):
             a = precomp * (2 << i)
             b = pj * (2 << i)
             self.assertEqual(precomp * (2 << i), pj * (2 << i))
+
+    def test_add_point(self):
+        j_g = PointJacobi.from_weierstrass(generator_256)

From 4c17c22a3b2f001564e4d44d40e4de2f23f11de8 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Sat, 12 Oct 2019 12:59:17 +0200
Subject: [PATCH 24/30] remove if in precomputed multiplication

(the ifs are already in the addition methods)
---
 src/ecdsa/ellipticcurve.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index 9550fa45..f5cbfc30 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -101,6 +101,8 @@ def __eq__(self, other):
       """Compare two points with each-other."""
       if (not self.__y or not self.__z) and other is INFINITY:
           return True
+      if self.__y and self.__z and other is INFINITY:
+          return False
       if isinstance(other, Point):
           x2, y2, z2 = other.x(), other.y(), 1
       elif isinstance(other, PointJacobi):
@@ -327,13 +329,10 @@ def __rmul__(self, other):
 
   def _mul_precompute(self, other):
       i = 1
-      result = None
+      result = INFINITY
       while i <= other:
           if i & other:
-              if result:
-                  result = result + self.__precompute[i]
-              else:
-                  result = self.__precompute[i]
+              result = result + self.__precompute[i]
           i *= 2
       return result
 

From ec6312deb113f81778d60b584a1d4e4b49e9b988 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Sat, 12 Oct 2019 12:59:52 +0200
Subject: [PATCH 25/30] just use PointJacobi when comparing with regular Point

---
 src/ecdsa/ellipticcurve.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index f5cbfc30..5abffd0d 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -384,8 +384,8 @@ def __init__(self, curve, x, y, order=None):
 
   def __eq__(self, other):
     """Return True if the points are identical, False otherwise."""
-    if isinstance(other, PointJacobi):
-        return self == other.to_weierstrass()
+    if not isinstance(other, Point):
+        return NotImplemented
     if self.__curve == other.__curve \
        and self.__x == other.__x \
        and self.__y == other.__y:

From 62616b7694284992d89c4b2c43e68cf6af51f9ae Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Sat, 12 Oct 2019 13:27:24 +0200
Subject: [PATCH 26/30] add method that does multiplication and addition of two
 points

---
 src/ecdsa/ellipticcurve.py | 49 ++++++++++++++++++++++++++++++
 src/ecdsa/test_jacobi.py   | 61 +++++++++++++++++++++++++++++++++++++-
 2 files changed, 109 insertions(+), 1 deletion(-)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index 5abffd0d..8df7b35c 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -366,6 +366,55 @@ def leftmost_bit(x):
               result = result + self
       return result
 
+  def _mul_add_precompute(self, self_mul, other, other_mul):
+      i = 1
+      result = INFINITY
+      m = max(self_mul, other_mul)
+      while i <= m:
+          if i & self_mul:
+              result = result + self.__precompute[i]
+          if i & other_mul:
+              result = result + other.__precompute[i]
+          i *= 2
+      return result
+
+  @staticmethod
+  def _leftmost_bit(x):
+    assert x > 0
+    result = 1
+    while result <= x:
+      result = 2 * result
+    return result // 2
+
+  def mul_add(self, self_mul, other, other_mul):
+      """
+      Do two multiplications at the same time, add results.
+
+      calculates self*self_mul + other*other_mul
+      """
+      if other is INFINITY or other_mul == 0:
+          return self * self_mul
+      if self_mul == 0:
+          return other * other_mul
+      if self.__precompute and other.__precompute:
+          return self._mul_add_precompute(self_mul, other, other_mul)
+
+      i = self._leftmost_bit(max(self_mul, other_mul))*2
+      result = INFINITY
+      self = self.scale()
+      other = other.scale()
+      both = (self + other).scale()
+      while i > 1:
+          result = result.double()
+          i = i // 2
+          if self_mul & i and other_mul & i:
+              result = result + both
+          elif self_mul & i:
+              result = result + self
+          elif other_mul & i:
+              result = result + other
+      return result
+
 
 class Point(object):
   """A point on an elliptic curve. Altering x and y is forbidding,
diff --git a/src/ecdsa/test_jacobi.py b/src/ecdsa/test_jacobi.py
index e6357e86..b2290967 100644
--- a/src/ecdsa/test_jacobi.py
+++ b/src/ecdsa/test_jacobi.py
@@ -4,7 +4,7 @@
 except ImportError:
     import unittest
 
-from .ellipticcurve import Point, PointJacobi
+from .ellipticcurve import Point, PointJacobi, INFINITY
 from .ecdsa import generator_256, curve_256
 
 class TestJacobi(unittest.TestCase):
@@ -71,3 +71,62 @@ def test_precompute(self):
 
     def test_add_point(self):
         j_g = PointJacobi.from_weierstrass(generator_256)
+
+        self.assertEqual(j_g.double(), j_g + j_g)
+
+    def test_add_point_3_times(self):
+        j_g = PointJacobi.from_weierstrass(generator_256)
+
+        self.assertEqual(j_g * 3, j_g + j_g + j_g)
+
+    def test_mul_add_inf(self):
+        j_g = PointJacobi.from_weierstrass(generator_256)
+
+        self.assertEqual(j_g, j_g.mul_add(1, INFINITY, 1))
+
+    def test_mul_add_same(self):
+        j_g = PointJacobi.from_weierstrass(generator_256)
+
+        self.assertEqual(j_g * 2, j_g.mul_add(1, j_g, 1))
+
+    def test_mul_add_precompute(self):
+        j_g = PointJacobi.from_weierstrass(generator_256, True)
+
+        b = PointJacobi.from_weierstrass(j_g * 255, True)
+
+        self.assertEqual(j_g * 256, j_g + b)
+
+        self.assertEqual(j_g * (5 + 255 * 7), j_g * 5 + b * 7)
+        self.assertEqual(j_g * (5 + 255 * 7), j_g.mul_add(5, b, 7))
+
+    def test_mul_add_precompute_large(self):
+        j_g = PointJacobi.from_weierstrass(generator_256, True)
+
+        b = PointJacobi.from_weierstrass(j_g * 255, True)
+
+        self.assertEqual(j_g * 256, j_g + b)
+
+        self.assertEqual(j_g * (0xff00 + 255 * 0xf0f0),
+                         j_g * 0xff00 + b * 0xf0f0)
+        self.assertEqual(j_g * (0xff00 + 255 * 0xf0f0),
+                         j_g.mul_add(0xff00, b, 0xf0f0))
+
+
+    def test_mul_add_to_mul(self):
+        j_g = PointJacobi.from_weierstrass(generator_256)
+
+        a = j_g * 3
+        b = j_g.mul_add(2, j_g, 1)
+        self.assertEqual(a, b)
+
+    def test_mul_add(self):
+        j_g = PointJacobi.from_weierstrass(generator_256)
+        w_a = generator_256 * 255
+        w_b = generator_256 * (0xa8*0xf0)
+
+        j_b = j_g * 0xa8
+
+        ret = j_g.mul_add(255, j_b, 0xf0)
+
+        self.assertEqual(ret.to_weierstrass(), w_a + w_b)
+

From e6807e3124977352175e748f4a623f8187ddba9b Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Sat, 12 Oct 2019 13:27:41 +0200
Subject: [PATCH 27/30] use combined multiplication and addition in
 verification

---
 src/ecdsa/ecdsa.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/ecdsa/ecdsa.py b/src/ecdsa/ecdsa.py
index 3cfc5d18..2add4366 100644
--- a/src/ecdsa/ecdsa.py
+++ b/src/ecdsa/ecdsa.py
@@ -151,7 +151,10 @@ def verifies(self, hash, signature):
     c = numbertheory.inverse_mod(s, n)
     u1 = (hash * c) % n
     u2 = (r * c) % n
-    xy = u1 * G + u2 * self.point
+    if hasattr(G, "mul_add"):
+        xy = G.mul_add(u1, self.point, u2)
+    else:
+        xy = u1 * G + u2 * self.point
     v = xy.x() % n
     return v == r
 

From 12ec5b00846d719656aa94d7a542606748ef7992 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Sat, 12 Oct 2019 13:28:04 +0200
Subject: [PATCH 28/30] more space for benchmark data

---
 speed.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/speed.py b/speed.py
index a74fcd7f..bcc8a085 100644
--- a/speed.py
+++ b/speed.py
@@ -32,7 +32,7 @@ def do(setup_statements, statement):
     import ecdsa
     c = getattr(ecdsa, curve)
     sig = ecdsa.SigningKey.generate(c).sign(six.b("msg"))
-    print("%9s: siglen=%3d, keygen=%.4fs, %5.1f/s, "
-          "sign=%.4fs, %5.1f/s, verify=%.4fs, %5.1f/s" \
+    print("%9s: siglen=%3d, keygen=%.4fs, %6.1f/s, "
+          "sign=%.4fs, %6.1f/s, verify=%.4fs, %6.1f/s" \
           % (curve, len(sig), keygen, 1.0/keygen, sign, 1.0/sign,
              verf, 1.0/verf))

From 3fbbd15ef6433371e82bb69000d5f98a446219a9 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Sat, 12 Oct 2019 13:55:12 +0200
Subject: [PATCH 29/30] use 2-ary NAF for multiplication with precomputation

---
 src/ecdsa/ellipticcurve.py | 33 ++++++++++++++++-----------------
 1 file changed, 16 insertions(+), 17 deletions(-)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index 8df7b35c..37dff372 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -90,6 +90,7 @@ def __init__(self, curve, x, y, z, order=None, generator=False):
           assert order
           i = 1
           doubler = PointJacobi(curve, x, y, z, order)
+          order *= 2
           self.__precompute[i] = doubler
 
           while i < order:
@@ -330,9 +331,15 @@ def __rmul__(self, other):
   def _mul_precompute(self, other):
       i = 1
       result = INFINITY
-      while i <= other:
-          if i & other:
-              result = result + self.__precompute[i]
+      while other:
+          if other % 2:
+              if other % 4 >= 2:
+                  result = result + (-self.__precompute[i])
+                  other -= -1
+              else:
+                  result = result + self.__precompute[i]
+                  other -= 1
+          other = other // 2
           i *= 2
       return result
 
@@ -366,18 +373,6 @@ def leftmost_bit(x):
               result = result + self
       return result
 
-  def _mul_add_precompute(self, self_mul, other, other_mul):
-      i = 1
-      result = INFINITY
-      m = max(self_mul, other_mul)
-      while i <= m:
-          if i & self_mul:
-              result = result + self.__precompute[i]
-          if i & other_mul:
-              result = result + other.__precompute[i]
-          i *= 2
-      return result
-
   @staticmethod
   def _leftmost_bit(x):
     assert x > 0
@@ -396,8 +391,8 @@ def mul_add(self, self_mul, other, other_mul):
           return self * self_mul
       if self_mul == 0:
           return other * other_mul
-      if self.__precompute and other.__precompute:
-          return self._mul_add_precompute(self_mul, other, other_mul)
+      if self.__precompute or other.__precompute:
+          return self * self_mul + other * other_mul
 
       i = self._leftmost_bit(max(self_mul, other_mul))*2
       result = INFINITY
@@ -415,6 +410,10 @@ def mul_add(self, self_mul, other, other_mul):
               result = result + other
       return result
 
+  def __neg__(self):
+      return PointJacobi(self.__curve, self.__x, -self.__y, self.__z,
+                         self.__order)
+
 
 class Point(object):
   """A point on an elliptic curve. Altering x and y is forbidding,

From 775535695394c66567f2591477d94f1badc38156 Mon Sep 17 00:00:00 2001
From: Hubert Kario <hubert@kario.pl>
Date: Sat, 12 Oct 2019 14:28:53 +0200
Subject: [PATCH 30/30] use 2-ary NAF for regular multiplication too

---
 src/ecdsa/ellipticcurve.py | 34 +++++++++++++++++++++-------------
 1 file changed, 21 insertions(+), 13 deletions(-)

diff --git a/src/ecdsa/ellipticcurve.py b/src/ecdsa/ellipticcurve.py
index 37dff372..00e2457b 100644
--- a/src/ecdsa/ellipticcurve.py
+++ b/src/ecdsa/ellipticcurve.py
@@ -343,6 +343,20 @@ def _mul_precompute(self, other):
           i *= 2
       return result
 
+  def _naf(self, mult):
+      ret = []
+      while mult:
+          if mult % 2:
+              nd = mult % 4
+              if nd >= 2:
+                  nd = nd - 4
+              ret += [nd]
+              mult -= nd
+          else:
+              ret += [0]
+          mult //= 2
+      return ret
+
   def __mul__(self, other):
       """Multiply point by an integer."""
       if self.__y == 0 or other == 0:
@@ -355,21 +369,15 @@ def __mul__(self, other):
       #if self.__order:
       #    other = other % self.__order
 
-      def leftmost_bit(x):
-        assert x > 0
-        result = 1
-        while result <= x:
-          result = 2 * result
-        return result // 2
-
-      e = other
-      i = leftmost_bit(e)
       self = self.scale()
-      result = self
-      while i > 1:
+      result = INFINITY
+      # since adding points when at least one of them is scaled
+      # is quicker, reverse the NAF order
+      for i in reversed(self._naf(other)):
           result = result.double()
-          i = i // 2
-          if e & i != 0:
+          if i < 0:
+              result = result + (-self)
+          elif i > 0:
               result = result + self
       return result