github.com/consensys/gnark-crypto@v0.14.0/internal/generator/ecc/template/point.go.tmpl

{{ $TAffine := print (toUpper .PointName) "Affine" }}
{{ $TJacobian := print (toUpper .PointName) "Jac" }}
{{ $TJacobianExtended := print (toLower .PointName) "JacExtended" }}
{{ $TProjective := print (toLower .PointName) "Proj" }}


import (
	{{- if eq .PointName "g2"}}
	"crypto/rand"
	{{- end}}
	"math/big"
	"runtime"

	{{- if .GLV}}
	"github.com/consensys/gnark-crypto/ecc"
	{{- end}}
	"github.com/consensys/gnark-crypto/internal/parallel"
	"github.com/consensys/gnark-crypto/ecc/{{.Name}}/fr"
	{{- if or (eq .CoordType "fptower.E2") (eq .CoordType "fptower.E4") }}
	"github.com/consensys/gnark-crypto/ecc/{{.Name}}/internal/fptower"
	{{else}}
	"github.com/consensys/gnark-crypto/ecc/{{.Name}}/fp"
	{{- end}}
)


// {{ $TAffine }} is a point in affine coordinates (x,y)
type {{ $TAffine }} struct {
	X, Y {{.CoordType}}
}

// {{ $TJacobian }} is a point in Jacobian coordinates (x=X/Z², y=Y/Z³)
type {{ $TJacobian }} struct {
	X, Y, Z {{.CoordType}}
}

// {{ $TJacobianExtended }} is a point in extended Jacobian coordinates (x=X/ZZ, y=Y/ZZZ, ZZ³=ZZZ²)
type {{ $TJacobianExtended }} struct {
	X, Y, ZZ, ZZZ {{.CoordType}}
}

{{- if .Projective}}
// {{ $TProjective }} point in projective coordinates
type {{ $TProjective }} struct {
	x, y, z {{.CoordType}}
}
{{- end}}



// -------------------------------------------------------------------------------------------------
// Affine coordinates

// Set sets p to a in affine coordinates.
func (p *{{ $TAffine }}) Set(a *{{ $TAffine }}) *{{ $TAffine }} {
       p.X, p.Y = a.X, a.Y
       return p
}

// setInfinity sets p to the infinity point, which is encoded as (0,0).
// N.B.: (0,0) is never on the curve for j=0 curves (Y²=X³+B).
func (p *{{ $TAffine }}) setInfinity() *{{ $TAffine }} {
	p.X.SetZero()
	p.Y.SetZero()
	return p
}

// ScalarMultiplication computes and returns p = [s]a
// where p and a are affine points.
func (p *{{ $TAffine }}) ScalarMultiplication(a *{{ $TAffine }}, s *big.Int) *{{ $TAffine }} {
	var _p {{ $TJacobian }}
	_p.FromAffine(a)
	{{- if .GLV}}
        _p.mulGLV(&_p, s)
	{{- else }}
        _p.mulWindowed(&_p, s)
	{{- end }}
	p.FromJacobian(&_p)
	return p
}

// ScalarMultiplicationBase computes and returns p = [s]g
// where g is the affine point generating the prime subgroup.
func (p *{{ $TAffine }}) ScalarMultiplicationBase(s *big.Int) *{{ $TAffine }} {
	var _p {{ $TJacobian }}
	{{- if .GLV}}
	_p.mulGLV(&{{ toLower .PointName}}Gen, s)
	{{- else }}
        _p.mulWindowed(&{{ toLower .PointName}}Gen, s)
	{{- end }}
	p.FromJacobian(&_p)
	return p
}


// Add adds two points in affine coordinates.
// It uses the Jacobian addition with a.Z=b.Z=1 and converts the result to affine coordinates.
//
// https://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-mmadd-2007-bl
func (p *{{ $TAffine }}) Add(a, b *{{ $TAffine }}) *{{ $TAffine }} {
	var q {{ $TJacobian }}
	// a is infinity, return b
	if a.IsInfinity() {
		p.Set(b)
		return p
	}
	// b is infinity, return a
	if b.IsInfinity() {
		p.Set(a)
		return p
	}
	if a.X.Equal(&b.X) {
		// if b == a, we double instead
		if a.Y.Equal(&b.Y) {
			q.DoubleMixed(a)
			return p.FromJacobian(&q)
		} else {
			// if b == -a, we return 0
			return p.setInfinity()
		}
	}
	var H, HH, I, J, r, V {{.CoordType}}
	H.Sub(&b.X, &a.X)
	HH.Square(&H)
	I.Double(&HH).Double(&I)
	J.Mul(&H, &I)
	r.Sub(&b.Y, &a.Y)
	r.Double(&r)
	V.Mul(&a.X, &I)
	q.X.Square(&r).
		Sub(&q.X, &J).
		Sub(&q.X, &V).
		Sub(&q.X, &V)
	q.Y.Sub(&V, &q.X).
		Mul(&q.Y, &r)
	J.Mul(&a.Y, &J).Double(&J)
	q.Y.Sub(&q.Y, &J)
	q.Z.Double(&H)

	return p.FromJacobian(&q)
}

// Double doubles a point in affine coordinates.
// It converts the point to Jacobian coordinates, doubles it using Jacobian
// addition with a.Z=1, and converts it back to affine coordinates.
//
// http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#doubling-mdbl-2007-bl
func (p *{{ $TAffine }}) Double(a *{{ $TAffine }}) *{{ $TAffine }} {
	var q {{ $TJacobian }}
	q.FromAffine(a)
	q.DoubleMixed(a)
	p.FromJacobian(&q)
	return p
}

// Sub subtracts two points in affine coordinates.
// It uses a similar approach to Add, but negates the second point before adding.
func (p *{{ $TAffine }}) Sub(a, b *{{ $TAffine }}) *{{ $TAffine }} {
	var bneg {{ $TAffine }}
	bneg.Neg(b)
	p.Add(a, &bneg)
	return p
}

// Equal tests if two points in affine coordinates are equal.
func (p *{{ $TAffine }}) Equal(a *{{ $TAffine }}) bool {
	return p.X.Equal(&a.X) && p.Y.Equal(&a.Y)
}


// Neg sets p to the affine negative point -a = (a.X, -a.Y).
func (p *{{ $TAffine }}) Neg(a *{{ $TAffine }}) *{{ $TAffine }} {
	p.X = a.X
	p.Y.Neg(&a.Y)
	return p
}

// FromJacobian converts a point p1 from Jacobian to affine coordinates.
func (p *{{ $TAffine }}) FromJacobian(p1 *{{ $TJacobian }}) *{{ $TAffine }} {

	var a, b {{.CoordType}}

	if p1.Z.IsZero() {
		p.X.SetZero()
		p.Y.SetZero()
		return p
	}

	a.Inverse(&p1.Z)
	b.Square(&a)
	p.X.Mul(&p1.X, &b)
	p.Y.Mul(&p1.Y, &b).Mul(&p.Y, &a)

	return p
}

// String returns the string representation E(x,y) of the affine point p or "O" if it is infinity.
func (p *{{ $TAffine }}) String() string {
	if p.IsInfinity() {
		return "O"
	}
	return "E([" + p.X.String() + "," + p.Y.String() + "])"
}

// IsInfinity checks if the affine point p is infinity, which is encoded as (0,0).
// N.B.: (0,0) is never on the curve for j=0 curves (Y²=X³+B).
func (p *{{ $TAffine }}) IsInfinity() bool {
	return p.X.IsZero() && p.Y.IsZero()
}

// IsOnCurve returns true if the affine point p in on the curve.
func (p *{{ $TAffine }}) IsOnCurve() bool {
	var point {{ $TJacobian }}
	point.FromAffine(p)
	return point.IsOnCurve() // call this function to handle infinity point
}

// IsInSubGroup returns true if the affine point p is in the correct subgroup, false otherwise.
func (p *{{ $TAffine }}) IsInSubGroup() bool {
	var _p {{ $TJacobian }}
	_p.FromAffine(p)
	return _p.IsInSubGroup()
}


// -------------------------------------------------------------------------------------------------
// Jacobian coordinates

// Set sets p to a in Jacobian coordinates.
func (p *{{ $TJacobian }}) Set(q *{{ $TJacobian }}) *{{ $TJacobian }} {
	p.X, p.Y, p.Z = q.X, q.Y, q.Z
	return p
}

