1// Copyright 2020 ConsenSys Software Inc.
2//
3// Licensed under the Apache License, Version 2.0 (the "License");
4// you may not use this file except in compliance with the License.
5// You may obtain a copy of the License at
6//
7//     http://www.apache.org/licenses/LICENSE-2.0
8//
9// Unless required by applicable law or agreed to in writing, software
10// distributed under the License is distributed on an "AS IS" BASIS,
11// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12// See the License for the specific language governing permissions and
13// limitations under the License.
14
15// Code generated by consensys/gnark-crypto DO NOT EDIT
16
17package bls12381
18
19import (
20	"math/big"
21
22	"github.com/consensys/gnark-crypto/ecc"
23	"github.com/consensys/gnark-crypto/ecc/bls12-381/fp"
24	"github.com/consensys/gnark-crypto/ecc/bls12-381/fr"
25	"github.com/consensys/gnark-crypto/internal/parallel"
26)
27
28// G1Affine point in affine coordinates
29type G1Affine struct {
30	X, Y fp.Element
31}
32
33// G1Jac is a point with fp.Element coordinates
34type G1Jac struct {
35	X, Y, Z fp.Element
36}
37
38//  g1JacExtended parameterized jacobian coordinates (x=X/ZZ, y=Y/ZZZ, ZZ**3=ZZZ**2)
39type g1JacExtended struct {
40	X, Y, ZZ, ZZZ fp.Element
41}
42
43// -------------------------------------------------------------------------------------------------
44// Affine
45
46// Set sets p to the provided point
47func (p *G1Affine) Set(a *G1Affine) *G1Affine {
48	p.X, p.Y = a.X, a.Y
49	return p
50}
51
52// ScalarMultiplication computes and returns p = a*s
53func (p *G1Affine) ScalarMultiplication(a *G1Affine, s *big.Int) *G1Affine {
54	var _p G1Jac
55	_p.FromAffine(a)
56	_p.mulGLV(&_p, s)
57	p.FromJacobian(&_p)
58	return p
59}
60
61// Equal tests if two points (in Affine coordinates) are equal
62func (p *G1Affine) Equal(a *G1Affine) bool {
63	return p.X.Equal(&a.X) && p.Y.Equal(&a.Y)
64}
65
66// Neg computes -G
67func (p *G1Affine) Neg(a *G1Affine) *G1Affine {
68	p.X = a.X
69	p.Y.Neg(&a.Y)
70	return p
71}
72
73// FromJacobian rescale a point in Jacobian coord in z=1 plane
74func (p *G1Affine) FromJacobian(p1 *G1Jac) *G1Affine {
75
76	var a, b fp.Element
77
78	if p1.Z.IsZero() {
79		p.X.SetZero()
80		p.Y.SetZero()
81		return p
82	}
83
84	a.Inverse(&p1.Z)
85	b.Square(&a)
86	p.X.Mul(&p1.X, &b)
87	p.Y.Mul(&p1.Y, &b).Mul(&p.Y, &a)
88
89	return p
90}
91
92func (p *G1Affine) String() string {
93	var x, y fp.Element
94	x.Set(&p.X)
95	y.Set(&p.Y)
96	return "E([" + x.String() + "," + y.String() + "]),"
97}
98
99// IsInfinity checks if the point is infinity (in affine, it's encoded as (0,0))
100func (p *G1Affine) IsInfinity() bool {
101	return p.X.IsZero() && p.Y.IsZero()
102}
103
104// IsOnCurve returns true if p in on the curve
105func (p *G1Affine) IsOnCurve() bool {
106	var point G1Jac
107	point.FromAffine(p)
108	return point.IsOnCurve() // call this function to handle infinity point
109}
110
111// IsInSubGroup returns true if p is in the correct subgroup, false otherwise
112func (p *G1Affine) IsInSubGroup() bool {
113	var _p G1Jac
114	_p.FromAffine(p)
115	return _p.IsOnCurve() && _p.IsInSubGroup()
116}
117
118// -------------------------------------------------------------------------------------------------
119// Jacobian
120
121// Set sets p to the provided point
122func (p *G1Jac) Set(a *G1Jac) *G1Jac {
123	p.X, p.Y, p.Z = a.X, a.Y, a.Z
124	return p
125}
126
127// Equal tests if two points (in Jacobian coordinates) are equal
128func (p *G1Jac) Equal(a *G1Jac) bool {
129
130	if p.Z.IsZero() && a.Z.IsZero() {
131		return true
132	}
133	_p := G1Affine{}
134	_p.FromJacobian(p)
135
136	_a := G1Affine{}
137	_a.FromJacobian(a)
138
139	return _p.X.Equal(&_a.X) && _p.Y.Equal(&_a.Y)
140}
141
142// Neg computes -G
143func (p *G1Jac) Neg(a *G1Jac) *G1Jac {
144	*p = *a
145	p.Y.Neg(&a.Y)
146	return p
147}
148
149// SubAssign substracts two points on the curve
150func (p *G1Jac) SubAssign(a *G1Jac) *G1Jac {
151	var tmp G1Jac
152	tmp.Set(a)
153	tmp.Y.Neg(&tmp.Y)
154	p.AddAssign(&tmp)
155	return p
156}
157
158// AddAssign point addition in montgomery form
159// https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-2007-bl
160func (p *G1Jac) AddAssign(a *G1Jac) *G1Jac {
161
162	// p is infinity, return a
163	if p.Z.IsZero() {
164		p.Set(a)
165		return p
166	}
167
168	// a is infinity, return p
169	if a.Z.IsZero() {
170		return p
171	}
172
173	var Z1Z1, Z2Z2, U1, U2, S1, S2, H, I, J, r, V fp.Element
174	Z1Z1.Square(&a.Z)
175	Z2Z2.Square(&p.Z)
176	U1.Mul(&a.X, &Z2Z2)
177	U2.Mul(&p.X, &Z1Z1)
178	S1.Mul(&a.Y, &p.Z).
179		Mul(&S1, &Z2Z2)
180	S2.Mul(&p.Y, &a.Z).
181		Mul(&S2, &Z1Z1)
182
183	// if p == a, we double instead
184	if U1.Equal(&U2) && S1.Equal(&S2) {
185		return p.DoubleAssign()
186	}
187
188	H.Sub(&U2, &U1)
189	I.Double(&H).
190		Square(&I)
191	J.Mul(&H, &I)
192	r.Sub(&S2, &S1).Double(&r)
193	V.Mul(&U1, &I)
194	p.X.Square(&r).
195		Sub(&p.X, &J).
196		Sub(&p.X, &V).
197		Sub(&p.X, &V)
198	p.Y.Sub(&V, &p.X).
199		Mul(&p.Y, &r)
200	S1.Mul(&S1, &J).Double(&S1)
201	p.Y.Sub(&p.Y, &S1)
202	p.Z.Add(&p.Z, &a.Z)
203	p.Z.Square(&p.Z).
204		Sub(&p.Z, &Z1Z1).
205		Sub(&p.Z, &Z2Z2).
206		Mul(&p.Z, &H)
207
208	return p
209}
210
211// AddMixed point addition
212// http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-madd-2007-bl
213func (p *G1Jac) AddMixed(a *G1Affine) *G1Jac {
214
215	//if a is infinity return p
216	if a.X.IsZero() && a.Y.IsZero() {
217		return p
218	}
219	// p is infinity, return a
220	if p.Z.IsZero() {
221		p.X = a.X
222		p.Y = a.Y
223		p.Z.SetOne()
224		return p
225	}
226
227	var Z1Z1, U2, S2, H, HH, I, J, r, V fp.Element
228	Z1Z1.Square(&p.Z)
229	U2.Mul(&a.X, &Z1Z1)
230	S2.Mul(&a.Y, &p.Z).
231		Mul(&S2, &Z1Z1)
232
233	// if p == a, we double instead
234	if U2.Equal(&p.X) && S2.Equal(&p.Y) {
235		return p.DoubleAssign()
236	}
237
238	H.Sub(&U2, &p.X)
239	HH.Square(&H)
240	I.Double(&HH).Double(&I)
241	J.Mul(&H, &I)
242	r.Sub(&S2, &p.Y).Double(&r)
243	V.Mul(&p.X, &I)
244	p.X.Square(&r).
245		Sub(&p.X, &J).
246		Sub(&p.X, &V).
247		Sub(&p.X, &V)
248	J.Mul(&J, &p.Y).Double(&J)
249	p.Y.Sub(&V, &p.X).
250		Mul(&p.Y, &r)
251	p.Y.Sub(&p.Y, &J)
252	p.Z.Add(&p.Z, &H)
253	p.Z.Square(&p.Z).
254		Sub(&p.Z, &Z1Z1).
255		Sub(&p.Z, &HH)
256
257	return p
258}
259
260// Double doubles a point in Jacobian coordinates
261// https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#doubling-dbl-2007-bl
262func (p *G1Jac) Double(q *G1Jac) *G1Jac {
263	p.Set(q)
264	p.DoubleAssign()
265	return p
266}
267
268// DoubleAssign doubles a point in Jacobian coordinates
269// https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#doubling-dbl-2007-bl
270func (p *G1Jac) DoubleAssign() *G1Jac {
271
272	var XX, YY, YYYY, ZZ, S, M, T fp.Element
273
274	XX.Square(&p.X)
275	YY.Square(&p.Y)
276	YYYY.Square(&YY)
277	ZZ.Square(&p.Z)
278	S.Add(&p.X, &YY)
279	S.Square(&S).
280		Sub(&S, &XX).
281		Sub(&S, &YYYY).
282		Double(&S)
283	M.Double(&XX).Add(&M, &XX)
284	p.Z.Add(&p.Z, &p.Y).
285		Square(&p.Z).
286		Sub(&p.Z, &YY).
287		Sub(&p.Z, &ZZ)
288	T.Square(&M)
289	p.X = T
290	T.Double(&S)
291	p.X.Sub(&p.X, &T)
292	p.Y.Sub(&S, &p.X).
293		Mul(&p.Y, &M)
294	YYYY.Double(&YYYY).Double(&YYYY).Double(&YYYY)
295	p.Y.Sub(&p.Y, &YYYY)
296
297	return p
298}
299
300// ScalarMultiplication computes and returns p = a*s
301// see https://www.iacr.org/archive/crypto2001/21390189.pdf
302func (p *G1Jac) ScalarMultiplication(a *G1Jac, s *big.Int) *G1Jac {
303	return p.mulGLV(a, s)
304}
305
306func (p *G1Jac) String() string {
307	if p.Z.IsZero() {
308		return "O"
309	}
310	_p := G1Affine{}
311	_p.FromJacobian(p)
312	return "E([" + _p.X.String() + "," + _p.Y.String() + "]),"
313}
314
315// FromAffine sets p = Q, p in Jacboian, Q in affine
316func (p *G1Jac) FromAffine(Q *G1Affine) *G1Jac {
317	if Q.X.IsZero() && Q.Y.IsZero() {
318		p.Z.SetZero()
319		p.X.SetOne()
320		p.Y.SetOne()
321		return p
322	}
323	p.Z.SetOne()
324	p.X.Set(&Q.X)
325	p.Y.Set(&Q.Y)
326	return p
327}
328
329// IsOnCurve returns true if p in on the curve
330func (p *G1Jac) IsOnCurve() bool {
331	var left, right, tmp fp.Element
332	left.Square(&p.Y)
333	right.Square(&p.X).Mul(&right, &p.X)
334	tmp.Square(&p.Z).
335		Square(&tmp).
336		Mul(&tmp, &p.Z).
337		Mul(&tmp, &p.Z).
338		Mul(&tmp, &bCurveCoeff)
339	right.Add(&right, &tmp)
340	return left.Equal(&right)
341}
342
343// IsInSubGroup returns true if p is on the r-torsion, false otherwise.
344// Z[r,0]+Z[-lambdaG1Affine, 1] is the kernel
345// of (u,v)->u+lambdaG1Affinev mod r. Expressing r, lambdaG1Affine as
346// polynomials in x, a short vector of this Zmodule is
347// 1, x**2. So we check that p+x**2*phi(p)
348// is the infinity.
349func (p *G1Jac) IsInSubGroup() bool {
350
351	var res G1Jac
352	res.phi(p).
353		ScalarMultiplication(&res, &xGen).
354		ScalarMultiplication(&res, &xGen).
355		AddAssign(p)
356
357	return res.IsOnCurve() && res.Z.IsZero()
358
359}
360
361// mulWindowed 2-bits windowed exponentiation
362func (p *G1Jac) mulWindowed(a *G1Jac, s *big.Int) *G1Jac {
363
364	var res G1Jac
365	var ops [3]G1Jac
366
367	res.Set(&g1Infinity)
368	ops[0].Set(a)
369	ops[1].Double(&ops[0])
370	ops[2].Set(&ops[0]).AddAssign(&ops[1])
371
372	b := s.Bytes()
373	for i := range b {
374		w := b[i]
375		mask := byte(0xc0)
376		for j := 0; j < 4; j++ {
377			res.DoubleAssign().DoubleAssign()
378			c := (w & mask) >> (6 - 2*j)
379			if c != 0 {
380				res.AddAssign(&ops[c-1])
381			}
382			mask = mask >> 2
383		}
384	}
385	p.Set(&res)
386
387	return p
388
389}
390
391// phi assigns p to phi(a) where phi: (x,y)->(ux,y), and returns p
392func (p *G1Jac) phi(a *G1Jac) *G1Jac {
393	p.Set(a)
394	p.X.Mul(&p.X, &thirdRootOneG1)
395	return p
396}
397
398// mulGLV performs scalar multiplication using GLV
399// see https://www.iacr.org/archive/crypto2001/21390189.pdf
400func (p *G1Jac) mulGLV(a *G1Jac, s *big.Int) *G1Jac {
401
402	var table [15]G1Jac
403	var zero big.Int
404	var res G1Jac
405	var k1, k2 fr.Element
406
407	res.Set(&g1Infinity)
408
409	// table[b3b2b1b0-1] = b3b2*phi(a) + b1b0*a
410	table[0].Set(a)
411	table[3].phi(a)
412
413	// split the scalar, modifies +-a, phi(a) accordingly
414	k := ecc.SplitScalar(s, &glvBasis)
415
416	if k[0].Cmp(&zero) == -1 {
417		k[0].Neg(&k[0])
418		table[0].Neg(&table[0])
419	}
420	if k[1].Cmp(&zero) == -1 {
421		k[1].Neg(&k[1])
422		table[3].Neg(&table[3])
423	}
424
425	// precompute table (2 bits sliding window)
426	// table[b3b2b1b0-1] = b3b2*phi(a) + b1b0*a if b3b2b1b0 != 0
427	table[1].Double(&table[0])
428	table[2].Set(&table[1]).AddAssign(&table[0])
429	table[4].Set(&table[3]).AddAssign(&table[0])
430	table[5].Set(&table[3]).AddAssign(&table[1])
431	table[6].Set(&table[3]).AddAssign(&table[2])
432	table[7].Double(&table[3])
433	table[8].Set(&table[7]).AddAssign(&table[0])
434	table[9].Set(&table[7]).AddAssign(&table[1])
435	table[10].Set(&table[7]).AddAssign(&table[2])
436	table[11].Set(&table[7]).AddAssign(&table[3])
437	table[12].Set(&table[11]).AddAssign(&table[0])
438	table[13].Set(&table[11]).AddAssign(&table[1])
439	table[14].Set(&table[11]).AddAssign(&table[2])
440
441	// bounds on the lattice base vectors guarantee that k1, k2 are len(r)/2 bits long max
442	k1.SetBigInt(&k[0]).FromMont()
443	k2.SetBigInt(&k[1]).FromMont()
444
445	// loop starts from len(k1)/2 due to the bounds
446	for i := len(k1)/2 - 1; i >= 0; i-- {
447		mask := uint64(3) << 62
448		for j := 0; j < 32; j++ {
449			res.Double(&res).Double(&res)
450			b1 := (k1[i] & mask) >> (62 - 2*j)
451			b2 := (k2[i] & mask) >> (62 - 2*j)
452			if b1|b2 != 0 {
453				s := (b2<<2 | b1)
454				res.AddAssign(&table[s-1])
455			}
456			mask = mask >> 2
457		}
458	}
459
460	p.Set(&res)
461	return p
462}
463
464// ClearCofactor ...
465func (p *G1Affine) ClearCofactor(a *G1Affine) *G1Affine {
466	var _p G1Jac
467	_p.FromAffine(a)
468	_p.ClearCofactor(&_p)
469	p.FromJacobian(&_p)
470	return p
471}
472
473// ClearCofactor maps a point in E(Fp) to E(Fp)[r]
474func (p *G1Jac) ClearCofactor(a *G1Jac) *G1Jac {
475	// cf https://eprint.iacr.org/2019/403.pdf, 5
476	var res G1Jac
477	res.ScalarMultiplication(a, &xGen).AddAssign(a)
478	p.Set(&res)
479	return p
480
481}
482
483// -------------------------------------------------------------------------------------------------
484// Jacobian extended
485
486// Set sets p to the provided point
487func (p *g1JacExtended) Set(a *g1JacExtended) *g1JacExtended {
488	p.X, p.Y, p.ZZ, p.ZZZ = a.X, a.Y, a.ZZ, a.ZZZ
489	return p
490}
491
492// setInfinity sets p to O
493func (p *g1JacExtended) setInfinity() *g1JacExtended {
494	p.X.SetOne()
495	p.Y.SetOne()
496	p.ZZ = fp.Element{}
497	p.ZZZ = fp.Element{}
498	return p
499}
500
501// fromJacExtended sets Q in affine coords
502func (p *G1Affine) fromJacExtended(Q *g1JacExtended) *G1Affine {
503	if Q.ZZ.IsZero() {
504		p.X = fp.Element{}
505		p.Y = fp.Element{}
506		return p
507	}
508	p.X.Inverse(&Q.ZZ).Mul(&p.X, &Q.X)
509	p.Y.Inverse(&Q.ZZZ).Mul(&p.Y, &Q.Y)
510	return p
511}
512
513// fromJacExtended sets Q in Jacobian coords
514func (p *G1Jac) fromJacExtended(Q *g1JacExtended) *G1Jac {
515	if Q.ZZ.IsZero() {
516		p.Set(&g1Infinity)
517		return p
518	}
519	p.X.Mul(&Q.ZZ, &Q.X).Mul(&p.X, &Q.ZZ)
520	p.Y.Mul(&Q.ZZZ, &Q.Y).Mul(&p.Y, &Q.ZZZ)
521	p.Z.Set(&Q.ZZZ)
522	return p
523}
524
525// unsafeFromJacExtended sets p in jacobian coords, but don't check for infinity
526func (p *G1Jac) unsafeFromJacExtended(Q *g1JacExtended) *G1Jac {
527	p.X.Square(&Q.ZZ).Mul(&p.X, &Q.X)
528	p.Y.Square(&Q.ZZZ).Mul(&p.Y, &Q.Y)
529	p.Z = Q.ZZZ
530	return p
531}
532
533// add point in ZZ coords
534// https://www.hyperelliptic.org/EFD/g1p/auto-shortw-xyzz.html#addition-add-2008-s
535func (p *g1JacExtended) add(q *g1JacExtended) *g1JacExtended {
536	//if q is infinity return p
537	if q.ZZ.IsZero() {
538		return p
539	}
540	// p is infinity, return q
541	if p.ZZ.IsZero() {
542		p.Set(q)
543		return p
544	}
545
546	var A, B, X1ZZ2, X2ZZ1, Y1ZZZ2, Y2ZZZ1 fp.Element
547
548	// p2: q, p1: p
549	X2ZZ1.Mul(&q.X, &p.ZZ)
550	X1ZZ2.Mul(&p.X, &q.ZZ)
551	A.Sub(&X2ZZ1, &X1ZZ2)
552	Y2ZZZ1.Mul(&q.Y, &p.ZZZ)
553	Y1ZZZ2.Mul(&p.Y, &q.ZZZ)
554	B.Sub(&Y2ZZZ1, &Y1ZZZ2)
555
556	if A.IsZero() {
557		if B.IsZero() {
558			return p.double(q)
559
560		}
561		p.ZZ = fp.Element{}
562		p.ZZZ = fp.Element{}
563		return p
564	}
565
566	var U1, U2, S1, S2, P, R, PP, PPP, Q, V fp.Element
567	U1.Mul(&p.X, &q.ZZ)
568	U2.Mul(&q.X, &p.ZZ)
569	S1.Mul(&p.Y, &q.ZZZ)
570	S2.Mul(&q.Y, &p.ZZZ)
571	P.Sub(&U2, &U1)
572	R.Sub(&S2, &S1)
573	PP.Square(&P)
574	PPP.Mul(&P, &PP)
575	Q.Mul(&U1, &PP)
576	V.Mul(&S1, &PPP)
577
578	p.X.Square(&R).
579		Sub(&p.X, &PPP).
580		Sub(&p.X, &Q).
581		Sub(&p.X, &Q)
582	p.Y.Sub(&Q, &p.X).
583		Mul(&p.Y, &R).
584		Sub(&p.Y, &V)
585	p.ZZ.Mul(&p.ZZ, &q.ZZ).
586		Mul(&p.ZZ, &PP)
587	p.ZZZ.Mul(&p.ZZZ, &q.ZZZ).
588		Mul(&p.ZZZ, &PPP)
589
590	return p
591}
592
593// double point in ZZ coords
594// http://www.hyperelliptic.org/EFD/g1p/auto-shortw-xyzz.html#doubling-dbl-2008-s-1
595func (p *g1JacExtended) double(q *g1JacExtended) *g1JacExtended {
596	var U, V, W, S, XX, M fp.Element
597
598	U.Double(&q.Y)
599	V.Square(&U)
600	W.Mul(&U, &V)
601	S.Mul(&q.X, &V)
602	XX.Square(&q.X)
603	M.Double(&XX).
604		Add(&M, &XX) // -> + a, but a=0 here
605	U.Mul(&W, &q.Y)
606
607	p.X.Square(&M).
608		Sub(&p.X, &S).
609		Sub(&p.X, &S)
610	p.Y.Sub(&S, &p.X).
611		Mul(&p.Y, &M).
612		Sub(&p.Y, &U)
613	p.ZZ.Mul(&V, &q.ZZ)
614	p.ZZZ.Mul(&W, &q.ZZZ)
615
616	return p
617}
618
619// subMixed same as addMixed, but will negate a.Y
620// http://www.hyperelliptic.org/EFD/g1p/auto-shortw-xyzz.html#addition-madd-2008-s
621func (p *g1JacExtended) subMixed(a *G1Affine) *g1JacExtended {
622
623	//if a is infinity return p
624	if a.X.IsZero() && a.Y.IsZero() {
625		return p
626	}
627	// p is infinity, return a
628	if p.ZZ.IsZero() {
629		p.X = a.X
630		p.Y.Neg(&a.Y)
631		p.ZZ.SetOne()
632		p.ZZZ.SetOne()
633		return p
634	}
635
636	var P, R fp.Element
637
638	// p2: a, p1: p
639	P.Mul(&a.X, &p.ZZ)
640	P.Sub(&P, &p.X)
641
642	R.Mul(&a.Y, &p.ZZZ)
643	R.Neg(&R)
644	R.Sub(&R, &p.Y)
645
646	if P.IsZero() {
647		if R.IsZero() {
648			return p.doubleNegMixed(a)
649
650		}
651		p.ZZ = fp.Element{}
652		p.ZZZ = fp.Element{}
653		return p
654	}
655
656	var PP, PPP, Q, Q2, RR, X3, Y3 fp.Element
657
658	PP.Square(&P)
659	PPP.Mul(&P, &PP)
660	Q.Mul(&p.X, &PP)
661	RR.Square(&R)
662	X3.Sub(&RR, &PPP)
663	Q2.Double(&Q)
664	p.X.Sub(&X3, &Q2)
665	Y3.Sub(&Q, &p.X).Mul(&Y3, &R)
666	R.Mul(&p.Y, &PPP)
667	p.Y.Sub(&Y3, &R)
668	p.ZZ.Mul(&p.ZZ, &PP)
669	p.ZZZ.Mul(&p.ZZZ, &PPP)
670
671	return p
672
673}
674
675// addMixed
676// http://www.hyperelliptic.org/EFD/g1p/auto-shortw-xyzz.html#addition-madd-2008-s
677func (p *g1JacExtended) addMixed(a *G1Affine) *g1JacExtended {
678
679	//if a is infinity return p
680	if a.X.IsZero() && a.Y.IsZero() {
681		return p
682	}
683	// p is infinity, return a
684	if p.ZZ.IsZero() {
685		p.X = a.X
686		p.Y = a.Y
687		p.ZZ.SetOne()
688		p.ZZZ.SetOne()
689		return p
690	}
691
692	var P, R fp.Element
693
694	// p2: a, p1: p
695	P.Mul(&a.X, &p.ZZ)
696	P.Sub(&P, &p.X)
697
698	R.Mul(&a.Y, &p.ZZZ)
699	R.Sub(&R, &p.Y)
700
701	if P.IsZero() {
702		if R.IsZero() {
703			return p.doubleMixed(a)
704
705		}
706		p.ZZ = fp.Element{}
707		p.ZZZ = fp.Element{}
708		return p
709	}
710
711	var PP, PPP, Q, Q2, RR, X3, Y3 fp.Element
712
713	PP.Square(&P)
714	PPP.Mul(&P, &PP)
715	Q.Mul(&p.X, &PP)
716	RR.Square(&R)
717	X3.Sub(&RR, &PPP)
718	Q2.Double(&Q)
719	p.X.Sub(&X3, &Q2)
720	Y3.Sub(&Q, &p.X).Mul(&Y3, &R)
721	R.Mul(&p.Y, &PPP)
722	p.Y.Sub(&Y3, &R)
723	p.ZZ.Mul(&p.ZZ, &PP)
724	p.ZZZ.Mul(&p.ZZZ, &PPP)
725
726	return p
727
728}
729
730// doubleNegMixed same as double, but will negate q.Y
731func (p *g1JacExtended) doubleNegMixed(q *G1Affine) *g1JacExtended {
732
733	var U, V, W, S, XX, M, S2, L fp.Element
734
735	U.Double(&q.Y)
736	U.Neg(&U)
737	V.Square(&U)
738	W.Mul(&U, &V)
739	S.Mul(&q.X, &V)
740	XX.Square(&q.X)
741	M.Double(&XX).
742		Add(&M, &XX) // -> + a, but a=0 here
743	S2.Double(&S)
744	L.Mul(&W, &q.Y)
745
746	p.X.Square(&M).
747		Sub(&p.X, &S2)
748	p.Y.Sub(&S, &p.X).
749		Mul(&p.Y, &M).
750		Add(&p.Y, &L)
751	p.ZZ.Set(&V)
752	p.ZZZ.Set(&W)
753
754	return p
755}
756
757// doubleMixed point in ZZ coords
758// http://www.hyperelliptic.org/EFD/g1p/auto-shortw-xyzz.html#doubling-dbl-2008-s-1
759func (p *g1JacExtended) doubleMixed(q *G1Affine) *g1JacExtended {
760
761	var U, V, W, S, XX, M, S2, L fp.Element
762
763	U.Double(&q.Y)
764	V.Square(&U)
765	W.Mul(&U, &V)
766	S.Mul(&q.X, &V)
767	XX.Square(&q.X)
768	M.Double(&XX).
769		Add(&M, &XX) // -> + a, but a=0 here
770	S2.Double(&S)
771	L.Mul(&W, &q.Y)
772
773	p.X.Square(&M).
774		Sub(&p.X, &S2)
775	p.Y.Sub(&S, &p.X).
776		Mul(&p.Y, &M).
777		Sub(&p.Y, &L)
778	p.ZZ.Set(&V)
779	p.ZZZ.Set(&W)
780
781	return p
782}
783
784// BatchJacobianToAffineG1 converts points in Jacobian coordinates to Affine coordinates
785// performing a single field inversion (Montgomery batch inversion trick)
786// result must be allocated with len(result) == len(points)
787func BatchJacobianToAffineG1(points []G1Jac, result []G1Affine) {
788	zeroes := make([]bool, len(points))
789	accumulator := fp.One()
790
791	// batch invert all points[].Z coordinates with Montgomery batch inversion trick
792	// (stores points[].Z^-1 in result[i].X to avoid allocating a slice of fr.Elements)
793	for i := 0; i < len(points); i++ {
794		if points[i].Z.IsZero() {
795			zeroes[i] = true
796			continue
797		}
798		result[i].X = accumulator
799		accumulator.Mul(&accumulator, &points[i].Z)
800	}
801
802	var accInverse fp.Element
803	accInverse.Inverse(&accumulator)
804
805	for i := len(points) - 1; i >= 0; i-- {
806		if zeroes[i] {
807			// do nothing, X and Y are zeroes in affine.
808			continue
809		}
810		result[i].X.Mul(&result[i].X, &accInverse)
811		accInverse.Mul(&accInverse, &points[i].Z)
812	}
813
814	// batch convert to affine.
815	parallel.Execute(len(points), func(start, end int) {
816		for i := start; i < end; i++ {
817			if zeroes[i] {
818				// do nothing, X and Y are zeroes in affine.
819				continue
820			}
821			var a, b fp.Element
822			a = result[i].X
823			b.Square(&a)
824			result[i].X.Mul(&points[i].X, &b)
825			result[i].Y.Mul(&points[i].Y, &b).
826				Mul(&result[i].Y, &a)
827		}
828	})
829
830}
831
832// BatchScalarMultiplicationG1 multiplies the same base (generator) by all scalars
833// and return resulting points in affine coordinates
834// uses a simple windowed-NAF like exponentiation algorithm
835func BatchScalarMultiplicationG1(base *G1Affine, scalars []fr.Element) []G1Affine {
836
837	// approximate cost in group ops is
838	// cost = 2^{c-1} + n(scalar.nbBits+nbChunks)
839
840	nbPoints := uint64(len(scalars))
841	min := ^uint64(0)
842	bestC := 0
843	for c := 2; c < 18; c++ {
844		cost := uint64(1 << (c - 1))
845		nbChunks := uint64(fr.Limbs * 64 / c)
846		if (fr.Limbs*64)%c != 0 {
847			nbChunks++
848		}
849		cost += nbPoints * ((fr.Limbs * 64) + nbChunks)
850		if cost < min {
851			min = cost
852			bestC = c
853		}
854	}
855	c := uint64(bestC) // window size
856	nbChunks := int(fr.Limbs * 64 / c)
857	if (fr.Limbs*64)%c != 0 {
858		nbChunks++
859	}
860	mask := uint64((1 << c) - 1) // low c bits are 1
861	msbWindow := uint64(1 << (c - 1))
862
863	// precompute all powers of base for our window
864	// note here that if performance is critical, we can implement as in the msmX methods
865	// this allocation to be on the stack
866	baseTable := make([]G1Jac, (1 << (c - 1)))
867	baseTable[0].Set(&g1Infinity)
868	baseTable[0].AddMixed(base)
869	for i := 1; i < len(baseTable); i++ {
870		baseTable[i] = baseTable[i-1]
871		baseTable[i].AddMixed(base)
872	}
873
874	pScalars := partitionScalars(scalars, c)
875
876	// compute offset and word selector / shift to select the right bits of our windows
877	selectors := make([]selector, nbChunks)
878	for chunk := 0; chunk < nbChunks; chunk++ {
879		jc := uint64(uint64(chunk) * c)
880		d := selector{}
881		d.index = jc / 64
882		d.shift = jc - (d.index * 64)
883		d.mask = mask << d.shift
884		d.multiWordSelect = (64%c) != 0 && d.shift > (64-c) && d.index < (fr.Limbs-1)
885		if d.multiWordSelect {
886			nbBitsHigh := d.shift - uint64(64-c)
887			d.maskHigh = (1 << nbBitsHigh) - 1
888			d.shiftHigh = (c - nbBitsHigh)
889		}
890		selectors[chunk] = d
891	}
892	// convert our base exp table into affine to use AddMixed
893	baseTableAff := make([]G1Affine, (1 << (c - 1)))
894	BatchJacobianToAffineG1(baseTable, baseTableAff)
895	toReturn := make([]G1Jac, len(scalars))
896
897	// for each digit, take value in the base table, double it c time, voila.
898	parallel.Execute(len(pScalars), func(start, end int) {
899		var p G1Jac
900		for i := start; i < end; i++ {
901			p.Set(&g1Infinity)
902			for chunk := nbChunks - 1; chunk >= 0; chunk-- {
903				s := selectors[chunk]
904				if chunk != nbChunks-1 {
905					for j := uint64(0); j < c; j++ {
906						p.DoubleAssign()
907					}
908				}
909
910				bits := (pScalars[i][s.index] & s.mask) >> s.shift
911				if s.multiWordSelect {
912					bits += (pScalars[i][s.index+1] & s.maskHigh) << s.shiftHigh
913				}
914
915				if bits == 0 {
916					continue
917				}
918
919				// if msbWindow bit is set, we need to substract
920				if bits&msbWindow == 0 {
921					// add
922					p.AddMixed(&baseTableAff[bits-1])
923				} else {
924					// sub
925					t := baseTableAff[bits & ^msbWindow]
926					t.Neg(&t)
927					p.AddMixed(&t)
928				}
929			}
930
931			// set our result point
932			toReturn[i] = p
933
934		}
935	})
936	toReturnAff := make([]G1Affine, len(scalars))
937	BatchJacobianToAffineG1(toReturn, toReturnAff)
938	return toReturnAff
939}
940