ginac/polynomial/newton_interpolate.h

   1 /** @file newton_interpolate.h
   2  *
   3  *  Functions for Newton interpolation. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2014 Johannes Gutenberg University Mainz, Germany
   7  *
   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.
  12  *
  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.
  17  *
  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
  21  */
  22
  23 #ifndef GINAC_PGCD_NEWTON_INTERPOLATE_H
  24 #define GINAC_PGCD_NEWTON_INTERPOLATE_H
  25
  26 #include "ex.h"
  27 #include "numeric.h"
  28 #include "smod_helpers.h"
  29
  30 namespace GiNaC {
  31
  32 /**
  33  * Newton interpolation -- incremental form.
  34  *
  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).
  38  *
  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.
  42  */
  43 ex newton_interp(const ex& e1, const long pt1,
  44                  const ex& prev, const ex& prevpts,
  45                  const ex& x, const long p)
  46 {
  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);
  55         return tmp;
  56 }
  57
  58 } // namespace GiNaC
  59
  60 #endif // ndef GINAC_PGCD_NEWTON_INTERPOLATE_H