// Equal tests if two points in Jacobian coordinates are equal.
func (p *{{ $TJacobian }}) Equal(q *{{ $TJacobian }}) bool {
	// If one point is infinity, the other must also be infinity.
	if p.Z.IsZero() {
		return q.Z.IsZero()
	}
	// If the other point is infinity, return false since we can't
	// the following checks would be incorrect.
	if q.Z.IsZero() {
		return false
	}

	var pZSquare, aZSquare {{.CoordType}}
	pZSquare.Square(&p.Z)
	aZSquare.Square(&q.Z)

	var lhs, rhs {{.CoordType}}
	lhs.Mul(&p.X, &aZSquare)
	rhs.Mul(&q.X, &pZSquare)
	if !lhs.Equal(&rhs) {
		return false
	}
	lhs.Mul(&p.Y, &aZSquare).Mul(&lhs, &q.Z)
	rhs.Mul(&q.Y, &pZSquare).Mul(&rhs, &p.Z)

	return lhs.Equal(&rhs)
}

// Neg sets p to the Jacobian negative point -q = (q.X, -q.Y, q.Z).
func (p *{{ $TJacobian }}) Neg(q *{{ $TJacobian }}) *{{ $TJacobian }} {
	*p = *q
	p.Y.Neg(&q.Y)
	return p
}

// AddAssign sets p to p+a in Jacobian coordinates.
//
// https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-2007-bl
func (p *{{ $TJacobian }}) AddAssign(q *{{ $TJacobian }}) *{{ $TJacobian }} {

	// p is infinity, return q
	if p.Z.IsZero() {
		p.Set(q)
		return p
	}

	// q is infinity, return p
	if q.Z.IsZero() {
		return p
	}

	var Z1Z1, Z2Z2, U1, U2, S1, S2, H, I, J, r, V {{.CoordType}}
	Z1Z1.Square(&q.Z)
	Z2Z2.Square(&p.Z)
	U1.Mul(&q.X, &Z2Z2)
	U2.Mul(&p.X, &Z1Z1)
	S1.Mul(&q.Y, &p.Z).
		Mul(&S1, &Z2Z2)
	S2.Mul(&p.Y, &q.Z).
		Mul(&S2, &Z1Z1)

	// if p == q, we double instead
	if U1.Equal(&U2) && S1.Equal(&S2) {
		return p.DoubleAssign()
	}

	H.Sub(&U2, &U1)
	I.Double(&H).
		Square(&I)
	J.Mul(&H, &I)
	r.Sub(&S2, &S1).Double(&r)
	V.Mul(&U1, &I)
	p.X.Square(&r).
		Sub(&p.X, &J).
		Sub(&p.X, &V).
		Sub(&p.X, &V)
	p.Y.Sub(&V, &p.X).
		Mul(&p.Y, &r)
	S1.Mul(&S1, &J).Double(&S1)
	p.Y.Sub(&p.Y, &S1)
	p.Z.Add(&p.Z, &q.Z)
	p.Z.Square(&p.Z).
		Sub(&p.Z, &Z1Z1).
		Sub(&p.Z, &Z2Z2).
		Mul(&p.Z, &H)

	return p
}

// SubAssign sets p to p-a in Jacobian coordinates.
// It uses a similar approach to AddAssign, but negates the point a before adding.
func (p *{{ $TJacobian }}) SubAssign(q *{{ $TJacobian }}) *{{ $TJacobian }} {
	var tmp {{ $TJacobian }}
	tmp.Set(q)
	tmp.Y.Neg(&tmp.Y)
	p.AddAssign(&tmp)
	return p
}

// Double sets p to [2]q in Jacobian coordinates.
//
// https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#doubling-dbl-2007-bl
func (p *{{ $TJacobian }}) DoubleMixed(a *{{ $TAffine }}) *{{ $TJacobian }} {
	var XX, YY, YYYY, S, M, T {{.CoordType}}
	XX.Square(&a.X)
	YY.Square(&a.Y)
	YYYY.Square(&YY)
	S.Add(&a.X, &YY).
		Square(&S).
		Sub(&S, &XX).
		Sub(&S, &YYYY).
		Double(&S)
	M.Double(&XX).
		Add(&M, &XX) // -> + A, but A=0 here
	T.Square(&M).
		Sub(&T, &S).
		Sub(&T, &S)
	p.X.Set(&T)
	p.Y.Sub(&S, &T).
		Mul(&p.Y, &M)
	YYYY.Double(&YYYY).
		Double(&YYYY).
		Double(&YYYY)
	p.Y.Sub(&p.Y, &YYYY)
	p.Z.Double(&a.Y)

	return p
}

// AddMixed sets p to p+a in Jacobian coordinates, where a.Z = 1.
//
// http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-madd-2007-bl
func (p *{{ $TJacobian }}) AddMixed(a *{{ $TAffine }}) *{{ $TJacobian }} {

	//if a is infinity return p
	if a.IsInfinity() {
		return p
	}
	// p is infinity, return a
	if p.Z.IsZero() {
		p.X = a.X
		p.Y = a.Y
		p.Z.SetOne()
		return p
	}

	var Z1Z1, U2, S2, H, HH, I, J, r, V {{.CoordType}}
	Z1Z1.Square(&p.Z)
	U2.Mul(&a.X, &Z1Z1)
	S2.Mul(&a.Y, &p.Z).
		Mul(&S2, &Z1Z1)

	// if p == a, we double instead
	if U2.Equal(&p.X) && S2.Equal(&p.Y) {
		return p.DoubleMixed(a)
	}

	H.Sub(&U2, &p.X)
	HH.Square(&H)
	I.Double(&HH).Double(&I)
	J.Mul(&H, &I)
	r.Sub(&S2, &p.Y).Double(&r)
	V.Mul(&p.X, &I)
	p.X.Square(&r).
		Sub(&p.X, &J).
		Sub(&p.X, &V).
		Sub(&p.X, &V)
	J.Mul(&J, &p.Y).Double(&J)
	p.Y.Sub(&V, &p.X).
		Mul(&p.Y, &r)
	p.Y.Sub(&p.Y, &J)
	p.Z.Add(&p.Z, &H)
	p.Z.Square(&p.Z).
		Sub(&p.Z, &Z1Z1).
		Sub(&p.Z, &HH)

	return p
}

// Double sets p to [2]q in Jacobian coordinates.
//
// https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#doubling-dbl-2007-bl
func (p *{{ $TJacobian }}) Double(q *{{ $TJacobian }}) *{{ $TJacobian }} {
	p.Set(q)
	p.DoubleAssign()
	return p
}

// DoubleAssign doubles p in Jacobian coordinates.
//
// https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#doubling-dbl-2007-bl
func (p *{{ $TJacobian }}) DoubleAssign() *{{ $TJacobian }} {

	var XX, YY, YYYY, ZZ, S, M, T {{.CoordType}}

	XX.Square(&p.X)
	YY.Square(&p.Y)
	YYYY.Square(&YY)
	ZZ.Square(&p.Z)
	S.Add(&p.X, &YY)
	S.Square(&S).
		Sub(&S, &XX).
		Sub(&S, &YYYY).
		Double(&S)
	M.Double(&XX).Add(&M, &XX)
	p.Z.Add(&p.Z, &p.Y).
		Square(&p.Z).
		Sub(&p.Z, &YY).
		Sub(&p.Z, &ZZ)
	T.Square(&M)
	p.X = T
	T.Double(&S)
	p.X.Sub(&p.X, &T)
	p.Y.Sub(&S, &p.X).
		Mul(&p.Y, &M)
	YYYY.Double(&YYYY).Double(&YYYY).Double(&YYYY)
	p.Y.Sub(&p.Y, &YYYY)

	return p
}

// ScalarMultiplication computes and returns p = [s]a
// where p and a are Jacobian points.
{{- if .GLV}}
// using the GLV technique.
// see https://www.iacr.org/archive/crypto2001/21390189.pdf
{{- else }}
// using a 2-bits windowed double-and-add method.
{{- end }}
func (p *{{ $TJacobian }}) ScalarMultiplication(q *{{ $TJacobian }}, s *big.Int) *{{ $TJacobian }} {
	{{- if .GLV}}
		return p.mulGLV(q, s)
	{{- else }}
		return p.mulWindowed(q, s)
	{{- end }}
}

// ScalarMultiplicationBase computes and returns p = [s]g
// where g is the prime subgroup generator.
func (p *{{ $TJacobian  }}) ScalarMultiplicationBase(s *big.Int) *{{ $TJacobian  }} {
    {{- if .GLV}}
        return p.mulGLV(&{{ toLower .PointName }}Gen, s)
    {{- else }}
        _p.mulWindowed(&{{ toLower .PointName }}Gen, s)
    {{- end }}

}

// String converts p to affine coordinates and returns its string representation E(x,y) or "O" if it is infinity.
func (p *{{ $TJacobian }}) String() string {
	_p := {{ $TAffine }}{}
	_p.FromJacobian(p)
	return _p.String()
}

