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

{{ $G1TAffine := print (toUpper .G1.PointName) "Affine" }}
{{ $G1TJacobian := print (toUpper .G1.PointName) "Jac" }}
{{ $G1TJacobianExtended := print (toLower .G1.PointName) "JacExtended" }}

{{ $G2TAffine := print (toUpper .G2.PointName) "Affine" }}
{{ $G2TJacobian := print (toUpper .G2.PointName) "Jac" }}
{{ $G2TJacobianExtended := print (toLower .G2.PointName) "JacExtended" }}


import (
	"fmt"
	"runtime"
    "math/rand/v2"
	"math/big"
	"testing"
    "math/bits"
	"sync"

	"github.com/consensys/gnark-crypto/ecc"
	"github.com/consensys/gnark-crypto/ecc/{{.Name}}/fr"
	"github.com/leanovate/gopter"
	"github.com/leanovate/gopter/prop"
)


{{- if ne .Name "secp256k1"}}
{{template "multiexp" dict "PointName" .G1.PointName "UPointName" (toUpper .G1.PointName) "TAffine" $G1TAffine "TJacobian" $G1TJacobian "TJacobianExtended" $G1TJacobianExtended "FrNbWords" .Fr.NbWords "CRange" .G1.CRange "cmax" 16}}
{{template "multiexp" dict "PointName" .G2.PointName "UPointName" (toUpper .G2.PointName) "TAffine" $G2TAffine "TJacobian" $G2TJacobian "TJacobianExtended" $G2TJacobianExtended "FrNbWords" .Fr.NbWords "CRange" .G2.CRange "cmax" 16}}
{{- else}}
{{template "multiexp" dict "PointName" .G1.PointName "UPointName" (toUpper .G1.PointName) "TAffine" $G1TAffine "TJacobian" $G1TJacobian "TJacobianExtended" $G1TJacobianExtended "FrNbWords" .Fr.NbWords "CRange" .G1.CRange "cmax" 15}}
{{- end}}

{{define "multiexp" }}

