00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include <mblad/bap.hpp>
00014 #include <bap.h>
00015
00017 namespace mmx {
00018 #define M vector<nat>
00019 #define Monomial monomial<M>
00020 #define Coordinate multivariate_coordinate<>
00021 #define Coordinates multivariate_coordinates<>
00022 #define Sparse_polynomial_Z sparse_polynomial<integer>
00023 #define Sparse_polynomial_Q sparse_polynomial<rational>
00024 #define Polynomial_Z multivariate<sparse_polynomial<integer> >
00025 #define Polynomial_Q multivariate<sparse_polynomial<rational> >
00026
00027 void*
00028 to_bap_polynomial_mpq (const blad_session* _blad, const Polynomial_Q& eqn) {
00029 vector<Coordinate> indeps= _blad -> indeps;
00030 vector<Coordinate> deps= _blad -> deps;
00031 vector<void*> bav_indeps= _blad -> bav_indeps;
00032 vector<void*> bav_deps= _blad -> bav_deps;
00033 Coordinates crds= eqn.coords;
00034 nat n= N(crds), j;
00035
00036 vector<bav_variable*> eqn_vars;
00037 for (nat i= 0; i < n; i++) {
00038 if ((j= find (indeps, crds[i])) < N(indeps))
00039 eqn_vars << (bav_variable*) bav_indeps[j];
00040 else if ((j= find (deps, crds[i])) < N(deps))
00041 eqn_vars << (bav_variable*) bav_deps[j];
00042 else {
00043 ASSERT (is<compound> (**crds[i]), "unexpected variable");
00044 vector<generic> v= compound_to_vector (**crds[i]);
00045 ASSERT (v[0] == GEN_CACHED_DERIVE, "unexpected variable");
00046 j= find (deps, Coordinate (v[1]));
00047 ASSERT (j < N(deps), "unexpected dependent variable");
00048 bav_variable* u= (bav_variable*) bav_deps[j];
00049 for (nat k= 2; k < N(v); k++) {
00050 j= find (indeps, Coordinate (v[k]));
00051 ASSERT (j < N(indeps), "unexpected independent variable");
00052 u= bav_R_derivative (u, ((bav_variable* ) bav_indeps[j]) -> root);
00053 }
00054 eqn_vars << u;
00055 }
00056 }
00057
00058 bap_polynom_mpq* pol;
00059 bap_polynom_mpq* acc= bap_new_polynom_mpq ();
00060 bap_polynom_mpq* aux= bap_new_polynom_mpq ();
00061 bav_term* term;
00062 bav_term* term_acc= bav_new_term ();
00063 bav_term* term_aux= bav_new_term ();
00064
00065 pol= bap_new_polynom_mpq ();
00066 for (nat i= 0; i < N(eqn); i++) {
00067 rational c= get_coefficient (eqn, i);
00068 M m= *get_monomial (eqn.rep, i);
00069 term= bav_new_term ();
00070 for (nat j= 0; j < N(m); j++) {
00071 nat e= m[j];
00072 bav_set_term_variable (term_aux, eqn_vars[j], m[j]);
00073 bav_mul_term (term_acc, term, term_aux);
00074 bav_set_term (term, term_acc);
00075 }
00076 mpq_t x; mpq_init (x); mpq_set (x, *c);
00077 bap_set_polynom_monom_mpq (aux, x, term);
00078 bap_add_polynom_mpq (acc, pol, aux);
00079 bap_set_polynom_mpq (pol, acc);
00080 }
00081 return (void*) pol; }
00082
00083 Polynomial_Z
00084 from_bap_polynomial_mpz (const blad_session* _blad, const void* pol) {
00085 vector<Coordinate> indeps= _blad -> indeps;
00086 vector<Coordinate> deps= _blad -> deps;
00087 nat r= N(indeps), s= N(deps);
00088 bav_tableof_variable* indets= (bav_tableof_variable*) ba0_new_table ();
00089 bav_R_mark_variables (false);
00090 bap_mark_indets_polynom_mpz ((bap_polynom_mpz*) pol);
00091 bav_R_marked_variables (indets, true);
00092 vector<Coordinate> coords;
00093 for (nat i= 0; i < indets -> size; i++) {
00094 bav_variable* aux= indets -> tab[i];
00095 nat k= aux -> root -> index;
00096 VERIFY (k >= r && k < N(deps) + r, "index out of range");
00097 Coordinate coord= deps [k - r];
00098 vector<Coordinate> ders;
00099 for (nat j= 0; j < r; j++)
00100 ders= append (ders, fill (indeps[j], (aux -> order).tab[j]));
00101 coords << derive (coord, ders);
00102 }
00103 struct bap_itermon_mpz iter;;
00104 bap_begin_itermon_mpz (&iter, (bap_polynom_mpz*) pol);
00105 bav_term* term= bav_new_term ();
00106 vector<integer> coeffs;
00107 vector<Monomial> monoms;
00108 while (! bap_outof_itermon_mpz (&iter)) {
00109 integer cf;
00110 mpz_set (*cf, *bap_coeff_itermon_mpz (&iter));
00111 coeffs << cf;
00112 bap_term_itermon_mpz (term, &iter);
00113 M monom;
00114 for (nat k= 0; k < N(coords); k++)
00115 monom << bav_degree_term (term, indets -> tab[k]);
00116 monoms << Monomial (monom);
00117
00118 bap_next_itermon_mpz (&iter);
00119 }
00120 bap_close_itermon_mpz (&iter);
00121 Sparse_polynomial_Z spol (coeffs, monoms);
00122 return Polynomial_Z (coords, spol); }
00123 }
00124
