Loading...
Searching...
No Matches
rings.cpp
1#include <symengine/add.h>
2#include <symengine/pow.h>
3#include <symengine/rings.h>
5#include <symengine/symengine_exception.h>
6
7namespace SymEngine
8{
9
10void expr2poly(const RCP<const Basic> &p, umap_basic_num &syms, umap_vec_mpz &P)
11{
12 if (is_a<Add>(*p)) {
13 auto n = syms.size();
14 const umap_basic_num &d = down_cast<const Add &>(*p).get_dict();
16 integer_class coef;
17 for (const auto &p : d) {
18 if (not is_a<Integer>(*p.second))
19 throw NotImplementedError("Not Implemented");
20 coef = down_cast<const Integer &>(*p.second).as_integer_class();
21 exp.assign(n, 0); // Initialize to [0]*n
22 if (is_a<Mul>(*p.first)) {
23 const map_basic_basic &term
24 = down_cast<const Mul &>(*p.first).get_dict();
25 for (const auto &q : term) {
26 RCP<const Basic> sym = q.first;
27 if (not is_a<Integer>(*syms.at(sym)))
28 throw NotImplementedError("Not Implemented");
29 int i = numeric_cast<int>(
30 down_cast<const Integer &>(*syms.at(sym)).as_int());
31 if (is_a<Integer>(*q.second)) {
32 exp[i] = numeric_cast<int>(
33 down_cast<const Integer &>(*q.second).as_int());
34 } else {
35 throw SymEngineException(
36 "Cannot convert symbolic exponents to sparse "
37 "polynomials with integer exponents.");
38 }
39 }
40 } else if (is_a<Pow>(*p.first)) {
41 RCP<const Basic> sym
42 = down_cast<const Pow &>(*p.first).get_base();
43 RCP<const Basic> exp_
44 = down_cast<const Pow &>(*p.first).get_exp();
45 if (not is_a<Integer>(*syms.at(sym)))
46 throw NotImplementedError("Not Implemented");
47 int i = numeric_cast<int>(
48 down_cast<const Integer &>(*syms.at(sym)).as_int());
49 if (not is_a<Integer>(*exp_))
50 throw NotImplementedError("Not Implemented");
51 exp[i] = numeric_cast<int>(
52 down_cast<const Integer &>(*exp_).as_int());
53 } else if (is_a<Symbol>(*p.first)) {
54 RCP<const Basic> sym = p.first;
55 if (not is_a<Integer>(*syms.at(sym)))
56 throw NotImplementedError("Not Implemented");
57 int i = numeric_cast<int>(
58 down_cast<const Integer &>(*syms.at(sym)).as_int());
59 exp[i] = 1;
60 } else {
61 throw NotImplementedError("Not Implemented");
62 }
63
64 P[exp] = coef;
65 }
66 } else {
67 throw NotImplementedError("Not Implemented");
68 }
69}
70
71void poly_mul(const umap_vec_mpz &A, const umap_vec_mpz &B, umap_vec_mpz &C)
72{
74 auto n = A.begin()->first.size();
75 exp.assign(n, 0); // Initialize to [0]*n
76 /*
77 std::cout << "A: " << A.load_factor() << " " << A.bucket_count() << " " <<
78 A.size() << " "
79 << A.max_bucket_count() << std::endl;
80 std::cout << "B: " << B.load_factor() << " " << B.bucket_count() << " " <<
81 B.size() << " "
82 << B.max_bucket_count() << std::endl;
83 std::cout << "C: " << C.load_factor() << " " << C.bucket_count() << " " <<
84 C.size() << " "
85 << C.max_bucket_count() << std::endl;
86 */
87 for (const auto &a : A) {
88 for (const auto &b : B) {
89 monomial_mul(a.first, b.first, exp);
90 mp_addmul(C[exp], a.second, b.second);
91 }
92 }
93 /*
94 std::cout << "C: " << C.load_factor() << " " << C.bucket_count() << " " <<
95 C.size() << " "
96 << C.max_bucket_count() << std::endl;
97 for (const std::size_t n=0; n < C.bucket_count(); n++) {
98 std::cout << n << ": " << C.bucket_size(n) << "|";
99 for (auto it = C.begin(n); it != C.end(n); ++it)
100 std::cout << " " << it->first << myhash2(it->first) %
101 C.bucket_count();
102 std::cout << std::endl;
103 }
104 */
105}
106
107} // namespace SymEngine
Classes and functions relating to the binary operation of addition.
T at(T... args)
T begin(T... args)
Main namespace for SymEngine package.
Definition: add.cpp:19
void monomial_mul(const vec_int &A, const vec_int &B, vec_int &C)
Monomial multiplication.
Definition: monomials.cpp:7
RCP< const Basic > exp(const RCP< const Basic > &x)
Returns the natural exponential function E**x = pow(E, x)
Definition: pow.cpp:271
void expr2poly(const RCP< const Basic > &p, umap_basic_num &syms, umap_vec_mpz &P)
Converts expression p into a polynomial P, with symbols sym
Definition: rings.cpp:10
void poly_mul(const umap_vec_mpz &A, const umap_vec_mpz &B, umap_vec_mpz &C)
Multiply two polynomials: C = A*B
Definition: rings.cpp:71
T size(T... args)