func TestMultiExp{{$.UPointName}}(t *testing.T) {

	parameters := gopter.DefaultTestParameters()
	if testing.Short() {
		parameters.MinSuccessfulTests = 3
	} else {
		parameters.MinSuccessfulTests = nbFuzzShort * 2
	}


	properties := gopter.NewProperties(parameters)

	genScalar := GenFr()

	// size of the multiExps
	const nbSamples = 73

	// multi exp points
	var samplePoints [nbSamples]{{ $.TAffine }}
	var g {{ $.TJacobian }}
	g.Set(&{{ toLower $.PointName }}Gen)
	for i := 1; i <= nbSamples; i++ {
		samplePoints[i-1].FromJacobian(&g)
		g.AddAssign(&{{ toLower $.PointName }}Gen)
	}

    // sprinkle some points at infinity
    samplePoints[rand.N(nbSamples)].setInfinity() //#nosec G404 weak rng is fine here
    samplePoints[rand.N(nbSamples)].setInfinity() //#nosec G404 weak rng is fine here
    samplePoints[rand.N(nbSamples)].setInfinity() //#nosec G404 weak rng is fine here
    samplePoints[rand.N(nbSamples)].setInfinity() //#nosec G404 weak rng is fine here

	// final scalar to use in double and add method (without mixer factor)
	// n(n+1)(2n+1)/6  (sum of the squares from 1 to n)
	var scalar big.Int
	scalar.SetInt64(nbSamples)
	scalar.Mul(&scalar, new(big.Int).SetInt64(nbSamples+1))
	scalar.Mul(&scalar, new(big.Int).SetInt64(2*nbSamples+1))
	scalar.Div(&scalar, new(big.Int).SetInt64(6))


	// ensure a multiexp that's splitted has the same result as a non-splitted one..
	properties.Property("[{{ $.UPointName }}] Multi exponentiation (cmax) should be consistent with splitted multiexp", prop.ForAll(
		func(mixer fr.Element) bool {
			var samplePointsLarge [nbSamples*13]{{ $.TAffine }}
			for i:=0; i<13; i++ {
				copy(samplePointsLarge[i*nbSamples:], samplePoints[:])
			}

			var rmax, splitted1, splitted2 {{ $.TJacobian }}

			// mixer ensures that all the words of a fpElement are set
			var sampleScalars [nbSamples*13]fr.Element

			for i := 1; i <= nbSamples; i++ {
				sampleScalars[i-1].SetUint64(uint64(i)).
					Mul(&sampleScalars[i-1], &mixer)
			}

			rmax.MultiExp(samplePointsLarge[:], sampleScalars[:], ecc.MultiExpConfig{})
			splitted1.MultiExp(samplePointsLarge[:], sampleScalars[:], ecc.MultiExpConfig{NbTasks: 128})
			splitted2.MultiExp(samplePointsLarge[:], sampleScalars[:], ecc.MultiExpConfig{NbTasks: 51})
			return rmax.Equal(&splitted1) && rmax.Equal(&splitted2)
		},
		genScalar,
	))

	// cRange is generated from template and contains the available parameters for the multiexp window size
	{{- if eq $.PointName "g1" }}
	cRange := []uint64{
		{{- range $c :=  $.CRange}}{{- if gt $c 1}}{{$c}},{{- end}}{{- end}}
	}
	if testing.Short() {
		// test only "odd" and "even" (ie windows size divide word size vs not)
		cRange = []uint64{5, 14}
	}
	{{- else }}
	// for g2, CI suffers with large c size since it needs to allocate a lot of memory for the buckets.
	// test only "odd" and "even" (ie windows size divide word size vs not)
	cRange := []uint64{5, 14}
	{{- end}}

	properties.Property(fmt.Sprintf("[{{ $.UPointName }}] Multi exponentiation (c in %v) should be consistent with sum of square", cRange), prop.ForAll(
		func(mixer fr.Element) bool {

			var expected {{ $.TJacobian }}

			// compute expected result with double and add
			var finalScalar,mixerBigInt big.Int
			finalScalar.Mul(&scalar, mixer.BigInt(&mixerBigInt))
			expected.ScalarMultiplication(&{{ toLower $.PointName }}Gen, &finalScalar)

			// mixer ensures that all the words of a fpElement are set
			var sampleScalars [nbSamples]fr.Element

			for i := 1; i <= nbSamples; i++ {
				sampleScalars[i-1].SetUint64(uint64(i)).
					Mul(&sampleScalars[i-1], &mixer)
			}


			results := make([]{{ $.TJacobian }}, len(cRange))
			for i, c := range cRange {
				_innerMsm{{ $.UPointName }}(&results[i], c, samplePoints[:], sampleScalars[:], ecc.MultiExpConfig{NbTasks: runtime.NumCPU()})
			}
			for i:=1; i < len(results);i++ {
				if !results[i].Equal(&results[i-1]) {
					t.Logf("result for c=%d != c=%d", cRange[i-1],cRange[i])
					return false
				}
			}
			return true
		},
		genScalar,
	))

	properties.Property(fmt.Sprintf("[{{ $.UPointName }}] Multi exponentiation (c in %v) of points at infinity should output a point at infinity", cRange), prop.ForAll(
		func(mixer fr.Element) bool {

			var samplePointsZero [nbSamples]{{ $.TAffine }}

			var expected {{ $.TJacobian }}

			// compute expected result with double and add
			var finalScalar, mixerBigInt big.Int
			finalScalar.Mul(&scalar, mixer.BigInt(&mixerBigInt))
			expected.ScalarMultiplication(&{{ toLower $.PointName }}Gen, &finalScalar)

			// mixer ensures that all the words of a fpElement are set
			var sampleScalars [nbSamples]fr.Element

			for i := 1; i <= nbSamples; i++ {
				sampleScalars[i-1].SetUint64(uint64(i)).
					Mul(&sampleScalars[i-1], &mixer)
				samplePointsZero[i-1].setInfinity()
			}

			results := make([]{{ $.TJacobian }}, len(cRange))
			for i, c := range cRange {
				_innerMsm{{ $.UPointName }}(&results[i], c, samplePointsZero[:], sampleScalars[:], ecc.MultiExpConfig{NbTasks: runtime.NumCPU()})
			}
			for i := 0; i < len(results); i++ {
				if !results[i].Z.IsZero() {
					t.Logf("result for c=%d is not infinity", cRange[i])
					return false
				}
			}
			return true
		},
		genScalar,
	))

	properties.Property(fmt.Sprintf("[{{ $.UPointName }}] Multi exponentiation (c in %v) with a vector of 0s as input should output a point at infinity", cRange), prop.ForAll(
		func(mixer fr.Element) bool {
			// mixer ensures that all the words of a fpElement are set
			var sampleScalars [nbSamples]fr.Element


			results := make([]{{ $.TJacobian }}, len(cRange))
			for i, c := range cRange {
				_innerMsm{{ $.UPointName }}(&results[i], c, samplePoints[:], sampleScalars[:], ecc.MultiExpConfig{NbTasks: runtime.NumCPU()})
			}
			for i := 0; i < len(results); i++ {
				if !results[i].Z.IsZero() {
					t.Logf("result for c=%d is not infinity", cRange[i])
					return false
				}
			}
			return true
		},
		genScalar,
	))


	// note : this test is here as we expect to have a different multiExp than the above bucket method
	// for small number of points
	properties.Property("[{{ $.UPointName }}] Multi exponentiation (<50points) should be consistent with sum of square", prop.ForAll(
		func(mixer fr.Element) bool {

			var g {{ $.TJacobian }}
			g.Set(&{{ toLower .PointName}}Gen)

			// mixer ensures that all the words of a fpElement are set
			samplePoints := make([]{{ $.TAffine }}, 30)
			sampleScalars := make([]fr.Element, 30)

			for i := 1; i <= 30; i++ {
				sampleScalars[i-1].SetUint64(uint64(i)).
					Mul(&sampleScalars[i-1], &mixer)
				samplePoints[i-1].FromJacobian(&g)
				g.AddAssign(&{{ toLower .PointName}}Gen)
			}

			var op1MultiExp {{ $.TAffine }}
			op1MultiExp.MultiExp(samplePoints, sampleScalars, ecc.MultiExpConfig{})

			var finalBigScalar fr.Element
			var finalBigScalarBi big.Int
			var op1ScalarMul {{ $.TAffine }}
			finalBigScalar.SetUint64(9455).Mul(&finalBigScalar, &mixer)
			finalBigScalar.BigInt(&finalBigScalarBi)
			op1ScalarMul.ScalarMultiplication(&{{ toLower .PointName}}GenAff, &finalBigScalarBi)

			return op1ScalarMul.Equal(&op1MultiExp)
		},
		genScalar,
	))


	properties.TestingRun(t, gopter.ConsoleReporter(false))
}


