ginac/polynomial/euclid_gcd_wrap.h

   1 /** @file euclid_gcd_wrap.h
   2  *
   3  *  Euclidean GCD and supporting functions. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2011 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_EUCLID_GCD_H
  24 #define GINAC_PGCD_EUCLID_GCD_H
  25
  26 #include "upoly.h"
  27 #include "gcd_euclid.h"
  28 #include "smod_helpers.h"
  29 #include "add.h"
  30 #include "ex.h"
  31 #include "operators.h"
  32 #include "power.h"
  33 #include "relational.h"
  34 #include "debug.h"
  35
  36 namespace GiNaC {
  37
  38 static void ex2upoly(umodpoly& u, ex e, const ex& var, const long p)
  39 {
  40         e = e.expand();
  41         cln::cl_modint_ring R = cln::find_modint_ring(cln::cl_I(p));
  42         u.resize(e.degree(var) + 1);
  43         for (int i = 0; i <= e.degree(var); ++i) {
  44                 ex ce = e.coeff(var, i);
  45                 bug_on(!is_a<numeric>(ce), "i = " << i << ", " <<
  46                         "coefficient is not a number: " << ce);
  47                 const cln::cl_I c = to_cl_I(ce);
  48                 u[i] = R->canonhom(c);
  49         }
  50 }
  51
  52 static ex umodpoly2ex(const umodpoly& a, const ex& var, const long p)
  53 {
  54         cln::cl_modint_ring R = cln::find_modint_ring(cln::cl_I(p));
  55         const numeric pnum(p);
  56         exvector ev(a.size());
  57         for (std::size_t i = a.size(); i-- != 0; ) {
  58                 const cln::cl_I c = smod(R->retract(a[i]), p);
  59                 const ex term = numeric(c)*power(var, i);
  60                 ev.push_back(term);
  61         }
  62         ex ret = (new add(ev))->setflag(status_flags::dynallocated);
  63         return ret;
  64 }
  65
  66 static ex euclid_gcd(ex A, ex B, const ex& var, const long p)
  67 {
  68         A = A.expand();
  69         B = B.expand();
  70
  71         umodpoly a, b;
  72         ex2upoly(a, A, var, p);
  73         ex2upoly(b, B, var, p);
  74         umodpoly g;
  75         gcd_euclid(g, a, b);
  76         ex ge = umodpoly2ex(g, var, p);
  77         return ge;
  78 }
  79
  80 } // namespace GiNaC
  81
  82 #endif // ndef GINAC_PGCD_EUCLID_GCD_H