ginac/polynomial/newton_interpolate.h

   1 #ifndef GINAC_PGCD_NEWTON_INTERPOLATE_H
   2 #define GINAC_PGCD_NEWTON_INTERPOLATE_H
   3 #include "ex.h"
   4 #include "numeric.h"
   5 #include "smod_helpers.h"
   6 namespace GiNaC
   7 {
   8
   9 /**
  10  * Newton interpolation -- incremental form.
  11  *
  12  * Given a polynomials e1 \in Z_p[Y], prev \in Z_p[x][Y], and evaluation
  13  * points pt1, prevpts \in Z_p compute a polynomial P \in Z_[x][Y] such
  14  * that P(x = pt1) = e1, P(x = pt \in prevpts) = prev(x = pt).
  15  *
  16  * @var{prevpts} enumerates previous evaluation points, should have a form
  17  * (x - pt_2) (x - pt_3) \ldots (x - pt_n).
  18  * @var{prev} encodes the result of previous interpolations.
  19  */
  20 ex newton_interp(const ex& e1, const long pt1,
  21                  const ex& prev, const ex& prevpts,
  22                  const ex& x, const long p)
  23 {
  24         const numeric pnum(p);
  25         const numeric nc = ex_to<numeric>(prevpts.subs(x == numeric(pt1)).smod(pnum));
  26         const numeric nc_1 = recip(nc, p);
  27         // result = prev + prevpts (e1 - prev|_{x = pt_1})/prevpts|_{x = pt_1)
  28         ex tmp = prev.subs(x == numeric(pt1)).smod(pnum);
  29         tmp = (((e1 - tmp).expand().smod(pnum))*nc_1).smod(pnum);
  30         tmp = (prevpts*tmp).expand().smod(pnum);
  31         tmp = (prev + tmp).expand().smod(pnum);
  32         return tmp;
  33 }
  34
  35 }
  36 #endif /* GINAC_PGCD_NEWTON_INTERPOLATE_H */