func TestCrossMultiExp{{ $.UPointName }}(t *testing.T) {
	const nbSamples = 1 << 14
	// multi exp points
	var samplePoints [nbSamples]{{ $.TAffine }}
	var g {{ $.TJacobian }}
	g.Set(&{{ toLower $.PointName }}Gen)
	for i := 1; i <= nbSamples; i++ {
		samplePoints[i-1].FromJacobian(&g)
		g.AddAssign(&{{ toLower $.PointName }}Gen)
	}

	// sprinkle some points at infinity
	samplePoints[rand.N(nbSamples)].setInfinity() //#nosec G404 weak rng is fine here
	samplePoints[rand.N(nbSamples)].setInfinity() //#nosec G404 weak rng is fine here
	samplePoints[rand.N(nbSamples)].setInfinity() //#nosec G404 weak rng is fine here
	samplePoints[rand.N(nbSamples)].setInfinity() //#nosec G404 weak rng is fine here


	var sampleScalars [nbSamples]fr.Element
	fillBenchScalars(sampleScalars[:])

	// sprinkle some doublings
	for i:=10; i < 100; i++ {
		samplePoints[i]  = samplePoints[0]
		sampleScalars[i] = sampleScalars[0]
	}

	// cRange is generated from template and contains the available parameters for the multiexp window size
	{{- if eq $.PointName "g1" }}
	cRange := []uint64{
		{{- range $c :=  $.CRange}}{{- if gt $c 1}}{{$c}},{{- end}}{{- end}}
	}
	if testing.Short() {
		// test only "odd" and "even" (ie windows size divide word size vs not)
		cRange = []uint64{5, 14}
	}
	{{- else }}
	// for g2, CI suffers with large c size since it needs to allocate a lot of memory for the buckets.
	// test only "odd" and "even" (ie windows size divide word size vs not)
	cRange := []uint64{5, 14}
	{{- end}}

	results := make([]{{ $.TJacobian }}, len(cRange))
	for i, c := range cRange {
		_innerMsm{{ $.UPointName }}(&results[i], c, samplePoints[:], sampleScalars[:], ecc.MultiExpConfig{NbTasks: runtime.NumCPU()})
	}

	var r {{ $.TJacobian }}
	_innerMsm{{ $.UPointName }}Reference(&r, samplePoints[:], sampleScalars[:], ecc.MultiExpConfig{NbTasks: runtime.NumCPU()})

	var expected, got {{ $.TAffine}}
	expected.FromJacobian(&r)

	for i:=0; i<len(results);i++ {
		got.FromJacobian(&results[i])
		if !expected.Equal(&got) {
			t.Fatalf("cross msm failed with c=%d", cRange[i])
		}
	}

}



