A small kernel of code for playing with Galois fields of arbitrary characteristic
Revision | a6fa33805b0aa13effbb5bd9202d50fb4e5d0035 (tree) |
---|---|
Time | 2017-06-17 02:12:59 |
Author | Eric Hopper <hopper@omni...> |
Commiter | Eric Hopper |
Add a function to do Lagrange interpolation.
@@ -168,3 +168,32 @@ | ||
168 | 168 | result += a_coef * xpower |
169 | 169 | xpower *= x |
170 | 170 | return result |
171 | + | |
172 | +def lagrange_interp(points): | |
173 | + from functools import reduce | |
174 | + there_is_no_one = object() | |
175 | + | |
176 | + def determine_mult_ident_and_numpts(points): | |
177 | + numpts = 0 | |
178 | + one = there_is_no_one | |
179 | + for x, y in points: | |
180 | + if (one is there_is_no_one) and (type(x)() != x): | |
181 | + one = x / x | |
182 | + numpts += 1 | |
183 | + return (one, numpts) | |
184 | + | |
185 | + (xone, numpts) = determine_mult_ident_and_numpts(points) | |
186 | + if xone is there_is_no_one: | |
187 | + if numpts > 1: | |
188 | + raise ValueError("Impossible to interpolate.") | |
189 | + else: | |
190 | + return (next(iter(points))[1],) | |
191 | + xzero = type(next(iter(points))[0])() | |
192 | + master_poly = reduce(polymul, | |
193 | + ((xzero - x, 1) for x, y in points), | |
194 | + (xone,)) | |
195 | + xpolys = ((polydiv(master_poly, (xzero - x, xone))[0], x, y) | |
196 | + for x, y in points) | |
197 | + xpolys = (polyscalarmul(xp, y / evalpoly(xp, x)) | |
198 | + for xp, x, y in xpolys) | |
199 | + return normalized(reduce(polyadd, xpolys, ())) |