Alien-libsecp256k1

 view release on metacpan or  search on metacpan

libsecp256k1/README.md  view on Meta::CPAN

  * Optimized implementation of arithmetic modulo the curve's field size (2^256 - 0x1000003D1).
    * Using 5 52-bit limbs
    * Using 10 26-bit limbs (including hand-optimized assembly for 32-bit ARM, by Wladimir J. van der Laan).
      * This is an experimental feature that has not received enough scrutiny to satisfy the standard of quality of this library but is made available for testing and review by the community.
* Scalar operations
  * Optimized implementation without data-dependent branches of arithmetic modulo the curve's order.
    * Using 4 64-bit limbs (relying on __int128 support in the compiler).
    * Using 8 32-bit limbs.
* Modular inverses (both field elements and scalars) based on [safegcd](https://gcd.cr.yp.to/index.html) with some modifications, and a variable-time variant (by Peter Dettman).
* Group operations
  * Point addition formula specifically simplified for the curve equation (y^2 = x^3 + 7).
  * Use addition between points in Jacobian and affine coordinates where possible.
  * Use a unified addition/doubling formula where necessary to avoid data-dependent branches.
  * Point/x comparison without a field inversion by comparison in the Jacobian coordinate space.
* Point multiplication for verification (a*P + b*G).
  * Use wNAF notation for point multiplicands.
  * Use a much larger window for multiples of G, using precomputed multiples.
  * Use Shamir's trick to do the multiplication with the public key and the generator simultaneously.
  * Use secp256k1's efficiently-computable endomorphism to split the P multiplicand into 2 half-sized ones.
* Point multiplication for signing
  * Use a precomputed table of multiples of powers of 16 multiplied with the generator, so general multiplication becomes a series of additions.
  * Intended to be completely free of timing sidechannels for secret-key operations (on reasonable hardware/toolchains)
    * Access the table with branch-free conditional moves so memory access is uniform.

libsecp256k1/doc/ellswift.md  view on Meta::CPAN


**Note**: in the paper these conditions result in $\infty$ as output, due to the use of projective coordinates there.
We wish to avoid the need for callers to deal with this special case.

This is implemented in `secp256k1_ellswift_xswiftec_frac_var` (which decodes to an x-coordinate represented as a fraction), and
in `secp256k1_ellswift_xswiftec_var` (which outputs the actual x-coordinate).

## 3. The encoding function

To implement $F_u^{-1}(x)$, the function to find the set of inverses $t$ for which $F_u(t) = x$, we have to reverse the process:
* Find all the $(X, Y) \in S_u$ that could have given rise to $x$, through the $x_1$, $x_2$, or $x_3$ formulas in $\psi_u.$
* Map those $(X, Y)$ solutions to $t$ values using $P_u^{-1}(X, Y).$
* For each of the found $t$ values, verify that $F_u(t) = x.$
* Return the remaining $t$ values.

The function $P_u^{-1}$, which finds $t$ given $(X, Y) \in S_u$, is significantly simpler than $P_u:$

$$
P_u^{-1}(X, Y) = \left\\{\begin{array}{ll}
Yu\sqrt{-3} - X & a = 0 \\
\dfrac{Y-Y_0(u)}{X-X_0(u)} & a \neq 0 \land X \neq X_0(u) \\

libsecp256k1/doc/ellswift.md  view on Meta::CPAN

before those values are computed.

It is not possible that an encoding found through the $x_1$ expression decodes to a different valid x-coordinate using $x_2$ (which would
take precedence), for the same reason: if both $x_1$ and $x_2$ decodings were valid, $x_3$ would be valid as well, and thus take
precedence over both. Because of this, the $g(-u-x)$ being square test for $x_1$ and $x_2$ is the only test necessary to guarantee the found $t$
values round-trip back to the input $x$ correctly. This is the reason for choosing the $(x_3, x_2, x_1)$ precedence order in the decoder;
any order which does not place $x_3$ first requires more complicated round-trip checks in the encoder.

### 3.1 Switching to *v, w* coordinates

Before working out the formulas for all this, we switch to different variables for $S_u.$ Let $v = (X/Y - u)/2$, and
$w = 2Y.$ Or in the other direction, $X = w(u/2 + v)$ and $Y = w/2:$
* $S_u'$ becomes the set of $(v, w)$ for which $w^2 (u^2 + uv + v^2 + a) = -g(u)$ and $w \neq 0.$
* For $a=0$ curves, $P_u^{-1}$ can be stated for $(v,w)$ as $P_u^{'-1}(v, w) = w\left(\frac{\sqrt{-3}-1}{2}u - v\right).$
* $\psi_u$ can be stated for $(v, w)$ as $\psi_u'(v, w) = (x_1, x_2, x_3, z)$, where

$$
\begin{array}{lcl}
  x_1 & = & v \\
  x_2 & = & -u - v \\
  x_3 & = & u + w^2 \\

libsecp256k1/doc/ellswift.md  view on Meta::CPAN

**Define** *ElligatorSwift(x)* as:
* Loop:
  * Pick a uniformly random field element $u.$
  * Compute the set $L = F_u^{-1}(x).$
  * Let $T$ be the 8-element vector consisting of the elements of $L$, plus $8 - \\#L$ times $\\{\bot\\}.$
  * Select a uniformly random $t \in T.$
  * If $t \neq \bot$, return $(u, t)$; restart loop otherwise.

Now notice that the order of elements in $T$ does not matter, as all we do is pick a uniformly
random element in it, so we do not need to have all $\bot$ values at the end.
As we have 8 distinct formulas for finding $(v, w)$ (taking the variants due to $\pm$ into account),
we can associate every index in $T$ with exactly one of those formulas, making sure that:
* Formulas that yield no solutions (due to division by zero or non-existing square roots) or invalid solutions are made to return $\bot.$
* For the $x_1$ and $x_2$ cases, if $g(-u-x)$ is a square, $\bot$ is returned instead (the round-trip check).
* In case multiple formulas would return the same non- $\bot$ result, all but one of those must be turned into $\bot$ to avoid biasing those.

The last condition above only occurs with negligible probability for cryptographically-sized curves, but is interesting
to take into account as it allows exhaustive testing in small groups. See [Section 3.4](#34-dealing-with-special-cases)
for an analysis of all the negligible cases.

If we define $T = (G_{0,u}(x), G_{1,u}(x), \ldots, G_{7,u}(x))$, with each $G_{i,u}$ matching one of the formulas,
the loop can be simplified to only compute one of the inverses instead of all of them:

**Define** *ElligatorSwift(x)* as:
* Loop:
  * Pick a uniformly random field element $u.$
  * Pick a uniformly random integer $c$ in $[0,8).$
  * Let $t = G_{c,u}(x).$
  * If $t \neq \bot$, return $(u, t)$; restart loop otherwise.

This is implemented in `secp256k1_ellswift_xelligatorswift_var`.

### 3.3 Finding the inverse

To implement $G_{c,u}$, we map $c=0$ to the $x_1$ formula, $c=1$ to the $x_2$ formula, and $c=2$ and $c=3$ to the $x_3$ formula.
Those are then repeated as $c=4$ through $c=7$ for the other sign of $w$ (noting that in each formula, $w$ is a square root of some expression).
Ignoring the negligible cases, we get:

**Define** $G_{c,u}(x)$ as:
* If $c \in \\{0, 1, 4, 5\\}$ (for $x_1$ and $x_2$ formulas):
  * If $g(-u-x)$ is square, return $\bot$ (as $x_3$ would be valid and take precedence).
  * If $c \in \\{0, 4\\}$ (the $x_1$ formula) let $v = x$, otherwise let $v = -u-x$ (the $x_2$ formula)
  * Let $s = -g(u)/(u^2 + uv + v^2 + a)$ (using $s = w^2$ in what follows).
* Otherwise, when $c \in \\{2, 3, 6, 7\\}$ (for $x_3$ formulas):
  * Let $s = x-u.$
  * Let $r = \sqrt{-s(4g(u) + sh(u))}.$
  * Let $v = (r/s - u)/2$ if $c \in \\{3, 7\\}$; $(-r/s - u)/2$ otherwise.
* Let $w = \sqrt{s}.$
* Depending on $c:$
  * If $c \in \\{0, 1, 2, 3\\}:$ return $P_u^{'-1}(v, w).$
  * If $c \in \\{4, 5, 6, 7\\}:$ return $P_u^{'-1}(v, -w).$

Whenever a square root of a non-square is taken, $\bot$ is returned; for both square roots this happens with roughly
50% on random inputs. Similarly, when a division by 0 would occur, $\bot$ is returned as well; this will only happen

libsecp256k1/doc/ellswift.md  view on Meta::CPAN

  * If $c \in \\{5, 7\\}$, let $t = P_u^{'-1}(-u-v, -w).$
* If $a=0$ and $t=0$, return $\bot$ (even curves only).
* If $a \neq 0$ and $h(u)t^2 = -1$, return $\bot.$
* Return $t.$

Given any $u$, using this algorithm over all $x$ and $c$ values, every $t$ value will be reached exactly once,
for an $x$ for which $F_u(t) = x$ holds, except for these cases that will not be reached:
* All cases where $P_u(t)$ is not defined:
  * For $a=0$ curves, when $u=0$, $t=0$, or $g(u) = -t^2.$
  * For $a \neq 0$ curves, when $h(u)t^2 = -1$, $X_0(u) = 0$, or $Y_0(u) (1 - h(u) t^2) = 2X_0(u)t.$
* When $g(u)=0$, the potentially many $t$ values that decode to an $x$ satisfying $g(x)=0$ using the $x_2$ formula. These were excluded by the $g(u)=0$ condition in the $c \in \\{0, 1, 4, 5\\}$ branch.

These cases form a negligible subset of all $(u, t)$ for cryptographically sized curves.

### 3.5 Encoding for `secp256k1`

Specialized for odd-ordered $a=0$ curves:

**Define** $G_{c,u}(x)$ as:
* If $u=0$, return $\bot.$
* If $c \in \\{0, 1, 4, 5\\}:$

libsecp256k1/doc/safegcd_implementation.md  view on Meta::CPAN

    # Compute w = -g/f mod 2**limit, using the table value for -1/f mod 2**4. Note that f is
    # always odd, so its inverse modulo a power of two always exists.
    w = (g * NEGINV16[(f & 15) // 2]) % (2**limit)
    # As w = -g/f mod (2**limit), g+w*f mod 2**limit = 0 mod 2**limit.
    g += w * f
    assert g % (2**limit) == 0
    # The next iteration will now shift out at least limit bottom zero bits from g.
```

By using a bigger table more bits can be cancelled at once. The table can also be implemented
as a formula. Several formulas are known for computing modular inverses modulo powers of two;
some can be found in Hacker's Delight second edition by Henry S. Warren, Jr. pages 245-247.
Here we need the negated modular inverse, which is a simple transformation of those:

- Instead of a 3-bit table:
  - *-f* or *f ^ 6*
- Instead of a 4-bit table:
  - *1 - f(f + 1)*
  - *-(f + (((f + 1) & 4) << 1))*
- For larger tables the following technique can be used: if *w=-1/f mod 2<sup>L</sup>*, then *w(w&thinsp;f+2)* is
  *-1/f mod 2<sup>2L</sup>*. This allows extending the previous formulas (or tables). In particular we
  have this 6-bit function (based on the 3-bit function above):
  - *f(f<sup>2</sup> - 2)*

This loop, again extended to also handle *u*, *v*, *q*, and *r* alongside *f* and *g*, placed in
`divsteps_n_matrix`, gives a significantly faster, but non-constant time version.


## 7. Final Python version

All together we need the following functions:

libsecp256k1/sage/group_prover.sage  view on Meta::CPAN

        if normalize_factor(zeroextra.reduce(f)) not in nonzero:
          ret.add("%s != 0" % normalize_factor(zeroextra.reduce(f)))
  return ", ".join(x for x in ret)


def check_symbolic(R, assumeLaw, assumeAssert, assumeBranch, require):
  """Check a set of zero and nonzero requirements, given a set of zero and nonzero assumptions"""
  assume = assumeLaw + assumeAssert + assumeBranch

  if conflicts(R, assume):
    # This formula does not apply
    return (True, None)

  describe = describe_extra(R, assumeLaw + assumeBranch, assumeAssert)
  if describe != "":
    describe = " (assuming " + describe + ")"

  ok, msg = prove_zero(R, require.zero, assume)
  if not ok:
    return (False, "FAIL, %s fails%s" % (str(msg), describe))

libsecp256k1/sage/prove_group_implementations.sage  view on Meta::CPAN

# Test libsecp256k1' group operation implementations using prover.sage

import sys

load("group_prover.sage")
load("weierstrass_prover.sage")

def formula_secp256k1_gej_double_var(a):
  """libsecp256k1's secp256k1_gej_double_var, used by various addition functions"""
  rz = a.Z * a.Y
  s = a.Y^2
  l = a.X^2
  l = l * 3
  l = l / 2
  t = -s
  t = t * a.X
  rx = l^2
  rx = rx + t
  rx = rx + t
  s = s^2
  t = t + rx
  ry = t * l
  ry = ry + s
  ry = -ry
  return jacobianpoint(rx, ry, rz)

def formula_secp256k1_gej_add_var(branch, a, b):
  """libsecp256k1's secp256k1_gej_add_var"""
  if branch == 0:
    return (constraints(), constraints(nonzero={a.Infinity : 'a_infinite'}), b)
  if branch == 1:
    return (constraints(), constraints(zero={a.Infinity : 'a_finite'}, nonzero={b.Infinity : 'b_infinite'}), a)
  z22 = b.Z^2
  z12 = a.Z^2
  u1 = a.X * z22
  u2 = b.X * z12
  s1 = a.Y * z22
  s1 = s1 * b.Z
  s2 = b.Y * z12
  s2 = s2 * a.Z
  h = -u1
  h = h + u2
  i = -s2
  i = i + s1
  if branch == 2:
    r = formula_secp256k1_gej_double_var(a)
    return (constraints(), constraints(zero={h : 'h=0', i : 'i=0', a.Infinity : 'a_finite', b.Infinity : 'b_finite'}), r)
  if branch == 3:
    return (constraints(), constraints(zero={h : 'h=0', a.Infinity : 'a_finite', b.Infinity : 'b_finite'}, nonzero={i : 'i!=0'}), point_at_infinity())
  t = h * b.Z
  rz = a.Z * t
  h2 = h^2
  h2 = -h2
  h3 = h2 * h
  t = u1 * h2
  rx = i^2
  rx = rx + h3
  rx = rx + t
  rx = rx + t
  t = t + rx
  ry = t * i
  h3 = h3 * s1
  ry = ry + h3
  return (constraints(), constraints(zero={a.Infinity : 'a_finite', b.Infinity : 'b_finite'}, nonzero={h : 'h!=0'}), jacobianpoint(rx, ry, rz))

def formula_secp256k1_gej_add_ge_var(branch, a, b):
  """libsecp256k1's secp256k1_gej_add_ge_var, which assume bz==1"""
  if branch == 0:
    return (constraints(zero={b.Z - 1 : 'b.z=1'}), constraints(nonzero={a.Infinity : 'a_infinite'}), b)
  if branch == 1:
    return (constraints(zero={b.Z - 1 : 'b.z=1'}), constraints(zero={a.Infinity : 'a_finite'}, nonzero={b.Infinity : 'b_infinite'}), a)
  z12 = a.Z^2
  u1 = a.X
  u2 = b.X * z12
  s1 = a.Y
  s2 = b.Y * z12
  s2 = s2 * a.Z
  h = -u1
  h = h + u2
  i = -s2
  i = i + s1
  if (branch == 2):
    r = formula_secp256k1_gej_double_var(a)
    return (constraints(zero={b.Z - 1 : 'b.z=1'}), constraints(zero={a.Infinity : 'a_finite', b.Infinity : 'b_finite', h : 'h=0', i : 'i=0'}), r)
  if (branch == 3):
    return (constraints(zero={b.Z - 1 : 'b.z=1'}), constraints(zero={a.Infinity : 'a_finite', b.Infinity : 'b_finite', h : 'h=0'}, nonzero={i : 'i!=0'}), point_at_infinity())
  rz = a.Z * h
  h2 = h^2
  h2 = -h2
  h3 = h2 * h
  t = u1 * h2
  rx = i^2
  rx = rx + h3
  rx = rx + t
  rx = rx + t
  t = t + rx
  ry = t * i
  h3 = h3 * s1
  ry = ry + h3
  return (constraints(zero={b.Z - 1 : 'b.z=1'}), constraints(zero={a.Infinity : 'a_finite', b.Infinity : 'b_finite'}, nonzero={h : 'h!=0'}), jacobianpoint(rx, ry, rz))

def formula_secp256k1_gej_add_zinv_var(branch, a, b):
  """libsecp256k1's secp256k1_gej_add_zinv_var"""
  bzinv = b.Z^(-1)
  if branch == 0:
    rinf = b.Infinity
    bzinv2 = bzinv^2
    bzinv3 = bzinv2 * bzinv
    rx = b.X * bzinv2
    ry = b.Y * bzinv3
    rz = 1
    return (constraints(), constraints(nonzero={a.Infinity : 'a_infinite'}), jacobianpoint(rx, ry, rz, rinf))

libsecp256k1/sage/prove_group_implementations.sage  view on Meta::CPAN

  u1 = a.X
  u2 = b.X * z12
  s1 = a.Y
  s2 = b.Y * z12
  s2 = s2 * azz
  h = -u1
  h = h + u2
  i = -s2
  i = i + s1
  if branch == 2:
    r = formula_secp256k1_gej_double_var(a)
    return (constraints(), constraints(zero={a.Infinity : 'a_finite', b.Infinity : 'b_finite', h : 'h=0', i : 'i=0'}), r)
  if branch == 3:
    return (constraints(), constraints(zero={a.Infinity : 'a_finite', b.Infinity : 'b_finite', h : 'h=0'}, nonzero={i : 'i!=0'}), point_at_infinity())
  rz = a.Z * h
  h2 = h^2
  h2 = -h2
  h3 = h2 * h
  t = u1 * h2
  rx = i^2
  rx = rx + h3
  rx = rx + t
  rx = rx + t
  t = t + rx
  ry = t * i
  h3 = h3 * s1
  ry = ry + h3
  return (constraints(), constraints(zero={a.Infinity : 'a_finite', b.Infinity : 'b_finite'}, nonzero={h : 'h!=0'}), jacobianpoint(rx, ry, rz))

def formula_secp256k1_gej_add_ge(branch, a, b):
  """libsecp256k1's secp256k1_gej_add_ge"""
  zeroes = {}
  nonzeroes = {}
  a_infinity = False
  if (branch & 2) != 0:
    nonzeroes.update({a.Infinity : 'a_infinite'})
    a_infinity = True
  else:
    zeroes.update({a.Infinity : 'a_finite'})
  zz = a.Z^2

libsecp256k1/sage/prove_group_implementations.sage  view on Meta::CPAN

    rx = b.X
    ry = b.Y
    rz = 1
  if (branch & 4) != 0:
    zeroes.update({rz : 'r.z = 0'})
    return (constraints(zero={b.Z - 1 : 'b.z=1', b.Infinity : 'b_finite'}), constraints(zero=zeroes, nonzero=nonzeroes), point_at_infinity())
  else:
    nonzeroes.update({rz : 'r.z != 0'})
  return (constraints(zero={b.Z - 1 : 'b.z=1', b.Infinity : 'b_finite'}), constraints(zero=zeroes, nonzero=nonzeroes), jacobianpoint(rx, ry, rz))

def formula_secp256k1_gej_add_ge_old(branch, a, b):
  """libsecp256k1's old secp256k1_gej_add_ge, which fails when ay+by=0 but ax!=bx"""
  a_infinity = (branch & 1) != 0
  zero = {}
  nonzero = {}
  if a_infinity:
    nonzero.update({a.Infinity : 'a_infinite'})
  else:
    zero.update({a.Infinity : 'a_finite'})
  zz = a.Z^2
  u1 = a.X

libsecp256k1/sage/prove_group_implementations.sage  view on Meta::CPAN

  t = t * (1 if a_infinity else 0)
  ry = ry + t
  t = (1 if a_infinity else 0)
  rz = rz + t
  if infinity:
    return (constraints(zero={b.Z - 1 : 'b.z=1', b.Infinity : 'b_finite'}), constraints(zero=zero, nonzero=nonzero), point_at_infinity())
  return (constraints(zero={b.Z - 1 : 'b.z=1', b.Infinity : 'b_finite'}), constraints(zero=zero, nonzero=nonzero), jacobianpoint(rx, ry, rz))

if __name__ == "__main__":
  success = True
  success = success & check_symbolic_jacobian_weierstrass("secp256k1_gej_add_var", 0, 7, 5, formula_secp256k1_gej_add_var)
  success = success & check_symbolic_jacobian_weierstrass("secp256k1_gej_add_ge_var", 0, 7, 5, formula_secp256k1_gej_add_ge_var)
  success = success & check_symbolic_jacobian_weierstrass("secp256k1_gej_add_zinv_var", 0, 7, 5, formula_secp256k1_gej_add_zinv_var)
  success = success & check_symbolic_jacobian_weierstrass("secp256k1_gej_add_ge", 0, 7, 8, formula_secp256k1_gej_add_ge)
  success = success & (not check_symbolic_jacobian_weierstrass("secp256k1_gej_add_ge_old [should fail]", 0, 7, 4, formula_secp256k1_gej_add_ge_old))

  if len(sys.argv) >= 2 and sys.argv[1] == "--exhaustive":
    success = success & check_exhaustive_jacobian_weierstrass("secp256k1_gej_add_var", 0, 7, 5, formula_secp256k1_gej_add_var, 43)
    success = success & check_exhaustive_jacobian_weierstrass("secp256k1_gej_add_ge_var", 0, 7, 5, formula_secp256k1_gej_add_ge_var, 43)
    success = success & check_exhaustive_jacobian_weierstrass("secp256k1_gej_add_zinv_var", 0, 7, 5, formula_secp256k1_gej_add_zinv_var, 43)
    success = success & check_exhaustive_jacobian_weierstrass("secp256k1_gej_add_ge", 0, 7, 8, formula_secp256k1_gej_add_ge, 43)
    success = success & (not check_exhaustive_jacobian_weierstrass("secp256k1_gej_add_ge_old [should fail]", 0, 7, 4, formula_secp256k1_gej_add_ge_old, 43))

  sys.exit(int(not success))

libsecp256k1/sage/weierstrass_prover.sage  view on Meta::CPAN

laws_jacobian_weierstrass = {
  'add': law_jacobian_weierstrass_add,
  'double': law_jacobian_weierstrass_double,
  'add_opposite': law_jacobian_weierstrass_add_opposites,
  'add_infinite_a': law_jacobian_weierstrass_add_infinite_a,
  'add_infinite_b': law_jacobian_weierstrass_add_infinite_b,
  'add_infinite_ab': law_jacobian_weierstrass_add_infinite_ab
}


def check_exhaustive_jacobian_weierstrass(name, A, B, branches, formula, p):
  """Verify an implementation of addition of Jacobian points on a Weierstrass curve, by executing and validating the result for every possible addition in a prime field"""
  F = Integers(p)
  print("Formula %s on Z%i:" % (name, p))
  points = []
  for x in range(0, p):
    for y in range(0, p):
      point = affinepoint(F(x), F(y))
      r, e = concrete_verify(on_weierstrass_curve(A, B, point))
      if r:
        points.append(point)

libsecp256k1/sage/weierstrass_prover.sage  view on Meta::CPAN

  ret = True
  for za in range(1, p):
    for zb in range(1, p):
      for pa in points:
        for pb in points:
          for ia in range(2):
            for ib in range(2):
              pA = jacobianpoint(pa.x * F(za)^2, pa.y * F(za)^3, F(za), ia)
              pB = jacobianpoint(pb.x * F(zb)^2, pb.y * F(zb)^3, F(zb), ib)
              for branch in range(0, branches):
                assumeAssert, assumeBranch, pC = formula(branch, pA, pB)
                pC.X = F(pC.X)
                pC.Y = F(pC.Y)
                pC.Z = F(pC.Z)
                pC.Infinity = F(pC.Infinity)
                r, e = concrete_verify(assumeAssert + assumeBranch)
                if r:
                  match = False
                  for key in laws_jacobian_weierstrass:
                    assumeLaw, require = laws_jacobian_weierstrass[key](A, B, pa, pb, pA, pB, pC)
                    r, e = concrete_verify(assumeLaw)

libsecp256k1/sage/weierstrass_prover.sage  view on Meta::CPAN

                        print("  failure in branch %i for (%s,%s,%s,%s) + (%s,%s,%s,%s) = (%s,%s,%s,%s): %s" % (branch, pA.X, pA.Y, pA.Z, pA.Infinity, pB.X, pB.Y, pB.Z, pB.Infinity, pC.X, pC.Y, pC.Z, pC.Infinity, e))

  print()
  return ret


def check_symbolic_function(R, assumeAssert, assumeBranch, f, A, B, pa, pb, pA, pB, pC):
  assumeLaw, require = f(A, B, pa, pb, pA, pB, pC)
  return check_symbolic(R, assumeLaw, assumeAssert, assumeBranch, require)

def check_symbolic_jacobian_weierstrass(name, A, B, branches, formula):
  """Verify an implementation of addition of Jacobian points on a Weierstrass curve symbolically"""
  R.<ax,bx,ay,by,Az,Bz,Ai,Bi> = PolynomialRing(QQ,8,order='invlex')
  lift = lambda x: fastfrac(R,x)
  ax = lift(ax)
  ay = lift(ay)
  Az = lift(Az)
  bx = lift(bx)
  by = lift(by)
  Bz = lift(Bz)
  Ai = lift(Ai)

libsecp256k1/sage/weierstrass_prover.sage  view on Meta::CPAN


  res = {}

  for key in laws_jacobian_weierstrass:
    res[key] = []

  print("Formula " + name + ":")
  count = 0
  ret = True
  for branch in range(branches):
    assumeFormula, assumeBranch, pC = formula(branch, pA, pB)
    assumeBranch = assumeBranch.map(lift)
    assumeFormula = assumeFormula.map(lift)
    pC.X = lift(pC.X)
    pC.Y = lift(pC.Y)
    pC.Z = lift(pC.Z)
    pC.Infinity = lift(pC.Infinity)

    for key in laws_jacobian_weierstrass:
      success, msg = check_symbolic_function(R, assumeFormula, assumeBranch, laws_jacobian_weierstrass[key], A, B, pa, pb, pA, pB, pC)
      if not success:

libsecp256k1/src/ecmult_const_impl.h  view on Meta::CPAN

         * the comment in ecmult_gen_impl.h for rationale. */ \
        secp256k1_fe_cmov(&(r)->x, &(pre)[m].x, m == index); \
        secp256k1_fe_cmov(&(r)->y, &(pre)[m].y, m == index); \
    } \
    (r)->infinity = 0; \
    secp256k1_fe_negate(&neg_y, &(r)->y, 1); \
    secp256k1_fe_cmov(&(r)->y, &neg_y, negative); \
} while(0)

/* For K as defined in the comment of secp256k1_ecmult_const, we have several precomputed
 * formulas/constants.
 * - in exhaustive test mode, we give an explicit expression to compute it at compile time: */
#ifdef EXHAUSTIVE_TEST_ORDER
static const secp256k1_scalar secp256k1_ecmult_const_K = ((SECP256K1_SCALAR_CONST(0, 0, 0, (1U << (ECMULT_CONST_BITS - 128)) - 2U, 0, 0, 0, 0) + EXHAUSTIVE_TEST_ORDER - 1U) * (1U + EXHAUSTIVE_TEST_LAMBDA)) % EXHAUSTIVE_TEST_ORDER;
/* - for the real secp256k1 group we have constants for various ECMULT_CONST_BITS values. */
#elif ECMULT_CONST_BITS == 129
/* For GROUP_SIZE = 1,3. */
static const secp256k1_scalar secp256k1_ecmult_const_K = SECP256K1_SCALAR_CONST(0xac9c52b3ul, 0x3fa3cf1ful, 0x5ad9e3fdul, 0x77ed9ba4ul, 0xa880b9fcul, 0x8ec739c2ul, 0xe0cfc810ul, 0xb51283ceul);
#elif ECMULT_CONST_BITS == 130
/* For GROUP_SIZE = 2,5. */
static const secp256k1_scalar secp256k1_ecmult_const_K = SECP256K1_SCALAR_CONST(0xa4e88a7dul, 0xcb13034eul, 0xc2bdd6bful, 0x7c118d6bul, 0x589ae848ul, 0x26ba29e4ul, 0xb5c2c1dcul, 0xde9798d9ul);

libsecp256k1/src/ecmult_const_impl.h  view on Meta::CPAN

    /* Verify that v1 and v2 are in range [0, 2^129-1]. */
    for (i = 129; i < 256; ++i) {
        VERIFY_CHECK(secp256k1_scalar_get_bits_limb32(&v1, i, 1) == 0);
        VERIFY_CHECK(secp256k1_scalar_get_bits_limb32(&v2, i, 1) == 0);
    }
#endif

    /* Calculate odd multiples of A and A*lambda.
     * All multiples are brought to the same Z 'denominator', which is stored
     * in global_z. Due to secp256k1' isomorphism we can do all operations pretending
     * that the Z coordinate was 1, use affine addition formulae, and correct
     * the Z coordinate of the result once at the end.
     */
    secp256k1_gej_set_ge(r, a);
    secp256k1_ecmult_const_odd_multiples_table_globalz(pre_a, &global_z, r);
    for (i = 0; i < ECMULT_CONST_TABLE_SIZE; i++) {
        secp256k1_ge_mul_lambda(&pre_a_lam[i], &pre_a[i]);
    }

    /* Next, we compute r = C_l(v1, A) + C_l(v2, lambda*A).
     *

libsecp256k1/src/ecmult_const_impl.h  view on Meta::CPAN

    /* This algorithm is a generalization of Peter Dettman's technique for
     * avoiding the square root in a random-basepoint x-only multiplication
     * on a Weierstrass curve:
     * https://mailarchive.ietf.org/arch/msg/cfrg/7DyYY6gg32wDgHAhgSb6XxMDlJA/
     *
     *
     * === Background: the effective affine technique ===
     *
     * Let phi_u be the isomorphism that maps (x, y) on secp256k1 curve y^2 = x^3 + 7 to
     * x' = u^2*x, y' = u^3*y on curve y'^2 = x'^3 + u^6*7. This new curve has the same order as
     * the original (it is isomorphic), but moreover, has the same addition/doubling formulas, as
     * the curve b=7 coefficient does not appear in those formulas (or at least does not appear in
     * the formulas implemented in this codebase, both affine and Jacobian). See also Example 9.5.2
     * in https://www.math.auckland.ac.nz/~sgal018/crypto-book/ch9.pdf.
     *
     * This means any linear combination of secp256k1 points can be computed by applying phi_u
     * (with non-zero u) on all input points (including the generator, if used), computing the
     * linear combination on the isomorphic curve (using the same group laws), and then applying
     * phi_u^{-1} to get back to secp256k1.
     *
     * Switching to Jacobian coordinates, note that phi_u applied to (X, Y, Z) is simply
     * (X, Y, Z/u). Thus, if we want to compute (X1, Y1, Z) + (X2, Y2, Z), with identical Z
     * coordinates, we can use phi_Z to transform it to (X1, Y1, 1) + (X2, Y2, 1) on an isomorphic
     * curve where the affine addition formula can be used instead.
     * If (X3, Y3, Z3) = (X1, Y1) + (X2, Y2) on that curve, then our answer on secp256k1 is
     * (X3, Y3, Z3*Z).
     *
     * This is the effective affine technique: if we have a linear combination of group elements
     * to compute, and all those group elements have the same Z coordinate, we can simply pretend
     * that all those Z coordinates are 1, perform the computation that way, and then multiply the
     * original Z coordinate back in.
     *
     * The technique works on any a=0 short Weierstrass curve. It is possible to generalize it to
     * other curves too, but there the isomorphic curves will have different 'a' coefficients,

libsecp256k1/src/ecmult_gen_impl.h  view on Meta::CPAN

     *   for comb_off in range(COMB_SPACING - 1, -1, -1):
     *     for block in range(COMB_BLOCKS):
     *       c += table(block, (d >> comb_off) & mask(block))
     *     if comb_off > 0:
     *       c = 2*c
     *   return c
     *
     * This computes c = comb(d, G/2), and thus finally R = c + ctx->ge_offset. Note that it would
     * be possible to apply an initial offset instead of a final offset (moving ge_offset to take
     * the place of infinity above), but the chosen approach allows using (in a future improvement)
     * an incomplete addition formula for most of the multiplication.
     *
     * The last question is how to implement the table(b, m) function. For any value of b,
     * m=(d & mask(b)) can only take on at most 2^COMB_TEETH possible values (the last one may have
     * fewer as there mask(b) may exceed the curve order). So we could create COMB_BLOCK tables
     * which contain a value for each such m value.
     *
     * Now note that if m=(d & mask(b)), then flipping the relevant bits of m results in negating
     * the result of table(b, m). This is because table(b,m XOR mask(b)) = table(b, mask(b) - m) =
     * (mask(b) - m - mask(b)/2)*G = (-m + mask(b)/2)*G = -(m - mask(b)/2)*G = -table(b, m).
     * Because of this it suffices to only store the first half of the m values for every b. If an

libsecp256k1/src/ecmult_impl.h  view on Meta::CPAN

        if (state->ps[no].bits_na_1 > bits) {
            bits = state->ps[no].bits_na_1;
        }
        if (state->ps[no].bits_na_lam > bits) {
            bits = state->ps[no].bits_na_lam;
        }

        /* Calculate odd multiples of a.
         * All multiples are brought to the same Z 'denominator', which is stored
         * in Z. Due to secp256k1' isomorphism we can do all operations pretending
         * that the Z coordinate was 1, use affine addition formulae, and correct
         * the Z coordinate of the result once at the end.
         * The exception is the precomputed G table points, which are actually
         * affine. Compared to the base used for other points, they have a Z ratio
         * of 1/Z, so we can use secp256k1_gej_add_zinv_var, which uses the same
         * isomorphism to efficiently add with a known Z inverse.
         */
        tmp = a[np];
        if (no) {
            secp256k1_gej_rescale(&tmp, &Z);
        }

libsecp256k1/src/group_impl.h  view on Meta::CPAN

    secp256k1_fe zz, u1, u2, s1, s2, t, tt, m, n, q, rr;
    secp256k1_fe m_alt, rr_alt;
    int degenerate;
    SECP256K1_GEJ_VERIFY(a);
    SECP256K1_GE_VERIFY(b);
    VERIFY_CHECK(!b->infinity);

    /*  In:
     *    Eric Brier and Marc Joye, Weierstrass Elliptic Curves and Side-Channel Attacks.
     *    In D. Naccache and P. Paillier, Eds., Public Key Cryptography, vol. 2274 of Lecture Notes in Computer Science, pages 335-345. Springer-Verlag, 2002.
     *  we find as solution for a unified addition/doubling formula:
     *    lambda = ((x1 + x2)^2 - x1 * x2 + a) / (y1 + y2), with a = 0 for secp256k1's curve equation.
     *    x3 = lambda^2 - (x1 + x2)
     *    2*y3 = lambda * (x1 + x2 - 2 * x3) - (y1 + y2).
     *
     *  Substituting x_i = Xi / Zi^2 and yi = Yi / Zi^3, for i=1,2,3, gives:
     *    U1 = X1*Z2^2, U2 = X2*Z1^2
     *    S1 = Y1*Z2^3, S2 = Y2*Z1^3
     *    Z = Z1*Z2
     *    T = U1+U2
     *    M = S1+S2
     *    Q = -T*M^2
     *    R = T^2-U1*U2
     *    X3 = R^2+Q
     *    Y3 = -(R*(2*X3+Q)+M^4)/2
     *    Z3 = M*Z
     *  (Note that the paper uses xi = Xi / Zi and yi = Yi / Zi instead.)
     *
     *  This formula has the benefit of being the same for both addition
     *  of distinct points and doubling. However, it breaks down in the
     *  case that either point is infinity, or that y1 = -y2. We handle
     *  these cases in the following ways:
     *
     *    - If b is infinity we simply bail by means of a VERIFY_CHECK.
     *
     *    - If a is infinity, we detect this, and at the end of the
     *      computation replace the result (which will be meaningless,
     *      but we compute to be constant-time) with b.x : b.y : 1.
     *

libsecp256k1/src/modinv64_impl.h  view on Meta::CPAN

        VERIFY_CHECK((q * f0 + r * g0) == g << (62 - i));
        /* Bounds on eta that follow from the bounds on iteration count (max 12*62 divsteps). */
        VERIFY_CHECK(eta >= -745 && eta <= 745);
        /* If eta is negative, negate it and replace f,g with g,-f. */
        if (eta < 0) {
            uint64_t tmp;
            eta = -eta;
            tmp = f; f = g; g = -tmp;
            tmp = u; u = q; q = -tmp;
            tmp = v; v = r; r = -tmp;
            /* Use a formula to cancel out up to 6 bits of g. Also, no more than i can be cancelled
             * out (as we'd be done before that point), and no more than eta+1 can be done as its
             * sign will flip again once that happens. */
            limit = ((int)eta + 1) > i ? i : ((int)eta + 1);
            VERIFY_CHECK(limit > 0 && limit <= 62);
            /* m is a mask for the bottom min(limit, 6) bits. */
            m = (UINT64_MAX >> (64 - limit)) & 63U;
            /* Find what multiple of f must be added to g to cancel its bottom min(limit, 6)
             * bits. */
            w = (f * g * (f * f - 2)) & m;
        } else {
            /* In this branch, use a simpler formula that only lets us cancel up to 4 bits of g, as
             * eta tends to be smaller here. */
            limit = ((int)eta + 1) > i ? i : ((int)eta + 1);
            VERIFY_CHECK(limit > 0 && limit <= 62);
            /* m is a mask for the bottom min(limit, 4) bits. */
            m = (UINT64_MAX >> (64 - limit)) & 15U;
            /* Find what multiple of f must be added to g to cancel its bottom min(limit, 4)
             * bits. */
            w = f + (((f + 1) & 4) << 1);
            w = (-w * g) & m;
        }

libsecp256k1/src/modinv64_impl.h  view on Meta::CPAN

        /* If eta is negative, negate it and replace f,g with g,f. */
        if (eta < 0) {
            uint64_t tmp;
            eta = -eta;
            tmp = f; f = g; g = tmp;
            tmp = u; u = q; q = tmp;
            tmp = v; v = r; r = tmp;
            /* Update bottom bit of jac: when swapping f and g, the Jacobi symbol changes sign
             * if both f and g are 3 mod 4. */
            jac ^= ((f & g) >> 1);
            /* Use a formula to cancel out up to 6 bits of g. Also, no more than i can be cancelled
             * out (as we'd be done before that point), and no more than eta+1 can be done as its
             * sign will flip again once that happens. */
            limit = ((int)eta + 1) > i ? i : ((int)eta + 1);
            VERIFY_CHECK(limit > 0 && limit <= 62);
            /* m is a mask for the bottom min(limit, 6) bits. */
            m = (UINT64_MAX >> (64 - limit)) & 63U;
            /* Find what multiple of f must be added to g to cancel its bottom min(limit, 6)
             * bits. */
            w = (f * g * (f * f - 2)) & m;
        } else {
            /* In this branch, use a simpler formula that only lets us cancel up to 4 bits of g, as
             * eta tends to be smaller here. */
            limit = ((int)eta + 1) > i ? i : ((int)eta + 1);
            VERIFY_CHECK(limit > 0 && limit <= 62);
            /* m is a mask for the bottom min(limit, 4) bits. */
            m = (UINT64_MAX >> (64 - limit)) & 15U;
            /* Find what multiple of f must be added to g to cancel its bottom min(limit, 4)
             * bits. */
            w = f + (((f + 1) & 4) << 1);
            w = (-w * g) & m;
        }

libsecp256k1/src/modules/ellswift/main_impl.h  view on Meta::CPAN

    int ret;

    secp256k1_fe_normalize_weak(&x);
    secp256k1_fe_normalize_weak(&u);

    VERIFY_CHECK(c >= 0 && c < 8);
    VERIFY_CHECK(secp256k1_ge_x_on_curve_var(&x));

    if (!(c & 2)) {
        /* c is in {0, 1, 4, 5}. In this case we look for an inverse under the x1 (if c=0 or
         * c=4) formula, or x2 (if c=1 or c=5) formula. */

        /* If -u-x is a valid X coordinate, fail. This would yield an encoding that roundtrips
         * back under the x3 formula instead (which has priority over x1 and x2, so the decoding
         * would not match x). */
        m = x;                                          /* m = x */
        secp256k1_fe_add(&m, &u);                       /* m = u+x */
        secp256k1_fe_negate(&m, &m, 2);                 /* m = -u-x */
        /* Test if (-u-x) is a valid X coordinate. If so, fail. */
        if (secp256k1_ge_x_on_curve_var(&m)) return 0;

        /* Let s = -(u^3 + 7)/(u^2 + u*x + x^2) [first part] */
        secp256k1_fe_sqr(&s, &m);                       /* s = (u+x)^2 */
        secp256k1_fe_negate(&s, &s, 1);                 /* s = -(u+x)^2 */

libsecp256k1/src/modules/ellswift/main_impl.h  view on Meta::CPAN

        secp256k1_fe_mul(&m, &s, &g);                   /* m = -(u^3 + 7)*(u^2 + u*x + x^2) */
        if (!secp256k1_fe_is_square_var(&m)) return 0;

        /* Let s = -(u^3 + 7)/(u^2 + u*x + x^2) [second part] */
        secp256k1_fe_inv_var(&s, &s);                   /* s = -1/(u^2 + u*x + x^2) [no div by 0] */
        secp256k1_fe_mul(&s, &s, &g);                   /* s = -(u^3 + 7)/(u^2 + u*x + x^2) */

        /* Let v = x. */
        v = x;
    } else {
        /* c is in {2, 3, 6, 7}. In this case we look for an inverse under the x3 formula. */

        /* Let s = x-u. */
        secp256k1_fe_negate(&m, &u, 1);                 /* m = -u */
        s = m;                                          /* s = -u */
        secp256k1_fe_add(&s, &x);                       /* s = x-u */

        /* If s is not square, fail. */
        if (!secp256k1_fe_is_square_var(&s)) return 0;

        /* Let r = sqrt(-s*(4*(u^3+7)+3*u^2*s)); fail if it doesn't exist. */

libsecp256k1/src/tests_exhaustive.c  view on Meta::CPAN

    uint64_t iter = 0;

    /* Sanity-check (and check infinity functions) */
    CHECK(secp256k1_ge_is_infinity(&group[0]));
    CHECK(secp256k1_gej_is_infinity(&groupj[0]));
    for (i = 1; i < EXHAUSTIVE_TEST_ORDER; i++) {
        CHECK(!secp256k1_ge_is_infinity(&group[i]));
        CHECK(!secp256k1_gej_is_infinity(&groupj[i]));
    }

    /* Check all addition formulae */
    for (j = 0; j < EXHAUSTIVE_TEST_ORDER; j++) {
        secp256k1_fe fe_inv;
        if (skip_section(&iter)) continue;
        secp256k1_fe_inv(&fe_inv, &groupj[j].z);
        for (i = 0; i < EXHAUSTIVE_TEST_ORDER; i++) {
            secp256k1_ge zless_gej;
            secp256k1_gej tmp;
            /* add_var */
            secp256k1_gej_add_var(&tmp, &groupj[i], &groupj[j], NULL);
            CHECK(secp256k1_gej_eq_ge_var(&tmp, &group[(i + j) % EXHAUSTIVE_TEST_ORDER]));



( run in 0.578 second using v1.01-cache-2.11-cpan-3cd7ad12f66 )