// _innerMsm{{ $.UPointName }}Reference always do ext jacobian with c == {{$.cmax}}
func _innerMsm{{ $.UPointName }}Reference(p *{{ $.TJacobian }}, points []{{ $.TAffine }}, scalars []fr.Element, config ecc.MultiExpConfig) *{{ $.TJacobian }} {
	// partition the scalars
	digits, _ := partitionScalars(scalars, {{$.cmax}},  config.NbTasks)

	nbChunks := computeNbChunks({{$.cmax}})

	// for each chunk, spawn one go routine that'll loop through all the scalars in the
	// corresponding bit-window
	// note that buckets is an array allocated on the stack and this is critical for performance

	// each go routine sends its result in chChunks[i] channel
	chChunks := make([]chan {{ $.TJacobianExtended }}, nbChunks)
	for i:=0; i < len(chChunks);i++ {
		chChunks[i] = make(chan {{ $.TJacobianExtended }}, 1)
	}

	// the last chunk may be processed with a different method than the rest, as it could be smaller.
	n := len(points)
	for j := int(nbChunks - 1); j >= 0; j-- {
		processChunk := processChunk{{ $.UPointName }}Jacobian[bucket{{ $.TJacobianExtended }}C{{$.cmax}}]
		go processChunk(uint64(j), chChunks[j], {{$.cmax}}, points, digits[j*n:(j+1)*n], nil)
	}

	return msmReduceChunk{{ $.TAffine }}(p, int({{$.cmax}}), chChunks[:])
}


func BenchmarkMultiExp{{ $.UPointName }}(b *testing.B) {

	const (
		pow = (bits.UintSize / 2) - (bits.UintSize / 8) // 24 on 64 bits arch, 12 on 32 bits
		nbSamples = 1 << pow
	)

	var (
		samplePoints [nbSamples]{{ $.TAffine }}
		sampleScalars [nbSamples]fr.Element
		sampleScalarsSmallValues [nbSamples]fr.Element
		sampleScalarsRedundant [nbSamples]fr.Element
	)

	fillBenchScalars(sampleScalars[:])
	copy(sampleScalarsSmallValues[:],sampleScalars[:])
	copy(sampleScalarsRedundant[:],sampleScalars[:])

	// this means first chunk is going to have more work to do and should be split into several go routines
	for i:=0; i < len(sampleScalarsSmallValues);i++ {
		if i % 5 == 0 {
			sampleScalarsSmallValues[i].SetZero()
			sampleScalarsSmallValues[i][0] = 1
		}
	}

	// bad case for batch affine because scalar distribution might look uniform
	// but over batchSize windows, we may hit a lot of conflicts and force the msm-affine
	// to process small batches of additions to flush its queue of conflicted points.
	for i:=0; i < len(sampleScalarsRedundant);i+=100 {
		for j:=i+1; j < i+100 && j < len(sampleScalarsRedundant);j++ {
			sampleScalarsRedundant[j] = sampleScalarsRedundant[i]
		}
	}

	fillBenchBases{{ $.UPointName }}(samplePoints[:])


	var testPoint {{ $.TAffine }}

	for i := 5; i <= pow; i++ {
		using := 1 << i

		b.Run(fmt.Sprintf("%d points", using), func(b *testing.B) {
			b.ResetTimer()
			for j := 0; j < b.N; j++ {
				testPoint.MultiExp(samplePoints[:using], sampleScalars[:using],ecc.MultiExpConfig{})
			}
		})

		b.Run(fmt.Sprintf("%d points-smallvalues", using), func(b *testing.B) {
			b.ResetTimer()
			for j := 0; j < b.N; j++ {
				testPoint.MultiExp(samplePoints[:using], sampleScalarsSmallValues[:using],ecc.MultiExpConfig{})
			}
		})

		b.Run(fmt.Sprintf("%d points-redundancy", using), func(b *testing.B) {
			b.ResetTimer()
			for j := 0; j < b.N; j++ {
				testPoint.MultiExp(samplePoints[:using], sampleScalarsRedundant[:using],ecc.MultiExpConfig{})
			}
		})
	}
}


