# BN128 curve --- C implementation
16, September, 2019
## Definition:
### BN
- Barreto, Naehrig et al. [Pairing-Friendly Elliptic Curves of Prime Order](https://cryptojedi.org/papers/pfcpo.pdf)
- 提出 embedding degree <font color=red>$k=12$</font>，存在有效率的演算法架構出橢圓曲線
- 其中：
- prime $p=36x^4+36x^3+24x^2+6x+1$ (架構出質數體的質數)
- order $n=36x^4+36x^3+18x^2+6x+1$ (橢圓曲線中的點個數)
- trace $t=6x^2+1$ (trace of the Frobenius)
- $n = p+ 1−t$
- 在Ethereum的`alt_bn128`曲線中
<font color=red>$x=4965661367192848881$</font>
(用[quickmath-solve polynomial](http://www.quickmath.com/webMathematica3/quickmath/equations/solve/intermediate.jsp#c=solve_basicsolveequation&v1=36z%255E4%2520%2B%252036z%255E3%2520%2B%252024z%255E2%2520%2B%25206z%2520%2B%25201%253D21888242871839275222246405745257275088696311157297823662689037894645226208583&v2=z)得出的解)
- 因此：
- $p=21888242871839275222246405745257275088696311157297823662689037894645226208583$
- $n=21888242871839275222246405745257275088548364400416034343698204186575808495617$
- $t=147946756881789318990833708069417712967$
- [*sage implementation*](https://cocalc.com):
```python=
sgae: x = 4965661367192848881
sgae: p = 36*x^4+36*x^3+24*x^2+6*x+1
sgae: n = 36*x^4+36*x^3+18*x^2+6*x+1
sage: t = 6*x^2+1
```
### Field
- [herumi/ate-pairing](https://github.com/herumi/ate-pairing)
- `include/bn.h`中的`BN_SUPPORT_SNARK`部分：
- $\mathbb{F}_{p^2} = \mathbb{F}_{p}[u]\ /\ (u^2 + 1)$
- $\mathbb{F}_{p^6} = \mathbb{F}_{p^2}[v]\ /\ (v^3 - \xi)\text{ where }\xi=u+9$
- $\mathbb{F}_{p^{12}} = \mathbb{F}_{p^6}[w]\ /\ (w^2 - v)$
- 根據定義：
- $u^2=1$
- $v^3=\xi=u+9$
- $w^2=v$
- *sage implementatoin*:
```python=
# GF(p) p的質數體，x為generator
sage: P.<x> = PolynomialRing(GF(p))
# 用GF(p) extension 建構Fp2，u為generator
sage: F2.<u> = GF(p).extension(x^2 + 1)
# Fp2的Polynomial ring P，t為generator
sage: P.<t> = F2[]
# 用Fp2 extension 建構Fp6，v為generator
sage: F6.<v> = F2.extension(t^3 - u-9)
# 若可以則執行下列：
# Fp6的Polynomial Ring P，y為generator
sage: P.<y> = F6[]
# 用Fp6 extension 建構Fp12，w為generator
sage: F12.<w> = F6.extension(y^2 - v)
```
> extension只能用同field的元素
> 如構成 $\mathbb{F}_{p^6}$ 的有: $\mathbb{F}_{p}$ extend而來的生成元 `a` 及 $\mathbb{F}_{p^2}$ 生成元`t`
> 構成 $\mathbb{F}_{p^12}$ 的有: $\mathbb{F}_{p^2}$ extend而來的生成元 `b` 及 $\mathbb{F}_{p^6}$ 生成元`y`
### Elliptic Curve
- [Ethereum Yellow Paper](https://ethereum.github.io/yellowpaper/paper.pdf) P.21
- curve: $y^2=x^3+3$
- $C_1\equiv \{(X,Y)\in F_p\times F_p | Y^2=X^3+3\}\bigcup\{(0,0)\}$
$P_1\equiv(1,2)\text{ on }C_1$
- $F_{p^2}$ be a field $F_p[i]\ /\ (i^2+1)$
$C_2\equiv \{(X,Y)\in F_{p^2}\times F_{p^2} | Y^2=X^3+3(i+9)^{-1}\}\bigcup\{(0,0)\}$
$P_2\equiv(11559732032986387107991004021392285783925812861821192530917403151452391805634\times i$
$+10857046999023057135944570762232829481370756359578518086990519993285655852781,$
$4082367875863433681332203403145435568316851327593401208105741076214120093531 \times i$
$+8495653923123431417604973247489272438418190587263600148770280649306958101930)$
- *sage implementation*:
```python=
# G1
sage: F1 = GF(21888242871839275222246405745257275088696311157297823662689037894645226208583)
sage: G1 = EllipticCurve(F1,[0,3])
sage: P1 = G1(1,2)
# G2
sage: F2 = GF(21888242871839275222246405745257275088696311157297823662689037894645226208583^2,"i",modulus=x^2 + 1)
sage: TwistB = 3*F2("9+i")^(-1)
sage: G2 = EllipticCurve(F2,[0,TwistB])
sage: P2x = F2("11559732032986387107991004021392285783925812861821192530917403151452391805634*i + 10857046999023057135944570762232829481370756359578518086990519993285655852781")
sage: P2y = F2("4082367875863433681332203403145435568316851327593401208105741076214120093531*i + 8495653923123431417604973247489272438418190587263600148770280649306958101930")
sage: P2 = G2(P2x,P2y)
```
### Pairing
- [EIP-197](https://github.com/ethereum/EIPs/blob/master/EIPS/eip-197.md)
- 定義：
```
Input: (a1, b1, a2, b2, ..., ak, bk) from (G_1 x G_2)^k
Output: If the length of the input is incorrect or any of the inputs are not elements of
the respective group or are not encoded correctly, the call fails.
Otherwise, return one if
log_P1(a1) * log_P2(b1) + ... + log_P1(ak) * log_P2(bk) = 0
(in F_q) and zero else.
```
- 因此目標為計算 $k$ 組 $G_1,G_2$ 點的pairing
- 目標： <font color=red>$e(a_1, b_1)\ *\ ...\ *\ e(a_k, b_k) = 1$</font>
- 其中 $log_P(x)$ 為滿足 $n\cdot P=x$ 的最小 $n$
- $log_{P_1}(a_1) * log_{P_2}(b_1) + ... + log_{P_1}(a_k) * log_{P_2}(b_k) = 0 \text{ (in }\mathbb{F}_q)$
$e(P_1, P_2)^{(log_{P_1}(a_1) * log_{P_2}(b_1) + ... + log_{P_1}(a_k) * log_{P_2}(b_k))} = 1 \text{ (in }\mathbb{G}_T)$
$=e(log_{P_1}(a1) * P_1, log_{P_2}(b_1) * P_2) * ... * e(log_{P_1}(a_k) * P_1, log_{P_2}(b_k) * P_2))$
$=e(a_1, b_1)\ \cdot\ ...\ \cdot\ e(a_k, b_k)$
## Packages
### [`ate-pairing`](https://github.com/herumi/ate-pairing)
- C++
- 建構[`libff`](https://github.com/scipr-lab/libff), [`libsnark`](https://github.com/scipr-lab/libsnark)的基礎
- default: $x = -(2^{62} + 2^{55} + 1)$
- SNARK: `$ make -j SUPPORT_SNARK=1`
- optimal ate-pairing
- loop count為 $29793968203157093288=(6z+2)$
- [sample code](https://github.com/herumi/ate-pairing/blob/master/test/sample.cpp)
- print出結果：
```cpp=
Fp12 e1;
opt_atePairing(e1, b1, a1);
std::cout<<"e1:\n"<<std::hex<<e1;
// std::hex 將輸出轉成10進位
```
### [`py_ecc`](https://github.com/ethereum/py_ecc)
- Python
- field及elliptic curve都如同ethereum用的bn128
- optimal ate-pairing
- sample:
```python
from py_ecc.bn128 import *
p = curve_order
x = randint(1, p-1) # out secret key
H_m = multiply(G1, randint(1, p-1)) # lets pretend it's HashToPoint
P = multiply(G2, x) # our public key in G2
S = multiply(H_m, x) # our signature in G1
a = pairing(P, H_m)
b = pairing(G2, S)
assert a == b # Verify signature
```
### [`bn256`](https://www.cryptojedi.org/crypto/index.shtml#bncurves)
- C
- download:
```bash=
wget http://cryptojedi.org/crypto/data/bn256-20090824.tar.bz2
tar xjvf bn256-20090824.tar.bz2
cd bn256-20090824
make
```
- default: bn256 curve
- ate-pairing,optimal ate-pairing, tate-pairing, eta-pairing
- compile: `gcc bn256.c -lgmp`
- **用此library改驗證contract**
### [`pbc`](https://crypto.stanford.edu/pbc/)
- C
- BN curve defined in [Type F pairing](https://crypto.stanford.edu/pbc/manual/ch08s08.html)
(參數`beta`,`alpha`尚未確認為何)
- sample:
```cpp=
$ pbc/pbc
$ init_pairing_f();
$ g := rnd(G1);
$ g;
[79852913720140618033238675189472612021317443593, 108845852893881268898046458273434478449744261041]
$ h := rnd(G2);
$ h;
[[66738494328564557862145166514672743296967736031, 138628962866914872444577256854665138260901670758], [44108488804575150857180778316233085344147224102, 73374821578866986221123150061800643945355771561]]
$ pairing(g,h);
[[127244296846151320229451765231937407345264275337, 194516214605916518197434454512307639138513347685], [72123392631556257776993010580007966809196163433, 103133235071738349318869312483762541080166756669], [115470402213399488913960527477665437360341563634, 200045067655991734256248847368109781819590012108], [165249791411356411563096446612379783127179161114, 90879321830240188846380135482009532401892893033], [14435220015022740410783909886514751711701242926, 60515085796840445289165108813993241043979134259], [49371634445738451451225200934854581888610639833, 127714555993719752867082276854684686305396949588]]
$ a := rnd(Zr);
$ b := rnd(Zr);
$ pairing(g^a,h^b);
[[88424061323103803164426643017917154342437690388, 17300053332458983889698342671869711927095355785], [105794066241143284372826714937348966757146275247, 98584892219681632258828536835126690818980273309], [183812827043409431154961438096794763381657047935, 5296345413835424569297807155298296475996049089], [119183393812191428328100091251793077351089043188, 22627592276984654933875258547844800793012505017], [107533753564271022761938415005984406018901618293, 149300949962179025507296347427720414700427056977], [154979478921340615280616064447497122078583184353,159326123140227494966318218788156290320376454546]]
$ pairing(g,h)^(a*b);
[[88424061323103803164426643017917154342437690388, 17300053332458983889698342671869711927095355785], [105794066241143284372826714937348966757146275247, 98584892219681632258828536835126690818980273309], [183812827043409431154961438096794763381657047935, 5296345413835424569297807155298296475996049089], [119183393812191428328100091251793077351089043188, 22627592276984654933875258547844800793012505017], [107533753564271022761938415005984406018901618293, 149300949962179025507296347427720414700427056977], [154979478921340615280616064447497122078583184353, 159326123140227494966318218788156290320376454546]]
```
## Pairing
- $$e:\mathbb{G}_2\times\mathbb{G}_1\longrightarrow\mathbb{G}_T$$
- $\mathbb{G}_T$ 即上述架構的 $F_{p^{12}}$
- 性質：
- $e(P, Q + R) = e(P, Q) * e(P, R)$
- $e(P + S, Q) = e(P, Q) * e(S, Q)$
- $e(a\cdot P,Q)=e(P,Q)^a$
- $e(a\cdot P,b\cdot Q)=e(P,Q)^{ab}$
- Reference:
- [The Eta Pairing Revisited](https://eprint.iacr.org/2006/110.pdf) -- Tate-pairing, Ate-pairing
- [High-Speed Software Implementation of the Optimal Ate Pairing over Barreto–Naehrig Curves](https://eprint.iacr.org/2010/354.pdf) -- Opt Ate-pairing
### Tate-pairing
- $$(Q,P)\longmapsto f_{r,Q}(P)^{\frac{q^k-1}{r}}$$
### Weil-pairing
- $$w_r(Q,P)=\frac{f_{r,Q}{(P)}}{f_{r,P}{(Q)}}$$
### Ate-pairing
- $$(Q,P)\longmapsto f_{t-1,Q}(P)^{\frac{q^k-1}{r}}$$
### Optimal Ate-pairing
- $$(Q,P)\longmapsto(f_{6t+2,Q}(P)\cdot l_{[6t+2]Q,\pi_p(Q)}(P)\cdot l_{[6t+2]Q+\pi_p(Q),-\pi_p^2(Q)}(P))^{\frac{q^k-1}{r}}$$
## Algorithm
### Frobenius
- 定義: Raising elements of $\mathbb{F}_{p^e}$ to there $p$-th powers gives a self-map of $\mathbb{F}_{p^e}$; it is called the Frobenius map and written as $F$.
$$F:\mathbb{F}_{p^e}→\mathbb{F}_{p^e},x↦x^p$$
- trace of Frobenius of the curve
```sage
sage: E=EllipticCurve(GF(101),[2,3])
sage: E.trace_of_frobenius()
6
```
- `fp6e_frobenius_p`:
- algorithm: [High-Speed Software Implementation of the Optimal Ate Pairing over Barreto–Naehrig Curves](https://eprint.iacr.org/2010/354.pdf)
- $\textbf{Algorithm 28}\text{ :Frobenius raised to }p\text{ of }f\in\mathbb{F}_{p^{12}}=\mathbb{F}_{p^6}[w]/(w^2-\gamma)$
- sage code:
```python=
# point1 is in F6
sage: F6.frobenius_endomorphism()(point1)
sage: F6.frobenius_endomorphism(1)(point1) # equal
# return another point in F6
```
- `fp6e_frobenius_p2`:
- algorithm:
- $\textbf{Algorithm 29}\text{ :Frobenius raised to }p^2\textit{ of }f\in\mathbb{F}_{p^{12}}=\mathbb{F}_{p^6}[w]/(w^2-\gamma)$
- sage code:
```python=
# point1 is in F6
sage: F6.frobenius_endomorphism(2)(point1)
# return another point in F6
```
- `fp6e_frobenius_p3`:
- 雖然有演算法，但也可用`frobenius_p1`和`frobenius_p2`建構出`frobenius_p3`
- rule: $f^{p^3}(x)=f^{p}(f^{p^2}(x))$
- sage code:
```python=
sage: f61 = F6.frobenius_endomorphism()(point1)
sage: f62 = F6.frobenius_endomorphism(2)(point1)
sage: f631 = F6.frobenius_endomorphism(3)(point1)
sage: f632 = F6.frobenius_endomorphism()(f62)
sage: f631 == f632
True
```
### Miller's Loop
- [$\text{Miller's algorithm:}$ (P.3-6)](https://www.math.u-bordeaux.fr/~damienrobert/csi2017/pairings.pdf)
- **Input**: $r\in\mathbb{N},I=[log r],P=(x_P,y_P)\in E[r](K),Q=(x_Q,y_Q)\in E(K)$
- **Output:** $f_{r,P}(Q)$
- 1. Compute the binary decomposition: $r:=\sum_{i=0}^Ib_i2^i$. Let $T=P,f=1$
- For $i$ in $[I-1..0]$:
(a) $f=f^2\mu_{T,T}(Q)$
(b) $T=2T$
(c ) if $b_i=1$:
i. $f=f\mu_{T,P}(Q)$
ii. $T=T+P$
- return $f$
- py_ecc code:
```python=85
def miller_loop(Q: Point2D[FQ12],
P: Point2D[FQ12]) -> FQ12:
...
R = Q # type: Point2D[FQ12]
f = FQ12.one()
for i in range(log_ate_loop_count, -1, -1):
f = f * f * linefunc(R, R, P)
R = double(R)
if ate_loop_count & (2**i):
f = f * linefunc(R, Q, P)
R = add(R, Q)
...
```
### Line function
- Linear function $l_{Q,Q}(P)\in \mathbb{F}_{p^{12}}$
- algorithm: [High-Speed Software Implementation of the Optimal Ate Pairing over Barreto–Naehrig Curves P. 28 $\text{Algorithm 26, 27}$](https://eprint.iacr.org/2010/354.pdf)
- *sage implementation (bn256)*:
```python=
# Algorithm 26 Point doubling and line evaluation
x_p = F1("19036326650351419380309844008588046887500192473820022916684059283212312810755")
y_p = F1("24416690099424885538680537874355857700418235418382750960697588420893490979447")
X_Q = F2("66020635303426444812813329252614299010925071937759307499201344123165041483630 * a + 19412881731638021337230777338399549073276632608516152708419727813812179213342")
Y_Q = F2("9496995991393763151014644875678504226651962851610752566278875877566242661313 * a + 64547729877332878386921358349209593549688364540149661508455335221364968200136")
Z_Q = F2("28265609773395302082323410388877418050796517628590810004689763505562099393371 * a + 11608911084161395568477488481039820448251922441323271942259565147644676056919")
tmp0 = X_Q^2
tmp1 = Y_Q^2
tmp2 = tmp1^2
tmp3 = (tmp1+X_Q)^2-tmp0-tmp2
tmp3 = 2*tmp3
tmp4 = 3*tmp0
tmp6 = X_Q +tmp4
tmp5 = tmp4^2
X_T = tmp5-2*tmp3
Z_T = (Y_Q+Z_Q)^2-tmp1-Z_Q^2
Y_T = (tmp3-X_T)*tmp4-8*tmp2
tmp3 = -2*(tmp4*Z_Q^2)
tmp3 = tmp3*x_p
tmp6 = tmp6^2-tmp0-tmp5-4*tmp1
tmp0 = 2*(Z_T*Z_Q^2)
tmp0 = tmp0*y_p
# return ((tmp6*b)+tmp3)*c+tmp0 and (X_T,Y_T,Z_T)
# tmp6 = 39467808497513321971418495975092757910573795172620731653284844774590128544949*a + 16064578912625871412495875021273158087538297781312732719918703103255685212978
# tmp3 = 26900272025442231202136810059174826466870068404883709542831107009130992101024*a + 36820176082059830702390334368466688997958373791692089380588742415783206487213
# tmp0 = 42540957778736336440218600222726529148774731933627463366204312833488321059843*a + 17279992842561602300512709260017417861197710148394396110364376447336168241241
# X_T = 59939262272742075282293750272562916432861452573095215973543955503737179845918*a + 71849103101821744092544346176918665246275637177231140546732043562364871274010
# Y_T = 55118971710910072378614458628513479971221383112282002676433416805306813327600*a + 22474579714957608756815244332228042411943570322630763049703532407058892183412
# Z_T = 58472854990350142177940989350497195394893722422701615335663755883768116801021*a + 15698510426968927131396167416642403602632680084454872605963655368201407093488
```
```python=
# Algorithm 27 Point addition and line evaluation
x_p = F1("19036326650351419380309844008588046887500192473820022916684059283212312810755")
y_p = F1("24416690099424885538680537874355857700418235418382750960697588420893490979447")
X_R = F2("60468327060250697695777504843760602312778240813685687476754545290340736660117 * a + 75207288700540388172087626426077159416239342674220757084925862172259324791798");
Y_R = F2("11747159676870954034830968911985740355803730873092166329012919399968203758626 * a + 77437491087823200008673050455309591478143842462808890777525559574909739529050");
Z_R = F2("48197990245400830303642899484867393594767420425921999847413039690362552391954 * a + 5774166803647311766191305052111952966489183189940295269365880795628328435859");
X_Q = F2("621974397985525318194063125170198577902540556498747469157675653746576620880 * a + 6251805185683980258298929625871846172079392273578387836099945035602364239135")
Y_Q = F2("49990908316789968724860573235913273901071484323829192314337721068681689250452 * a + 19258290530482130955008764894053126575339703184772660063788583927864324464626")
Z_Q = F2("0 *a + 1")
t0 = X_Q*Z_R^2
t1 = (Y_Q+Z_R)^2-Y_Q^2-Z_R^2
t1 = t1*Z_R^2
t2 = t0-X_R
t3 = t2^2
t4 = 4*t3
t5 = t4*t2
t6 = t1-2*Y_R
t9 = t6*X_Q
t7 = X_R*t4
X_T = t6^2-t5-2*t7
Z_T = (Z_R+t2)^2-Z_R^2-t3
t10 = Y_Q+Z_T
t8 = (t7-X_T)*t6
t0 = 2*(Y_R*t5)
Y_T = t8-t0
t10 = t10^2-Y_Q^2-Z_T^2
t9 = 2*t9-t10
t10 = 2*(Z_T*y_p)
t6 = -1*t6
t1 = 2*(t6*x_p)
t9
t1
t10
# return ((t9*b)+t1)*c+t10 and (X_T,Y_T,Z_T)
# t9 = 76479445506672858954637202806298445894378948504602958825436274313311101562445*a + 64652353141305220762250925929760150731422897775111546669321115421764662312129
# t1 = 31824713071612796591909199573998992099967176584471593726257655335602275248320*a + 59393640435344032559925061282433789629607126201787190693026971066582909965704
# t10 = 6526929262762995979386081461373064697408294715716098178260463132377538058359*a + 70402627597488067633993521117962910027448315214449308123592677402636620499031
# X_T = 66020635303426444812813329252614299010925071937759307499201344123165041483630*a + 19412881731638021337230777338399549073276632608516152708419727813812179213342
# Y_T = 9496995991393763151014644875678504226651962851610752566278875877566242661313*a + 64547729877332878386921358349209593549688364540149661508455335221364968200136
# Z_T = 28265609773395302082323410388877418050796517628590810004689763505562099393371*a + 11608911084161395568477488481039820448251922441323271942259565147644676056919
```
### Final Exponential
- $f$ raised to $f^{\frac{p^{12}-1}{n}}$
- $\frac{p^{12}-1}{r}=(p^6-1)\cdot (p^2+1)\cdot \frac{p^4-p^2+1}{r}$
- [High-Speed Software Implementation of the Optimal Ate Pairing over Barreto–Naehrig Curves P. 32 $\text{Algorithm 31}$](https://eprint.iacr.org/2010/354.pdf)
- $f^{(p^6-1)}=\bar{f}\cdot f^{-1}$
- $\bar{f}$: 共軛
- $f^{-1}\cdot f=1$
- [Implementing Cryptographic Pairings over Barreto-Naehrig Curves P.7 $\text{Algorithm 3}$](https://eprint.iacr.org/2007/390.pdf)
- $p$ 和 $n$ 為 $x$ generate而來
- $\frac{p^4-p^2+1}{r}=p^3 + (6x^2 + 1)p^2 + (36x^3 − 18x^2 + 12x + 1)p + (36x^3 − 30x^2 + 18x − 2)$
## Parameters
### BN_P
- prime $p=36x^4+36x^3+24x^2+6x+1$
- ```python=
sgae: x = 4965661367192848881
sgae: p = 36*x^4+36*x^3+24*x^2+6*x+1
```
- `#define BN_P "21888242871839275222246405745257275088696311157297823662689037894645226208583"`
### BN_PINV32
- `p^{-1} mod 2^{GMP_LIMB_BITS}` used in Montgomery reduction
`GMP_LIMB_BITS = 32`
- ```python=
sgae: (-p)^(-1) % 2^(32)
```
- `#define BN_PINV32 3834012553UL`
### BN_PINV64
- `p^{-1} mod 2^{GMP_LIMB_BITS}` used in Montgomery reduction
`GMP_LIMB_BITS = 32`
- ```python=
sgae: (-p)^(-1) % 2^(64)
```
- `#define BN_PINV64 9786893198990664585UL`
### ALPHA
- $\alpha:\mathbb{F}_{p^2} = \mathbb{F}_{p}[u]\ /\ (u^2 - \alpha)$
- $\alpha=-1$ by definition
- `#define ALPHA (-1) // constant coefficient in the irreducible polynomial x^2 - alpha, used to construct F_{p^2}`
### BN_X
- $x=4965661367192848881$
(用[quickmath-solve polynomial](http://www.quickmath.com/webMathematica3/quickmath/equations/solve/intermediate.jsp#c=solve_basicsolveequation&v1=36z%255E4%2520%2B%252036z%255E3%2520%2B%252024z%255E2%2520%2B%25206z%2520%2B%25201%253D21888242871839275222246405745257275088696311157297823662689037894645226208583&v2=z)得出的解)
- `#define BN_X "4965661367192848881" // parameter x used to generate the curve (see "Pairing-Friendly Elliptic Curves of Prime Order")`
### BN_N
- order $p=36x^4+36x^3+18x^2+6x+1$
- ```python=
sgae: x = 4965661367192848881
sgae: p = 36*x^4+36*x^3+18*x^2+6*x+1
```
- `#define BN_N "21888242871839275222246405745257275088548364400416034343698204186575808495617" // prime order of E(F_p)`
### BN_TRACE
- trace $t=6x^2+1$ (trace of the Frobenius)
- ```python=
sage: G1.trace_of_frobenius()
147946756881789318990833708069417712967
```
- `#define BN_TRACE "147946756881789318990833708069417712967" // trace of Frobenius of the curve`
### BN_XI
- $\mathbb{F}_{p^6} = \mathbb{F}_{p^2}[v]\ /\ (v^3 - \xi)\text{ where }\xi=u+9$
- $\xi$ : Xi
- `#define BN_XI "1", "9"`
### BN_YPMINUS1
- $Y^{p-1}$ lies in $\mathbb{F}_{p^2}$
where $\mathbb{F}_{p^6}[Y]=\mathbb{F}_{p^2}[Y]/(Y^3-\xi)$
-
```python=
# construct Fp6 = Fp[v] / (v^3 - xi)
sage: P.<t> = F2[]
sage: F6.<b> = F2.extension(t^3 - a-9)
sage: b^(p-1)
10307601595873709700152284273816112264069230130616436755625194854815875713954*a + 21575463638280843010398324269430826099269044274347216827212613867836435027261
```
- `#define BN_YPMINUS1 "10307601595873709700152284273816112264069230130616436755625194854815875713954", "21575463638280843010398324269430826099269044274347216827212613867836435027261"`
> `b`是 $\mathbb{F}_{p^6}[Y]=\mathbb{F}_{p^2}[Y]/(Y^3-\xi)$ 裡的生成元，即 $Y$，因此取 $Y$ 的 $p-1$ 次方
> $\mathbb{F}_{p^6}$ 其實是由`a`的二次方+`b`的三次方組合而成，例如:
>
```sage
sage: for i in range(10):
sage: print b^i
1
b
b^2
a + 1
(a + 1)*b
(a + 1)*b^2
2*a + 82434016654300679721217353503190038836571781811386228921167322412819029493182
(2*a + 82434016654300679721217353503190038836571781811386228921167322412819029493182)*b
(2*a + 82434016654300679721217353503190038836571781811386228921167322412819029493182)*b^2
a + 82434016654300679721217353503190038836571781811386228921167322412819029493178
```
> 剛好 $(p-1)\ mod\ 3=0$，沒有`b`次方項只有`a`次方項
>
### BN_ZETA
- Third root of unity in $F_p$ fulfilling $Z^{p^2} = -\zeta * Z$
- $\zeta=Z^{p^2-1}$
- sage:
```python=
# 若能執行則執行：
sage: P.<y> = F6[]
sage: F12.<c> = F6.extension(y^2 - b)
sage: c^((p^2)-1)
```
- 無法用sage generate出F12，因此`c^((p^2)-1)`無法使用
- 由`extension(y^2-b)`得知：$y^2=b$
- 原來的 $c^{p^2-1}$ 可變成 $b^{\frac{p^2-1}{2}}$
- 可算出 $\zeta=-b^{\frac{p^2-1}{2}}$，為 $2203960485148121921418603742825762020974279258880205651966$
- `#define BN_ZETA "2203960485148121921418603742825762020974279258880205651966"`
### BN_TAU
- def: $\tau$ = `"0", "0", "0", "1", "0", "0" // constant tau used to construct F_p^12 as F_p^6[Z]/ (Z^2 - tau)`
- `fp6e_multau(a) = a*tau`
- where `tau=(0*a+0)*b^2+(0*a+1)*b+(0*a+0)`
$$(0*a+0)*b^2+(0*a+1)*b+(0*a+0)\\=b$$
- `a` is in $\mathbb{F}_{p^6}$
- sage code:
```python=
# point1 is in F6
pointx = (0*a+0)*b^2+(0*a+1)*b+(0*a+0)
point1 * pointx
# return fp6e_multau(point1)
```
- ref: [Multiplication and Squaring on Pairing-Friendly Fields](https://pdfs.semanticscholar.org/3e01/de88d7428076b2547b60072088507d881bf1.pdf)
- `#define BN_TAU "0", "0", "0", "1", "0", "0" // constant tau used to construct F_p^12 as F_p^6[Z]/ (Z^2 - tau)`
### BN_ZPMINUS1
- `BN_ZPMINUS1`: $c^{p-1}=b^{\frac{p-1}{2}}$
-
```python=
sage: # c^(p-1)
sage: # (p-1)/2 = 10944121435919637611123202872628637544348155578648911831344518947322613104291
sage: b^10944121435919637611123202872628637544348155578648911831344518947322613104291 # c^2 = b
```
- `#define BN_ZPMINUS1 "16469823323077808223889137241176536799009286646108169935659301613961712198316", "8376118865763821496583973867626364092589906065868298776909617916018768340080" // Z^(p-1)`
### BN_ZPMINUS1INV
- `BN_ZPMINUS1`: $c^{1-p}=b^{\frac{1-p}{2}}$
-
```python=
sage: # c^(1-p)
sage: # (1-p)/2 = -10944121435919637611123202872628637544348155578648911831344518947322613104291
sage: b^(-10944121435919637611123202872628637544348155578648911831344518947322613104291) # c^2 = b
```
- `#define BN_ZPMINUS1INV "5722266937896532885780051958958348231143373700109372999374820235121374419868", "18566938241244942414004596690298913868373833782006617400804628704885040364344" // Z^(1-p)`
### BN_CURVEGEN
- base point of $G_1$，定義在[ethereum yellow paper](https://ethereum.github.io/yellowpaper/paper.pdf)中
- $P_1=(1,2)$
- <font color=red>需轉換為**Jacobian coordinates**</font>
- Affine to Jacobian:
$(X,Y)\Rightarrow(X,Y,1)$
- Jacobian to Affine:
$(X,Y,Z)\Rightarrow(\frac{X}{Z^2},\frac{Y}{Z^3})$
- ref: [How can convert affine to Jacobian coordinates?](https://crypto.stackexchange.com/questions/19598/how-can-convert-affine-to-jacobian-coordinates)
- `#define BN_CURVEGEN "1", "2", "1"`
### BN_TWISTGEN_X
- base point of $G_2$，定義在[ethereum yellow paper](https://ethereum.github.io/yellowpaper/paper.pdf)中
- $P_2=(11559732032986387107991004021392285783925812861821192530917403151452391805634\times i$
$+10857046999023057135944570762232829481370756359578518086990519993285655852781,$
$4082367875863433681332203403145435568316851327593401208105741076214120093531 \times i$
$+8495653923123431417604973247489272438418190587263600148770280649306958101930)$
- `#define BN_TWISTGEN_X "11559732032986387107991004021392285783925812861821192530917403151452391805634", "10857046999023057135944570762232829481370756359578518086990519993285655852781"`
### BN_TWISTGEN_Y
- 同上
- `#define BN_TWISTGEN_Y "4082367875863433681332203403145435568316851327593401208105741076214120093531", "8495653923123431417604973247489272438418190587263600148770280649306958101930"`
### BN_ZETA2
- $\zeta^2$
- $\zeta=2203960485148121921418603742825762020974279258880205651966$
- ```python=
sage: P("2203960485148121921418603742825762020974279258880205651966^2")
```
- `#define BN_ZETA2 "21888242871839275220042445260109153167277707414472061641714758635765020556616" // zeta^2`
### BN_Z2P
- $c^{2p}=b^p$
- sage:
```python=
sage: print(" BN_Z2P // Z^(2p)")
sage: # c^(2p)
sage: b^21888242871839275222246405745257275088696311157297823662689037894645226208583
BN_Z2P // Z^(2p)
(10307601595873709700152284273816112264069230130616436755625194854815875713954*a + 21575463638280843010398324269430826099269044274347216827212613867836435027261)*b
```
- code:
```cpp=80
#define BN_Z2P "10307601595873709700152284273816112264069230130616436755625194854815875713954", "21575463638280843010398324269430826099269044274347216827212613867836435027261" // Z^(2p)
```
### BN_Z3P
- $c^{3p}=b^\frac{3p-1}{2}\times c$
- sgae:
```python=
sage: print(" BN_Z3P // Z^(3p)")
sage: # c^(3p)
sage: # ((3*p-1)/2)
sage: b^32832364307758912833369608617885912633044466735946735494033556841967839312874
BN_Z3P // Z^(3p)
(3505843767911556378687030309984248845540243509899259641013678093033130930403*a + 2821565182194536844548159561693502659359617185244120367078079554186484126554)*b
```
- code:
```cpp=81
#define BN_Z3P "3505843767911556378687030309984248845540243509899259641013678093033130930403", "2821565182194536844548159561693502659359617185244120367078079554186484126554" // Z^(3p)
```
## Notes
- `py_ecc`及`ate-pairing`是用<font color=red>**Opt ate-pairing**</font>，loop count為29793968203157093288(6z+2)
- 原來有`fp2e_mulxi`這個函數，但因此函數還沒實作 $\alpha$ 為其他數的可能性，因此用`fp2e_mul(rop,op,xi)`代替，效果一樣
- 注意哪些函數的input需要符合<font color=red>**affine coordinate**</font>(ex. `curvepoint_fp_mixadd`)
在呼叫函數前可以用`curvepoint_fp_makeaffine(op)` 轉成affine coordinate
- 在原verifier.sol有此行code:
```javascript=
Pairing.G1Point memory vk_x = Pairing.G1Point(0, 0);
for (uint i = 0; i < input.length; i++)
vk_x = Pairing.addition(vk_x, Pairing.scalar_mul(vk.IC[i + 1], input[i]));
```
$(0,0)$ 不在 $Y^2=X^3+3$ 上，為ethereum yellow paper定義的<font color=red>無窮原點</font>
在此code中若是無窮原點則jacobian的 $Z$ 座標為零(根據橢圓曲線定義$g_1^n=(0,0)，$用`curve_gen`的 $n$ 次方檢查無窮原點為何)
在此code`vk_x=(x,y,0)`在第一個loop不能用加法，會出錯。
改寫成
```cpp=
if(fpe_iszero(vk_x->m_z)){
curvepoint_fp_set(vk_x, tmp0);
}
else{
curvepoint_fp_mixadd(vk_x, vk_x, tmp0);
curvepoint_fp_makeaffine(vk_x);
}
```
若非零才加法，零則直接指定
## Appendix -- \<gmp\>
### \<gmp\> on macOS
0. Install
- [Install gmp on Mac OSX](http://macappstore.org/gmp/)
- `ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)" < /dev/null 2> /dev/null`
- `brew install gmp`
1. Compile
```
gcc code.c -L. -lpbc -lgmp
```
2. Execution
```
./a.out
```
```
Enter Elliptic Curve Parameters i.e. a,b and p
Enter Choice of Operation
Enter 1 For Point Addition Operation
Enter 2 For Scalar Multiplication Operation
1
Enter Points P(x,y) and/or Q(x,y) to beAdded
Resultant Point is 53245742458311851613963450768604695125869175507241963909214915539587838403999,86019658436333621282556383871045654836888327576374321095942401925234007301028
```
```
Enter Elliptic Curve Parameters i.e. a,b and p
Enter Choice of Operation
Enter 1 For Point Addition Operation
Enter 2 For Scalar Multiplication Operation
2
Enter Points P(x,y) and/or Q(x,y) to beAdded
Enter m to Find mP
1
Resultant Point is 55066263022277343669578718895168534326250603453777594175500187360389116729240,32670510020758816978083085130507043184471273380659243275938904335757337482424
```
- code
```cpp=
#include<stdio.h>
#include<stdlib.h>
#include<gmp.h>
struct Elliptic_Curve
{
mpz_t a;
mpz_t b;
mpz_t p;
};
struct Point
{
mpz_t x;
mpz_t y;
};
struct Elliptic_Curve EC;
void Select_EC();
void Point_Addition(struct Point P,struct Point Q,struct Point* R);
void Point_Doubling(struct Point P,struct Point *R);
void Scalar_Multiplication(struct Point P,struct Point* R, mpz_t m);
int main(void)
{
int choice;
mpz_init(EC.a); mpz_init(EC.b); mpz_init(EC.p);
Select_EC();
printf("\n Enter Choice of Operation\n");
printf("\n Enter 1 For Point Addition Operation\n");
printf("\n Enter 2 For Scalar Multiplication Operation\n");
scanf("%d",&choice);
struct Point P,R;
mpz_init(P.x);
mpz_init(P.y);
mpz_init_set_ui(R.x,0);
mpz_init_set_ui(R.y,0);
printf("\n Enter Points P(x,y) and/or Q(x,y) to beAdded\n");
mpz_init_set_str (P.x, "55066263022277343669578718895168534326250603453777594175500187360389116729240", 10);
mpz_init_set_str (P.y, "32670510020758816978083085130507043184471273380659243275938904335757337482424", 10);
//gmp_scanf("%Zd",&P.x);
//gmp_scanf("%Zd",&P.y);
if(choice==1)
{
struct Point Q;
mpz_init(Q.x);
mpz_init(Q.y);
mpz_init_set_str (Q.x, "48840125481190545212233038815866595080052446190700286140915905752173158787212", 10);
mpz_init_set_str (Q.y, "55258565891714835096390285893857201264507235115974921462591453996678450950230", 10);
// gmp_scanf("%Zd",&Q.x);
// gmp_scanf("%Zd",&Q.y);
Point_Addition(P,Q,&R);
}
else
{
printf("\n Enter m to Find mP\n");
mpz_t m;
mpz_init(m);
gmp_scanf("%Zd",&m);
Scalar_Multiplication(P,&R,m);
}
gmp_printf("\n Resultant Point is %Zd,%Zd",R.x,R.y);
}
void Select_EC()
{
printf("\n Enter Elliptic Curve Parameters i.e. a,b and p");
mpz_init_set_str (EC.a, "0", 10);
mpz_init_set_str (EC.b, "7", 10);
mpz_init_set_str (EC.p, "115792089237316195423570985008687907853269984665640564039457584007908834671663", 10);
//gmp_scanf("%Zd",&EC.a);
//gmp_scanf("%Zd",&EC.b);
//gmp_scanf("%Zd",&EC.p);
}
void Point_Addition(struct Point P,struct Point Q,struct Point* R)
{
mpz_mod(P.x,P.x,EC.p);
mpz_mod(P.y,P.y,EC.p);
mpz_mod(Q.x,Q.x,EC.p);
mpz_mod(Q.y,Q.y,EC.p);
mpz_t temp,slope;
mpz_init(temp);
mpz_init_set_ui(slope,0);
if(mpz_cmp_ui(P.x,0)==0 && mpz_cmp_ui(P.y,0)==0)
{ mpz_set(R->x,Q.x); mpz_set(R->y,Q.y); return;}
if(mpz_cmp_ui(Q.x,0)==0 && mpz_cmp_ui(Q.y,0)==0)
{ mpz_set(R->x,P.x); mpz_set(R->y,P.y); return;}
if(mpz_cmp_ui(Q.y,0)!=0)
{mpz_sub(temp,EC.p,Q.y);mpz_mod(temp,temp,EC
.p);}
else
mpz_set_ui(temp,0);
// gmp_printf("\n temp=%Zd\n",temp);
if(mpz_cmp(P.y,temp)==0 &&
mpz_cmp(P.x,Q.x)==0)
{ mpz_set_ui(R->x,0); mpz_set_ui(R->y,0); return;}
if(mpz_cmp(P.x,Q.x)==0 &&
mpz_cmp(P.y,Q.y)==0)
{
Point_Doubling(P,R);
return;
}
else
{
mpz_sub(temp,P.x,Q.x);
mpz_mod(temp,temp,EC.p);
mpz_invert(temp,temp,EC.p);
mpz_sub(slope,P.y,Q.y);
mpz_mul(slope,slope,temp);
mpz_mod(slope,slope,EC.p);
mpz_mul(R->x,slope,slope);
mpz_sub(R->x,R->x,P.x);
mpz_sub(R->x,R->x,Q.x);
mpz_mod(R->x,R->x,EC.p);
mpz_sub(temp,P.x,R->x);
mpz_mul(R->y,slope,temp);
mpz_sub(R->y,R->y,P.y);
mpz_mod(R->y,R->y,EC.p);
return;
}
}
void Point_Doubling(struct Point P,struct Point *R)
{
mpz_t slope,temp;
mpz_init(temp);
mpz_init(slope);
if(mpz_cmp_ui(P.y,0)!=0)
{
mpz_mul_ui(temp,P.y,2);
mpz_invert(temp,temp,EC.p);
mpz_mul(slope,P.x,P.x);
mpz_mul_ui(slope,slope,3);
mpz_add(slope,slope,EC.a);
mpz_mul(slope,slope,temp);
mpz_mod(slope,slope,EC.p);
mpz_mul(R->x,slope,slope);
mpz_sub(R->x,R->x,P.x);
mpz_sub(R->x,R->x,P.x);
mpz_mod(R->x,R->x,EC.p);
mpz_sub(temp,P.x,R->x);
mpz_mul(R->y,slope,temp);
mpz_sub(R->y,R->y,P.y);
mpz_mod(R->y,R->y,EC.p);
}
else
{
mpz_set_ui(R->x,0);
mpz_set_ui(R->y,0);
}
}
void Scalar_Multiplication(struct Point P,struct Point* R, mpz_t m)
{
struct Point Q,T;
mpz_init(Q.x); mpz_init(Q.y);
mpz_init(T.x); mpz_init(T.y);
long no_of_bits,loop;
no_of_bits=mpz_sizeinbase(m,2);
mpz_set_ui(R->x,0);mpz_set_ui(R->y,0);
if(mpz_cmp_ui(m,0)==0)
return;
mpz_set(Q.x,P.x);
mpz_set(Q.y,P.y);
if(mpz_tstbit(m,0)==1)
{mpz_set(R->x,P.x);mpz_set(R->y,P.y);}
for(loop=1;loop<no_of_bits;loop++)
{
mpz_set_ui(T.x,0);
mpz_set_ui(T.y,0);
Point_Doubling(Q,&T);
gmp_printf("\n %Zd %Zd %Zd %Zd ",Q.x,Q.y,T.x,T.y);
mpz_set(Q.x,T.x);
mpz_set(Q.y,T.y);
mpz_set(T.x,R->x);
mpz_set(T.y,R->y);
if(mpz_tstbit(m,loop))
Point_Addition(T,Q,R);
}
}
// int contract_main(int argc, char **argv){
// mpz_init(EC.a); mpz_init(EC.b); mpz_init(EC.p);
// mpz_init_set_str (EC.a, "0", 10);
// mpz_init_set_str (EC.b, "7", 10);
// mpz_init_set_str (EC.p, "115792089237316195423570985008687907853269984665640564039457584007908834671663", 10);
// struct Point P,R;
// mpz_init(P.x);
// mpz_init(P.y);
// mpz_init_set_ui(R.x,0);
// mpz_init_set_ui(R.y,0);
// mpz_init_set_str (P.x, "55066263022277343669578718895168534326250603453777594175500187360389116729240", 10);
// mpz_init_set_str (P.y, "32670510020758816978083085130507043184471273380659243275938904335757337482424", 10);
// //gmp_scanf("%Zd",&P.x);
// //gmp_scanf("%Zd",&P.y);
// struct Point Q;
// mpz_init(Q.x);
// mpz_init(Q.y);
// mpz_init_set_str (Q.x, "48840125481190545212233038815866595080052446190700286140915905752173158787212", 10);
// mpz_init_set_str (Q.y, "55258565891714835096390285893857201264507235115974921462591453996678450950230", 10);
// //gmp_scanf("%Zd",&Q.x);
// //gmp_scanf("%Zd",&Q.y);
// Point_Addition(P,Q,&R);
// gmp_printf("\n Resultant Point is %Zd,%Zd",R.x,R.y);
// //printf("%d, %d",R.x,R.y);
// return 0;
// }
```
### \<gmp\> on linux
1. download from [gmplib](https://gmplib.org/)
2. download `lzip`:
```bash
sudo apt install lzip
```
3. Conver `gmp-6.1.2.tar.lz` to `gmp-6.1.2.tar`
```bash
sudo lzip -d gmp-6.1.2.tar.lz
```
4. Unzip `gmp-6.1.2.tar`
```bash
tar -xvf gmp-6.1.2.tar
```
5. Install GMP
```bash
cd gmp-6.1.2/
sudo ./configure
sudo make
sudo make check
sudo make install
```
### \<gmp\> tutorial
#### Headers and Libraries
- header: `#include <gmp.h>`
- C compile: `gcc code.c -lgmp`
- C++ compile `g++ code.cc -lgmpxx -lgmp`
- 若GMP被安裝在非標準的地方，則用`-I`及`-L`使編譯時指向正確的資料夾
#### Nomenclature and Types
- `mpz_t`: 一般整數(Integers)
- 函數開頭為 `mpz_`
- `mpq_t`: 有理數(Rational number)
- 函數開頭為 `mpq_`
- `mpf_t`: 浮點數(Floating point number)
- 函數開頭為 `mpf_`
#### Variable Conventions
- input和output可以吃相同的變數
- 如整數乘法`mpz_mul`，若想要對`x`平方
則用`mpz_mul(x,x,x);`
- 在assign到一個GMP變數前，必須呼叫一個**初始化函數**
在用完該變數後，要呼叫清除函數
- 範例：
```cpp=
mpz_t n;
mpz_init(n); // Initializing
mpz_mul(n,...);
mpz_clear(n); //clearing
```
#### Parameter Conventions
- 是call-by-reference
- 函數呼叫只返回座標，所以用法是將回傳結果`result`放入函數參數中
- 範例：
```cpp=
void foo(mpz_t result, const mpz_t param, unsigned long n){
mpz_mul_ui(result, param, n);
}
int main(void){
mpz_t r,n;
// init
mpz_init(r);
mpz_init_set_str (n,"123456",0);
// call function foo
foo(r,n,20L);
// print
gmp_printf("%Zd\n",r);
return 0;
}
```
#### Demonstration programs
1. `demos/` 資料夾中有範例程式
```bash
cd demos/
```
2. make
```bash
sudo make pexpr
```
3. 執行
```bash
./pexpr 68^975+10
```
#### Integer functions
- `void mpz_init(mpz_t x)`: Initialize x, and set its value to 0
- `void mpz_clear(mpz_t x)`: free the space
- `void mpz_set_str(mpz_t rop, const char *str, int base)`: set the value of $rop$ from $str$
- `int mpz_init_set_str(mpz_t rop, const char *str, int base)`: initialize $rop$ and set value like `mpz_set_str` (-1 if an error occurs)
- `void mpz_add(mpz_t rop, const mpz_t op1, const mpz_t op2)`: $rop:=op1+op2$
- `void mpz_powm(mpz_t rop, const mpz_t base, const mpz_t exp, const mpz_t mod)`: $rop:=base^{exp}\ \text{mod}\ mod$
- formatted output
```cpp=
mpz_t z;
gmp_printf("%s is an mpz %Zd\n","here",z);
```
- formatted input
```cpp=
// to read say "a(5)=1234"
int n;
mpz_t z;
gmp_scanf("a(%d)=%Zd\n",&n,z);
```
Sign in with Wallet
Connect another wallet