text: # --- ELLIPTIC CURVE MATH ------------------------------------------------------
  #
  #   curve definition:   y^2 = x^3 - p*x - q
  #   over finite field:  Z/nZ* (prime residue classes modulo a prime number n)
  #
  #
  #   COPYRIGHT (c) 2010 by Toni Mattis <solaris@live.de>
  # ------------------------------------------------------------------------------
  
  """
  Module for elliptic curve arithmetic over a prime field GF(n).
  E(GF(n)) takes the form y**2 == x**3 - p*x - q (mod n) for a prime n.
  
  0. Structures used by this module
  
      PARAMETERS and SCALARS are non-negative (long) integers.
  
      A POINT (x, y), usually denoted p1, p2, ...
      is a pair of (long) integers where 0 <= x < n and 0 <= y < n
  
      A POINT in PROJECTIVE COORDINATES, usually denoted jp1, jp2, ...
      takes the form (X, Y, Z, Z**2, Z**3) where x = X / Z**2
      and y = Y / z**3. This form is called Jacobian coordinates.
  
      The NEUTRAL element "0" or "O" is represented by None
      in both coordinate systems.
  
  1. Basic Functions
  
      euclid()            Is the Extended Euclidean Algorithm.
      inv()               Computes the multiplicative inversion modulo n.
      curve_q()           Finds the curve parameter q (mod n)
                          when p and a point are given.
      element()           Tests whether a point (x, y) is on the curve.
  
  2. Point transformations
  
      to_projective()     Converts a point (x, y) to projective coordinates.
      from_projective()   Converts a point from projective coordinates
                          to (x, y) using the transformation described above.
      neg()               Computes the inverse point -P in both coordinate
                          systems.
  
  3. Slow point arithmetic
  
      These algorithms make use of basic geometry and modular arithmetic
      thus being suitable for small numbers and academic study.
  
      add()               Computes the sum of two (x, y)-points
      mul()               Perform scalar multiplication using "double & add"
  
  4. Fast point arithmetic
  
      These algorithms make use of projective coordinates, signed binary
      expansion and a JSP-like approach (joint sparse form).
  
      The following functions consume and return projective coordinates:
  
      addf()              Optimized point addition.
      doublef()           Optimized point doubling.
      mulf()              Highly optimized scalar multiplication.
      muladdf()           Highly optimized addition of two products.
      
      The following functions use the optimized ones above but consume
      and output (x, y)-coordinates for a more convenient usage:
  
      mulp()              Encapsulates mulf()
      muladdp()           Encapsulates muladdf()
  
      For single additions add() is generally faster than an encapsulation of
      addf() which would involve expensive coordinate transformations.
      Hence there is no addp() and doublep().
  """
  from __future__ import division
  try:
      from builtins import zip
      from builtins import range
  except ImportError:
      pass
  #from past.utils import old_div
  
  # BASIC MATH -------------------------------------------------------------------
  
  
  def euclid(a, b):
      """Solve x*a + y*b = ggt(a, b) and return (x, y, ggt(a, b))"""
      # Non-recursive approach hence suitable for large numbers
      x = yy = 0
      y = xx = 1
      while b:
          q = a // b
          a, b = b, a % b
          x, xx = xx - q * x, x
          y, yy = yy - q * y, y
      return xx, yy, a
  
  
  def inv(a, n):
      """Perform inversion 1/a modulo n. a and n should be COPRIME."""
      # coprimality is not checked here in favour of performance
      i = euclid(a, n)[0]
      while i < 0:
          i += n
      return i
  
  
  def curve_q(x, y, p, n):
      """Find curve parameter q mod n having point (x, y) and parameter p"""
      return ((x * x - p) * x - y * y) % n
  
  
  def element(point, p, q, n):
      """Test, whether the given point is on the curve (p, q, n)"""
      if point:
          x, y = point
          return (x * x * x - p * x - q) % n == (y * y) % n
      else:
          return True
  
  
  def to_projective(p):
      """Transform point p given as (x, y) to projective coordinates"""
      if p:
          return (p[0], p[1], 1, 1, 1)
      else:
          return None     # Identity point (0)
  
  
  def from_projective(jp, n):
      """Transform a point from projective coordinates to (x, y) mod n"""
      if jp:
          return (jp[0] * inv(jp[3], n)) % n, (jp[1] * inv(jp[4], n)) % n
      else:
          return None     # Identity point (0)
  
  
  def neg(p, n):
      """Compute the inverse point to p in any coordinate system"""
      return (p[0], (n - p[1]) % n) + p[2:] if p else None
  
  
  # POINT ADDITION ---------------------------------------------------------------
  
  # addition of points in y**2 = x**3 - p*x - q over <Z/nZ*; +>
  def add(p, q, n, p1, p2):
      """Add points p1 and p2 over curve (p, q, n)"""
      if p1 and p2:
          x1, y1 = p1
          x2, y2 = p2
          if (x1 - x2) % n:
              s = ((y1 - y2) * inv(x1 - x2, n)) % n   # slope
              x = (s * s - x1 - x2) % n               # intersection with curve
              return x, n - (y1 + s * (x - x1)) % n
          else:
              if (y1 + y2) % n:       # slope s calculated by derivation
                  s = ((3 * x1 * x1 - p) * inv(2 * y1, n)) % n
                  x = (s * s - 2 * x1) % n            # intersection with curve
                  return x, n - (y1 + s * (x - x1)) % n
              else:
                  return None
      else:   # either p1 is not none -> ret. p1, otherwiese p2, which may be
          return p1 if p1 else p2     # none too.
  
  
  # faster addition: redundancy in projective coordinates eliminates
  # expensive inversions mod n.
  def addf(p, q, n, jp1, jp2):
      """Add jp1 and jp2 in projective (jacobian) coordinates."""
      if jp1 and jp2:
  
          x1, y1, z1, z1s, z1c = jp1
          x2, y2, z2, z2s, z2c = jp2
  
          s1 = (y1 * z2c) % n
          s2 = (y2 * z1c) % n
  
          u1 = (x1 * z2s) % n
          u2 = (x2 * z1s) % n
  
          if (u1 - u2) % n:
  
              h = (u2 - u1) % n
              r = (s2 - s1) % n
  
              hs = (h * h) % n
              hc = (hs * h) % n
  
              x3 = (-hc - 2 * u1 * hs + r * r) % n
              y3 = (-s1 * hc + r * (u1 * hs - x3)) % n
              z3 = (z1 * z2 * h) % n
  
              z3s = (z3 * z3) % n
              z3c = (z3s * z3) % n
  
              return x3, y3, z3, z3s, z3c
  
          else:
              if (s1 + s2) % n:
                  return doublef(p, q, n, jp1)
              else:
                  return None
      else:
          return jp1 if jp1 else jp2
  
  # explicit point doubling using redundant coordinates
  def doublef(p, q, n, jp):
      """Double jp in projective (jacobian) coordinates"""
      if not jp:
          return None
      x1, y1, z1, z1p2, z1p3 = jp
  
      y1p2 = (y1 * y1) % n
      a = (4 * x1 * y1p2) % n
      b = (3 * x1 * x1 - p * z1p3 * z1) % n
      x3 = (b * b - 2 * a) % n
      y3 = (b * (a - x3) - 8 * y1p2 * y1p2) % n
      z3 = (2 * y1 * z1) % n
      z3p2 = (z3 * z3) % n
  
      return x3, y3, z3, z3p2, (z3p2 * z3) % n
  
  
  # SCALAR MULTIPLICATION --------------------------------------------------------
  
  # scalar multiplication p1 * c = p1 + p1 + ... + p1 (c times) in O(log(n))
  def mul(p, q, n, p1, c):
      """multiply point p1 by scalar c over curve (p, q, n)"""
      res = None
      while c > 0:
          if c & 1:
              res = add(p, q, n, res, p1)
          c >>= 1                     # c = c / 2
          p1 = add(p, q, n, p1, p1)   # p1 = p1 * 2
      return res
  
  
  # this method allows _signed_bin() to choose between 1 and -1. It will select
  # the sign which leaves the higher number of zeroes in the binary
  # representation (the higher GDB).
  def _gbd(n):
      """Compute second greatest base-2 divisor"""
      i = 1
      if n <= 0: return 0
      while not n % i:
          i <<= 1
      return i >> 2
  
  
  # This method transforms n into a binary representation having signed bits.
  # A signed binary expansion contains more zero-bits hence reducing the number
  # of additions required by a multiplication algorithm.
  #
  # Example:  15 ( 0b1111 ) can be written as 16 - 1, resulting in (1,0,0,0,-1)
  #           and saving 2 additions. Subtraction can be performed as
  #           efficiently as addition.
  def _signed_bin(n):
      """Transform n into an optimized signed binary representation"""
      r = []
      while n > 1:
          if n & 1:
              cp = _gbd(n + 1)
              cn = _gbd(n - 1)
              if cp > cn:         # -1 leaves more zeroes -> subtract -1 (= +1)
                  r.append(-1)
                  n += 1
              else:               # +1 leaves more zeroes -> subtract +1 (= -1)
                  r.append(+1)
                  n -= 1
          else:
              r.append(0)         # be glad about one more zero
          n >>= 1
      r.append(n)
      return r[::-1]
  
  
  # This multiplication algorithm combines signed binary expansion and
  # fast addition using projective coordinates resulting in 5 to 10 times
  # faster multiplication.
  def mulf(p, q, n, jp1, c):
      """Multiply point jp1 by c in projective coordinates"""
      sb = _signed_bin(c)
      res = None
      jp0 = neg(jp1, n)  # additive inverse of jp1 to be used fot bit -1
      for s in sb:
          res = doublef(p, q, n, res)
          if s:
              res = addf(p, q, n, res, jp1) if s > 0 else \
                  addf(p, q, n, res, jp0)
      return res
  
  
  # Encapsulates mulf() in order to enable flat coordinates (x, y)
  def mulp(p, q, n, p1, c):
      """Multiply point p by c using fast multiplication"""
      return from_projective(mulf(p, q, n, to_projective(p1), c), n)
  
  
  # Sum of two products using Shamir's trick and signed binary expansion
  def muladdf(p, q, n, jp1, c1, jp2, c2):
      """Efficiently compute c1 * jp1 + c2 * jp2 in projective coordinates"""
      s1 = _signed_bin(c1)
      s2 = _signed_bin(c2)
      diff = len(s2) - len(s1)
      if diff > 0:
          s1 = [0] * diff + s1
      elif diff < 0:
          s2 = [0] * -diff + s2
  
      jp1p2 = addf(p, q, n, jp1, jp2)
      jp1n2 = addf(p, q, n, jp1, neg(jp2, n))
  
      precomp = ((None, jp2, neg(jp2, n)),
                 (jp1, jp1p2, jp1n2),
                 (neg(jp1, n), neg(jp1n2, n), neg(jp1p2, n)))
      res = None
  
      for i, j in zip(s1, s2):
          res = doublef(p, q, n, res)
          if i or j:
              res = addf(p, q, n, res, precomp[i][j])
      return res
  
  
  # Encapsulate muladdf()
  def muladdp(p, q, n, p1, c1, p2, c2):
      """Efficiently compute c1 * p1 + c2 * p2 in (x, y)-coordinates"""
      return from_projective(muladdf(p, q, n,
                                     to_projective(p1), c1,
                                     to_projective(p2), c2), n)
  
  # POINT COMPRESSION ------------------------------------------------------------
  
  # Compute the square root modulo n
  
  
  # Determine the sign-bit of a point allowing to reconstruct y-coordinates
  # when x and the sign-bit are given:
  def sign_bit(p1):
      """Return the signedness of a point p1"""
      return p1[1] % 2 if p1 else 0
  
  # Reconstruct the y-coordinate when curve parameters, x and the sign-bit of
  # the y coordinate are given:
  def y_from_x(x, p, q, n, sign):
      """Return the y coordinate over curve (p, q, n) for given (x, sign)"""
  
      # optimized form of (x**3 - p*x - q) % n
      a = (((x * x) % n - p) * x - q) % n
  
  
  if __name__ == "__main__":
      from Cryptodome.Random.random import randint
      from Cryptodome.Util.number import getPrime
      import time
  
      t = time.time()
      n = getPrime(int(256/8))
      #n = rsa.get_prime(256 / 8, 20)
      tp = time.time() - t
      p = randint(1, n)
      p1 = (randint(1, n), randint(1, n))
      q = curve_q(p1[0], p1[1], p, n)
      r1 = randint(1, n)
      r2 = randint(1, n)
      q1 = mulp(p, q, n, p1, r1)
      q2 = mulp(p, q, n, p1, r2)
      s1 = mulp(p, q, n, q1, r2)
      s2 = mulp(p, q, n, q2, r1)
      s1 == s2
      tt = time.time() - t
  
      def test(tcount, bits=256):
          n = getPrime(int(bits/8))
          #n = rsa.get_prime(bits / 8, 20)
          p = randint(1, n)
          p1 = (randint(1, n), randint(1, n))
          q = curve_q(p1[0], p1[1], p, n)
          p2 = mulp(p, q, n, p1, randint(1, n))
  
          c1 = [randint(1, n) for i in range(tcount)]
          c2 = [randint(1, n) for i in range(tcount)]
          c = list(zip(c1, c2))
  
          t = time.time()
          for i, j in c:
              from_projective(addf(p, q, n,
                                   mulf(p, q, n, to_projective(p1), i),
                                   mulf(p, q, n, to_projective(p2), j)), n)
          t1 = time.time() - t
          t = time.time()
          for i, j in c:
              muladdp(p, q, n, p1, i, p2, j)
          t2 = time.time() - t
  
          return tcount, t1, t2
          
  
          
  