func BenchmarkMultiExp{{ $.UPointName }}Reference(b *testing.B) {
	const nbSamples = 1 << 20

	var (
		samplePoints [nbSamples]{{ $.TAffine }}
		sampleScalars [nbSamples]fr.Element
	)

	fillBenchScalars(sampleScalars[:])
	fillBenchBases{{ $.UPointName }}(samplePoints[:])

	var testPoint {{ $.TAffine }}

	b.ResetTimer()
	for j := 0; j < b.N; j++ {
		testPoint.MultiExp(samplePoints[:], sampleScalars[:],ecc.MultiExpConfig{})
	}
}


func BenchmarkManyMultiExp{{ $.UPointName }}Reference(b *testing.B) {
	const nbSamples = 1 << 20

	var (
		samplePoints [nbSamples]{{ $.TAffine }}
		sampleScalars [nbSamples]fr.Element
	)

	fillBenchScalars(sampleScalars[:])
	fillBenchBases{{ $.UPointName }}(samplePoints[:])


	var t1, t2, t3 {{ $.TAffine }}
	b.ResetTimer()
	for j := 0; j < b.N; j++ {
		var wg sync.WaitGroup
		wg.Add(3)
		go func() {
			t1.MultiExp(samplePoints[:], sampleScalars[:],ecc.MultiExpConfig{})
			wg.Done()
		}()
		go func() {
			t2.MultiExp(samplePoints[:], sampleScalars[:],ecc.MultiExpConfig{})
			wg.Done()
		}()
		go func() {
			t3.MultiExp(samplePoints[:], sampleScalars[:],ecc.MultiExpConfig{})
			wg.Done()
		}()
		wg.Wait()
	}
}

// WARNING: this return points that are NOT on the curve and is meant to be use for benchmarking
// purposes only. We don't check that the result is valid but just measure "computational complexity".
//
// Rationale for generating points that are not on the curve is that for large benchmarks, generating
// a vector of different points can take minutes. Using the same point or subset will bias the benchmark result
// since bucket additions in extended jacobian coordinates will hit doubling algorithm instead of add.
func fillBenchBases{{ $.UPointName }}(samplePoints []{{ $.TAffine }}) {
	var r big.Int
	r.SetString("340444420969191673093399857471996460938405", 10)
	samplePoints[0].ScalarMultiplication(&samplePoints[0], &r)

	one := samplePoints[0].X
	one.SetOne()

	for i := 1; i < len(samplePoints); i++ {
		samplePoints[i].X.Add(&samplePoints[i-1].X, &one)
		samplePoints[i].Y.Sub(&samplePoints[i-1].Y, &one)
	}
}


{{end }}


func fillBenchScalars(sampleScalars []fr.Element) {
	// ensure every words of the scalars are filled
	for i := 0; i < len(sampleScalars); i++ {
		sampleScalars[i].SetRandom()
	}
}