carl  24.04
Computer ARithmetic Library
OldGinacConverter.tpp
Go to the documentation of this file.
1 namespace carl
2 {
3  template<typename Poly>
4  Poly OldGinacConverter<Poly>::convertToCarl(const GiNaC::ex& _toConvert, const std::map<GiNaC::ex, carl::Variable, GiNaC::ex_is_less>& vars)
5  {
6  std::lock_guard<std::recursive_mutex> lock( mMutex );
7  Poly result;
8  GiNaC::ex ginacPoly = _toConvert.expand();
9  if(GiNaC::is_exactly_a<GiNaC::add>(ginacPoly))
10  {
11  result = Poly(typename Poly::CoeffType(0));
12  for(auto summand = ginacPoly.begin(); summand != ginacPoly.end(); ++summand)
13  {
14  const GiNaC::ex summandEx = *summand;
15  if(GiNaC::is_exactly_a<GiNaC::mul>(summandEx))
16  {
17  Poly carlSummand = Poly(typename Poly::CoeffType(1));
18  for(auto factor = summandEx.begin(); factor != summandEx.end(); ++factor)
19  {
20  const GiNaC::ex factorEx = *factor;
21  if(GiNaC::is_exactly_a<GiNaC::symbol>(factorEx))
22  {
23  auto iter = vars.find(factorEx);
24  assert(iter != vars.end());
25  carlSummand *= createPolynomial( iter->second );
26  }
27  else if(GiNaC::is_exactly_a<GiNaC::numeric>(factorEx))
28  {
29  carlSummand *= carl::convert<cln::cl_RA, typename Poly::CoeffType>(cln::rationalize(cln::realpart(GiNaC::ex_to<GiNaC::numeric>(factorEx).to_cl_N())));
30  }
31  else if(GiNaC::is_exactly_a<GiNaC::power>(factorEx))
32  {
33  assert(factorEx.nops() == 2);
34  GiNaC::ex exponent = *(++(factorEx.begin()));
35  assert(!exponent.info(GiNaC::info_flags::negative));
36  unsigned exp = static_cast<unsigned>(exponent.integer_content().to_int());
37  GiNaC::ex subterm = *factorEx.begin();
38  assert(GiNaC::is_exactly_a<GiNaC::symbol>(subterm));
39  auto iter = vars.find(subterm);
40  assert(iter != vars.end());
41  carlSummand *= carl::pow(createPolynomial(iter->second), exp);
42  }
43  else assert(false);
44  }
45  result += carlSummand;
46  }
47  else if(GiNaC::is_exactly_a<GiNaC::symbol>(summandEx))
48  {
49  auto iter = vars.find(summandEx);
50  assert(iter != vars.end());
51  result += createPolynomial(iter->second);
52  }
53  else if(GiNaC::is_exactly_a<GiNaC::numeric>(summandEx))
54  {
55  result += carl::convert<cln::cl_RA, typename Poly::CoeffType>(cln::rationalize(cln::realpart(GiNaC::ex_to<GiNaC::numeric>(summandEx).to_cl_N())));
56  }
57  else if(GiNaC::is_exactly_a<GiNaC::power>(summandEx))
58  {
59  assert(summandEx.nops() == 2);
60  GiNaC::ex exponent = *(++(summandEx.begin()));
61  assert(!exponent.info( GiNaC::info_flags::negative));
62  unsigned exp = static_cast<unsigned>(exponent.integer_content().to_int());
63  GiNaC::ex subterm = *summandEx.begin();
64  assert(GiNaC::is_exactly_a<GiNaC::symbol>(subterm));
65  auto iter = vars.find(subterm);
66  assert(iter != vars.end());
67  result += carl::pow(createPolynomial(iter->second), exp);
68  }
69  else assert(false);
70  }
71  }
72  else if(GiNaC::is_exactly_a<GiNaC::mul>(ginacPoly))
73  {
74  result = Poly(typename Poly::CoeffType(1));
75  for(auto factor = ginacPoly.begin(); factor != ginacPoly.end(); ++factor)
76  {
77  const GiNaC::ex factorEx = *factor;
78  if(GiNaC::is_exactly_a<GiNaC::symbol>(factorEx))
79  {
80  auto iter = vars.find(factorEx);
81  assert(iter != vars.end());
82  result *= createPolynomial(iter->second);
83  }
84  else if(GiNaC::is_exactly_a<GiNaC::numeric>(factorEx))
85  {
86  result *= carl::convert<cln::cl_RA, typename Poly::CoeffType>(cln::rationalize(cln::realpart(GiNaC::ex_to<GiNaC::numeric>(factorEx).to_cl_N())));
87  }
88  else if(GiNaC::is_exactly_a<GiNaC::power>(factorEx))
89  {
90  assert(factorEx.nops() == 2);
91  GiNaC::ex exponent = *(++(factorEx.begin()));
92  assert(!exponent.info(GiNaC::info_flags::negative));
93  unsigned exp = static_cast<unsigned>(exponent.integer_content().to_int());
94  GiNaC::ex subterm = *factorEx.begin();
95  assert(GiNaC::is_exactly_a<GiNaC::symbol>(subterm));
96  auto iter = vars.find(subterm);
97  assert(iter != vars.end());
98  result *= carl::pow(createPolynomial(iter->second), exp);
99  }
100  else assert( false );
101  }
102  }
103  else if(GiNaC::is_exactly_a<GiNaC::symbol>(ginacPoly))
104  {
105  auto iter = vars.find(ginacPoly);
106  assert(iter != vars.end());
107  result = createPolynomial(iter->second);
108  }
109  else if(GiNaC::is_exactly_a<GiNaC::numeric>(ginacPoly))
110  {
111  result = Poly(carl::convert<cln::cl_RA, typename Poly::CoeffType>(cln::rationalize(cln::realpart( GiNaC::ex_to<GiNaC::numeric>(ginacPoly).to_cl_N()))));
112  }
113  else if(GiNaC::is_exactly_a<GiNaC::power>(ginacPoly))
114  {
115  assert(ginacPoly.nops() == 2);
116  GiNaC::ex exponent = *(++ginacPoly.begin());
117  assert(!exponent.info(GiNaC::info_flags::negative));
118  unsigned exp = static_cast<unsigned>(exponent.integer_content().to_int());
119  GiNaC::ex subterm = *ginacPoly.begin();
120  assert(GiNaC::is_exactly_a<GiNaC::symbol>(subterm));
121  auto iter = vars.find(subterm);
122  assert(iter != vars.end());
123  result = carl::pow(createPolynomial(iter->second), exp);
124  }
125  else assert( false );
126  return result;
127  }
128 
129  template<typename Poly>
130  Poly OldGinacConverter<Poly>::ginacGcd(const Poly& polyA, const Poly& polyB)
131  {
132  std::lock_guard<std::recursive_mutex> lock( mMutex );
133  Poly result;
134  std::map<Variable, GiNaC::ex> carlToGinacVarMap;
135  std::map<GiNaC::ex, Variable, GiNaC::ex_is_less> ginacToCarlVarMap;
136  gatherVariables(polyA, carlToGinacVarMap, ginacToCarlVarMap);
137  gatherVariables(polyB, carlToGinacVarMap, ginacToCarlVarMap);
138  GiNaC::ex ginacResult = GiNaC::gcd(convertToGinac(polyA, carlToGinacVarMap), convertToGinac(polyB, carlToGinacVarMap));
139  result = convertToCarl(ginacResult, ginacToCarlVarMap);
140  if( !carl::is_zero(result) && result.lcoeff() < carl::constant_zero<typename Poly::CoeffType>().get() )
141  return -result;
142  return result;
143  }
144 
145  template<typename Poly>
146  bool OldGinacConverter<Poly>::checkConversion(const Poly& polyA)
147  {
148  std::lock_guard<std::recursive_mutex> lock( mMutex );
149  std::map<Variable, GiNaC::ex> carlToGinacVarMap;
150  std::map<GiNaC::ex, Variable, GiNaC::ex_is_less> ginacToCarlVarMap;
151  gatherVariables(polyA, carlToGinacVarMap, ginacToCarlVarMap);
152  Poly result = convertToCarl(convertToGinac(polyA, carlToGinacVarMap), ginacToCarlVarMap);
153  return polyA == result;
154  }
155 
156  template<typename Poly>
157  bool OldGinacConverter<Poly>::ginacDivide(const Poly& polyA, const Poly& polyB, Poly& result)
158  {
159  std::lock_guard<std::recursive_mutex> lock( mMutex );
160  std::map<Variable, GiNaC::ex> carlToGinacVarMap;
161  std::map<GiNaC::ex, Variable, GiNaC::ex_is_less> ginacToCarlVarMap;
162  gatherVariables(polyA, carlToGinacVarMap, ginacToCarlVarMap);
163  gatherVariables(polyB, carlToGinacVarMap, ginacToCarlVarMap);
164  GiNaC::ex ginacResult;
165  bool divided = GiNaC::divide(convertToGinac(polyA, carlToGinacVarMap), convertToGinac(polyB, carlToGinacVarMap), ginacResult);
166  if(divided)
167  result = convertToCarl(ginacResult, ginacToCarlVarMap);
168  return divided;
169  }
170 
171  template<typename Poly>
172  Factors<Poly> OldGinacConverter<Poly>::ginacFactorization(const Poly& poly)
173  {
174  std::lock_guard<std::recursive_mutex> lock( mMutex );
175  Factors<Poly> result;
176  std::map<Variable, GiNaC::ex> carlToGinacVarMap;
177  std::map<GiNaC::ex, Variable, GiNaC::ex_is_less> ginacToCarlVarMap;
178  gatherVariables(poly, carlToGinacVarMap, ginacToCarlVarMap);
179  GiNaC::ex ginacResult = GiNaC::factor(convertToGinac(poly, carlToGinacVarMap));
180  if(GiNaC::is_exactly_a<GiNaC::mul>(ginacResult))
181  {
182  for(auto factor = ginacResult.begin(); factor != ginacResult.end(); ++factor)
183  {
184  const GiNaC::ex& factorEx = *factor;
185  if(GiNaC::is_exactly_a<GiNaC::power>(factorEx))
186  {
187  assert(factorEx.nops() == 2);
188  GiNaC::ex exponent = *(++factorEx.begin());
189  assert(!exponent.info(GiNaC::info_flags::negative));
190  unsigned exp = static_cast<unsigned>(exponent.integer_content().to_int());
191  GiNaC::ex subterm = *factorEx.begin();
192  Poly carlFactor = convertToCarl(subterm, ginacToCarlVarMap);
193  assert(result.find(carlFactor) == result.end());
194  result.insert(std::pair<Poly, unsigned>(carlFactor, exp));
195  }
196  else
197  {
198  Poly carlFactor = convertToCarl(factorEx, ginacToCarlVarMap);
199  assert(result.find(carlFactor) == result.end());
200  result.insert(std::pair<Poly, unsigned>(carlFactor, 1));
201  }
202  }
203  }
204  else if(GiNaC::is_exactly_a<GiNaC::power>(ginacResult))
205  {
206  assert(ginacResult.nops() == 2);
207  GiNaC::ex exponent = *(++ginacResult.begin());
208  assert(!exponent.info(GiNaC::info_flags::negative));
209  unsigned exp = static_cast<unsigned>(exponent.integer_content().to_int());
210  GiNaC::ex subterm = *ginacResult.begin();
211  result.insert(std::pair<Poly, unsigned>(convertToCarl(subterm, ginacToCarlVarMap), exp));
212  }
213  else
214  {
215  result.insert(std::pair<Poly, unsigned>(convertToCarl(ginacResult, ginacToCarlVarMap), 1));
216  }
217  return result;
218  }
219 } // end namespace carl