// FromAffine converts a point a from affine to Jacobian coordinates.
func (p *{{ $TJacobian }}) FromAffine(a *{{ $TAffine }}) *{{ $TJacobian }} {
	if a.IsInfinity() {
		p.Z.SetZero()
		p.X.SetOne()
		p.Y.SetOne()
		return p
	}
	p.Z.SetOne()
	p.X.Set(&a.X)
	p.Y.Set(&a.Y)
	return p
}

// IsOnCurve returns true if the Jacobian point p in on the curve.
func (p *{{ $TJacobian }}) IsOnCurve() bool {
	var left, right, tmp, ZZ  {{.CoordType}}
	left.Square(&p.Y)
	right.Square(&p.X).Mul(&right, &p.X)
	ZZ.Square(&p.Z)
    tmp.Square(&ZZ).Mul(&tmp, &ZZ)
		{{- if eq .PointName "g1"}}
            {{- if or (eq .Name "bls12-381") (eq .Name "bls24-317") (eq .Name "bw6-633")}}
                // Mul tmp by bCurveCoeff=4
                tmp.Double(&tmp).Double(&tmp)
            {{- else if eq .Name "bn254"}}
                // Mul tmp by bCurveCoeff=3
                fp.MulBy3(&tmp)
            {{- else if eq .Name "bw6-761"}}
                // Mul tmp by bCurveCoeff=-1
                tmp.Neg(&tmp)
            {{- else if or (eq .Name "bls12-377") (eq .Name "bls12-378") (eq .Name "bls24-315") (eq .Name "bw6-756")}}
                // Mul tmp by bCurveCoeff=1 (nothing to do)
            {{- else}}
                tmp.Mul(&tmp, &bCurveCoeff)
            {{- end}}
		{{- else}}
            {{- if or (eq .Name "bls12-377") (eq .Name "bls12-381") (eq .Name "bn254") (eq .Name "bls24-315") (eq .Name "bls24-317")}}
                tmp.MulBybTwistCurveCoeff(&tmp)
            {{- else if eq .Name "bw6-761"}}
                // Mul tmp by bTwistCurveCoeff=4
                tmp.Double(&tmp).Double(&tmp)
            {{- else if eq .Name "bw6-633"}}
                // Mul tmp by bTwistCurveCoeff=8
                tmp.Double(&tmp).Double(&tmp).Double(&tmp)
            {{- else}}
                tmp.Mul(&tmp, &bTwistCurveCoeff)
            {{- end}}
		{{- end}}
	right.Add(&right, &tmp)
	return left.Equal(&right)
}



