carl  24.04
Computer ARithmetic Library
SqrtEx.tpp
Go to the documentation of this file.
1 /**
2  * Class to create a square root expression object.
3  * @author Florian Corzilius
4  * @since 2011-05-26
5  * @version 2013-10-07
6  */
7 
8 #include <carl-arith/poly/umvpoly/functions/GCD.h>
9 
10 namespace carl
11 {
12 
13  template<typename Poly>
14  SqrtEx<Poly>::SqrtEx():
15  m_constant_part(),
16  m_factor(),
17  m_denominator(1),
18  m_radicand()
19  {}
20 
21  template<typename Poly>
22  SqrtEx<Poly>::SqrtEx( Poly&& _poly ):
23  m_constant_part( std::move( _poly ) ),
24  m_factor(),
25  m_denominator(1),
26  m_radicand()
27  {
28  normalize();
29  }
30 
31  template<typename Poly>
32  SqrtEx<Poly>::SqrtEx( Poly&& _constantPart, Poly&& _factor, Poly&& _denominator, Poly&& _radicand ):
33  m_constant_part( std::move( _constantPart ) ),
34  m_factor( is_zero(_radicand) ? std::move( _radicand ) : std::move( _factor ) ),
35  m_denominator( (is_zero(m_factor) && is_zero(m_constant_part)) ? constant_one<Poly>::get() : std::move( _denominator ) ),
36  m_radicand( is_zero(m_factor) ? m_factor : std::move( _radicand ) )
37  {
38  assert( !is_zero(m_denominator) );
39  assert( !m_radicand.is_constant() || is_zero(m_radicand) || constant_zero<Rational>::get() <= m_radicand.trailingTerm().coeff() );
40  normalize();
41  }
42 
43  template<typename Poly>
44  void SqrtEx<Poly>::normalize()
45  {
46 // std::cout << *this << std::endl;
47  Poly gcdA;
48  if( is_zero(m_factor) )
49  {
50  gcdA = m_constant_part;
51  }
52  else
53  {
54  Poly sqrtOfRadicand;
55  if( m_radicand.sqrt( sqrtOfRadicand ) )
56  {
57  m_constant_part += m_factor * sqrtOfRadicand;
58  m_factor = constant_zero<Poly>::get();
59  m_radicand = constant_zero<Poly>::get();
60  }
61  else
62  {
63  assert( !is_zero(radicand()) );
64  Rational absOfLCoeff = abs( radicand().coprime_factor() );
65  Rational sqrtResult = 0;
66  if( carl::sqrt_exact( absOfLCoeff, sqrtResult ) )
67  {
68  m_factor *= constant_one<Rational>::get()/sqrtResult;
69  m_radicand *= absOfLCoeff;
70  }
71  }
72  if( is_zero(m_factor) )
73  {
74  gcdA = m_constant_part;
75  }
76  else
77  {
78  if( is_zero(m_constant_part) )
79  {
80  gcdA = m_factor;
81  }
82  else
83  {
84  Rational ccConstantPart = m_constant_part.coprime_factor();
85  Poly cpConstantPart = m_constant_part * ccConstantPart;
86  Rational ccFactor = m_factor.coprime_factor();
87  Poly cpFactor = m_factor * ccFactor;
88  gcdA = carl::gcd( cpConstantPart, cpFactor )*carl::gcd(ccConstantPart,ccFactor);
89  }
90  }
91  }
92  if( is_zero(gcdA) ) return;
93  Rational ccGcdA = gcdA.coprime_factor();
94  Poly cpGcdA = gcdA * ccGcdA;
95  Rational ccDenominator = m_denominator.coprime_factor();
96  Poly cpDenominator = m_denominator * ccDenominator;
97  gcdA = carl::gcd( cpGcdA, cpDenominator )*carl::gcd(ccGcdA,ccDenominator);
98  // Make sure that the polynomial to divide by cannot be negative, otherwise the sign of the square root expression could change.
99  if( !(gcdA == constant_one<Poly>::get()) && carl::definiteness(gcdA) == carl::Definiteness::POSITIVE_SEMI )
100  {
101  if( !is_zero(m_constant_part) )
102  {
103  carl::try_divide(m_constant_part, gcdA, m_constant_part );
104  }
105  if( !is_zero(m_factor) )
106  {
107  carl::try_divide(m_factor, gcdA, m_factor );
108  }
109  carl::try_divide(m_denominator, gcdA, m_denominator );
110  }
111  Rational numGcd = constant_zero<Rational>::get();
112  Rational denomLcm = constant_one<Rational>::get();
113  if( is_zero(factor()) )
114  {
115  if( !is_zero(constant_part()) )
116  {
117  Rational cpOfConstantPart = constant_part().coprime_factor();
118  numGcd = carl::get_num( cpOfConstantPart );
119  denomLcm = carl::get_denom( cpOfConstantPart );
120  }
121  else
122  {
123 // std::cout << " to " << *this << std::endl;
124  return; // Sqrt expression corresponds to 0.
125  }
126  }
127  else
128  {
129  Rational cpOfFactorPart = factor().coprime_factor();
130  if( is_zero(constant_part()) )
131  {
132  numGcd = carl::get_num( cpOfFactorPart );
133  denomLcm = carl::get_denom( cpOfFactorPart );
134  }
135  else
136  {
137  Rational cpOfConstantPart = constant_part().coprime_factor();
138  numGcd = carl::gcd( carl::get_num( cpOfConstantPart ), carl::get_num( cpOfFactorPart ) );
139  denomLcm = carl::lcm( carl::get_denom( cpOfConstantPart ), carl::get_denom( cpOfFactorPart ) );
140  }
141  }
142  assert( numGcd != constant_zero<Rational>::get() );
143  Rational cpFactor = numGcd/denomLcm;
144  m_constant_part *= cpFactor;
145  m_factor *= cpFactor;
146  Rational cpOfDenominator = denominator().coprime_factor();
147  m_denominator *= cpOfDenominator;
148  Rational sqrtExFactor = (denomLcm*carl::get_num( cpOfDenominator ))/(numGcd*carl::get_denom( cpOfDenominator ));
149  m_constant_part *= carl::get_num( sqrtExFactor );
150  m_factor *= carl::get_num( sqrtExFactor );
151  m_denominator *= carl::get_denom( sqrtExFactor );
152 // std::cout << " to " << *this << std::endl;
153  //TODO: implement this method further
154  }
155 
156  template<typename Poly>
157  bool SqrtEx<Poly>::operator==( const SqrtEx& _toCompareWith ) const
158  {
159  return m_radicand == _toCompareWith.radicand() && m_denominator == _toCompareWith.denominator()
160  && m_factor == _toCompareWith.factor() && m_constant_part == _toCompareWith.constant_part();
161  }
162 
163  /*
164  template<typename Poly>
165  SqrtEx<Poly>& SqrtEx<Poly>::operator=( const SqrtEx<Poly>& _sqrtEx )
166  {
167  if( this != &_sqrtEx )
168  {
169  m_constant_part = _sqrtEx.constant_part();
170  m_factor = _sqrtEx.factor();
171  m_denominator = _sqrtEx.denominator();
172  if (is_zero(factor()))
173  m_radicand = constant_zero<Poly>::get();
174  else
175  m_radicand = _sqrtEx.radicand();
176  }
177  return *this;
178  }
179  */
180 
181  template<typename Poly>
182  SqrtEx<Poly>& SqrtEx<Poly>::operator=( const Poly& _poly )
183  {
184  m_constant_part = _poly;
185  m_factor = constant_zero<Poly>::get();
186  m_denominator = constant_one<Poly>::get();
187  m_radicand = constant_zero<Poly>::get();
188  return *this;
189  }
190 
191  template<typename Poly>
192  SqrtEx<Poly> SqrtEx<Poly>::operator+( const SqrtEx<Poly>& rhs ) const
193  {
194  assert( !has_sqrt() ||!rhs.has_sqrt() || radicand() == rhs.radicand() );
195  SqrtEx<Poly> result = SqrtEx<Poly>( rhs.denominator() * constant_part() + rhs.constant_part() * denominator(),
196  rhs.denominator() * factor() + rhs.factor() * denominator(),
197  denominator() * rhs.denominator(), radicand() );
198  return result;
199  }
200 
201  template<typename Poly>
202  SqrtEx<Poly> SqrtEx<Poly>::operator-( const SqrtEx<Poly>& rhs ) const
203  {
204  assert( !has_sqrt() || !rhs.has_sqrt() || radicand() == rhs.radicand() );
205  SqrtEx<Poly> result = SqrtEx<Poly>( rhs.denominator() * constant_part() - rhs.constant_part() * denominator(),
206  rhs.denominator() * factor() - rhs.factor() * denominator(),
207  denominator() * rhs.denominator(), radicand() );
208  return result;
209  }
210 
211  template<typename Poly>
212  SqrtEx<Poly> SqrtEx<Poly>::operator*( const SqrtEx<Poly>& rhs ) const
213  {
214  assert( !has_sqrt() || !rhs.has_sqrt() || radicand() == rhs.radicand() );
215  SqrtEx<Poly> result = SqrtEx<Poly>( rhs.constant_part() * constant_part() + rhs.factor() * factor() * radicand(),
216  rhs.constant_part() * factor() + rhs.factor() * constant_part(),
217  denominator() * rhs.denominator(), radicand() );
218  return result;
219  }
220 
221  template<typename Poly>
222  SqrtEx<Poly> SqrtEx<Poly>::operator/( const SqrtEx<Poly>& rhs ) const
223  {
224  assert( !rhs.has_sqrt() );
225  SqrtEx<Poly> result = SqrtEx<Poly>( constant_part() * rhs.denominator(), factor() * rhs.denominator(),
226  denominator() * rhs.factor(), radicand() );
227  return result;
228  }
229 
230  template<typename Poly>
231  inline std::ostream& operator<<( std::ostream& _out, const SqrtEx<Poly>& _substitution )
232  {
233  return (_out << _substitution.toString( true ) );
234  }
235 
236  template<typename Poly>
237  std::string SqrtEx<Poly>::toString( bool _infix, bool /*_friendlyNames*/ ) const
238  {
239  if( _infix )
240  {
241  bool complexNum = has_sqrt() && !m_constant_part.is_constant();
242  std::stringstream result;
243  if( complexNum && !carl::is_one(m_denominator) )
244  result << "(";
245  if( has_sqrt() )
246  {
247  if( m_constant_part.is_constant() )
248  result << m_constant_part;
249  else
250  result << "(" << m_constant_part << ")";
251  result << "+";
252  if( m_factor.is_constant() )
253  result << m_factor;
254  else
255  result << "(" << m_factor << ")";
256  result << "*sqrt(" << m_radicand << ")";
257  }
258  else
259  {
260  if( m_constant_part.is_constant() || carl::is_one(m_denominator))
261  result << m_constant_part;
262  else
263  result << "(" << m_constant_part << ")";
264  }
265  if (!carl::is_one(m_denominator))
266  {
267  if( complexNum )
268  result << ")";
269  result << "/";
270  if( m_denominator.is_constant() )
271  result << m_denominator;
272  else
273  result << "(" << m_denominator << ")";
274  }
275  return result.str();
276  }
277  else
278  {
279  std::stringstream result;
280  result << "(/ (+ ";
281  result << m_constant_part;
282  result << " (* ";
283  result << m_factor;
284  result << " ";
285  result << "(sqrt ";
286  result << m_radicand;
287  result << "))) ";
288  result << m_denominator;
289  result << ")";
290  return result.str();
291  }
292  }
293 } // end namspace vs