diff options
Diffstat (limited to 'vendor/github.com/tangerine-network/mcl/include/mcl/lagrange.hpp')
-rw-r--r-- | vendor/github.com/tangerine-network/mcl/include/mcl/lagrange.hpp | 107 |
1 files changed, 107 insertions, 0 deletions
diff --git a/vendor/github.com/tangerine-network/mcl/include/mcl/lagrange.hpp b/vendor/github.com/tangerine-network/mcl/include/mcl/lagrange.hpp new file mode 100644 index 000000000..18e0597ec --- /dev/null +++ b/vendor/github.com/tangerine-network/mcl/include/mcl/lagrange.hpp @@ -0,0 +1,107 @@ +#pragma once +/** + @file + @brief Lagrange Interpolation + @author MITSUNARI Shigeo(@herumi) + @license modified new BSD license + http://opensource.org/licenses/BSD-3-Clause +*/ +namespace mcl { + +/* + recover out = f(0) by { (x, y) | x = S[i], y = f(x) = vec[i] } + @retval 0 if succeed else -1 +*/ +template<class G, class F> +void LagrangeInterpolation(bool *pb, G& out, const F *S, const G *vec, size_t k) +{ + if (k == 0) { + *pb = false; + return; + } + if (k == 1) { + out = vec[0]; + *pb = true; + return; + } + /* + delta_{i,S}(0) = prod_{j != i} S[j] / (S[j] - S[i]) = a / b + where a = prod S[j], b = S[i] * prod_{j != i} (S[j] - S[i]) + */ + F a = S[0]; + for (size_t i = 1; i < k; i++) { + a *= S[i]; + } + if (a.isZero()) { + *pb = false; + return; + } + /* + f(0) = sum_i f(S[i]) delta_{i,S}(0) + */ + G r; + r.clear(); + for (size_t i = 0; i < k; i++) { + F b = S[i]; + for (size_t j = 0; j < k; j++) { + if (j != i) { + F v = S[j] - S[i]; + if (v.isZero()) { + *pb = false; + return; + } + b *= v; + } + } + G t; + G::mul(t, vec[i], a / b); + r += t; + } + out = r; + *pb = true; +} + +/* + out = f(x) = c[0] + c[1] * x + c[2] * x^2 + ... + c[cSize - 1] * x^(cSize - 1) + @retval 0 if succeed else -1 (if cSize == 0) +*/ +template<class G, class T> +void evaluatePolynomial(bool *pb, G& out, const G *c, size_t cSize, const T& x) +{ + if (cSize == 0) { + *pb = false; + return; + } + if (cSize == 1) { + out = c[0]; + *pb = true; + return; + } + G y = c[cSize - 1]; + for (int i = (int)cSize - 2; i >= 0; i--) { + G::mul(y, y, x); + G::add(y, y, c[i]); + } + out = y; + *pb = true; +} + +#ifndef CYBOZU_DONT_USE_EXCEPTION +template<class G, class F> +void LagrangeInterpolation(G& out, const F *S, const G *vec, size_t k) +{ + bool b; + LagrangeInterpolation(&b, out, S, vec, k); + if (!b) throw cybozu::Exception("LagrangeInterpolation"); +} + +template<class G, class T> +void evaluatePolynomial(G& out, const G *c, size_t cSize, const T& x) +{ + bool b; + evaluatePolynomial(&b, out, c, cSize, x); + if (!b) throw cybozu::Exception("evaluatePolynomial"); +} +#endif + +} // mcl |