1#! /usr/bin/env python
2# -*- coding: utf-8 -*-
3#
4# Implementation of elliptic curves, for cryptographic applications.
5#
6# This module doesn't provide any way to choose a random elliptic
7# curve, nor to verify that an elliptic curve was chosen randomly,
8# because one can simply use NIST's standard curves.
9#
10# Notes from X9.62-1998 (draft):
11#   Nomenclature:
12#     - Q is a public key.
13#     The "Elliptic Curve Domain Parameters" include:
14#     - q is the "field size", which in our case equals p.
15#     - p is a big prime.
16#     - G is a point of prime order (5.1.1.1).
17#     - n is the order of G (5.1.1.1).
18#   Public-key validation (5.2.2):
19#     - Verify that Q is not the point at infinity.
20#     - Verify that X_Q and Y_Q are in [0,p-1].
21#     - Verify that Q is on the curve.
22#     - Verify that nQ is the point at infinity.
23#   Signature generation (5.3):
24#     - Pick random k from [1,n-1].
25#   Signature checking (5.4.2):
26#     - Verify that r and s are in [1,n-1].
27#
28# Version of 2008.11.25.
29#
30# Revision history:
31#    2005.12.31 - Initial version.
32#    2008.11.25 - Change CurveFp.is_on to contains_point.
33#
34# Written in 2005 by Peter Pearson and placed in the public domain.
35
36from __future__ import division
37
38try:
39    from gmpy2 import mpz
40    GMPY = True
41except ImportError:
42    try:
43        from gmpy import mpz
44        GMPY = True
45    except ImportError:
46        GMPY = False
47
48
49from six import python_2_unicode_compatible
50from . import numbertheory
51from ._rwlock import RWLock
52
53
54@python_2_unicode_compatible
55class CurveFp(object):
56  """Elliptic Curve over the field of integers modulo a prime."""
57
58  if GMPY:
59      def __init__(self, p, a, b, h=None):
60        """
61        The curve of points satisfying y^2 = x^3 + a*x + b (mod p).
62
63        h is an integer that is the cofactor of the elliptic curve domain
64        parameters; it is the number of points satisfying the elliptic curve
65        equation divided by the order of the base point. It is used for selection
66        of efficient algorithm for public point verification.
67        """
68        self.__p = mpz(p)
69        self.__a = mpz(a)
70        self.__b = mpz(b)
71        # h is not used in calculations and it can be None, so don't use
72        # gmpy with it
73        self.__h = h
74  else:
75      def __init__(self, p, a, b, h=None):
76        """
77        The curve of points satisfying y^2 = x^3 + a*x + b (mod p).
78
79        h is an integer that is the cofactor of the elliptic curve domain
80        parameters; it is the number of points satisfying the elliptic curve
81        equation divided by the order of the base point. It is used for selection
82        of efficient algorithm for public point verification.
83        """
84        self.__p = p
85        self.__a = a
86        self.__b = b
87        self.__h = h
88
89  def __eq__(self, other):
90    if isinstance(other, CurveFp):
91      """Return True if the curves are identical, False otherwise."""
92      return self.__p == other.__p \
93          and self.__a == other.__a \
94          and self.__b == other.__b
95    return NotImplemented
96
97  def __hash__(self):
98    return hash((self.__p, self.__a, self.__b))
99
100  def p(self):
101    return self.__p
102
103  def a(self):
104    return self.__a
105
106  def b(self):
107    return self.__b
108
109  def cofactor(self):
110    return self.__h
111
112  def contains_point(self, x, y):
113    """Is the point (x,y) on this curve?"""
114    return (y * y - ((x * x + self.__a) * x + self.__b)) % self.__p == 0
115
116  def __str__(self):
117    return "CurveFp(p=%d, a=%d, b=%d, h=%d)" % (
118        self.__p, self.__a, self.__b, self.__h)
119
120
121class PointJacobi(object):
122  """
123  Point on an elliptic curve. Uses Jacobi coordinates.
124
125  In Jacobian coordinates, there are three parameters, X, Y and Z.
126  They correspond to affine parameters 'x' and 'y' like so:
127
128  x = X / Z²
129  y = Y / Z³
130  """
131  def __init__(self, curve, x, y, z, order=None, generator=False):
132      """
133      Initialise a point that uses Jacobi representation internally.
134
135      :param CurveFp curve: curve on which the point resides
136      :param int x: the X parameter of Jacobi representation (equal to x when
137        converting from affine coordinates
138      :param int y: the Y parameter of Jacobi representation (equal to y when
139        converting from affine coordinates
140      :param int z: the Z parameter of Jacobi representation (equal to 1 when
141        converting from affine coordinates
142      :param int order: the point order, must be non zero when using
143        generator=True
144      :param bool generator: the point provided is a curve generator, as
145        such, it will be commonly used with scalar multiplication. This will
146        cause to precompute multiplication table for it
147      """
148      self.__curve = curve
149      # since it's generally better (faster) to use scaled points vs unscaled
150      # ones, use writer-biased RWLock for locking:
151      self._scale_lock = RWLock()
152      if GMPY:
153          self.__x = mpz(x)
154          self.__y = mpz(y)
155          self.__z = mpz(z)
156          self.__order = order and mpz(order)
157      else:
158          self.__x = x
159          self.__y = y
160          self.__z = z
161          self.__order = order
162      self.__precompute = []
163      if generator:
164          assert order
165          i = 1
166          order *= 2
167          doubler = PointJacobi(curve, x, y, z, order)
168          order *= 2
169          self.__precompute.append((doubler.x(), doubler.y()))
170
171          while i < order:
172              i *= 2
173              doubler = doubler.double().scale()
174              self.__precompute.append((doubler.x(), doubler.y()))
175
176  def __eq__(self, other):
177      """Compare two points with each-other."""
178      try:
179          self._scale_lock.reader_acquire()
180          if other is INFINITY:
181              return not self.__y or not self.__z
182          x1, y1, z1 = self.__x, self.__y, self.__z
183      finally:
184          self._scale_lock.reader_release()
185      if isinstance(other, Point):
186          x2, y2, z2 = other.x(), other.y(), 1
187      elif isinstance(other, PointJacobi):
188          try:
189              other._scale_lock.reader_acquire()
190              x2, y2, z2 = other.__x, other.__y, other.__z
191          finally:
192              other._scale_lock.reader_release()
193      else:
194          return NotImplemented
195      if self.__curve != other.curve():
196          return False
197      p = self.__curve.p()
198
199      zz1 = z1 * z1 % p
200      zz2 = z2 * z2 % p
201
202      # compare the fractions by bringing them to the same denominator
203      # depend on short-circuit to save 4 multiplications in case of inequality
204      return (x1 * zz2 - x2 * zz1) % p == 0 and \
205             (y1 * zz2 * z2 - y2 * zz1 * z1) % p == 0
206
207  def order(self):
208      """Return the order of the point.
209
210      None if it is undefined.
211      """
212      return self.__order
213
214  def curve(self):
215      """Return curve over which the point is defined."""
216      return self.__curve
217
218  def x(self):
219      """
220      Return affine x coordinate.
221
222      This method should be used only when the 'y' coordinate is not needed.
223      It's computationally more efficient to use `to_affine()` and then
224      call x() and y() on the returned instance. Or call `scale()`
225      and then x() and y() on the returned instance.
226      """
227      try:
228          self._scale_lock.reader_acquire()
229          if self.__z == 1:
230              return self.__x
231          x = self.__x
232          z = self.__z
233      finally:
234          self._scale_lock.reader_release()
235      p = self.__curve.p()
236      z = numbertheory.inverse_mod(z, p)
237      return x * z**2 % p
238
239  def y(self):
240      """
241      Return affine y coordinate.
242
243      This method should be used only when the 'x' coordinate is not needed.
244      It's computationally more efficient to use `to_affine()` and then
245      call x() and y() on the returned instance. Or call `scale()`
246      and then x() and y() on the returned instance.
247      """
248      try:
249          self._scale_lock.reader_acquire()
250          if self.__z == 1:
251              return self.__y
252          y = self.__y
253          z = self.__z
254      finally:
255          self._scale_lock.reader_release()
256      p = self.__curve.p()
257      z = numbertheory.inverse_mod(z, p)
258      return y * z**3 % p
259
260  def scale(self):
261      """
262      Return point scaled so that z == 1.
263
264      Modifies point in place, returns self.
265      """
266      try:
267          self._scale_lock.reader_acquire()
268          if self.__z == 1:
269              return self
270      finally:
271          self._scale_lock.reader_release()
272
273      try:
274          self._scale_lock.writer_acquire()
275          # scaling already scaled point is safe (as inverse of 1 is 1) and
276          # quick so we don't need to optimise for the unlikely event when
277          # two threads hit the lock at the same time
278          p = self.__curve.p()
279          z_inv = numbertheory.inverse_mod(self.__z, p)
280          zz_inv = z_inv * z_inv % p
281          self.__x = self.__x * zz_inv % p
282          self.__y = self.__y * zz_inv * z_inv % p
283          # we are setting the z last so that the check above will return true
284          # only after all values were already updated
285          self.__z = 1
286      finally:
287          self._scale_lock.writer_release()
288      return self
289
290  def to_affine(self):
291      """Return point in affine form."""
292      if not self.__y or not self.__z:
293          return INFINITY
294      self.scale()
295      # after point is scaled, it's immutable, so no need to perform locking
296      return Point(self.__curve, self.__x,
297                   self.__y, self.__order)
298
299  @staticmethod
300  def from_affine(point, generator=False):
301      """Create from an affine point.
302
303      :param bool generator: set to True to make the point to precalculate
304        multiplication table - useful for public point when verifying many
305        signatures (around 100 or so) or for generator points of a curve.
306      """
307      return PointJacobi(point.curve(), point.x(), point.y(), 1,
308                         point.order(), generator)
309
310  # plese note that all the methods that use the equations from hyperelliptic
311  # are formatted in a way to maximise performance.
312  # Things that make code faster: multiplying instead of taking to the power
313  # (`xx = x * x; xxxx = xx * xx % p` is faster than `xxxx = x**4 % p` and
314  # `pow(x, 4, p)`),
315  # multiple assignments at the same time (`x1, x2 = self.x1, self.x2` is
316  # faster than `x1 = self.x1; x2 = self.x2`),
317  # similarly, sometimes the `% p` is skipped if it makes the calculation
318  # faster and the result of calculation is later reduced modulo `p`
319
320  def _double_with_z_1(self, X1, Y1, p, a):
321      """Add a point to itself with z == 1."""
322      # after:
323      # http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#doubling-mdbl-2007-bl
324      XX, YY = X1 * X1 % p, Y1 * Y1 % p
325      if not YY:
326          return 0, 0, 1
327      YYYY = YY * YY % p
328      S = 2 * ((X1 + YY)**2 - XX - YYYY) % p
329      M = 3 * XX + a
330      T = (M * M - 2 * S) % p
331      # X3 = T
332      Y3 = (M * (S - T) - 8 * YYYY) % p
333      Z3 = 2 * Y1 % p
334      return T, Y3, Z3
335
336  def _double(self, X1, Y1, Z1, p, a):
337      """Add a point to itself, arbitrary z."""
338      if Z1 == 1:
339          return self._double_with_z_1(X1, Y1, p, a)
340      if not Z1:
341          return 0, 0, 1
342      # after:
343      # http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#doubling-dbl-2007-bl
344      XX, YY = X1 * X1 % p, Y1 * Y1 % p
345      if not YY:
346          return 0, 0, 1
347      YYYY = YY * YY % p
348      ZZ = Z1 * Z1 % p
349      S = 2 * ((X1 + YY)**2 - XX - YYYY) % p
350      M = (3 * XX + a * ZZ * ZZ) % p
351      T = (M * M - 2 * S) % p
352      # X3 = T
353      Y3 = (M * (S - T) - 8 * YYYY) % p
354      Z3 = ((Y1 + Z1)**2 - YY - ZZ) % p
355
356      return T, Y3, Z3
357
358  def double(self):
359      """Add a point to itself."""
360      if not self.__y:
361          return INFINITY
362
363      p, a = self.__curve.p(), self.__curve.a()
364
365      try:
366          self._scale_lock.reader_acquire()
367          X1, Y1, Z1 = self.__x, self.__y, self.__z
368      finally:
369          self._scale_lock.reader_release()
370
371      X3, Y3, Z3 = self._double(X1, Y1, Z1, p, a)
372
373      if not Y3 or not Z3:
374          return INFINITY
375      return PointJacobi(self.__curve, X3, Y3, Z3, self.__order)
376
377  def _add_with_z_1(self, X1, Y1, X2, Y2, p):
378      """add points when both Z1 and Z2 equal 1"""
379      # after:
380      # http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-mmadd-2007-bl
381      H = X2 - X1
382      HH = H * H
383      I = 4 * HH % p
384      J = H * I
385      r = 2 * (Y2 - Y1)
386      if not H and not r:
387          return self._double_with_z_1(X1, Y1, p, self.__curve.a())
388      V = X1 * I
389      X3 = (r**2 - J - 2 * V) % p
390      Y3 = (r * (V - X3) - 2 * Y1 * J) % p
391      Z3 = 2 * H % p
392      return X3, Y3, Z3
393
394  def _add_with_z_eq(self, X1, Y1, Z1, X2, Y2, p):
395      """add points when Z1 == Z2"""
396      # after:
397      # http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-zadd-2007-m
398      A = (X2 - X1)**2 % p
399      B = X1 * A % p
400      C = X2 * A
401      D = (Y2 - Y1)**2 % p
402      if not A and not D:
403          return self._double(X1, Y1, Z1, p, self.__curve.a())
404      X3 = (D - B - C) % p
405      Y3 = ((Y2 - Y1) * (B - X3) - Y1 * (C - B)) % p
406      Z3 = Z1 * (X2 - X1) % p
407      return X3, Y3, Z3
408
409  def _add_with_z2_1(self, X1, Y1, Z1, X2, Y2, p):
410      """add points when Z2 == 1"""
411      # after:
412      # http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-madd-2007-bl
413      Z1Z1 = Z1 * Z1 % p
414      U2, S2 = X2 * Z1Z1 % p, Y2 * Z1 * Z1Z1 % p
415      H = (U2 - X1) % p
416      HH = H * H % p
417      I = 4 * HH % p
418      J = H * I
419      r = 2 * (S2 - Y1) % p
420      if not r and not H:
421          return self._double_with_z_1(X2, Y2, p, self.__curve.a())
422      V = X1 * I
423      X3 = (r * r - J - 2 * V) % p
424      Y3 = (r * (V - X3) - 2 * Y1 * J) % p
425      Z3 = ((Z1 + H)**2 - Z1Z1 - HH) % p
426      return X3, Y3, Z3
427
428  def _add_with_z_ne(self, X1, Y1, Z1, X2, Y2, Z2, p):
429      """add points with arbitrary z"""
430      # after:
431      # http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-add-2007-bl
432      Z1Z1 = Z1 * Z1 % p
433      Z2Z2 = Z2 * Z2 % p
434      U1 = X1 * Z2Z2 % p
435      U2 = X2 * Z1Z1 % p
436      S1 = Y1 * Z2 * Z2Z2 % p
437      S2 = Y2 * Z1 * Z1Z1 % p
438      H = U2 - U1
439      I = 4 * H * H % p
440      J = H * I % p
441      r = 2 * (S2 - S1) % p
442      if not H and not r:
443          return self._double(X1, Y1, Z1, p, self.__curve.a())
444      V = U1 * I
445      X3 = (r * r - J - 2 * V) % p
446      Y3 = (r * (V - X3) - 2 * S1 * J) % p
447      Z3 = ((Z1 + Z2)**2 - Z1Z1 - Z2Z2) * H % p
448
449      return X3, Y3, Z3
450
451  def __radd__(self, other):
452      """Add other to self."""
453      return self + other
454
455  def _add(self, X1, Y1, Z1, X2, Y2, Z2, p):
456      """add two points, select fastest method."""
457      if not Y1 or not Z1:
458          return X2, Y2, Z2
459      if not Y2 or not Z2:
460          return X1, Y1, Z1
461      if Z1 == Z2:
462          if Z1 == 1:
463              return self._add_with_z_1(X1, Y1, X2, Y2, p)
464          return self._add_with_z_eq(X1, Y1, Z1, X2, Y2, p)
465      if Z1 == 1:
466          return self._add_with_z2_1(X2, Y2, Z2, X1, Y1, p)
467      if Z2 == 1:
468          return self._add_with_z2_1(X1, Y1, Z1, X2, Y2, p)
469      return self._add_with_z_ne(X1, Y1, Z1, X2, Y2, Z2, p)
470
471  def __add__(self, other):
472      """Add two points on elliptic curve."""
473      if self == INFINITY:
474          return other
475      if other == INFINITY:
476          return self
477      if isinstance(other, Point):
478          other = PointJacobi.from_affine(other)
479      if self.__curve != other.__curve:
480          raise ValueError("The other point is on different curve")
481
482      p = self.__curve.p()
483      try:
484          self._scale_lock.reader_acquire()
485          X1, Y1, Z1 = self.__x, self.__y, self.__z
486      finally:
487          self._scale_lock.reader_release()
488      try:
489          other._scale_lock.reader_acquire()
490          X2, Y2, Z2 = other.__x, other.__y, other.__z
491      finally:
492          other._scale_lock.reader_release()
493      X3, Y3, Z3 = self._add(X1, Y1, Z1, X2, Y2, Z2, p)
494
495      if not Y3 or not Z3:
496          return INFINITY
497      return PointJacobi(self.__curve, X3, Y3, Z3, self.__order)
498
499  def __rmul__(self, other):
500      """Multiply point by an integer."""
501      return self * other
502
503  def _mul_precompute(self, other):
504      """Multiply point by integer with precomputation table."""
505      X3, Y3, Z3, p = 0, 0, 1, self.__curve.p()
506      _add = self._add
507      for X2, Y2 in self.__precompute:
508          if other % 2:
509              if other % 4 >= 2:
510                  other = (other + 1)//2
511                  X3, Y3, Z3 = _add(X3, Y3, Z3, X2, -Y2, 1, p)
512              else:
513                  other = (other - 1)//2
514                  X3, Y3, Z3 = _add(X3, Y3, Z3, X2, Y2, 1, p)
515          else:
516              other //= 2
517
518      if not Y3 or not Z3:
519          return INFINITY
520      return PointJacobi(self.__curve, X3, Y3, Z3, self.__order)
521
522  @staticmethod
523  def _naf(mult):
524      """Calculate non-adjacent form of number."""
525      ret = []
526      while mult:
527          if mult % 2:
528              nd = mult % 4
529              if nd >= 2:
530                  nd = nd - 4
531              ret += [nd]
532              mult -= nd
533          else:
534              ret += [0]
535          mult //= 2
536      return ret
537
538  def __mul__(self, other):
539      """Multiply point by an integer."""
540      if not self.__y or not other:
541          return INFINITY
542      if other == 1:
543          return self
544      if self.__order:
545          # order*2 as a protection for Minerva
546          other = other % (self.__order*2)
547      if self.__precompute:
548          return self._mul_precompute(other)
549
550      self = self.scale()
551      # once scaled, point is immutable, not need to lock
552      X2, Y2 = self.__x, self.__y
553      X3, Y3, Z3 = 0, 0, 1
554      p, a = self.__curve.p(), self.__curve.a()
555      _double = self._double
556      _add = self._add
557      # since adding points when at least one of them is scaled
558      # is quicker, reverse the NAF order
559      for i in reversed(self._naf(other)):
560          X3, Y3, Z3 = _double(X3, Y3, Z3, p, a)
561          if i < 0:
562              X3, Y3, Z3 = _add(X3, Y3, Z3, X2, -Y2, 1, p)
563          elif i > 0:
564              X3, Y3, Z3 = _add(X3, Y3, Z3, X2, Y2, 1, p)
565
566      if not Y3 or not Z3:
567          return INFINITY
568
569      return PointJacobi(self.__curve, X3, Y3, Z3, self.__order)
570
571  @staticmethod
572  def _leftmost_bit(x):
573      """Return integer with the same magnitude as x but hamming weight of 1"""
574      assert x > 0
575      result = 1
576      while result <= x:
577        result = 2 * result
578      return result // 2
579
580  def mul_add(self, self_mul, other, other_mul):
581      """
582      Do two multiplications at the same time, add results.
583
584      calculates self*self_mul + other*other_mul
585      """
586      if other is INFINITY or other_mul == 0:
587          return self * self_mul
588      if self_mul == 0:
589          return other * other_mul
590      if not isinstance(other, PointJacobi):
591          other = PointJacobi.from_affine(other)
592      # when the points have precomputed answers, then multiplying them alone
593      # is faster (as it uses NAF)
594      if self.__precompute and other.__precompute:
595          return self * self_mul + other * other_mul
596
597      if self.__order:
598          self_mul = self_mul % self.__order
599          other_mul = other_mul % self.__order
600
601      i = self._leftmost_bit(max(self_mul, other_mul))*2
602      X3, Y3, Z3 = 0, 0, 1
603      p, a = self.__curve.p(), self.__curve.a()
604      self = self.scale()
605      # after scaling, point is immutable, no need for locking
606      X1, Y1 = self.__x, self.__y
607      other = other.scale()
608      X2, Y2 = other.__x, other.__y
609      both = (self + other).scale()
610      X4, Y4 = both.__x, both.__y
611      _double = self._double
612      _add = self._add
613      while i > 1:
614          X3, Y3, Z3 = _double(X3, Y3, Z3, p, a)
615          i = i // 2
616
617          if self_mul & i and other_mul & i:
618              X3, Y3, Z3 = _add(X3, Y3, Z3, X4, Y4, 1, p)
619          elif self_mul & i:
620              X3, Y3, Z3 = _add(X3, Y3, Z3, X1, Y1, 1, p)
621          elif other_mul & i:
622              X3, Y3, Z3 = _add(X3, Y3, Z3, X2, Y2, 1, p)
623
624      if not Y3 or not Z3:
625          return INFINITY
626
627      return PointJacobi(self.__curve, X3, Y3, Z3, self.__order)
628
629  def __neg__(self):
630      """Return negated point."""
631      try:
632          self._scale_lock.reader_acquire()
633          return PointJacobi(self.__curve, self.__x, -self.__y, self.__z,
634                             self.__order)
635      finally:
636          self._scale_lock.reader_release()
637
638
639class Point(object):
640  """A point on an elliptic curve. Altering x and y is forbidding,
641     but they can be read by the x() and y() methods."""
642  def __init__(self, curve, x, y, order=None):
643    """curve, x, y, order; order (optional) is the order of this point."""
644    self.__curve = curve
645    if GMPY:
646        self.__x = x and mpz(x)
647        self.__y = y and mpz(y)
648        self.__order = order and mpz(order)
649    else:
650        self.__x = x
651        self.__y = y
652        self.__order = order
653    # self.curve is allowed to be None only for INFINITY:
654    if self.__curve:
655      assert self.__curve.contains_point(x, y)
656    # for curves with cofactor 1, all points that are on the curve are scalar
657    # multiples of the base point, so performing multiplication is not
658    # necessary to verify that. See Section 3.2.2.1 of SEC 1 v2
659    if curve and curve.cofactor() != 1 and order:
660      assert self * order == INFINITY
661
662  def __eq__(self, other):
663    """Return True if the points are identical, False otherwise."""
664    if isinstance(other, Point):
665      return self.__curve == other.__curve \
666          and self.__x == other.__x \
667          and self.__y == other.__y
668    return NotImplemented
669
670  def __neg__(self):
671    return Point(self.__curve, self.__x, self.__curve.p() - self.__y)
672
673  def __add__(self, other):
674    """Add one point to another point."""
675
676    # X9.62 B.3:
677
678    if not isinstance(other, Point):
679        return NotImplemented
680    if other == INFINITY:
681      return self
682    if self == INFINITY:
683      return other
684    assert self.__curve == other.__curve
685    if self.__x == other.__x:
686      if (self.__y + other.__y) % self.__curve.p() == 0:
687        return INFINITY
688      else:
689        return self.double()
690
691    p = self.__curve.p()
692
693    l = ((other.__y - self.__y) * \
694         numbertheory.inverse_mod(other.__x - self.__x, p)) % p
695
696    x3 = (l * l - self.__x - other.__x) % p
697    y3 = (l * (self.__x - x3) - self.__y) % p
698
699    return Point(self.__curve, x3, y3)
700
701  def __mul__(self, other):
702    """Multiply a point by an integer."""
703
704    def leftmost_bit(x):
705      assert x > 0
706      result = 1
707      while result <= x:
708        result = 2 * result
709      return result // 2
710
711    e = other
712    if e == 0 or (self.__order and e % self.__order == 0):
713      return INFINITY
714    if self == INFINITY:
715      return INFINITY
716    if e < 0:
717      return (-self) * (-e)
718
719    # From X9.62 D.3.2:
720
721    e3 = 3 * e
722    negative_self = Point(self.__curve, self.__x, -self.__y, self.__order)
723    i = leftmost_bit(e3) // 2
724    result = self
725    # print_("Multiplying %s by %d (e3 = %d):" % (self, other, e3))
726    while i > 1:
727      result = result.double()
728      if (e3 & i) != 0 and (e & i) == 0:
729        result = result + self
730      if (e3 & i) == 0 and (e & i) != 0:
731        result = result + negative_self
732      # print_(". . . i = %d, result = %s" % ( i, result ))
733      i = i // 2
734
735    return result
736
737  def __rmul__(self, other):
738    """Multiply a point by an integer."""
739
740    return self * other
741
742  def __str__(self):
743    if self == INFINITY:
744      return "infinity"
745    return "(%d,%d)" % (self.__x, self.__y)
746
747  def double(self):
748    """Return a new point that is twice the old."""
749
750    if self == INFINITY:
751      return INFINITY
752
753    # X9.62 B.3:
754
755    p = self.__curve.p()
756    a = self.__curve.a()
757
758    l = ((3 * self.__x * self.__x + a) * \
759         numbertheory.inverse_mod(2 * self.__y, p)) % p
760
761    x3 = (l * l - 2 * self.__x) % p
762    y3 = (l * (self.__x - x3) - self.__y) % p
763
764    return Point(self.__curve, x3, y3)
765
766  def x(self):
767    return self.__x
768
769  def y(self):
770    return self.__y
771
772  def curve(self):
773    return self.__curve
774
775  def order(self):
776    return self.__order
777
778
779# This one point is the Point At Infinity for all purposes:
780INFINITY = Point(None, None, None)
781