1 /** @file newton_interpolate.h
3 * Functions for Newton interpolation. */
6 * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 #ifndef GINAC_PGCD_NEWTON_INTERPOLATE_H
24 #define GINAC_PGCD_NEWTON_INTERPOLATE_H
28 #include "smod_helpers.h"
33 * Newton interpolation -- incremental form.
35 * Given a polynomials e1 \in Z_p[Y], prev \in Z_p[x][Y], and evaluation
36 * points pt1, prevpts \in Z_p compute a polynomial P \in Z_[x][Y] such
37 * that P(x = pt1) = e1, P(x = pt \in prevpts) = prev(x = pt).
39 * @var{prevpts} enumerates previous evaluation points, should have a form
40 * (x - pt_2) (x - pt_3) \ldots (x - pt_n).
41 * @var{prev} encodes the result of previous interpolations.
43 ex newton_interp(const ex& e1, const long pt1,
44 const ex& prev, const ex& prevpts,
45 const ex& x, const long p)
47 const numeric pnum(p);
48 const numeric nc = ex_to<numeric>(prevpts.subs(x == numeric(pt1)).smod(pnum));
49 const numeric nc_1 = recip(nc, p);
50 // result = prev + prevpts (e1 - prev|_{x = pt_1})/prevpts|_{x = pt_1)
51 ex tmp = prev.subs(x == numeric(pt1)).smod(pnum);
52 tmp = (((e1 - tmp).expand().smod(pnum))*nc_1).smod(pnum);
53 tmp = (prevpts*tmp).expand().smod(pnum);
54 tmp = (prev + tmp).expand().smod(pnum);
60 #endif // ndef GINAC_PGCD_NEWTON_INTERPOLATE_H