Try   HackMD

Computing the Optimal Ate Pairing over the BLS12-381 Curve

The post BLS12-381 For The Rest Of Us constitutes the base point for this article, so I highly recommend reading it to get a solid background of the BLS1-381 curve. This post will reuse some parts with an emphasy on the computational aspect of the optimal Ate pairing over such curve.

For those with a poor background on Pairing-Based Cryptography (PBC) but with a decent background on Elliptic-Curve Cryptoraphy (ECC) I recommend reading the following resources in order:

  1. Pairings For Beginners.
  2. Constructing Elliptic Curves with Prescribed Embedding Degrees.
  3. BLS12-381: New zk-SNARK Elliptic Curve Construction.

Only after (and not before) you feel sufficiently comfortable with both PBC and ECC come back to this post and start to make sense out of it. The full list of consulted resources can be found at the end of this post.

The implementation in zkASM has been done in this PR.

The Curve

The BLS12-381 is the elliptic curve

E of the form:
y2=x3+4

defined over a finite field
Fp
, where the prime
p
is the following
381
-bit number:

p = 4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129015664037894272559787

The number of points

#E(Fp)=p+1t can be factorized as follows:

3 * 11^2 * 10177^2 * 859267^2 * 52437899^2 * 52435875175126190479447740508185965837690552500527637822603658699938581184513

where t = x + 1 = -15132376222941642751 is known as the trace of Frobenius and x is the BLS-family parameter that generates this specific curve.

Therefore, its

255-bit prime factor:

r = 52435875175126190479447740508185965837690552500527637822603658699938581184513

is the only one with sufficiently size for cryptographic purposes.

The point

PE(Fp) with coordinates:

x = 3685416753713387016781088315183077757961620795782546409894578378688607592378376318836054947676345821548104185464507
y = 1339506544944476473020471379941921221584933875938349620426543736416511423956333506472724655353366534992391756441569

is a point of order

r, i.e. a generator of the subgroup of size
r
of the curve.

The embedding degree

k of this curve is
12
.

Tower of Extension Fields

To represent

Fp12 we use the following tower of extension fields:
fields:
Fp2=Fp[u]/(u2β),Fp6=Fp2[v]/(v3ξ),Fp12=Fp6[w]/(w2v),

where
β
is not a quadratic residue in
Fp
and
ξ
is neither a quadratic residue or a cubic residue in
Fp2
.

In the particular case of the BLS12-381, we have

β=1 and
ξ=1+u
and therefore:
Fp2=Fp[u]/(u2+1),Fp6=Fp2[v]/(v3(1+u)),Fp12=Fp6[w]/(w2v),

from which we can extract the identities

u2=1,
w2=v
,
v3=1+u
and consequently the idenity
w6=1+u
.

Using the latter identity, we can also write:

Fp12=Fp2[w]/(w6(1+u))

These two representations of

Fp12 will allow us to apply some optimizations:

  1. Writing any element
    fFp12
    as
    f=g+hw
    , with
    g,hFp6
    implies
    fp6=f¯=ghw
    . Thus, computing the
    p6
    -th power of any element in
    Fp12
    is computationally "free".
  2. Now, if we write
    g=g0+g1v+g2v2
    and
    h=h0+h1v+h2v2
    where
    g0
    ,
    g1
    ,
    g2
    ,
    h0
    ,
    h1
    ,
    h2Fp2
    we have:
    f=g0+h0w+g1w2+h1w3+g2w4+h2w5,

    which will allow us to compute
    fp
    and
    fp2
    (needed for the final exponentiation) very efficiently.

Twist

BLS curves admits a sextic twist

E defined over
Fp2
of the form:
y2=x3+4·(1+u).

This means that pairing computations can be restricted to points

PE(Fp) and
QE(Fp2)
and avoid computations over
Fp12
entirely.

The point

PE(Fp2) with coordinates:

x = 352701069587466618187139116011060144890029952792775240219908644239793785735715026873347600343865175952761926303160 + 3059144344244213709971259814753781636986470325476647558659373206291635324768958432433509563104347017837885763365758·u
y = 1985150602287291935568054521177171638300868978215655730859378665066344726373823718423869104263333984641494340347905 + 927553665492332455747201965776037880757740193453592970025027978793976877002675564980949289727957565575433344219582·u

is a point of order

r, i.e. a generator of the subgroup of size
r
of the curve.

The following mapping, known as the untwist:

φ(x,y)=(x/w2,y/w3)
with:

1/w^2 = (2001204777610833696708894912867952078278441409969503942666029068062015825245418932221343814564507832018947136279894 + 2001204777610833696708894912867952078278441409969503942666029068062015825245418932221343814564507832018947136279893·u)·w^4
1/w^3 = (2001204777610833696708894912867952078278441409969503942666029068062015825245418932221343814564507832018947136279894 + 2001204777610833696708894912867952078278441409969503942666029068062015825245418932221343814564507832018947136279893·u)·w^3

defines a group isomorphism from

E(Fp2)[r] to
E(Fp12)[r]Ker(πp[p])
and we are going to use it to convert points from the twisted curve
E
to the curve
E
only whenever required.

The inverse mapping, the twist, is defined as:

φ1(x,y)=(xw2,yw3).

The Pairing

The optimal Ate pairing

e:G1×G2GT over the BLS12-381 curve is defined by:
e(P,Q)=f|x|,Q(P)p121r

where:

  • G1=E(Fp)[r] is the only subgroup of the
    r
    -torsion of
    E
    over
    Fp
    .

  • G2=E(Fp12)[r]Ker(πp[p]), represented for point arithmetic efficiency by
    E(Fp2)[r]
    . The former is the trace zero subgroup of the
    r
    -torsion of
    E
    over
    Fp12
    , whereas the latter is the only subgroup of the
    r
    -torsion of
    E
    over
    Fp2
    .

  • GT=μr is the set of
    r
    -th roots of unity over
    Fp12
    .

  • x = -15132376222941642752, which is a number of

    64 bits that can be expressed in binary as:

    ​​2^63 + 2^62 + 2^60 + 2^57 + 2^48 + 2^16
    ​​[1,1,0,1,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
    ​​0b1101001000000001000000000000000000000000000000010000000000000000
    

    and has a Hamming weight of 6.

  • fi,Q, for any
    iN
    and
    QG2
    , is a
    Fp12
    rational function on
    E
    that is computed using Miller's "double-and-add"-style algorithm:
    fi+j,Q=fi,Qfj,Q[i]φ(Q),[j]φ(Q)v[i+j]φ(Q),fi+1,Q=fm,Q[i]φ(Q),φ(Q)v[i+1]φ(Q),

    starting with
    f1,Q=1
    .

The conjugate of an element in

f=a+bw in
Fp12
is
f¯=abw
. Hence, if we write
f=g0+h0w+g1w2+h1w3+g2w4+h2w5
we have:
f¯=g0h0w+g1w2h1w3+g2w4h2w5.

Therefore, what we want to compute is the following:

def e(P,Q): assert(subgroup_check_G1(P)) assert(subgroup_check_G2(Q)) if ((P == 0) or (Q == 0)) return 1 # Here, 0 represents the # point at infinity R = Q f = 1 for i in range(len(x)-2,-1,-1): f = f * f * line(untwist(R),untwist(R),P) R = 2 * R if x[i] == 1: f = f * line(untwist(R),untwist(Q),P) R = R + Q return final_exponentiation(conjugate(f))

In what follows we explain how to compute each of the components of the previous pseudocode snippet.

Subgroup Checks

To summary 2019/814, we have:

  • A point

    PE(Fp) belongs to
    G1=E(Fp)[r]
    if and only if:
    [x213]([2]σ(P)Pσ2(P))=σ2(P)

    where
    σ:E(Fp)E(Fp)
    is the well-known endomorphism[1] defined by
    (a,b)(γ·a,b)
    , where
    γFp
    is of order
    3
    .

    ​​​​𝛄 = 4002409555221667392624310435006688643935503118305586438271171395842971157480381377015405980053539358417135540939436
    

    This is clearly much cheaper than the naive scalar multiplication by

    r, since its main cost is the scalar multiplication by
    (x21)/3
    , which is half the size of
    r
    and has low Hamming weight.

  • A point

    QE(Fp2) belongs to
    G2=E(Fp2)[r]
    if and only if:
    [x]ψ3(Q)+Q=ψ2(Q)

    where
    ψ:E(Fp2)E(Fp2)
    is the untwist-Frobenius-twist endomorphism, i.e.,
    ψ:=φ1πpφ
    .

    Again, this is much cheaper than the naive scalar multiplication by

    r, since its main cost is the scalar multiplication by
    x
    , which is roughly a fourth the size of
    r
    .

Line Equations

Different Points

Given two points

Q1=(x1,y1) and
Q2=(x2,y2)
in the twisted curve
E(Fp2)
such that
x1y1
and
x2y2
and a point
P=(x,y)
in the curve
E(Fp)
, the evaluation at
P
of the line passing through
φ(Q1)
and
φ(Q2)
is given by:
φ(Q1),φ(Q2)(P)=(x1y2x2y1)+w2(y1y2)x+w3(x2x1)y

The Same Points

Given a point

Q=(x,y) in the twisted curve
E(Fp2)
and a point
P=(x,y)
in the curve
E(Fp)
, the evaluation at
P
of the line passing through
φ(Q)
is given by:
φ(Q),φ(Q)(P)=w(3x32y2)+w3(3xx2)+w4(2yy)

Final Exponentiation

In the final exponentiation, we raise an element

fFp12 to the power
(p121)/r
.

First, divide the exponent as follows:

p121r=(p61)(p2+1)easy partp4p2+1rhard part

The Easy Part

Let's start with the easy part of the easy part, when we first attempt to compute

fp61. But since,
fp6=f¯
, we have
fp61=f¯f1
, just one conjugate, one inversion and one multiplication in
Fp12
.

To finish the easy part, we first compute

(f(p61))p2 which is the Frobenius operator applied to
fp61
and then multiply the result by
fp61
to finally obtain
m:=f(p61)(p2+1)Fp12
.

Importantly, after the computation of the easy part, the resulting field element becomes a member of the cyclotomic subgroup

Gϕ6(p2)=Gϕ12(p)={aFp12aϕ12(p)=ap4p2+1=ap6+1=1}, which means that inversion from now on can be computed by simply conjugation. Moreover, when an element
a
belongs to
Gϕ6(p2)
, one can compute
a2
using fewer multiplications and therefore the overall costs of performing exponentiations (see this section).

The Hard Part

Following Section 5 of this 2020/875, we can write:

p4p2+1r=(|x|+1)2(p|x|)(x2+p21)3+1.

To compute it, we proceed as follows:

def hard_part(m): x = abs(x) # 1] m^{(x+1)/3} y1 = m**((x+1)/3) # 2] m^{(x+1)^2/3} y2 = y1**(x+1) # 3.1] m^{(x+1)^2/3*-x} y31 = conjugate(y2)**x # 3.2] m^{(x+1)^2/3*p} y32 = Frobenius1(y2) # 4] m^{(x+1)^2/3*(p-x)} y4 = y31*y32 # 5.1] m^{(x+1)^2/3*(p-x)*x^2} y51 = (y4**x)**x # 5.2] m^{(x+1)^2/3*(p-x)*p^2} y52 = Frobenius2(y4) # 5.3] m^{(x+1)^2/3*(p-x)*-1} y53 = conjugate(y4) # 6] m^{(x+1)^2/3*(p-x)*(x^2+p^2)} y6 = y51*y52 # 7] m^{(x+1)^2/3*(p-x)*(x^2+p^2-1)} y7 = y6*y53 # 8] m^{(x+1)^2/3*(p-x)*(x^2+p^2-1)+1} y8 = y7*m return y8

(recall that, at this point, computing the conjugate is the same as computing the inverse)

Frobenius Operator

As we have shown, in the final exponentation we must compute the first and second Frobenius operator

πp,πp2 over an element in
Fp12
. Recalling that
πpi(x)=xpi
, a naive implementation would take
O(logp)
Fp12
-operations, but we will do it much better.

Recall that we can write any element

fFp12 as:
f=g+hw=g0+h0w+g1w2+h1w3+g2w4+h2w5

where
g=g0+g1v+g2v2
and
h=h0+h1v+h2v2
with
g,hFp6
and
gi,hiFp2
.

Then, we have:

fp=g¯0+h¯0γ1,1w+g¯1γ1,2w2+h¯1γ1,3w3+g¯2γ1,4w4+h¯2γ1,5w5,fp2=g0+h0γ2,1w+g1γ2,2w2+h1γ2,3w3+g2γ2,4w4+h2γ2,5w5,
with
γ1,i=(1+u)ip16
and
γ2,i=(1+u)ip216=γ1,iγ¯1,i
.

Expand to see the proof

We will now show that raising an

Fp2-element over any power of
p
is almost cost free and therefore so is raising
f
.

  • Take
    b=b0+b1uFp2
    with
    b0,b1Fp
    . Then,
    bp2i=b0+b1(u(p2)i)=b0+b1u=b
    from Lagrange theorem. Similarly,
    bp2i1=bp2ibp2bp=bbbp=1bp=bp=b¯
    .
  • Notice
    wp=(w6)p16w=(1+u)p16w
    and therefore
    (wi)p=γ1,iwi
    , where
    γ1,i=(1+u)ip16
    .

Using all the previous observations, we have:

fp=(g0+h0w+g1w2+h1w3+g2w4+h2w5)p=g¯0+h¯0wp+g¯1w2p+h¯1w3p+g¯2w4p+h¯2w5=g¯0+h¯0γ1,1w+g¯1γ1,2w2+h¯1γ1,3w3+g¯2γ1,4w4+h¯2γ1,5w5,
which can be computed using a few multiplications and a few conjugations, assuming each
γ1,i
is pre-computed for
i[5]
.

Similarly, notice

wp2=(1+u)p216w=(γ1,1)1+pw=γ1,1γ¯1,1w and therefore
(wi)p2=γ2,iwi
, where
γ2,i=γ1,iγ¯1,i
.

Hence, assuming the precomputation of

γ1,i,γ2,i for
i[5]
, we can compute
fp
and
fp2
almost for free.

Speeding up Exponentiations in the Cyclotomic Subgroup

I follow 2010/542 closely for this section.

Recall that:

Fp2=Fp[u]/(u2+1),Fp6=Fp2[v]/(v3(1+u)),Fp12=Fp6[w]/(w2v).
and also:
Fp12=Fp2[w]/(w6(1+u))

For this optimitzation, it will be useful to introduce

Fp4 and represent
Fp12
as a degree-3 extension of
Fp4
:
Fp4=Fp2[V]/(V2(1+u)),Fp12=Fp4[w]/(w3V),

where
V=w3
; so that an element
f
in
Fp12
can be written as
f=(g0+g1V)+(g2+g3V)w+(g4+g5V)w2
. We can equivalently write
f=g0+g2w+g4w2+g1w3+g3w4+g5w5
.

These optimizations are due to "compression" and "decompression" techniques, which exploits the algebraic relationships arising from elements belonging to the cyclotomic subgroup

Gϕ6(p2).

Compression and Decompression

If

f=g0+g2w+g4w2+g1w3+g3w4+g5w5Gϕ6(p2) then the compression function
C
is:
C(f)=[g2,g3,g4,g5].

and the decompression function
D
is:
D([g~2,g~3,g~4,g~5])=g~0+g~2w+g~4w2+g~1w3+g~3w4+g~5w5,

where:
{g~1=g~52(1+u)+3g~422g~34g~2,g~0=(2g~12+g~2g~53g~3g~4)(1+u)+1,if g~20g~1=2g~4g~5g~3,g~0=(2g~123g~3g~4)(1+u)+1,if g~2=0

Compressed Squaring

Now, on input

fGϕ6(p2) we will be able to compute its square on compressed form. That is,
C(f2)=[h2,h3,h4,h5]
, where
h=f2=h0+h2w+h4w2+h1w3+h3w4+h5w5
:
h2=2(g2+3(1+u)B4,5),h3=3(A4,5(2+u)B4,5)2g3,h4=3(A2,3(2+u)B2,3)2g4,h5=2(g5+3B2,3),Ai,j=(gi+gj)(gi+(1+u)gj),Bi,j=gigj.

Exponentiation in the Cyclotomic Subgroup

Given

fGϕ6(p2)Fp12 and
e1
, the idea is to compute
fe
using the squaring formula in compressed form.

Represent

e in binary as
e=e1e2e1e0
with
e1=1
and define
He={i:1i1 and ei=1}
. Then:
fe=i=01fei2i=i=01(f2i)ei=fe0iHeD(C(f)2i).

Therefore, we can use the square-and-multiply algorithm to compute such product. Since we are going to be squaring in compressed form, we start at

C(f) rather than the standard
f
in such algorithm.

def exp_cyclo(f,e): Cf = comp(f); result = 1; for i in range(len(e)): if (e[i] == 1) { result = decomp(Cf) * result; } Cf = square_comp(Cf); } return result;

Resources


  1. The endomorphism

    σ acts as a multiplication map
    [λ]P
    , where
    λ=t21
    . ↩︎