{{- if or (eq .Name "bn254") (eq .Name "secp256k1")}}
	{{- if eq .PointName "g1"}}
		// IsInSubGroup returns true if p is on the r-torsion, false otherwise.
        // the curve is of prime order i.e. E(𝔽p) is the full group
        // so we just check that the point is on the curve.
		func (p *{{ $TJacobian }}) IsInSubGroup() bool {

			return p.IsOnCurve()

		}
	{{else if eq .PointName "g2"}}
		// IsInSubGroup returns true if p is on the r-torsion, false otherwise.
        // https://eprint.iacr.org/2022/348.pdf, sec. 3 and 5.1
        // [r]P == 0 <==> [x₀+1]P + ψ([x₀]P) + ψ²([x₀]P) = ψ³([2x₀]P)
		func (p *{{ $TJacobian }}) IsInSubGroup() bool {
            var a, b, c, res G2Jac
            a.ScalarMultiplication(p, &xGen)
            b.psi(&a)
            a.AddAssign(p)
            res.psi(&b)
            c.Set(&res).
                AddAssign(&b).
                AddAssign(&a)
            res.psi(&res).
                Double(&res).
            SubAssign(&c)

            return res.IsOnCurve() && res.Z.IsZero()
		}
	{{- end}}
{{else if or (eq .Name "bw6-761") (eq .Name "bw6-756")}}
	// IsInSubGroup returns true if p is on the r-torsion, false otherwise.
    {{ if .GLV}}
	// Z[r,0]+Z[-lambda{{ $TAffine }}, 1] is the kernel
	// of (u,v)->u+lambda{{ $TAffine }}v mod r. Expressing r, lambda{{ $TAffine }} as
	// polynomials in x, a short vector of this Zmodule is
	// (x+1), (x³-x²+1). So we check that (x+1)p+(x³-x²+1)ϕ(p)
	// is the infinity.
	func (p *{{ $TJacobian }}) IsInSubGroup() bool {

		var res, phip {{ $TJacobian }}
		phip.phi(p)
		res.ScalarMultiplication(&phip, &xGen).
			SubAssign(&phip).
			ScalarMultiplication(&res, &xGen).
			ScalarMultiplication(&res, &xGen).
			AddAssign(&phip)

		phip.ScalarMultiplication(p, &xGen).AddAssign(p).AddAssign(&res)

		return phip.IsOnCurve() && phip.Z.IsZero()
    {{ else}}
	func (p *{{ $TJacobian }}) IsInSubGroup() bool {

		var res {{ $TJacobian }}
        res.ScalarMultiplication(p, fr.Modulus())
		return res.IsOnCurve() && res.Z.IsZero()
    {{ end}}

	}
{{else if eq .Name "bw6-633"}}
    // IsInSubGroup returns true if p is on the r-torsion, false otherwise.
    {{ if .GLV}}
    // 3r P = (x+1)ϕ(P) + (-x^5 + x⁴ + x)P
	func (p *{{ $TJacobian }}) IsInSubGroup() bool {

        var uP, u4P, u5P, q, r {{ $TJacobian }}
        uP.ScalarMultiplication(p, &xGen)
        u4P.ScalarMultiplication(&uP, &xGen).
            ScalarMultiplication(&u4P, &xGen).
            ScalarMultiplication(&u4P, &xGen)
        u5P.ScalarMultiplication(&u4P, &xGen)
        q.Set(p).SubAssign(&uP)
        r.phi(&q).SubAssign(&uP).
            AddAssign(&u4P).
            AddAssign(&u5P)

    {{ else}}
	func (p *{{ $TJacobian }}) IsInSubGroup() bool {

        var r {{ $TJacobian }}
        r.ScalarMultiplication(p, fr.Modulus())
    {{ end}}
		return r.IsOnCurve() && r.Z.IsZero()
	}
{{else if or (eq .Name "bls24-315") (eq .Name "bls24-317")}}
	{{- if eq .PointName "g1"}}
        // IsInSubGroup returns true if p is on the r-torsion, false otherwise.
        {{ if .GLV}}
        // Z[r,0]+Z[-lambda{{ $TAffine }}, 1] is the kernel
        // of (u,v)->u+lambda{{ $TAffine }}v mod r. Expressing r, lambda{{ $TAffine }} as
        // polynomials in x, a short vector of this Zmodule is
        // 1, x⁴. So we check that p+x⁴ϕ(p)
        // is the infinity.
        func (p *{{ $TJacobian }}) IsInSubGroup() bool {

            var res {{ $TJacobian }}
            res.phi(p).
                ScalarMultiplication(&res, &xGen).
                ScalarMultiplication(&res, &xGen).
                ScalarMultiplication(&res, &xGen).
                ScalarMultiplication(&res, &xGen).
                AddAssign(p)
    {{ else}}
	func (p *{{ $TJacobian }}) IsInSubGroup() bool {

        var res {{ $TJacobian }}
        res.ScalarMultiplication(p, fr.Modulus())
		return res.IsOnCurve() && res.Z.IsZero()
    {{ end}}

            return res.IsOnCurve() && res.Z.IsZero()

        }
	{{else if eq .PointName "g2"}}
        // IsInSubGroup returns true if p is on the r-torsion, false otherwise.
        // https://eprint.iacr.org/2021/1130.pdf, sec.4
        // and https://eprint.iacr.org/2022/352.pdf, sec. 4.2
        // ψ(p) = [x₀]P
        func (p *{{ $TJacobian }}) IsInSubGroup() bool {
            var res, tmp {{ $TJacobian }}
            tmp.psi(p)
            res.ScalarMultiplication(p, &xGen).
            {{ if eq .Name "bls24-315"}}
                AddAssign(&tmp)
            {{ else }}
                SubAssign(&tmp)
            {{ end }}

            return res.IsOnCurve() && res.Z.IsZero()

        }
    {{- end}}
{{else}}
	{{- if eq .PointName "g1"}}
        // IsInSubGroup returns true if p is on the r-torsion, false otherwise.
        // Z[r,0]+Z[-lambda{{ $TAffine }}, 1] is the kernel
        // of (u,v)->u+lambda{{ $TAffine }}v mod r. Expressing r, lambda{{ $TAffine }} as
        // polynomials in x, a short vector of this Zmodule is
        // 1, x². So we check that p+x²ϕ(p)
        // is the infinity.
        func (p *{{ $TJacobian }}) IsInSubGroup() bool {

            var res {{ $TJacobian }}
            {{ if .GLV}}
            res.phi(p).
                ScalarMultiplication(&res, &xGen).
                ScalarMultiplication(&res, &xGen).
                AddAssign(p)
            {{ else}}
            res.ScalarMultiplication(p, fr.Modulus())
            {{ end}}

            return res.IsOnCurve() && res.Z.IsZero()

        }
	{{else if eq .PointName "g2"}}
        {{if eq .Name "bls12-381"}}
            // IsInSubGroup returns true if p is on the r-torsion, false otherwise.
            // https://eprint.iacr.org/2021/1130.pdf, sec.4
            // and https://eprint.iacr.org/2022/352.pdf, sec. 4.2
            // ψ(p) = [x₀]P
            func (p *{{ $TJacobian }}) IsInSubGroup() bool {
                var res, tmp {{ $TJacobian }}
                tmp.psi(p)
                res.ScalarMultiplication(p, &xGen).
                    AddAssign(&tmp)

                return res.IsOnCurve() && res.Z.IsZero()
        }
        {{else if or (eq .Name "bls12-377") (eq .Name "bls12-378")}}
            // https://eprint.iacr.org/2021/1130.pdf, sec.4
            // and https://eprint.iacr.org/2022/352.pdf, sec. 4.2
            // ψ(p) = [x₀]P
            func (p *{{ $TJacobian }}) IsInSubGroup() bool {
                var res, tmp {{ $TJacobian }}
                tmp.psi(p)
                res.ScalarMultiplication(p, &xGen).
                    SubAssign(&tmp)

                return res.IsOnCurve() && res.Z.IsZero()
        }
        {{- end}}
	{{- end}}
{{- end}}


// mulWindowed computes the 2-bits windowed double-and-add scalar
// multiplication p=[s]q in Jacobian coordinates.
func (p *{{ $TJacobian }}) mulWindowed(q *{{ $TJacobian }}, s *big.Int) *{{ $TJacobian }} {

	var res {{ $TJacobian }}
	var ops [3]{{ $TJacobian }}

	ops[0].Set(q)
	if s.Sign() == -1 {
		ops[0].Neg(&ops[0])
	}
	res.Set(&{{ toLower .PointName}}Infinity)
	ops[1].Double(&ops[0])
	ops[2].Set(&ops[0]).AddAssign(&ops[1])

	b := s.Bytes()
	for i := range b {
		w := b[i]
		mask := byte(0xc0)
		for j := 0; j < 4; j++ {
			res.DoubleAssign().DoubleAssign()
			c := (w & mask) >> (6 - 2*j)
			if c != 0 {
				res.AddAssign(&ops[c-1])
			}
			mask = mask >> 2
		}
	}
	p.Set(&res)

	return p

}

{{ if eq .CoordType "fptower.E2"  }}
	// psi sets p to ψ(q) = u o π o u⁻¹ where u:E'→E is the isomorphism from the twist to the curve E and π is the Frobenius map.
	func (p *{{ $TJacobian }}) psi(q *{{ $TJacobian }}) *{{ $TJacobian }} {
		p.Set(q)
		p.X.Conjugate(&p.X).Mul(&p.X, &endo.u)
		p.Y.Conjugate(&p.Y).Mul(&p.Y, &endo.v)
		p.Z.Conjugate(&p.Z)
		return p
	}
{{ else if eq .CoordType "fptower.E4"}}
	// psi sets p to ψ(q) = u o π o u⁻¹ where u:E'→E is the isomorphism from the twist to the curve E and π is the Frobenius map.
	func (p *{{ $TJacobian }}) psi(q *{{ $TJacobian }}) *{{ $TJacobian }} {
		p.Set(q)
		p.X.Frobenius(&p.X).Mul(&p.X, &endo.u)
		p.Y.Frobenius(&p.Y).Mul(&p.Y, &endo.v)
		p.Z.Frobenius(&p.Z)
		return p
	}
{{ end }}

{{ if .GLV}}

// phi sets p to ϕ(a) where ϕ: (x,y) → (w x,y),
// where w is a third root of unity.
func (p *{{ $TJacobian }}) phi(q *{{ $TJacobian }}) *{{ $TJacobian }} {
	p.Set(q)
	{{- if or (eq .CoordType "fptower.E2" ) (eq .CoordType "fptower.E4" )}}
		p.X.MulByElement(&p.X, &thirdRootOne{{toUpper .PointName}})
	{{- else}}
		p.X.Mul(&p.X, &thirdRootOne{{toUpper .PointName}})
	{{- end}}
	return p
}

// mulGLV computes the scalar multiplication using a windowed-GLV method
//
// see https://www.iacr.org/archive/crypto2001/21390189.pdf
func (p *{{ $TJacobian }}) mulGLV(q *{{ $TJacobian }}, s *big.Int) *{{ $TJacobian }} {

	var table [15]{{ $TJacobian }}
	var res {{ $TJacobian }}
	var k1, k2 fr.Element

	res.Set(&{{ toLower .PointName}}Infinity)

	// table[b3b2b1b0-1] = b3b2 ⋅ ϕ(q) + b1b0*q
	table[0].Set(q)
	table[3].phi(q)

	// split the scalar, modifies ±q, ϕ(q) accordingly
	k := ecc.SplitScalar(s, &glvBasis)

	if k[0].Sign() == -1 {
		k[0].Neg(&k[0])
		table[0].Neg(&table[0])
	}
	if k[1].Sign() == -1 {
		k[1].Neg(&k[1])
		table[3].Neg(&table[3])
	}

	// precompute table (2 bits sliding window)
	// table[b3b2b1b0-1] = b3b2 ⋅ ϕ(q) + b1b0 ⋅ q if b3b2b1b0 != 0
	table[1].Double(&table[0])
	table[2].Set(&table[1]).AddAssign(&table[0])
	table[4].Set(&table[3]).AddAssign(&table[0])
	table[5].Set(&table[3]).AddAssign(&table[1])
	table[6].Set(&table[3]).AddAssign(&table[2])
	table[7].Double(&table[3])
	table[8].Set(&table[7]).AddAssign(&table[0])
	table[9].Set(&table[7]).AddAssign(&table[1])
	table[10].Set(&table[7]).AddAssign(&table[2])
	table[11].Set(&table[7]).AddAssign(&table[3])
	table[12].Set(&table[11]).AddAssign(&table[0])
	table[13].Set(&table[11]).AddAssign(&table[1])
	table[14].Set(&table[11]).AddAssign(&table[2])

	// bounds on the lattice base vectors guarantee that k1, k2 are len(r)/2 or len(r)/2+1 bits long max
	// this is because we use a probabilistic scalar decomposition that replaces a division by a right-shift
	k1 = k1.SetBigInt(&k[0]).Bits()
	k2 = k2.SetBigInt(&k[1]).Bits()

    // we don't target constant-timeness so we check first if we increase the bounds or not
	maxBit := k1.BitLen()
	if k2.BitLen() > maxBit {
		maxBit = k2.BitLen()
	}
	hiWordIndex := (maxBit - 1) / 64

	// loop starts from len(k1)/2 or len(k1)/2+1 due to the bounds
	for i := hiWordIndex; i >= 0; i-- {
		mask := uint64(3) << 62
		for j := 0; j < 32; j++ {
			res.Double(&res).Double(&res)
			b1 := (k1[i] & mask) >> (62 - 2*j)
			b2 := (k2[i] & mask) >> (62 - 2*j)
			if b1|b2 != 0 {
				s := (b2<<2 | b1)
				res.AddAssign(&table[s-1])
			}
			mask = mask >> 2
		}
	}

	p.Set(&res)
	return p
}


{{ end }}

{{ if .CofactorCleaning}}

// ClearCofactor maps a point in curve to r-torsion
func (p *{{ $TAffine }}) ClearCofactor(a *{{ $TAffine }}) *{{ $TAffine }} {
	var _p {{$TJacobian}}
	_p.FromAffine(a)
	_p.ClearCofactor(&_p)
	p.FromJacobian(&_p)
	return p
}

{{- if eq .PointName "g1"}}

// ClearCofactor maps a point in E(Fp) to E(Fp)[r]
func (p *{{$TJacobian}}) ClearCofactor(q *{{$TJacobian}}) *{{$TJacobian}} {
{{- if or (eq .Name "bls12-381") (eq .Name "bls24-315")}}
	// cf https://eprint.iacr.org/2019/403.pdf, 5
	var res {{$TJacobian}}
	res.ScalarMultiplication(q, &xGen).AddAssign(q)
	p.Set(&res)
	return p
{{else if or (eq .Name "bls12-377") (eq .Name "bls12-378") (eq .Name "bls24-317")}}
	// cf https://eprint.iacr.org/2019/403.pdf, 5
	var res {{$TJacobian}}
	res.ScalarMultiplication(q, &xGen).Neg(&res).AddAssign(q)
	p.Set(&res)
	return p
{{else if eq .Name "bw6-761"}}
{{ if .GLV}}
	// https://eprint.iacr.org/2020/351.pdf
	var points [4]{{$TJacobian}}
	points[0].Set(q)
	points[1].ScalarMultiplication(q, &xGen)
	points[2].ScalarMultiplication(&points[1], &xGen)
	points[3].ScalarMultiplication(&points[2], &xGen)

	var scalars [7]big.Int
	scalars[0].SetInt64(103)
	scalars[1].SetInt64(83)
	scalars[2].SetInt64(40)
	scalars[3].SetInt64(136)

	scalars[4].SetInt64(7)
	scalars[5].SetInt64(89)
	scalars[6].SetInt64(130)

	var p1, p2, tmp {{$TJacobian}}
	p1.ScalarMultiplication(&points[3], &scalars[0])
	tmp.ScalarMultiplication(&points[2], &scalars[1]).Neg(&tmp)
	p1.AddAssign(&tmp)
	tmp.ScalarMultiplication(&points[1], &scalars[2]).Neg(&tmp)
	p1.AddAssign(&tmp)
	tmp.ScalarMultiplication(&points[0], &scalars[3])
	p1.AddAssign(&tmp)

	p2.ScalarMultiplication(&points[2], &scalars[4])
	tmp.ScalarMultiplication(&points[1], &scalars[5])
	p2.AddAssign(&tmp)
	tmp.ScalarMultiplication(&points[0], &scalars[6])
	p2.AddAssign(&tmp)
	p2.phi(&p2)

	p.Set(&p1).AddAssign(&p2)
{{ else}}
    var c1 big.Int
    c1.SetString("26642435879335816683987677701488073867751118270052650655942102502312977592501693353047140953112195348280268661194876", 10)
    p.ScalarMultiplication(q, &c1)
{{ end}}

	return p
{{else if eq .Name "bw6-633"}}
{{ if .GLV}}
	var uP, vP, wP, L0, L1, tmp {{$TJacobian}}
	var v, one, uPlusOne, uMinusOne, d1, d2, ht big.Int
	one.SetInt64(1)
	uPlusOne.Add(&one, &xGen)
	uMinusOne.Sub(&xGen, &one)
	d1.SetInt64(13)
	d2.SetInt64(5)
	ht.SetInt64(7)
	v.Mul(&xGen, &xGen).Add(&v, &one).Mul(&v, &uPlusOne)

	uP.ScalarMultiplication(q, &xGen).Neg(&uP)
	vP.Set(q).SubAssign(&uP).
        ScalarMultiplication(&vP, &v)
	wP.ScalarMultiplication(&vP, &uMinusOne).Neg(&wP).
        AddAssign(&uP)
	L0.ScalarMultiplication(&wP, &d1)
	tmp.ScalarMultiplication(&vP, &ht)
	L0.AddAssign(&tmp)
	tmp.Double(q)
	L0.AddAssign(&tmp)
	L1.Set(&uP).AddAssign(q).ScalarMultiplication(&L1, &d1)
	tmp.ScalarMultiplication(&vP, &d2)
	L1.AddAssign(&tmp)
	tmp.ScalarMultiplication(q, &ht)
	L1.AddAssign(&tmp)

	p.phi(&L1).AddAssign(&L0)
{{ else}}
    var c1 big.Int
    c1.SetString("516166855112631370346774477030598579858367278343565509012644853411927535599366632765988905418773", 10)
    p.ScalarMultiplication(q, &c1)
{{ end}}

    return p
{{else if eq .Name "bw6-756"}}
{{ if .GLV}}
	var L0, L1, uP, u2P, u3P, tmp G1Jac

	uP.ScalarMultiplication(q, &xGen)
	u2P.ScalarMultiplication(&uP, &xGen)
	u3P.ScalarMultiplication(&u2P, &xGen)

	L0.Set(q).AddAssign(&u3P).
		SubAssign(&u2P)
	tmp.Set(q).AddAssign(&u2P).
		SubAssign(&uP).
		SubAssign(&uP).
		Double(&tmp)
	L0.SubAssign(&tmp).
		SubAssign(q)

	L1.Set(q).AddAssign(&uP)
	tmp.Set(&uP).SubAssign(q).
		Double(&tmp).
		SubAssign(&u2P)
	L1.AddAssign(&tmp).
		SubAssign(q)

	p.phi(&L1).
		AddAssign(&L0)
{{ else}}
    var c1 big.Int
    c1.SetString("605248206075306171568857128027361794400937215108643640003009340657451546212610770151705515081537938829431808196608", 10)
    p.ScalarMultiplication(q, &c1)
{{ end}}

	return p
{{- end}}
}
{{ else }}

// ClearCofactor maps a point in curve to r-torsion
func (p *{{$TJacobian}}) ClearCofactor(q *{{$TJacobian}}) *{{$TJacobian}} {
{{- if eq .Name "bn254"}}
	// cf http://cacr.uwaterloo.ca/techreports/2011/cacr2011-26.pdf, 6.1
	var points [4]{{$TJacobian}}

	points[0].ScalarMultiplication(q, &xGen)

	points[1].Double(&points[0]).
		AddAssign(&points[0]).
		psi(&points[1])

	points[2].psi(&points[0]).
		psi(&points[2])

	points[3].psi(q).psi(&points[3]).psi(&points[3])

	var res {{$TJacobian}}
	res.Set(&g2Infinity)
	for i := 0; i < 4; i++ {
		res.AddAssign(&points[i])
	}
	p.Set(&res)
	return p
{{else if eq .Name "bls12-381"}}
	// https://eprint.iacr.org/2017/419.pdf, 4.1
	var xg, xxg, res, t G2Jac
	xg.ScalarMultiplication(q, &xGen).Neg(&xg)
	xxg.ScalarMultiplication(&xg, &xGen).Neg(&xxg)

	res.Set(&xxg).
		SubAssign(&xg).
		SubAssign(q)

	t.Set(&xg).
		SubAssign(q).
		psi(&t)

	res.AddAssign(&t)

	t.Double(q)
	t.X.MulByElement(&t.X, &thirdRootOneG1)

	res.SubAssign(&t)

	p.Set(&res)

	return p

{{else if or (eq .Name "bls12-377") (eq .Name "bls12-378")}}
    // https://eprint.iacr.org/2017/419.pdf, 4.1
	var xg, xxg, res, t G2Jac
	xg.ScalarMultiplication(q, &xGen)
	xxg.ScalarMultiplication(&xg, &xGen)

	res.Set(&xxg).
		SubAssign(&xg).
		SubAssign(q)

	t.Set(&xg).
		SubAssign(q).
		psi(&t)

	res.AddAssign(&t)

	t.Double(q)
	t.X.MulByElement(&t.X, &thirdRootOneG1)

	res.SubAssign(&t)

	p.Set(&res)

	return p

{{else if or (eq .Name "bls24-315") (eq .Name "bls24-317")}}
	// https://eprint.iacr.org/2017/419.pdf, section 4.2
	// multiply by (3x⁴-3)*cofacor
    {{ if eq .Name "bls24-315"}}
	var xg, xxg, xxxg, xxxxg, res, t G2Jac
	xg.ScalarMultiplication(q, &xGen).Neg(&xg).SubAssign(q)
	xxg.ScalarMultiplication(&xg, &xGen).Neg(&xxg)
	xxxg.ScalarMultiplication(&xxg, &xGen).Neg(&xxxg)
	xxxxg.ScalarMultiplication(&xxxg, &xGen).Neg(&xxxxg)
    {{ else }}
	var xg, xxg, xxxg, xxxxg, res, t G2Jac
	xg.ScalarMultiplication(q, &xGen).SubAssign(q)
	xxg.ScalarMultiplication(&xg, &xGen)
	xxxg.ScalarMultiplication(&xxg, &xGen)
	xxxxg.ScalarMultiplication(&xxxg, &xGen)
    {{ end }}

	res.Set(&xxxxg).
		SubAssign(q)

	t.Set(&xxxg).
		psi(&t).
		AddAssign(&res)

	res.Set(&xxg).
		psi(&res).
		psi(&res).
		AddAssign(&t)

	t.Set(&xg).
		psi(&t).
		psi(&t).
		psi(&t).
		AddAssign(&res)

	res.Double(q).
		psi(&res).
		psi(&res).
		psi(&res).
		psi(&res).
		AddAssign(&t)

	p.Set(&res)

	return p
{{else if eq .Name "bw6-761"}}

{{- if .GLV}}
	var points [4]{{$TJacobian}}
	points[0].Set(q)
	points[1].ScalarMultiplication(q, &xGen)
	points[2].ScalarMultiplication(&points[1], &xGen)
	points[3].ScalarMultiplication(&points[2], &xGen)

	var scalars [7]big.Int
	scalars[0].SetInt64(103)
	scalars[1].SetInt64(83)
	scalars[2].SetInt64(143)
	scalars[3].SetInt64(27)

	scalars[4].SetInt64(7)
	scalars[5].SetInt64(117)
	scalars[6].SetInt64(109)

	var p1, p2, tmp {{$TJacobian}}
	p1.ScalarMultiplication(&points[3], &scalars[0])
	tmp.ScalarMultiplication(&points[2], &scalars[1]).Neg(&tmp)
	p1.AddAssign(&tmp)
	tmp.ScalarMultiplication(&points[1], &scalars[2]).Neg(&tmp)
	p1.AddAssign(&tmp)
	tmp.ScalarMultiplication(&points[0], &scalars[3])
	p1.AddAssign(&tmp)

	p2.ScalarMultiplication(&points[2], &scalars[4])
	tmp.ScalarMultiplication(&points[1], &scalars[5]).Neg(&tmp)
	p2.AddAssign(&tmp)
	tmp.ScalarMultiplication(&points[0], &scalars[6]).Neg(&tmp)
	p2.AddAssign(&tmp)
	p2.phi(&p2).phi(&p2)

	p.Set(&p1).AddAssign(&p2)
{{else}}
    var c2 big.Int
    c2.SetString("26642435879335816683987677701488073867751118270052650655942102502312977592501693353047140953112195348280268661194869", 10)
    p.ScalarMultiplication(q, &c2)
{{end}}

	return p
{{else if eq .Name "bw6-633"}}
{{- if .GLV}}
	var uP, u2P, u3P, u4P, u5P, xP, vP, wP, L0, L1, tmp {{$TJacobian}}
	var ht, d1, d3 big.Int
	ht.SetInt64(7) // negative
	d1.SetInt64(13)
	d3.SetInt64(5) // negative

	uP.ScalarMultiplication(q, &xGen) // negative
	u2P.ScalarMultiplication(&uP, &xGen)
	u3P.ScalarMultiplication(&u2P, &xGen) // negative
	u4P.ScalarMultiplication(&u3P, &xGen)
	u5P.ScalarMultiplication(&u4P, &xGen) // negative
	vP.Set(&u2P).AddAssign(&uP).
	    AddAssign(&u3P).
	    Double(&vP).
	    AddAssign(&u4P).
	    AddAssign(q)
	wP.Set(&uP).SubAssign(&u4P).SubAssign(&u5P)
    xP.Set(q).AddAssign(&vP)
	L0.Set(&uP).SubAssign(q).ScalarMultiplication(&L0, &d1)
	tmp.ScalarMultiplication(&xP, &d3)
	L0.AddAssign(&tmp)
    tmp.ScalarMultiplication(q, &ht) // negative
	L0.SubAssign(&tmp)
	L1.ScalarMultiplication(&wP, &d1)
	tmp.ScalarMultiplication(&vP, &ht)
	L1.AddAssign(&tmp)
	tmp.ScalarMultiplication(q, &d3)
	L1.AddAssign(&tmp)

	p.phi(&L1).AddAssign(&L0)
{{else}}
    var c2 big.Int
    c2.SetString("516166855112631370346774477030598579858367278343565509012644853411927535599366632765988905418768", 10)
    p.ScalarMultiplication(q, &c2)
{{end}}

	return p
{{else if eq .Name "bw6-756"}}

{{- if .GLV}}
	var L0, L1, uP, u2P, u3P, tmp G2Jac

	uP.ScalarMultiplication(q, &xGen)
	u2P.ScalarMultiplication(&uP, &xGen)
	u3P.ScalarMultiplication(&u2P, &xGen)
	// ht=-2, hy=0
	// d1=1, d2=-1, d3=-1

	L0.Set(q).
		AddAssign(&u2P).
		SubAssign(&uP)
	tmp.Set(&u2P).
		AddAssign(q).
		SubAssign(&uP).
		Double(&tmp)
	L1.Set(&u3P).
		SubAssign(&tmp)

	p.phi(&L0).
		AddAssign(&L1)
{{else}}
    var c2 big.Int
    c2.SetString("605248206075306171568857128027361794400937215108643640003009340657451546212610770151705515081537938829431808196609", 10)
    p.ScalarMultiplication(q, &c2)
{{end}}

	return p

{{- end}}
}
{{- end}}


{{ end }}

{{ if eq .PointName "g1" }}
// JointScalarMultiplication computes [s1]a1+[s2]a2 using Strauss-Shamir technique
// where a1 and a2 are affine points.
func (p *{{$TJacobian}}) JointScalarMultiplication(a1, a2 *G1Affine, s1, s2 *big.Int) *{{$TJacobian}} {

	var res, p1, p2 {{$TJacobian}}
	res.Set(&{{ toLower .PointName }}Infinity)
	p1.FromAffine(a1)
	p2.FromAffine(a2)

	var table [15]{{$TJacobian}}

	var k1, k2 big.Int
	if s1.Sign() == -1 {
		k1.Neg(s1)
		table[0].Neg(&p1)
	} else {
		k1.Set(s1)
		table[0].Set(&p1)
	}
	if s2.Sign() == -1 {
		k2.Neg(s2)
		table[3].Neg(&p2)
	} else {
		k2.Set(s2)
		table[3].Set(&p2)
	}

	// precompute table (2 bits sliding window)
	table[1].Double(&table[0])
	table[2].Set(&table[1]).AddAssign(&table[0])
	table[4].Set(&table[3]).AddAssign(&table[0])
	table[5].Set(&table[3]).AddAssign(&table[1])
	table[6].Set(&table[3]).AddAssign(&table[2])
	table[7].Double(&table[3])
	table[8].Set(&table[7]).AddAssign(&table[0])
	table[9].Set(&table[7]).AddAssign(&table[1])
	table[10].Set(&table[7]).AddAssign(&table[2])
	table[11].Set(&table[7]).AddAssign(&table[3])
	table[12].Set(&table[11]).AddAssign(&table[0])
	table[13].Set(&table[11]).AddAssign(&table[1])
	table[14].Set(&table[11]).AddAssign(&table[2])

	var s [2]fr.Element
	s[0] = s[0].SetBigInt(&k1).Bits()
	s[1] = s[1].SetBigInt(&k2).Bits()

	maxBit := k1.BitLen()
	if k2.BitLen() > maxBit {
		maxBit = k2.BitLen()
	}
	hiWordIndex := (maxBit - 1) / 64

	for i := hiWordIndex; i >= 0; i-- {
		mask := uint64(3) << 62
		for j := 0; j < 32; j++ {
			res.Double(&res).Double(&res)
			b1 := (s[0][i] & mask) >> (62 - 2*j)
			b2 := (s[1][i] & mask) >> (62 - 2*j)
			if b1|b2 != 0 {
				s := (b2<<2 | b1)
				res.AddAssign(&table[s-1])
			}
			mask = mask >> 2
		}
	}

	p.Set(&res)
	return p

}

// JointScalarMultiplicationBase computes [s1]g+[s2]a using Straus-Shamir technique
// where g is the prime subgroup generator.
func (p *{{$TJacobian}}) JointScalarMultiplicationBase(a *G1Affine, s1, s2 *big.Int) *{{$TJacobian}} {
    return p.JointScalarMultiplication(&g1GenAff, a, s1, s2)

}

{{ end }}


// -------------------------------------------------------------------------------------------------
// extended Jacobian coordinates

// Set sets p to a in extended Jacobian coordinates.
func (p *{{ $TJacobianExtended }}) Set(q *{{ $TJacobianExtended }}) *{{ $TJacobianExtended }} {
	p.X, p.Y, p.ZZ, p.ZZZ = q.X, q.Y, q.ZZ, q.ZZZ
	return p
}

// setInfinity sets p to the infinity point (1,1,0,0).
func (p *{{ $TJacobianExtended }}) setInfinity() *{{ $TJacobianExtended }} {
	p.X.SetOne()
	p.Y.SetOne()
	p.ZZ = {{.CoordType}}{}
	p.ZZZ = {{.CoordType}}{}
	return p
}

// IsInfinity checks if the p is infinity, i.e. p.ZZ=0.
func (p *{{ $TJacobianExtended }}) IsInfinity() bool {
	return p.ZZ.IsZero()
}

// fromJacExtended converts an extended Jacobian point to an affine point.
func (p *{{ $TAffine }})  fromJacExtended (q *{{ $TJacobianExtended }}) *{{ $TAffine }} {
	if q.ZZ.IsZero() {
		p.X = {{.CoordType}}{}
		p.Y = {{.CoordType}}{}
		return p
	}
	p.X.Inverse(&q.ZZ).Mul(&p.X, &q.X)
	p.Y.Inverse(&q.ZZZ).Mul(&p.Y, &q.Y)
	return p
}

// fromJacExtended converts an extended Jacobian point to a Jacobian point.
func (p *{{ $TJacobian }}) fromJacExtended(q *{{ $TJacobianExtended }}) *{{ $TJacobian }} {
	if q.ZZ.IsZero() {
		p.Set(&{{ toLower .PointName }}Infinity)
		return p
	}
	p.X.Mul(&q.ZZ, &q.X).Mul(&p.X, &q.ZZ)
	p.Y.Mul(&q.ZZZ, &q.Y).Mul(&p.Y, &q.ZZZ)
	p.Z.Set(&q.ZZZ)
	return p
}

// unsafeFromJacExtended converts an extended Jacobian point, distinct from Infinity, to a Jacobian point.
func (p *{{ $TJacobian }}) unsafeFromJacExtended(q *{{ $TJacobianExtended }}) *{{ $TJacobian }} {
	p.X.Square(&q.ZZ).Mul(&p.X, &q.X)
	p.Y.Square(&q.ZZZ).Mul(&p.Y, &q.Y)
	p.Z = q.ZZZ
	return p
}

// add sets p to p+q in extended Jacobian coordinates.
//
// https://www.hyperelliptic.org/EFD/g1p/auto-shortw-xyzz.html#addition-add-2008-s
func (p *{{ $TJacobianExtended }}) add(q *{{ $TJacobianExtended }}) *{{ $TJacobianExtended }} {
	//if q is infinity return p
	if q.ZZ.IsZero() {
		return p
	}
	// p is infinity, return q
	if p.ZZ.IsZero() {
		p.Set(q)
		return p
	}

	var A, B, U1, U2, S1, S2 {{.CoordType}}

	// p2: q, p1: p
	U2.Mul(&q.X, &p.ZZ)
	U1.Mul(&p.X, &q.ZZ)
	A.Sub(&U2, &U1)
	S2.Mul(&q.Y, &p.ZZZ)
	S1.Mul(&p.Y, &q.ZZZ)
	B.Sub(&S2, &S1)

 	if A.IsZero() {
		if B.IsZero() {
				return p.double(q)

		}
		p.ZZ = {{.CoordType}}{}
		p.ZZZ = {{.CoordType}}{}
		return p
	}


	var P, R, PP, PPP, Q, V {{.CoordType}}
    P.Sub(&U2, &U1)
    R.Sub(&S2, &S1)
    PP.Square(&P)
    PPP.Mul(&P, &PP)
    Q.Mul(&U1, &PP)
    V.Mul(&S1, &PPP)

    p.X.Square(&R).
        Sub(&p.X, &PPP).
        Sub(&p.X, &Q).
        Sub(&p.X, &Q)
    p.Y.Sub(&Q, &p.X).
        Mul(&p.Y, &R).
        Sub(&p.Y, &V)
    p.ZZ.Mul(&p.ZZ, &q.ZZ).
        Mul(&p.ZZ, &PP)
    p.ZZZ.Mul(&p.ZZZ, &q.ZZZ).
        Mul(&p.ZZZ, &PPP)

    return p
}

// double sets p to [2]q in Jacobian extended coordinates.
//
// http://www.hyperelliptic.org/EFD/g1p/auto-shortw-xyzz.html#doubling-dbl-2008-s-1
// N.B.: since we consider any point on Z=0 as the point at infinity
// this doubling formula works for infinity points as well.
func (p *{{ $TJacobianExtended }}) double(q *{{ $TJacobianExtended }}) *{{ $TJacobianExtended }} {
	var U, V, W, S, XX, M {{.CoordType}}

	U.Double(&q.Y)
	V.Square(&U)
	W.Mul(&U, &V)
	S.Mul(&q.X, &V)
	XX.Square(&q.X)
	M.Double(&XX).
		Add(&M, &XX) // -> + A, but A=0 here
	U.Mul(&W, &q.Y)

	p.X.Square(&M).
		Sub(&p.X, &S).
		Sub(&p.X, &S)
	p.Y.Sub(&S, &p.X).
		Mul(&p.Y, &M).
		Sub(&p.Y, &U)
	p.ZZ.Mul(&V, &q.ZZ)
	p.ZZZ.Mul(&W, &q.ZZZ)

	return p
}

// addMixed sets p to p+q in extended Jacobian coordinates, where a.ZZ=1.
//
// http://www.hyperelliptic.org/EFD/g1p/auto-shortw-xyzz.html#addition-madd-2008-s
func (p *{{ $TJacobianExtended }}) addMixed(a *{{ $TAffine }}) *{{ $TJacobianExtended }} {
	{{ template "mAdd" dict "all" . "negate" false}}
}

// subMixed works the same as addMixed, but negates a.Y.
//
// http://www.hyperelliptic.org/EFD/g1p/auto-shortw-xyzz.html#addition-madd-2008-s
func (p *{{ $TJacobianExtended }}) subMixed(a *{{ $TAffine }}) *{{ $TJacobianExtended }} {
	{{ template "mAdd" dict "all" . "negate" true}}
}

// doubleNegMixed works the same as double, but negates q.Y.
func (p *{{ $TJacobianExtended }}) doubleNegMixed(a *{{ $TAffine }}) *{{ $TJacobianExtended }} {
	{{ template "mDouble" dict "all" . "negate" true}}
}

// doubleMixed sets p to [2]a in Jacobian extended coordinates, where a.ZZ=1.
//
// http://www.hyperelliptic.org/EFD/g1p/auto-shortw-xyzz.html#doubling-dbl-2008-s-1
func (p *{{ $TJacobianExtended }}) doubleMixed(a *{{ $TAffine }}) *{{ $TJacobianExtended }} {
	{{ template "mDouble" dict "all" . "negate" false}}
}

{{define "mDouble" }}
	var U, V, W, S, XX, M, S2, L {{.all.CoordType}}

	U.Double(&a.Y)
	{{- if .negate}}
		U.Neg(&U)
	{{- end}}
	V.Square(&U)
    W.Mul(&U, &V)
    S.Mul(&a.X, &V)
    XX.Square(&a.X)
    M.Double(&XX).
        Add(&M, &XX) // -> + A, but A=0 here
    S2.Double(&S)
    L.Mul(&W, &a.Y)

    p.X.Square(&M).
        Sub(&p.X, &S2)
    p.Y.Sub(&S, &p.X).
        Mul(&p.Y, &M).
	{{- if .negate}}
        Add(&p.Y, &L)
	{{- else}}
        Sub(&p.Y, &L)
	{{- end}}
    p.ZZ.Set(&V)
    p.ZZZ.Set(&W)

	return p
{{- end}}

{{define "mAdd" }}
	//if a is infinity return p
	if a.IsInfinity() {
		return p
	}
	// p is infinity, return a
	if p.ZZ.IsZero() {
		p.X = a.X
		{{- if .negate}}
			p.Y.Neg(&a.Y)
		{{- else }}
			p.Y = a.Y
		{{- end}}
		p.ZZ.SetOne()
		p.ZZZ.SetOne()
		return p
	}

	var P, R {{.all.CoordType}}

	// p2: a, p1: p
	P.Mul(&a.X, &p.ZZ)
	P.Sub(&P, &p.X)

	R.Mul(&a.Y, &p.ZZZ)
	{{- if .negate}}
		R.Neg(&R)
	{{- end}}
	R.Sub(&R, &p.Y)

	if P.IsZero() {
		if R.IsZero() {
			{{- if .negate}}
				return p.doubleNegMixed(a)
			{{ else }}
				return p.doubleMixed(a)
			{{- end}}

		}
		p.ZZ = {{.all.CoordType}}{}
		p.ZZZ = {{.all.CoordType}}{}
		return p
	}

	var PP, PPP, Q, Q2, RR, X3, Y3 {{.all.CoordType}}


	PP.Square(&P)
	PPP.Mul(&P, &PP)
	Q.Mul(&p.X, &PP)
	RR.Square(&R)
	X3.Sub(&RR, &PPP)
	Q2.Double(&Q)
	p.X.Sub(&X3, &Q2)
	Y3.Sub(&Q, &p.X).Mul(&Y3, &R)
	R.Mul(&p.Y, &PPP)
	p.Y.Sub(&Y3, &R)
	p.ZZ.Mul(&p.ZZ, &PP)
	p.ZZZ.Mul(&p.ZZZ, &PPP)

	return p
{{ end }}


{{- if .Projective }}
// -------------------------------------------------------------------------------------------------
// Homogenous projective coordinates

// Set sets p to a in projective coordinates.
func (p *{{ $TProjective }}) Set(q *{{ $TProjective }}) *{{ $TProjective }} {
	p.x, p.y, p.z = q.x, q.y, q.z
	return p
}

// Neg sets p to the projective negative point -q = (q.X, -q.Y).
func (p *{{ $TProjective }}) Neg(q *{{ $TProjective }}) *{{ $TProjective }} {
	*p = *q
	p.y.Neg(&q.y)
	return p
}

// FromAffine converts q in affine to p in projective coordinates.
func (p *{{ $TProjective }}) FromAffine(a *{{ $TAffine }}) *{{ $TProjective }} {
	if a.X.IsZero() && a.Y.IsZero() {
		p.z.SetZero()
		p.x.SetOne()
		p.y.SetOne()
		return p
	}
	p.z.SetOne()
	p.x.Set(&a.X)
	p.y.Set(&a.Y)
	return p
}

{{end }}


{{/* note batch inversion for g2 elements with E2 that is curve specific is a bit more troublesome to implement */}}
{{- if eq .PointName "g1"}}

// BatchJacobianToAffine{{ toUpper .PointName }} converts points in Jacobian coordinates to Affine coordinates
// performing a single field inversion using the Montgomery batch inversion trick.
func BatchJacobianToAffine{{ toUpper .PointName }}(points []{{ $TJacobian }}) []{{ $TAffine }} {
	result := make([]{{ $TAffine }}, len(points))
	zeroes := make([]bool, len(points))
	accumulator := fp.One()

	// batch invert all points[].Z coordinates with Montgomery batch inversion trick
	// (stores points[].Z^-1 in result[i].X to avoid allocating a slice of fr.Elements)
	for i:=0; i < len(points); i++ {
		if points[i].Z.IsZero() {
			zeroes[i] = true
			continue
		}
		result[i].X = accumulator
		accumulator.Mul(&accumulator, &points[i].Z)
	}

	var accInverse fp.Element
	accInverse.Inverse(&accumulator)

	for i := len(points) - 1; i >= 0; i-- {
		if zeroes[i] {
			// do nothing, (X=0, Y=0) is infinity point in affine
			continue
		}
		result[i].X.Mul(&result[i].X, &accInverse)
		accInverse.Mul(&accInverse, &points[i].Z)
	}

	// batch convert to affine.
	parallel.Execute( len(points), func(start, end int) {
		for i:=start; i < end; i++ {
			if zeroes[i] {
				// do nothing, (X=0, Y=0) is infinity point in affine
				continue
			}
			var a, b fp.Element
			a = result[i].X
			b.Square(&a)
			result[i].X.Mul(&points[i].X, &b)
			result[i].Y.Mul(&points[i].Y, &b).
				Mul(&result[i].Y, &a)
		}
	})

    return result
}
{{- end}}


// BatchScalarMultiplication{{ toUpper .PointName }} multiplies the same base by all scalars
// and return resulting points in affine coordinates
// uses a simple windowed-NAF-like multiplication algorithm.
func BatchScalarMultiplication{{ toUpper .PointName }}(base *{{ $TAffine }}, scalars []fr.Element) []{{ $TAffine }} {
	// approximate cost in group ops is
	// cost = 2^{c-1} + n(scalar.nbBits+nbChunks)

	nbPoints := uint64(len(scalars))
	min := ^uint64(0)
	bestC := 0
	for c := 2; c <= 16; c++  {
		cost := uint64(1 << (c-1)) // pre compute the table
		nbChunks := computeNbChunks(uint64(c))
		cost += nbPoints * (uint64(c) + 1) * nbChunks // doublings + point add
		if cost < min {
			min = cost
			bestC = c
		}
	}
	c := uint64(bestC)  // window size
	nbChunks := int(computeNbChunks(c))

	// last window may be slightly larger than c; in which case we need to compute one
	// extra element in the baseTable
	maxC := lastC(c)
	if c > maxC {
		maxC = c
	}

	// precompute all powers of base for our window
	// note here that if performance is critical, we can implement as in the msmX methods
	// this allocation to be on the stack
	baseTable := make([]{{ $TJacobian }}, (1<<(maxC-1)))
	baseTable[0].FromAffine(base)
	for i:=1;i<len(baseTable);i++ {
		baseTable[i] = baseTable[i-1]
		baseTable[i].AddMixed(base)
	}

	{{- if eq .PointName "g1"}}
		// convert our base exp table into affine to use AddMixed
		baseTableAff := BatchJacobianToAffine{{ toUpper .PointName}}(baseTable)
		toReturn := make([]{{ $TJacobian }}, len(scalars))
	{{- else}}
		toReturn := make([]{{ $TAffine }}, len(scalars))
	{{- end}}

	// partition the scalars into digits
	digits, _ := partitionScalars(scalars, c, runtime.NumCPU())

	// for each digit, take value in the base table, double it c time, voilà.
	parallel.Execute( len(scalars), func(start, end int) {
		var p {{ $TJacobian }}
		for i:=start; i < end; i++ {
			p.Set(&{{ toLower .PointName}}Infinity)
			for chunk := nbChunks - 1; chunk >=0; chunk-- {
				if chunk != nbChunks -1 {
					for j:=uint64(0); j<c; j++ {
						p.DoubleAssign()
					}
				}
				offset := chunk*len(scalars)
				digit := digits[i+offset]

				if digit == 0 {
					continue
				}

				// if msbWindow bit is set, we need to subtract
				if digit & 1 == 0 {
					// add
					{{- if eq .PointName "g1"}}
						p.AddMixed(&baseTableAff[(digit >> 1)-1])
					{{- else}}
						p.AddAssign(&baseTable[(digit >> 1)-1])
					{{- end}}
				} else {
					// sub
					{{- if eq .PointName "g1"}}
						t := baseTableAff[digit >> 1]
						t.Neg(&t)
						p.AddMixed(&t)
					{{- else}}
						t := baseTable[digit >> 1]
						t.Neg(&t)
						p.AddAssign(&t)
					{{- end}}
				}
			}

			// set our result point
			{{- if eq .PointName "g1"}}
				toReturn[i] = p
			{{- else}}
				toReturn[i].FromJacobian(&p)
			{{- end}}

		}
	})

	{{- if eq .PointName "g1"}}
		toReturnAff := BatchJacobianToAffine{{ toUpper .PointName}}(toReturn)
		return toReturnAff
	{{- else}}
		return toReturn
	{{- end}}
}

// batchAdd{{ $TAffine }} adds affine points using the Montgomery batch inversion trick.
// Special cases (doubling, infinity) must be filtered out before this call.
func batchAdd{{ $TAffine }}[TP p{{ $TAffine }}, TPP pp{{ $TAffine }}, TC c{{ $TAffine }}](R *TPP,P *TP, batchSize int) {
	var lambda, lambdain TC


	// add part
	for j := 0; j < batchSize; j++ {
		lambdain[j].Sub(&(*P)[j].X, &(*R)[j].X)
	}

	// invert denominator using montgomery batch invert technique
	{
		var accumulator {{.CoordType}}
		lambda[0].SetOne()
		accumulator.Set(&lambdain[0])

		for i := 1; i < batchSize; i++ {
			lambda[i] = accumulator
			accumulator.Mul(&accumulator, &lambdain[i])
		}

		accumulator.Inverse(&accumulator)

		for i := batchSize - 1; i > 0; i-- {
			lambda[i].Mul(&lambda[i], &accumulator)
			accumulator.Mul(&accumulator, &lambdain[i])
		}
		lambda[0].Set(&accumulator)
	}

	var d {{.CoordType}}
	var rr {{ $TAffine }}

	// add part
	for j := 0; j < batchSize; j++ {
		// computa lambda
		d.Sub(&(*P)[j].Y, &(*R)[j].Y)
		lambda[j].Mul(&lambda[j], &d)

		// compute X, Y
		rr.X.Square(&lambda[j])
		rr.X.Sub(&rr.X, &(*R)[j].X)
		rr.X.Sub(&rr.X, &(*P)[j].X)
		d.Sub(&(*R)[j].X, &rr.X)
		rr.Y.Mul(&lambda[j], &d)
		rr.Y.Sub(&rr.Y, &(*R)[j].Y)
		(*R)[j].Set(&rr)
	}
}

{{ if eq .PointName "g2"}}
// RandomOnG2 produces a random point in G2
// using standard map-to-curve methods, which means the relative discrete log
// of the generated point with respect to the canonical generator is not known.
func RandomOnG2() (G2Affine, error) {
	if gBytes, err := randomFrSizedBytes(); err != nil {
		return G2Affine{}, err
	} else {
		return HashToG2(gBytes, []byte("random on g2"))
	}
}

func randomFrSizedBytes() ([]byte, error) {
	res := make([]byte, fr.Bytes)
	_, err := rand.Read(res)
	return res, err
}
{{- end}}