github.com/consensys/gnark-crypto@v0.14.0/ecc/bls12-377/multiexp.go (about)

     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  
    17  package bls12377
    18  
    19  import (
    20  	"errors"
    21  	"github.com/consensys/gnark-crypto/ecc"
    22  	"github.com/consensys/gnark-crypto/ecc/bls12-377/fr"
    23  	"github.com/consensys/gnark-crypto/internal/parallel"
    24  	"math"
    25  	"runtime"
    26  )
    27  
    28  // MultiExp implements section 4 of https://eprint.iacr.org/2012/549.pdf
    29  //
    30  // This call return an error if len(scalars) != len(points) or if provided config is invalid.
    31  func (p *G1Affine) MultiExp(points []G1Affine, scalars []fr.Element, config ecc.MultiExpConfig) (*G1Affine, error) {
    32  	var _p G1Jac
    33  	if _, err := _p.MultiExp(points, scalars, config); err != nil {
    34  		return nil, err
    35  	}
    36  	p.FromJacobian(&_p)
    37  	return p, nil
    38  }
    39  
    40  // MultiExp implements section 4 of https://eprint.iacr.org/2012/549.pdf
    41  //
    42  // This call return an error if len(scalars) != len(points) or if provided config is invalid.
    43  func (p *G1Jac) MultiExp(points []G1Affine, scalars []fr.Element, config ecc.MultiExpConfig) (*G1Jac, error) {
    44  	// TODO @gbotrel replace the ecc.MultiExpConfig by a Option pattern for maintainability.
    45  	// note:
    46  	// each of the msmCX method is the same, except for the c constant it declares
    47  	// duplicating (through template generation) these methods allows to declare the buckets on the stack
    48  	// the choice of c needs to be improved:
    49  	// there is a theoretical value that gives optimal asymptotics
    50  	// but in practice, other factors come into play, including:
    51  	// * if c doesn't divide 64, the word size, then we're bound to select bits over 2 words of our scalars, instead of 1
    52  	// * number of CPUs
    53  	// * cache friendliness (which depends on the host, G1 or G2... )
    54  	//	--> for example, on BN254, a G1 point fits into one cache line of 64bytes, but a G2 point don't.
    55  
    56  	// for each msmCX
    57  	// step 1
    58  	// we compute, for each scalars over c-bit wide windows, nbChunk digits
    59  	// if the digit is larger than 2^{c-1}, then, we borrow 2^c from the next window and subtract
    60  	// 2^{c} to the current digit, making it negative.
    61  	// negative digits will be processed in the next step as adding -G into the bucket instead of G
    62  	// (computing -G is cheap, and this saves us half of the buckets)
    63  	// step 2
    64  	// buckets are declared on the stack
    65  	// notice that we have 2^{c-1} buckets instead of 2^{c} (see step1)
    66  	// we use jacobian extended formulas here as they are faster than mixed addition
    67  	// msmProcessChunk places points into buckets base on their selector and return the weighted bucket sum in given channel
    68  	// step 3
    69  	// reduce the buckets weighed sums into our result (msmReduceChunk)
    70  
    71  	// ensure len(points) == len(scalars)
    72  	nbPoints := len(points)
    73  	if nbPoints != len(scalars) {
    74  		return nil, errors.New("len(points) != len(scalars)")
    75  	}
    76  
    77  	// if nbTasks is not set, use all available CPUs
    78  	if config.NbTasks <= 0 {
    79  		config.NbTasks = runtime.NumCPU() * 2
    80  	} else if config.NbTasks > 1024 {
    81  		return nil, errors.New("invalid config: config.NbTasks > 1024")
    82  	}
    83  
    84  	// here, we compute the best C for nbPoints
    85  	// we split recursively until nbChunks(c) >= nbTasks,
    86  	bestC := func(nbPoints int) uint64 {
    87  		// implemented msmC methods (the c we use must be in this slice)
    88  		implementedCs := []uint64{4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}
    89  		var C uint64
    90  		// approximate cost (in group operations)
    91  		// cost = bits/c * (nbPoints + 2^{c})
    92  		// this needs to be verified empirically.
    93  		// for example, on a MBP 2016, for G2 MultiExp > 8M points, hand picking c gives better results
    94  		min := math.MaxFloat64
    95  		for _, c := range implementedCs {
    96  			cc := (fr.Bits + 1) * (nbPoints + (1 << c))
    97  			cost := float64(cc) / float64(c)
    98  			if cost < min {
    99  				min = cost
   100  				C = c
   101  			}
   102  		}
   103  		return C
   104  	}
   105  
   106  	C := bestC(nbPoints)
   107  	nbChunks := int(computeNbChunks(C))
   108  
   109  	// should we recursively split the msm in half? (see below)
   110  	// we want to minimize the execution time of the algorithm;
   111  	// splitting the msm will **add** operations, but if it allows to use more CPU, it might be worth it.
   112  
   113  	// costFunction returns a metric that represent the "wall time" of the algorithm
   114  	costFunction := func(nbTasks, nbCpus, costPerTask int) int {
   115  		// cost for the reduction of all tasks (msmReduceChunk)
   116  		totalCost := nbTasks
   117  
   118  		// cost for the computation of each task (msmProcessChunk)
   119  		for nbTasks >= nbCpus {
   120  			nbTasks -= nbCpus
   121  			totalCost += costPerTask
   122  		}
   123  		if nbTasks > 0 {
   124  			totalCost += costPerTask
   125  		}
   126  		return totalCost
   127  	}
   128  
   129  	// costPerTask is the approximate number of group ops per task
   130  	costPerTask := func(c uint64, nbPoints int) int { return (nbPoints + int((1 << c))) }
   131  
   132  	costPreSplit := costFunction(nbChunks, config.NbTasks, costPerTask(C, nbPoints))
   133  
   134  	cPostSplit := bestC(nbPoints / 2)
   135  	nbChunksPostSplit := int(computeNbChunks(cPostSplit))
   136  	costPostSplit := costFunction(nbChunksPostSplit*2, config.NbTasks, costPerTask(cPostSplit, nbPoints/2))
   137  
   138  	// if the cost of the split msm is lower than the cost of the non split msm, we split
   139  	if costPostSplit < costPreSplit {
   140  		config.NbTasks = int(math.Ceil(float64(config.NbTasks) / 2.0))
   141  		var _p G1Jac
   142  		chDone := make(chan struct{}, 1)
   143  		go func() {
   144  			_p.MultiExp(points[:nbPoints/2], scalars[:nbPoints/2], config)
   145  			close(chDone)
   146  		}()
   147  		p.MultiExp(points[nbPoints/2:], scalars[nbPoints/2:], config)
   148  		<-chDone
   149  		p.AddAssign(&_p)
   150  		return p, nil
   151  	}
   152  
   153  	// if we don't split, we use the best C we found
   154  	_innerMsmG1(p, C, points, scalars, config)
   155  
   156  	return p, nil
   157  }
   158  
   159  func _innerMsmG1(p *G1Jac, c uint64, points []G1Affine, scalars []fr.Element, config ecc.MultiExpConfig) *G1Jac {
   160  	// partition the scalars
   161  	digits, chunkStats := partitionScalars(scalars, c, config.NbTasks)
   162  
   163  	nbChunks := computeNbChunks(c)
   164  
   165  	// for each chunk, spawn one go routine that'll loop through all the scalars in the
   166  	// corresponding bit-window
   167  	// note that buckets is an array allocated on the stack and this is critical for performance
   168  
   169  	// each go routine sends its result in chChunks[i] channel
   170  	chChunks := make([]chan g1JacExtended, nbChunks)
   171  	for i := 0; i < len(chChunks); i++ {
   172  		chChunks[i] = make(chan g1JacExtended, 1)
   173  	}
   174  
   175  	// we use a semaphore to limit the number of go routines running concurrently
   176  	// (only if nbTasks < nbCPU)
   177  	var sem chan struct{}
   178  	if config.NbTasks < runtime.NumCPU() {
   179  		// we add nbChunks because if chunk is overweight we split it in two
   180  		sem = make(chan struct{}, config.NbTasks+int(nbChunks))
   181  		for i := 0; i < config.NbTasks; i++ {
   182  			sem <- struct{}{}
   183  		}
   184  		defer func() {
   185  			close(sem)
   186  		}()
   187  	}
   188  
   189  	// the last chunk may be processed with a different method than the rest, as it could be smaller.
   190  	n := len(points)
   191  	for j := int(nbChunks - 1); j >= 0; j-- {
   192  		processChunk := getChunkProcessorG1(c, chunkStats[j])
   193  		if j == int(nbChunks-1) {
   194  			processChunk = getChunkProcessorG1(lastC(c), chunkStats[j])
   195  		}
   196  		if chunkStats[j].weight >= 115 {
   197  			// we split this in more go routines since this chunk has more work to do than the others.
   198  			// else what would happen is this go routine would finish much later than the others.
   199  			chSplit := make(chan g1JacExtended, 2)
   200  			split := n / 2
   201  
   202  			if sem != nil {
   203  				sem <- struct{}{} // add another token to the semaphore, since we split in two.
   204  			}
   205  			go processChunk(uint64(j), chSplit, c, points[:split], digits[j*n:(j*n)+split], sem)
   206  			go processChunk(uint64(j), chSplit, c, points[split:], digits[(j*n)+split:(j+1)*n], sem)
   207  			go func(chunkID int) {
   208  				s1 := <-chSplit
   209  				s2 := <-chSplit
   210  				close(chSplit)
   211  				s1.add(&s2)
   212  				chChunks[chunkID] <- s1
   213  			}(j)
   214  			continue
   215  		}
   216  		go processChunk(uint64(j), chChunks[j], c, points, digits[j*n:(j+1)*n], sem)
   217  	}
   218  
   219  	return msmReduceChunkG1Affine(p, int(c), chChunks[:])
   220  }
   221  
   222  // getChunkProcessorG1 decides, depending on c window size and statistics for the chunk
   223  // to return the best algorithm to process the chunk.
   224  func getChunkProcessorG1(c uint64, stat chunkStat) func(chunkID uint64, chRes chan<- g1JacExtended, c uint64, points []G1Affine, digits []uint16, sem chan struct{}) {
   225  	switch c {
   226  
   227  	case 2:
   228  		return processChunkG1Jacobian[bucketg1JacExtendedC2]
   229  	case 4:
   230  		return processChunkG1Jacobian[bucketg1JacExtendedC4]
   231  	case 5:
   232  		return processChunkG1Jacobian[bucketg1JacExtendedC5]
   233  	case 6:
   234  		return processChunkG1Jacobian[bucketg1JacExtendedC6]
   235  	case 7:
   236  		return processChunkG1Jacobian[bucketg1JacExtendedC7]
   237  	case 8:
   238  		return processChunkG1Jacobian[bucketg1JacExtendedC8]
   239  	case 9:
   240  		return processChunkG1Jacobian[bucketg1JacExtendedC9]
   241  	case 10:
   242  		const batchSize = 80
   243  		// here we could check some chunk statistic (deviation, ...) to determine if calling
   244  		// the batch affine version is worth it.
   245  		if stat.nbBucketFilled < batchSize {
   246  			// clear indicator that batch affine method is not appropriate here.
   247  			return processChunkG1Jacobian[bucketg1JacExtendedC10]
   248  		}
   249  		return processChunkG1BatchAffine[bucketg1JacExtendedC10, bucketG1AffineC10, bitSetC10, pG1AffineC10, ppG1AffineC10, qG1AffineC10, cG1AffineC10]
   250  	case 11:
   251  		const batchSize = 150
   252  		// here we could check some chunk statistic (deviation, ...) to determine if calling
   253  		// the batch affine version is worth it.
   254  		if stat.nbBucketFilled < batchSize {
   255  			// clear indicator that batch affine method is not appropriate here.
   256  			return processChunkG1Jacobian[bucketg1JacExtendedC11]
   257  		}
   258  		return processChunkG1BatchAffine[bucketg1JacExtendedC11, bucketG1AffineC11, bitSetC11, pG1AffineC11, ppG1AffineC11, qG1AffineC11, cG1AffineC11]
   259  	case 12:
   260  		const batchSize = 200
   261  		// here we could check some chunk statistic (deviation, ...) to determine if calling
   262  		// the batch affine version is worth it.
   263  		if stat.nbBucketFilled < batchSize {
   264  			// clear indicator that batch affine method is not appropriate here.
   265  			return processChunkG1Jacobian[bucketg1JacExtendedC12]
   266  		}
   267  		return processChunkG1BatchAffine[bucketg1JacExtendedC12, bucketG1AffineC12, bitSetC12, pG1AffineC12, ppG1AffineC12, qG1AffineC12, cG1AffineC12]
   268  	case 13:
   269  		const batchSize = 350
   270  		// here we could check some chunk statistic (deviation, ...) to determine if calling
   271  		// the batch affine version is worth it.
   272  		if stat.nbBucketFilled < batchSize {
   273  			// clear indicator that batch affine method is not appropriate here.
   274  			return processChunkG1Jacobian[bucketg1JacExtendedC13]
   275  		}
   276  		return processChunkG1BatchAffine[bucketg1JacExtendedC13, bucketG1AffineC13, bitSetC13, pG1AffineC13, ppG1AffineC13, qG1AffineC13, cG1AffineC13]
   277  	case 14:
   278  		const batchSize = 400
   279  		// here we could check some chunk statistic (deviation, ...) to determine if calling
   280  		// the batch affine version is worth it.
   281  		if stat.nbBucketFilled < batchSize {
   282  			// clear indicator that batch affine method is not appropriate here.
   283  			return processChunkG1Jacobian[bucketg1JacExtendedC14]
   284  		}
   285  		return processChunkG1BatchAffine[bucketg1JacExtendedC14, bucketG1AffineC14, bitSetC14, pG1AffineC14, ppG1AffineC14, qG1AffineC14, cG1AffineC14]
   286  	case 15:
   287  		const batchSize = 500
   288  		// here we could check some chunk statistic (deviation, ...) to determine if calling
   289  		// the batch affine version is worth it.
   290  		if stat.nbBucketFilled < batchSize {
   291  			// clear indicator that batch affine method is not appropriate here.
   292  			return processChunkG1Jacobian[bucketg1JacExtendedC15]
   293  		}
   294  		return processChunkG1BatchAffine[bucketg1JacExtendedC15, bucketG1AffineC15, bitSetC15, pG1AffineC15, ppG1AffineC15, qG1AffineC15, cG1AffineC15]
   295  	case 16:
   296  		const batchSize = 640
   297  		// here we could check some chunk statistic (deviation, ...) to determine if calling
   298  		// the batch affine version is worth it.
   299  		if stat.nbBucketFilled < batchSize {
   300  			// clear indicator that batch affine method is not appropriate here.
   301  			return processChunkG1Jacobian[bucketg1JacExtendedC16]
   302  		}
   303  		return processChunkG1BatchAffine[bucketg1JacExtendedC16, bucketG1AffineC16, bitSetC16, pG1AffineC16, ppG1AffineC16, qG1AffineC16, cG1AffineC16]
   304  	default:
   305  		// panic("will not happen c != previous values is not generated by templates")
   306  		return processChunkG1Jacobian[bucketg1JacExtendedC16]
   307  	}
   308  }
   309  
   310  // msmReduceChunkG1Affine reduces the weighted sum of the buckets into the result of the multiExp
   311  func msmReduceChunkG1Affine(p *G1Jac, c int, chChunks []chan g1JacExtended) *G1Jac {
   312  	var _p g1JacExtended
   313  	totalj := <-chChunks[len(chChunks)-1]
   314  	_p.Set(&totalj)
   315  	for j := len(chChunks) - 2; j >= 0; j-- {
   316  		for l := 0; l < c; l++ {
   317  			_p.double(&_p)
   318  		}
   319  		totalj := <-chChunks[j]
   320  		_p.add(&totalj)
   321  	}
   322  
   323  	return p.unsafeFromJacExtended(&_p)
   324  }
   325  
   326  // Fold computes the multi-exponentiation \sum_{i=0}^{len(points)-1} points[i] *
   327  // combinationCoeff^i and stores the result in p. It returns error in case
   328  // configuration is invalid.
   329  func (p *G1Affine) Fold(points []G1Affine, combinationCoeff fr.Element, config ecc.MultiExpConfig) (*G1Affine, error) {
   330  	var _p G1Jac
   331  	if _, err := _p.Fold(points, combinationCoeff, config); err != nil {
   332  		return nil, err
   333  	}
   334  	p.FromJacobian(&_p)
   335  	return p, nil
   336  }
   337  
   338  // Fold computes the multi-exponentiation \sum_{i=0}^{len(points)-1} points[i] *
   339  // combinationCoeff^i and stores the result in p. It returns error in case
   340  // configuration is invalid.
   341  func (p *G1Jac) Fold(points []G1Affine, combinationCoeff fr.Element, config ecc.MultiExpConfig) (*G1Jac, error) {
   342  	scalars := make([]fr.Element, len(points))
   343  	scalar := fr.NewElement(1)
   344  	for i := 0; i < len(points); i++ {
   345  		scalars[i].Set(&scalar)
   346  		scalar.Mul(&scalar, &combinationCoeff)
   347  	}
   348  	return p.MultiExp(points, scalars, config)
   349  }
   350  
   351  // MultiExp implements section 4 of https://eprint.iacr.org/2012/549.pdf
   352  //
   353  // This call return an error if len(scalars) != len(points) or if provided config is invalid.
   354  func (p *G2Affine) MultiExp(points []G2Affine, scalars []fr.Element, config ecc.MultiExpConfig) (*G2Affine, error) {
   355  	var _p G2Jac
   356  	if _, err := _p.MultiExp(points, scalars, config); err != nil {
   357  		return nil, err
   358  	}
   359  	p.FromJacobian(&_p)
   360  	return p, nil
   361  }
   362  
   363  // MultiExp implements section 4 of https://eprint.iacr.org/2012/549.pdf
   364  //
   365  // This call return an error if len(scalars) != len(points) or if provided config is invalid.
   366  func (p *G2Jac) MultiExp(points []G2Affine, scalars []fr.Element, config ecc.MultiExpConfig) (*G2Jac, error) {
   367  	// TODO @gbotrel replace the ecc.MultiExpConfig by a Option pattern for maintainability.
   368  	// note:
   369  	// each of the msmCX method is the same, except for the c constant it declares
   370  	// duplicating (through template generation) these methods allows to declare the buckets on the stack
   371  	// the choice of c needs to be improved:
   372  	// there is a theoretical value that gives optimal asymptotics
   373  	// but in practice, other factors come into play, including:
   374  	// * if c doesn't divide 64, the word size, then we're bound to select bits over 2 words of our scalars, instead of 1
   375  	// * number of CPUs
   376  	// * cache friendliness (which depends on the host, G1 or G2... )
   377  	//	--> for example, on BN254, a G1 point fits into one cache line of 64bytes, but a G2 point don't.
   378  
   379  	// for each msmCX
   380  	// step 1
   381  	// we compute, for each scalars over c-bit wide windows, nbChunk digits
   382  	// if the digit is larger than 2^{c-1}, then, we borrow 2^c from the next window and subtract
   383  	// 2^{c} to the current digit, making it negative.
   384  	// negative digits will be processed in the next step as adding -G into the bucket instead of G
   385  	// (computing -G is cheap, and this saves us half of the buckets)
   386  	// step 2
   387  	// buckets are declared on the stack
   388  	// notice that we have 2^{c-1} buckets instead of 2^{c} (see step1)
   389  	// we use jacobian extended formulas here as they are faster than mixed addition
   390  	// msmProcessChunk places points into buckets base on their selector and return the weighted bucket sum in given channel
   391  	// step 3
   392  	// reduce the buckets weighed sums into our result (msmReduceChunk)
   393  
   394  	// ensure len(points) == len(scalars)
   395  	nbPoints := len(points)
   396  	if nbPoints != len(scalars) {
   397  		return nil, errors.New("len(points) != len(scalars)")
   398  	}
   399  
   400  	// if nbTasks is not set, use all available CPUs
   401  	if config.NbTasks <= 0 {
   402  		config.NbTasks = runtime.NumCPU() * 2
   403  	} else if config.NbTasks > 1024 {
   404  		return nil, errors.New("invalid config: config.NbTasks > 1024")
   405  	}
   406  
   407  	// here, we compute the best C for nbPoints
   408  	// we split recursively until nbChunks(c) >= nbTasks,
   409  	bestC := func(nbPoints int) uint64 {
   410  		// implemented msmC methods (the c we use must be in this slice)
   411  		implementedCs := []uint64{4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}
   412  		var C uint64
   413  		// approximate cost (in group operations)
   414  		// cost = bits/c * (nbPoints + 2^{c})
   415  		// this needs to be verified empirically.
   416  		// for example, on a MBP 2016, for G2 MultiExp > 8M points, hand picking c gives better results
   417  		min := math.MaxFloat64
   418  		for _, c := range implementedCs {
   419  			cc := (fr.Bits + 1) * (nbPoints + (1 << c))
   420  			cost := float64(cc) / float64(c)
   421  			if cost < min {
   422  				min = cost
   423  				C = c
   424  			}
   425  		}
   426  		return C
   427  	}
   428  
   429  	C := bestC(nbPoints)
   430  	nbChunks := int(computeNbChunks(C))
   431  
   432  	// should we recursively split the msm in half? (see below)
   433  	// we want to minimize the execution time of the algorithm;
   434  	// splitting the msm will **add** operations, but if it allows to use more CPU, it might be worth it.
   435  
   436  	// costFunction returns a metric that represent the "wall time" of the algorithm
   437  	costFunction := func(nbTasks, nbCpus, costPerTask int) int {
   438  		// cost for the reduction of all tasks (msmReduceChunk)
   439  		totalCost := nbTasks
   440  
   441  		// cost for the computation of each task (msmProcessChunk)
   442  		for nbTasks >= nbCpus {
   443  			nbTasks -= nbCpus
   444  			totalCost += costPerTask
   445  		}
   446  		if nbTasks > 0 {
   447  			totalCost += costPerTask
   448  		}
   449  		return totalCost
   450  	}
   451  
   452  	// costPerTask is the approximate number of group ops per task
   453  	costPerTask := func(c uint64, nbPoints int) int { return (nbPoints + int((1 << c))) }
   454  
   455  	costPreSplit := costFunction(nbChunks, config.NbTasks, costPerTask(C, nbPoints))
   456  
   457  	cPostSplit := bestC(nbPoints / 2)
   458  	nbChunksPostSplit := int(computeNbChunks(cPostSplit))
   459  	costPostSplit := costFunction(nbChunksPostSplit*2, config.NbTasks, costPerTask(cPostSplit, nbPoints/2))
   460  
   461  	// if the cost of the split msm is lower than the cost of the non split msm, we split
   462  	if costPostSplit < costPreSplit {
   463  		config.NbTasks = int(math.Ceil(float64(config.NbTasks) / 2.0))
   464  		var _p G2Jac
   465  		chDone := make(chan struct{}, 1)
   466  		go func() {
   467  			_p.MultiExp(points[:nbPoints/2], scalars[:nbPoints/2], config)
   468  			close(chDone)
   469  		}()
   470  		p.MultiExp(points[nbPoints/2:], scalars[nbPoints/2:], config)
   471  		<-chDone
   472  		p.AddAssign(&_p)
   473  		return p, nil
   474  	}
   475  
   476  	// if we don't split, we use the best C we found
   477  	_innerMsmG2(p, C, points, scalars, config)
   478  
   479  	return p, nil
   480  }
   481  
   482  func _innerMsmG2(p *G2Jac, c uint64, points []G2Affine, scalars []fr.Element, config ecc.MultiExpConfig) *G2Jac {
   483  	// partition the scalars
   484  	digits, chunkStats := partitionScalars(scalars, c, config.NbTasks)
   485  
   486  	nbChunks := computeNbChunks(c)
   487  
   488  	// for each chunk, spawn one go routine that'll loop through all the scalars in the
   489  	// corresponding bit-window
   490  	// note that buckets is an array allocated on the stack and this is critical for performance
   491  
   492  	// each go routine sends its result in chChunks[i] channel
   493  	chChunks := make([]chan g2JacExtended, nbChunks)
   494  	for i := 0; i < len(chChunks); i++ {
   495  		chChunks[i] = make(chan g2JacExtended, 1)
   496  	}
   497  
   498  	// we use a semaphore to limit the number of go routines running concurrently
   499  	// (only if nbTasks < nbCPU)
   500  	var sem chan struct{}
   501  	if config.NbTasks < runtime.NumCPU() {
   502  		// we add nbChunks because if chunk is overweight we split it in two
   503  		sem = make(chan struct{}, config.NbTasks+int(nbChunks))
   504  		for i := 0; i < config.NbTasks; i++ {
   505  			sem <- struct{}{}
   506  		}
   507  		defer func() {
   508  			close(sem)
   509  		}()
   510  	}
   511  
   512  	// the last chunk may be processed with a different method than the rest, as it could be smaller.
   513  	n := len(points)
   514  	for j := int(nbChunks - 1); j >= 0; j-- {
   515  		processChunk := getChunkProcessorG2(c, chunkStats[j])
   516  		if j == int(nbChunks-1) {
   517  			processChunk = getChunkProcessorG2(lastC(c), chunkStats[j])
   518  		}
   519  		if chunkStats[j].weight >= 115 {
   520  			// we split this in more go routines since this chunk has more work to do than the others.
   521  			// else what would happen is this go routine would finish much later than the others.
   522  			chSplit := make(chan g2JacExtended, 2)
   523  			split := n / 2
   524  
   525  			if sem != nil {
   526  				sem <- struct{}{} // add another token to the semaphore, since we split in two.
   527  			}
   528  			go processChunk(uint64(j), chSplit, c, points[:split], digits[j*n:(j*n)+split], sem)
   529  			go processChunk(uint64(j), chSplit, c, points[split:], digits[(j*n)+split:(j+1)*n], sem)
   530  			go func(chunkID int) {
   531  				s1 := <-chSplit
   532  				s2 := <-chSplit
   533  				close(chSplit)
   534  				s1.add(&s2)
   535  				chChunks[chunkID] <- s1
   536  			}(j)
   537  			continue
   538  		}
   539  		go processChunk(uint64(j), chChunks[j], c, points, digits[j*n:(j+1)*n], sem)
   540  	}
   541  
   542  	return msmReduceChunkG2Affine(p, int(c), chChunks[:])
   543  }
   544  
   545  // getChunkProcessorG2 decides, depending on c window size and statistics for the chunk
   546  // to return the best algorithm to process the chunk.
   547  func getChunkProcessorG2(c uint64, stat chunkStat) func(chunkID uint64, chRes chan<- g2JacExtended, c uint64, points []G2Affine, digits []uint16, sem chan struct{}) {
   548  	switch c {
   549  
   550  	case 2:
   551  		return processChunkG2Jacobian[bucketg2JacExtendedC2]
   552  	case 4:
   553  		return processChunkG2Jacobian[bucketg2JacExtendedC4]
   554  	case 5:
   555  		return processChunkG2Jacobian[bucketg2JacExtendedC5]
   556  	case 6:
   557  		return processChunkG2Jacobian[bucketg2JacExtendedC6]
   558  	case 7:
   559  		return processChunkG2Jacobian[bucketg2JacExtendedC7]
   560  	case 8:
   561  		return processChunkG2Jacobian[bucketg2JacExtendedC8]
   562  	case 9:
   563  		return processChunkG2Jacobian[bucketg2JacExtendedC9]
   564  	case 10:
   565  		const batchSize = 80
   566  		// here we could check some chunk statistic (deviation, ...) to determine if calling
   567  		// the batch affine version is worth it.
   568  		if stat.nbBucketFilled < batchSize {
   569  			// clear indicator that batch affine method is not appropriate here.
   570  			return processChunkG2Jacobian[bucketg2JacExtendedC10]
   571  		}
   572  		return processChunkG2BatchAffine[bucketg2JacExtendedC10, bucketG2AffineC10, bitSetC10, pG2AffineC10, ppG2AffineC10, qG2AffineC10, cG2AffineC10]
   573  	case 11:
   574  		const batchSize = 150
   575  		// here we could check some chunk statistic (deviation, ...) to determine if calling
   576  		// the batch affine version is worth it.
   577  		if stat.nbBucketFilled < batchSize {
   578  			// clear indicator that batch affine method is not appropriate here.
   579  			return processChunkG2Jacobian[bucketg2JacExtendedC11]
   580  		}
   581  		return processChunkG2BatchAffine[bucketg2JacExtendedC11, bucketG2AffineC11, bitSetC11, pG2AffineC11, ppG2AffineC11, qG2AffineC11, cG2AffineC11]
   582  	case 12:
   583  		const batchSize = 200
   584  		// here we could check some chunk statistic (deviation, ...) to determine if calling
   585  		// the batch affine version is worth it.
   586  		if stat.nbBucketFilled < batchSize {
   587  			// clear indicator that batch affine method is not appropriate here.
   588  			return processChunkG2Jacobian[bucketg2JacExtendedC12]
   589  		}
   590  		return processChunkG2BatchAffine[bucketg2JacExtendedC12, bucketG2AffineC12, bitSetC12, pG2AffineC12, ppG2AffineC12, qG2AffineC12, cG2AffineC12]
   591  	case 13:
   592  		const batchSize = 350
   593  		// here we could check some chunk statistic (deviation, ...) to determine if calling
   594  		// the batch affine version is worth it.
   595  		if stat.nbBucketFilled < batchSize {
   596  			// clear indicator that batch affine method is not appropriate here.
   597  			return processChunkG2Jacobian[bucketg2JacExtendedC13]
   598  		}
   599  		return processChunkG2BatchAffine[bucketg2JacExtendedC13, bucketG2AffineC13, bitSetC13, pG2AffineC13, ppG2AffineC13, qG2AffineC13, cG2AffineC13]
   600  	case 14:
   601  		const batchSize = 400
   602  		// here we could check some chunk statistic (deviation, ...) to determine if calling
   603  		// the batch affine version is worth it.
   604  		if stat.nbBucketFilled < batchSize {
   605  			// clear indicator that batch affine method is not appropriate here.
   606  			return processChunkG2Jacobian[bucketg2JacExtendedC14]
   607  		}
   608  		return processChunkG2BatchAffine[bucketg2JacExtendedC14, bucketG2AffineC14, bitSetC14, pG2AffineC14, ppG2AffineC14, qG2AffineC14, cG2AffineC14]
   609  	case 15:
   610  		const batchSize = 500
   611  		// here we could check some chunk statistic (deviation, ...) to determine if calling
   612  		// the batch affine version is worth it.
   613  		if stat.nbBucketFilled < batchSize {
   614  			// clear indicator that batch affine method is not appropriate here.
   615  			return processChunkG2Jacobian[bucketg2JacExtendedC15]
   616  		}
   617  		return processChunkG2BatchAffine[bucketg2JacExtendedC15, bucketG2AffineC15, bitSetC15, pG2AffineC15, ppG2AffineC15, qG2AffineC15, cG2AffineC15]
   618  	case 16:
   619  		const batchSize = 640
   620  		// here we could check some chunk statistic (deviation, ...) to determine if calling
   621  		// the batch affine version is worth it.
   622  		if stat.nbBucketFilled < batchSize {
   623  			// clear indicator that batch affine method is not appropriate here.
   624  			return processChunkG2Jacobian[bucketg2JacExtendedC16]
   625  		}
   626  		return processChunkG2BatchAffine[bucketg2JacExtendedC16, bucketG2AffineC16, bitSetC16, pG2AffineC16, ppG2AffineC16, qG2AffineC16, cG2AffineC16]
   627  	default:
   628  		// panic("will not happen c != previous values is not generated by templates")
   629  		return processChunkG2Jacobian[bucketg2JacExtendedC16]
   630  	}
   631  }
   632  
   633  // msmReduceChunkG2Affine reduces the weighted sum of the buckets into the result of the multiExp
   634  func msmReduceChunkG2Affine(p *G2Jac, c int, chChunks []chan g2JacExtended) *G2Jac {
   635  	var _p g2JacExtended
   636  	totalj := <-chChunks[len(chChunks)-1]
   637  	_p.Set(&totalj)
   638  	for j := len(chChunks) - 2; j >= 0; j-- {
   639  		for l := 0; l < c; l++ {
   640  			_p.double(&_p)
   641  		}
   642  		totalj := <-chChunks[j]
   643  		_p.add(&totalj)
   644  	}
   645  
   646  	return p.unsafeFromJacExtended(&_p)
   647  }
   648  
   649  // Fold computes the multi-exponentiation \sum_{i=0}^{len(points)-1} points[i] *
   650  // combinationCoeff^i and stores the result in p. It returns error in case
   651  // configuration is invalid.
   652  func (p *G2Affine) Fold(points []G2Affine, combinationCoeff fr.Element, config ecc.MultiExpConfig) (*G2Affine, error) {
   653  	var _p G2Jac
   654  	if _, err := _p.Fold(points, combinationCoeff, config); err != nil {
   655  		return nil, err
   656  	}
   657  	p.FromJacobian(&_p)
   658  	return p, nil
   659  }
   660  
   661  // Fold computes the multi-exponentiation \sum_{i=0}^{len(points)-1} points[i] *
   662  // combinationCoeff^i and stores the result in p. It returns error in case
   663  // configuration is invalid.
   664  func (p *G2Jac) Fold(points []G2Affine, combinationCoeff fr.Element, config ecc.MultiExpConfig) (*G2Jac, error) {
   665  	scalars := make([]fr.Element, len(points))
   666  	scalar := fr.NewElement(1)
   667  	for i := 0; i < len(points); i++ {
   668  		scalars[i].Set(&scalar)
   669  		scalar.Mul(&scalar, &combinationCoeff)
   670  	}
   671  	return p.MultiExp(points, scalars, config)
   672  }
   673  
   674  // selector stores the index, mask and shifts needed to select bits from a scalar
   675  // it is used during the multiExp algorithm or the batch scalar multiplication
   676  type selector struct {
   677  	index uint64 // index in the multi-word scalar to select bits from
   678  	mask  uint64 // mask (c-bit wide)
   679  	shift uint64 // shift needed to get our bits on low positions
   680  
   681  	multiWordSelect bool   // set to true if we need to select bits from 2 words (case where c doesn't divide 64)
   682  	maskHigh        uint64 // same than mask, for index+1
   683  	shiftHigh       uint64 // same than shift, for index+1
   684  }
   685  
   686  // return number of chunks for a given window size c
   687  // the last chunk may be bigger to accommodate a potential carry from the NAF decomposition
   688  func computeNbChunks(c uint64) uint64 {
   689  	return (fr.Bits + c - 1) / c
   690  }
   691  
   692  // return the last window size for a scalar;
   693  // this last window should accommodate a carry (from the NAF decomposition)
   694  // it can be == c if we have 1 available bit
   695  // it can be > c if we have 0 available bit
   696  // it can be < c if we have 2+ available bits
   697  func lastC(c uint64) uint64 {
   698  	nbAvailableBits := (computeNbChunks(c) * c) - fr.Bits
   699  	return c + 1 - nbAvailableBits
   700  }
   701  
   702  type chunkStat struct {
   703  	// relative weight of work compared to other chunks. 100.0 -> nominal weight.
   704  	weight float32
   705  
   706  	// percentage of bucket filled in the window;
   707  	ppBucketFilled float32
   708  	nbBucketFilled int
   709  }
   710  
   711  // partitionScalars  compute, for each scalars over c-bit wide windows, nbChunk digits
   712  // if the digit is larger than 2^{c-1}, then, we borrow 2^c from the next window and subtract
   713  // 2^{c} to the current digit, making it negative.
   714  // negative digits can be processed in a later step as adding -G into the bucket instead of G
   715  // (computing -G is cheap, and this saves us half of the buckets in the MultiExp or BatchScalarMultiplication)
   716  func partitionScalars(scalars []fr.Element, c uint64, nbTasks int) ([]uint16, []chunkStat) {
   717  	// no benefit here to have more tasks than CPUs
   718  	if nbTasks > runtime.NumCPU() {
   719  		nbTasks = runtime.NumCPU()
   720  	}
   721  
   722  	// number of c-bit radixes in a scalar
   723  	nbChunks := computeNbChunks(c)
   724  
   725  	digits := make([]uint16, len(scalars)*int(nbChunks))
   726  
   727  	mask := uint64((1 << c) - 1) // low c bits are 1
   728  	max := int(1<<(c-1)) - 1     // max value (inclusive) we want for our digits
   729  	cDivides64 := (64 % c) == 0  // if c doesn't divide 64, we may need to select over multiple words
   730  
   731  	// compute offset and word selector / shift to select the right bits of our windows
   732  	selectors := make([]selector, nbChunks)
   733  	for chunk := uint64(0); chunk < nbChunks; chunk++ {
   734  		jc := uint64(chunk * c)
   735  		d := selector{}
   736  		d.index = jc / 64
   737  		d.shift = jc - (d.index * 64)
   738  		d.mask = mask << d.shift
   739  		d.multiWordSelect = !cDivides64 && d.shift > (64-c) && d.index < (fr.Limbs-1)
   740  		if d.multiWordSelect {
   741  			nbBitsHigh := d.shift - uint64(64-c)
   742  			d.maskHigh = (1 << nbBitsHigh) - 1
   743  			d.shiftHigh = (c - nbBitsHigh)
   744  		}
   745  		selectors[chunk] = d
   746  	}
   747  
   748  	parallel.Execute(len(scalars), func(start, end int) {
   749  		for i := start; i < end; i++ {
   750  			if scalars[i].IsZero() {
   751  				// everything is 0, no need to process this scalar
   752  				continue
   753  			}
   754  			scalar := scalars[i].Bits()
   755  
   756  			var carry int
   757  
   758  			// for each chunk in the scalar, compute the current digit, and an eventual carry
   759  			for chunk := uint64(0); chunk < nbChunks-1; chunk++ {
   760  				s := selectors[chunk]
   761  
   762  				// init with carry if any
   763  				digit := carry
   764  				carry = 0
   765  
   766  				// digit = value of the c-bit window
   767  				digit += int((scalar[s.index] & s.mask) >> s.shift)
   768  
   769  				if s.multiWordSelect {
   770  					// we are selecting bits over 2 words
   771  					digit += int(scalar[s.index+1]&s.maskHigh) << s.shiftHigh
   772  				}
   773  
   774  				// if the digit is larger than 2^{c-1}, then, we borrow 2^c from the next window and subtract
   775  				// 2^{c} to the current digit, making it negative.
   776  				if digit > max {
   777  					digit -= (1 << c)
   778  					carry = 1
   779  				}
   780  
   781  				// if digit is zero, no impact on result
   782  				if digit == 0 {
   783  					continue
   784  				}
   785  
   786  				var bits uint16
   787  				if digit > 0 {
   788  					bits = uint16(digit) << 1
   789  				} else {
   790  					bits = (uint16(-digit-1) << 1) + 1
   791  				}
   792  				digits[int(chunk)*len(scalars)+i] = bits
   793  			}
   794  
   795  			// for the last chunk, we don't want to borrow from a next window
   796  			// (but may have a larger max value)
   797  			chunk := nbChunks - 1
   798  			s := selectors[chunk]
   799  			// init with carry if any
   800  			digit := carry
   801  			// digit = value of the c-bit window
   802  			digit += int((scalar[s.index] & s.mask) >> s.shift)
   803  			if s.multiWordSelect {
   804  				// we are selecting bits over 2 words
   805  				digit += int(scalar[s.index+1]&s.maskHigh) << s.shiftHigh
   806  			}
   807  			digits[int(chunk)*len(scalars)+i] = uint16(digit) << 1
   808  		}
   809  
   810  	}, nbTasks)
   811  
   812  	// aggregate  chunk stats
   813  	chunkStats := make([]chunkStat, nbChunks)
   814  	if c <= 9 {
   815  		// no need to compute stats for small window sizes
   816  		return digits, chunkStats
   817  	}
   818  	parallel.Execute(len(chunkStats), func(start, end int) {
   819  		// for each chunk compute the statistics
   820  		for chunkID := start; chunkID < end; chunkID++ {
   821  			// indicates if a bucket is hit.
   822  			var b bitSetC16
   823  
   824  			// digits for the chunk
   825  			chunkDigits := digits[chunkID*len(scalars) : (chunkID+1)*len(scalars)]
   826  
   827  			totalOps := 0
   828  			nz := 0 // non zero buckets count
   829  			for _, digit := range chunkDigits {
   830  				if digit == 0 {
   831  					continue
   832  				}
   833  				totalOps++
   834  				bucketID := digit >> 1
   835  				if digit&1 == 0 {
   836  					bucketID -= 1
   837  				}
   838  				if !b[bucketID] {
   839  					nz++
   840  					b[bucketID] = true
   841  				}
   842  			}
   843  			chunkStats[chunkID].weight = float32(totalOps) // count number of ops for now, we will compute the weight after
   844  			chunkStats[chunkID].ppBucketFilled = (float32(nz) * 100.0) / float32(int(1<<(c-1)))
   845  			chunkStats[chunkID].nbBucketFilled = nz
   846  		}
   847  	}, nbTasks)
   848  
   849  	totalOps := float32(0.0)
   850  	for _, stat := range chunkStats {
   851  		totalOps += stat.weight
   852  	}
   853  
   854  	target := totalOps / float32(nbChunks)
   855  	if target != 0.0 {
   856  		// if target == 0, it means all the scalars are 0 everywhere, there is no work to be done.
   857  		for i := 0; i < len(chunkStats); i++ {
   858  			chunkStats[i].weight = (chunkStats[i].weight * 100.0) / target
   859  		}
   860  	}
   861  
   862  	return digits, chunkStats
   863  }