ginac/polynomial/euclid_gcd_wrap.h

   1 #ifndef GINAC_PGCD_EUCLID_GCD_H
   2 #define GINAC_PGCD_EUCLID_GCD_H
   3 #include "upoly.hpp"
   4 #include "gcd_euclid.tcc"
   5 #include "smod_helpers.h"
   6 #include "add.h"
   7 #include "ex.h"
   8 #include "operators.h"
   9 #include "power.h"
  10 #include "relational.h"
  11 #include "debug.hpp"
  12
  13 namespace GiNaC
  14 {
  15
  16 static void ex2upoly(umodpoly& u, ex e, const ex& var, const long p)
  17 {
  18         e = e.expand();
  19         cln::cl_modint_ring R = cln::find_modint_ring(cln::cl_I(p));
  20         u.resize(e.degree(var) + 1);
  21         for (int i = 0; i <= e.degree(var); ++i) {
  22                 ex ce = e.coeff(var, i);
  23                 bug_on(!is_a<numeric>(ce), "i = " << i << ", " <<
  24                         "coefficient is not a number: " << ce);
  25                 const cln::cl_I c = to_cl_I(ce);
  26                 u[i] = R->canonhom(c);
  27         }
  28 }
  29
  30 static ex umodpoly2ex(const umodpoly& a, const ex& var, const long p)
  31 {
  32         cln::cl_modint_ring R = cln::find_modint_ring(cln::cl_I(p));
  33         const numeric pnum(p);
  34         exvector ev(a.size());
  35         for (std::size_t i = a.size(); i-- != 0; ) {
  36                 const cln::cl_I c = smod(R->retract(a[i]), p);
  37                 const ex term = numeric(c)*power(var, i);
  38                 ev.push_back(term);
  39         }
  40         ex ret = (new add(ev))->setflag(status_flags::dynallocated);
  41         return ret;
  42 }
  43
  44 static ex euclid_gcd(ex A, ex B, const ex& var, const long p)
  45 {
  46         A = A.expand();
  47         B = B.expand();
  48
  49         umodpoly a, b;
  50         ex2upoly(a, A, var, p);
  51         ex2upoly(b, B, var, p);
  52         umodpoly g;
  53         gcd_euclid(g, a, b);
  54         ex ge = umodpoly2ex(g, var, p);
  55         return ge;
  56 }
  57
  58 } // namespace GiNaC
  59
  60 #endif