1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
|
#pragma once
/**
@file
@brief Lagrange Interpolation
@author MITSUNARI Shigeo(@herumi)
@license modified new BSD license
http://opensource.org/licenses/BSD-3-Clause
*/
#include <vector>
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(G& out, const F *S, const G *vec, size_t k)
{
/*
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])
*/
if (k < 2) throw cybozu::Exception("LagrangeInterpolation:smalll k") << k;
std::vector<F> delta(k);
F a = S[0];
for (size_t i = 1; i < k; i++) {
a *= S[i];
}
if (a.isZero()) throw cybozu::Exception("LagrangeInterpolation:S has zero");
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()) throw cybozu::Exception("LagrangeInterpolation:same S") << i << j;
b *= v;
}
}
delta[i] = a / b;
}
/*
f(0) = sum_i f(S[i]) delta_{i,S}(0)
*/
G r, t;
r.clear();
for (size_t i = 0; i < delta.size(); i++) {
G::mul(t, vec[i], delta[i]);
r += t;
}
out = r;
}
/*
out = f(x) = c[0] + c[1] * x + c[2] * x^2 + ... + c[cSize - 1] * x^(cSize - 1)
@retval 0 if succeed else -1
*/
template<class G, class T>
void evaluatePolynomial(G& out, const G *c, size_t cSize, const T& x)
{
if (cSize < 2) throw cybozu::Exception("evaluatePolynomial:small cSize") << cSize;
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;
}
} // mcl
|
