github.com/platonnetwork/platon-go@v0.7.6/cases/tool/win/bls_win/include/mcl/lagrange.hpp (about)

     1  #pragma once
     2  /**
     3  	@file
     4  	@brief Lagrange Interpolation
     5  	@author MITSUNARI Shigeo(@herumi)
     6  	@license modified new BSD license
     7  	http://opensource.org/licenses/BSD-3-Clause
     8  */
     9  namespace mcl {
    10  
    11  /*
    12  	recover out = f(0) by { (x, y) | x = S[i], y = f(x) = vec[i] }
    13  	@retval 0 if succeed else -1
    14  */
    15  template<class G, class F>
    16  void LagrangeInterpolation(bool *pb, G& out, const F *S, const G *vec, size_t k)
    17  {
    18  	if (k == 0) {
    19  		*pb = false;
    20  		return;
    21  	}
    22  	if (k == 1) {
    23  		out = vec[0];
    24  		*pb = true;
    25  		return;
    26  	}
    27  	/*
    28  		delta_{i,S}(0) = prod_{j != i} S[j] / (S[j] - S[i]) = a / b
    29  		where a = prod S[j], b = S[i] * prod_{j != i} (S[j] - S[i])
    30  	*/
    31  	F a = S[0];
    32  	for (size_t i = 1; i < k; i++) {
    33  		a *= S[i];
    34  	}
    35  	if (a.isZero()) {
    36  		*pb = false;
    37  		return;
    38  	}
    39  	/*
    40  		f(0) = sum_i f(S[i]) delta_{i,S}(0)
    41  	*/
    42  	G r;
    43  	r.clear();
    44  	for (size_t i = 0; i < k; i++) {
    45  		F b = S[i];
    46  		for (size_t j = 0; j < k; j++) {
    47  			if (j != i) {
    48  				F v = S[j] - S[i];
    49  				if (v.isZero()) {
    50  					*pb = false;
    51  					return;
    52  				}
    53  				b *= v;
    54  			}
    55  		}
    56  		G t;
    57  		G::mul(t, vec[i], a / b);
    58  		r += t;
    59  	}
    60  	out = r;
    61  	*pb = true;
    62  }
    63  
    64  /*
    65  	out = f(x) = c[0] + c[1] * x + c[2] * x^2 + ... + c[cSize - 1] * x^(cSize - 1)
    66  	@retval 0 if succeed else -1 (if cSize == 0)
    67  */
    68  template<class G, class T>
    69  void evaluatePolynomial(bool *pb, G& out, const G *c, size_t cSize, const T& x)
    70  {
    71  	if (cSize == 0) {
    72  		*pb = false;
    73  		return;
    74  	}
    75  	if (cSize == 1) {
    76  		out = c[0];
    77  		*pb = true;
    78  		return;
    79  	}
    80  	G y = c[cSize - 1];
    81  	for (int i = (int)cSize - 2; i >= 0; i--) {
    82  		G::mul(y, y, x);
    83  		G::add(y, y, c[i]);
    84  	}
    85  	out = y;
    86  	*pb = true;
    87  }
    88  
    89  #ifndef CYBOZU_DONT_USE_EXCEPTION
    90  template<class G, class F>
    91  void LagrangeInterpolation(G& out, const F *S, const G *vec, size_t k)
    92  {
    93  	bool b;
    94  	LagrangeInterpolation(&b, out, S, vec, k);
    95  	if (!b) throw cybozu::Exception("LagrangeInterpolation");
    96  }
    97  
    98  template<class G, class T>
    99  void evaluatePolynomial(G& out, const G *c, size_t cSize, const T& x)
   100  {
   101  	bool b;
   102  	evaluatePolynomial(&b, out, c, cSize, x);
   103  	if (!b) throw cybozu::Exception("evaluatePolynomial");
   104  }
   105  #endif
   106  
   107  } // mcl