carl  24.04
Computer ARithmetic Library
FactorizedPolynomial.tpp
Go to the documentation of this file.
1 /*
2  * File: FactorizedPolynomial.tpp
3  * Author: Florian Corzilius
4  *
5  * Created on September 4, 2014, 10:55 AM
6  */
7 
8 #include "FactorizedPolynomial.h"
9 
10 #include <carl-arith/poly/umvpoly/functions/Division.h>
11 #include <carl-arith/poly/umvpoly/functions/Substitution.h>
12 #include <carl-arith/poly/umvpoly/UnivariatePolynomial.h>
13 
14 #pragma once
15 
16 namespace carl
17 {
18  template<typename P>
19  FactorizedPolynomial<P>::FactorizedPolynomial():
20  mCacheRef( CACHE::NO_REF ),
21  mpCache( nullptr ),
22  mCoefficient( 0 )
23  {
24  assert( mpCache == nullptr || mCacheRef != CACHE::NO_REF );
25  }
26 
27  template<typename P>
28  FactorizedPolynomial<P>::FactorizedPolynomial( const CoeffType& _coefficient ):
29  mCacheRef( CACHE::NO_REF ),
30  mpCache( nullptr ),
31  mCoefficient( _coefficient )
32  {
33  }
34 
35  template<typename P>
36  FactorizedPolynomial<P>::FactorizedPolynomial( const P& _polynomial, const std::shared_ptr<CACHE>& _pCache, bool _polyNormalized ):
37  mCacheRef( CACHE::NO_REF ),
38  mpCache( carl::is_zero(_polynomial) ? nullptr : _pCache ),
39  mCoefficient( carl::is_zero(_polynomial) ? CoeffType(0) : (_polyNormalized ? CoeffType(1) : CoeffType(1)/_polynomial.coprime_factor()) )
40  {
41  assert( !_polyNormalized || (_polynomial.coprime_factor() == CoeffType(1)) );
42  if ( _polynomial.is_constant() )
43  {
44  mpCache = nullptr;
45  }
46  else
47  {
48  P poly = _polyNormalized ? _polynomial : P(_polynomial * (CoeffType(1) / mCoefficient));
49  if ( !_polyNormalized && poly.lcoeff() < 0 )
50  {
51  poly *= CoeffType(-1);
52  mCoefficient *= CoeffType(-1);
53  }
54  /*
55  * The following is not very nice, but we know, that the hash won't change, once the polynomial
56  * representation is fixed, so we can add the factorizations content belatedly. It is necessary to do so
57  * as otherwise the factorized polynomial (this) being the only factor, is not yet cached which leads to an assertion.
58  */
59  if ( _polyNormalized || carl::is_one(mCoefficient) )
60  {
61  assert( mpCache != nullptr );
62  Factorization<P> factorization;
63  PolynomialFactorizationPair<P>* pfPair = new PolynomialFactorizationPair<P>( std::move( factorization), new P(poly) );
64  //Factorization is not set yet
65  auto ret = mpCache->cache( pfPair );//, &carl::canBeUpdated, &carl::update );
66  mCacheRef = ret.first;
67  mpCache->reg( mCacheRef );
68  if( ret.second )
69  {
70  assert( content().mFactorization.empty() );
71  CARL_LOG_DEBUG("carl.core.factorizedpolynomial", "Adding single factor ( " << poly << " )^1");
72  content().mFactorization.insert( std::make_pair( *this, 1 ) );
73  }
74  else
75  {
76  delete pfPair;
77  }
78  }
79  else
80  {
81  // Factor is polynomial without coefficient
82  FactorizedPolynomial factor( poly, mpCache, true );
83  mCacheRef = factor.mCacheRef;
84  mpCache->reg( mCacheRef );
85  }
86  //We can not check the factorization yet, but as we have set it, it should be correct.
87  //pfPair->assertFactorization();
88  }
89  ASSERT_CACHE_REF_LEGAL( (*this) );
90  assert(computePolynomial(*this) == _polynomial);
91  }
92 
93  template<typename P>
94  FactorizedPolynomial<P>::FactorizedPolynomial( Factorization<P>&& _factorization, const CoeffType& _coefficient, const std::shared_ptr<CACHE>& _pCache ):
95  mCacheRef( CACHE::NO_REF ),
96  mpCache( _pCache ),
97  mCoefficient( _coefficient )
98  {
99  assert( !carl::is_zero(_coefficient) );
100  if ( _factorization.empty() )
101  mpCache = nullptr;
102  else
103  {
104  assert( mpCache != nullptr );
105  // TODO expensive
106  for ( auto factor = _factorization.begin(); factor != _factorization.end(); factor++ )
107  assert( carl::is_one(factor->first.coefficient()) );
108  PolynomialFactorizationPair<P>* pfPair = new PolynomialFactorizationPair<P>( std::move( _factorization ) );
109  auto ret = mpCache->cache( pfPair );//, &carl::canBeUpdated, &carl::update );
110  mCacheRef = ret.first;
111  mpCache->reg( mCacheRef );
112  if( !ret.second )
113  {
114  delete pfPair;
115  }
116  }
117  ASSERT_CACHE_REF_LEGAL( (*this) );
118  }
119 
120  template<typename P>
121  FactorizedPolynomial<P>::FactorizedPolynomial( const FactorizedPolynomial<P>& _toCopy ):
122  mCacheRef( _toCopy.mCacheRef ),
123  mpCache( _toCopy.mpCache ),
124  mCoefficient( _toCopy.mCoefficient )
125  {
126  if ( mpCache != nullptr )
127  {
128  mpCache->reg( mCacheRef );
129  }
130  ASSERT_CACHE_REF_LEGAL( (*this) );
131  }
132 
133  template<typename P>
134  FactorizedPolynomial<P>::FactorizedPolynomial( FactorizedPolynomial<P>&& rhs ):
135  mCacheRef( rhs.mCacheRef ),
136  mpCache( rhs.mpCache ),
137  mCoefficient( rhs.mCoefficient )
138  {
139  ASSERT_CACHE_REF_LEGAL( (*this) );
140  rhs.mCacheRef = CACHE::NO_REF;
141  rhs.mpCache = nullptr;
142  rhs.mCoefficient = 0;
143  }
144 
145  template<typename P>
146  FactorizedPolynomial<P>::FactorizedPolynomial(const std::pair<ConstructorOperation, std::vector<FactorizedPolynomial>>& _pair):
147  FactorizedPolynomial<P>::FactorizedPolynomial()
148  {
149  auto op = _pair.first;
150  auto sub = _pair.second;
151  assert(!sub.empty());
152  auto it = sub.begin();
153 
154  *this = *it;
155  if ((op == ConstructorOperation::SUB) && (sub.size() == 1)) {
156  // special treatment of unary minus
157  *this = -(*this);
158  return;
159  }
160  if( op == ConstructorOperation::DIV )
161  {
162  // division shall have at least two arguments
163  assert(sub.size() >= 2);
164  }
165  for( ++it; it != sub.end(); ++it)
166  {
167  switch( op )
168  {
169  case ConstructorOperation::ADD:
170  *this += *it;
171  break;
172  case ConstructorOperation::SUB:
173  *this -= *it;
174  break;
175  case ConstructorOperation::MUL:
176  *this *= *it;
177  break;
178  case ConstructorOperation::DIV:
179  assert(it->is_constant());
180  *this /= it->constant_part();
181  break;
182  }
183  }
184  }
185 
186  template<typename P>
187  FactorizedPolynomial<P>::~FactorizedPolynomial()
188  {
189  if ( mpCache != nullptr )
190  {
191  mpCache->dereg( mCacheRef );
192  mpCache = nullptr;
193  }
194  }
195 
196  template<typename P>
197  FactorizedPolynomial<P>& FactorizedPolynomial<P>::operator=( const FactorizedPolynomial<P>& _fpoly )
198  {
199  CARL_LOG_DEBUG("carl.core.factorizedpolynomial", "Copying " << _fpoly);
200  ASSERT_CACHE_EQUAL( mpCache, _fpoly.pCache() );
201  mCoefficient = _fpoly.mCoefficient;
202  if( mCacheRef != _fpoly.cacheRef() )
203  {
204  if( mpCache != nullptr )
205  {
206  mpCache->dereg( mCacheRef );
207  mCacheRef = _fpoly.cacheRef();
208  assert( _fpoly.pCache() == nullptr || mCacheRef != CACHE::NO_REF );
209  if( _fpoly.pCache() != nullptr )
210  mpCache->reg( mCacheRef );
211  else
212  mpCache = nullptr;
213  }
214  else if( _fpoly.mpCache != nullptr )
215  {
216  mpCache = _fpoly.mpCache;
217  mCacheRef = _fpoly.cacheRef();
218  mpCache->reg( mCacheRef );
219  }
220  }
221  else if( mpCache != nullptr )
222  {
223  mpCache->dereg(mCacheRef);
224  mCacheRef = _fpoly.cacheRef();
225  mpCache->reg( mCacheRef );
226  }
227  ASSERT_CACHE_REF_LEGAL( (*this) );
228  assert(computePolynomial(*this) == computePolynomial(_fpoly));
229  CARL_LOG_DEBUG("carl.core.factorizedpolynomial", "Done.");
230  return *this;
231  }
232 
233  template<typename P>
234  template<typename C, EnableIf<is_subset_of_rationals_type<C>>>
235  typename FactorizedPolynomial<P>::CoeffType FactorizedPolynomial<P>::coprime_factor_without_constant() const
236  {
237  if( existsFactorization( *this ) )
238  {
239  if( factorizedTrivially() )
240  {
241  assert( computePolynomial(*this).coprime_factor_without_constant() == polynomial().coprime_factor_without_constant() / mCoefficient );
242  return polynomial().coprime_factor_without_constant() / mCoefficient;
243  }
244  auto factor = content().factorization().begin();
245  CoeffType cf = factor->first.coprime_factor_without_constant();
246  auto resultNum = carl::get_num( cf );
247  auto resultDenom = carl::get_denom( cf );
248  ++factor;
249  while( factor != content().factorization().end() )
250  {
251  cf = factor->first.coprime_factor_without_constant();
252  resultNum = carl::lcm( resultNum, carl::get_num( cf ) );
253  resultDenom = carl::gcd( resultDenom, carl::get_denom( cf ) );
254  ++factor;
255  }
256  assert( computePolynomial(*this).coprime_factor_without_constant() == (resultNum / (mCoefficient*resultDenom)) );
257  return resultNum / (mCoefficient*resultDenom);
258  }
259  return constant_zero<CoeffType>::get();
260  }
261 
262  template<typename P>
263  typename FactorizedPolynomial<P>::CoeffType FactorizedPolynomial<P>::constant_part() const
264  {
265  if( existsFactorization( *this ) )
266  {
267  if( factorizedTrivially() )
268  return polynomial().constant_part() * mCoefficient;
269  CoeffType result( mCoefficient );
270  for( const auto& factor : content().factorization() )
271  {
272  result *= carl::pow( factor.first.constant_part(), factor.second );
273  }
274  assert( computePolynomial(*this).constant_part() == result );
275  return result;
276  }
277  return mCoefficient;
278  }
279 
280  template<typename P>
281  typename FactorizedPolynomial<P>::CoeffType FactorizedPolynomial<P>::lcoeff() const
282  {
283  if( existsFactorization( *this ) )
284  {
285  if( factorizedTrivially() )
286  return polynomial().lcoeff() * mCoefficient;
287  CoeffType result( mCoefficient );
288  for( const auto& factor : content().factorization() )
289  {
290  result *= carl::pow( factor.first.lcoeff(), factor.second );
291  }
292  assert( computePolynomial(*this).lcoeff() == result );
293  return result;
294  }
295  return mCoefficient;
296  }
297 
298  template<typename P>
299  size_t FactorizedPolynomial<P>::total_degree() const
300  {
301  if( existsFactorization( *this ) )
302  {
303  if( factorizedTrivially() )
304  return polynomial().total_degree();
305  exponent result = 0;
306  for( const auto& factor : content().factorization() )
307  {
308  result += factor.first.total_degree() * factor.second;
309  }
310  assert( computePolynomial(*this).total_degree() == result );
311  return result;
312  }
313  assert( !carl::is_zero(mCoefficient) );
314  return 0;
315  }
316 
317  template<typename P>
318  typename FactorizedPolynomial<P>::TermType FactorizedPolynomial<P>::lterm() const
319  {
320  if( existsFactorization( *this ) )
321  {
322  if( factorizedTrivially() )
323  return polynomial().lterm() * mCoefficient;
324  TermType result( mCoefficient );
325  for( const auto& factor : content().factorization() )
326  {
327  result *= factor.first.lterm().pow( factor.second );
328  }
329  assert( computePolynomial(*this).lterm() == result );
330  return result;
331  }
332  return TermType( mCoefficient );
333  }
334 
335  template<typename P>
336  typename FactorizedPolynomial<P>::TermType FactorizedPolynomial<P>::trailingTerm() const
337  {
338  if( existsFactorization( *this ) )
339  {
340  if( factorizedTrivially() )
341  {
342  return polynomial().trailingTerm() * mCoefficient;
343  }
344  TermType result( mCoefficient );
345  for( const auto& factor : content().factorization() )
346  {
347  result *= factor.first.trailingTerm().pow( factor.second );
348  }
349  assert( computePolynomial(*this).trailingTerm() == result );
350  return result;
351  }
352  return TermType( mCoefficient );
353  }
354 
355  template<typename P>
356  bool FactorizedPolynomial<P>::is_univariate() const
357  {
358  if( existsFactorization( *this ) )
359  {
360  if( factorizedTrivially() )
361  {
362  return polynomial().is_univariate();
363  }
364  Variable foundVar = Variable::NO_VARIABLE;
365  for( const auto& factor : content().factorization() )
366  {
367  if( factor.first.is_univariate() )
368  {
369  if( foundVar == Variable::NO_VARIABLE )
370  {
371  foundVar = factor.first.single_variable();
372  }
373  else if( foundVar != factor.first.single_variable() )
374  {
375  return false;
376  }
377  }
378  else
379  return false;
380  }
381  return true;
382  }
383  return true;
384  }
385 
386  template<typename P>
387  UnivariatePolynomial<FactorizedPolynomial<P>> FactorizedPolynomial<P>::toUnivariatePolynomial( Variable _var ) const
388  {
389  // TODO: This should maybe done directly on the factorization.
390  UnivariatePolynomial<P> univPol = polynomial().toUnivariatePolynomial( _var );
391  std::vector<FactorizedPolynomial<P>> resultCoeffs;
392  for( const P& coeff : univPol.coefficients() )
393  {
394  resultCoeffs.push_back( std::move(FactorizedPolynomial<P>( coeff, mpCache ) *= mCoefficient) );
395  }
396  return UnivariatePolynomial<FactorizedPolynomial<P>>( _var, std::move( resultCoeffs ) );
397  }
398 
399  template<typename P>
400  bool FactorizedPolynomial<P>::has_constant_term() const
401  {
402  if( existsFactorization( *this ) )
403  {
404  assert( !carl::is_zero( mCoefficient ) );
405  if( factorizedTrivially() )
406  return polynomial().has_constant_term();
407  TermType result( mCoefficient );
408  for( const auto& factor : content().factorization() )
409  {
410  if( !factor.first.has_constant_term() )
411  {
412  assert( !computePolynomial(*this).has_constant_term() );
413  return false;
414  }
415  }
416  assert( computePolynomial(*this).has_constant_term() );
417  return true;
418  }
419  return !carl::is_zero( mCoefficient );
420  }
421 
422  template<typename P>
423  bool FactorizedPolynomial<P>::has( Variable _var ) const
424  {
425  if( existsFactorization( *this ) )
426  {
427  if( factorizedTrivially() )
428  return polynomial().has( _var );
429  for( const auto& factor : content().factorization() )
430  {
431  if( factor.first.has( _var ) )
432  return true;
433  }
434  return false;
435  }
436  return false;
437  }
438 
439  template<typename P>
440  template<bool gatherCoeff>
441  VarInfo<FactorizedPolynomial<P>> FactorizedPolynomial<P>::var_info(Variable var) const
442  {
443  // TODO: Maybe we should use the factorization for collecting degrees and coefficients.
444  VarInfo<P> vi = carl::var_info(polynomial(), var, gatherCoeff);
445  if( gatherCoeff )
446  {
447  std::map<unsigned,FactorizedPolynomial<P>> coeffs;
448  for( auto const& expCoeffPair : vi.coeffs() )
449  {
450  if( expCoeffPair.second.is_constant() )
451  {
452  coeffs.insert( coeffs.end(), std::make_pair( expCoeffPair.first, FactorizedPolynomial<P>( expCoeffPair.second.constant_part() * mCoefficient ) ) );
453  }
454  else
455  {
456  coeffs.insert( coeffs.end(), std::make_pair( expCoeffPair.first, FactorizedPolynomial<P>( expCoeffPair.second, mpCache ) * mCoefficient ) );
457  }
458  }
459  return VarInfo<FactorizedPolynomial<P>>( vi.max_degree(), vi.min_degree(), vi.num_occurences(), std::move( coeffs ) );
460  }
461  return VarInfo<FactorizedPolynomial<P>>( vi.max_degree(), vi.min_degree(), vi.num_occurences() );
462 
463  }
464 
465  template<typename P>
466  Definiteness FactorizedPolynomial<P>::definiteness( bool /*_fullEffort*/ ) const
467  {
468  if( existsFactorization( *this ) )
469  {
470  if( factorizedTrivially() )
471  return polynomial().definiteness();
472  Definiteness result = Definiteness::POSITIVE;
473  for( const auto& factor : content().factorization() )
474  {
475  Definiteness factorDefiniteness = factor.first.definiteness();
476  if( factorDefiniteness == Definiteness::NON )
477  return Definiteness::NON;
478  if( factor.second % 2 == 0 )
479  {
480  if( factorDefiniteness == Definiteness::NEGATIVE_SEMI )
481  factorDefiniteness = Definiteness::POSITIVE_SEMI;
482  if( factorDefiniteness == Definiteness::NEGATIVE )
483  factorDefiniteness = Definiteness::POSITIVE;
484  }
485  switch( result )
486  {
487  case Definiteness::POSITIVE:
488  {
489  result = factorDefiniteness;
490  break;
491  }
492  case Definiteness::POSITIVE_SEMI:
493  {
494  if( result == Definiteness::NEGATIVE_SEMI && factorDefiniteness == Definiteness::NEGATIVE )
495  {
496  result = Definiteness::NEGATIVE_SEMI;
497  }
498  break;
499  }
500  case Definiteness::NEGATIVE:
501  {
502  switch( factorDefiniteness )
503  {
504  case Definiteness::NEGATIVE_SEMI:
505  {
506  result = Definiteness::POSITIVE_SEMI;
507  break;
508  }
509  case Definiteness::NEGATIVE:
510  {
511  result = Definiteness::POSITIVE;
512  break;
513  }
514  case Definiteness::POSITIVE:
515  {
516  result = Definiteness::NEGATIVE;
517  break;
518  }
519  default:
520  {
521  assert( factorDefiniteness == Definiteness::POSITIVE_SEMI );
522  result = Definiteness::NEGATIVE_SEMI;
523  break;
524  }
525  }
526  break;
527  }
528  default:
529  {
530  assert( result == Definiteness::NEGATIVE_SEMI );
531  if( result == Definiteness::NEGATIVE_SEMI && factorDefiniteness == Definiteness::NEGATIVE )
532  {
533  result = Definiteness::POSITIVE_SEMI;
534  }
535  break;
536  }
537  }
538  }
539  assert( result == Definiteness::NON || computePolynomial(*this).definiteness() == Definiteness::NON || result == computePolynomial(*this).definiteness() );
540  return result;
541  }
542  return mCoefficient < 0 ? Definiteness::NEGATIVE : (mCoefficient > 0 ? Definiteness::POSITIVE : Definiteness::POSITIVE_SEMI);
543  }
544 
545  template<typename P>
546  FactorizedPolynomial<P> FactorizedPolynomial<P>::derivative( const carl::Variable& _var, unsigned _nth ) const
547  {
548  assert(_nth == 1);
549  if (this->is_constant()) {
550  return FactorizedPolynomial<P>( constant_zero<CoeffType>::get() );
551  }
552  // TODO VERY NAIVE
553  FactorizedPolynomial<P> result(carl::derivative(polynomial(), _var), mpCache);
554  result *= coefficient();
555  return result;
556  }
557 
558  template<typename P>
559  FactorizedPolynomial<P> FactorizedPolynomial<P>::pow(unsigned _exp) const
560  {
561  if( _exp == 0 )
562  return FactorizedPolynomial<P>( constant_one<CoeffType>::get() );
563  if( is_zero() )
564  return FactorizedPolynomial<P>( constant_zero<CoeffType>::get() );
565  FactorizedPolynomial<P> res( constant_one<CoeffType>::get() );
566  FactorizedPolynomial<P> mult( *this );
567  while( _exp > 0 )
568  {
569  if( (_exp & 1) != 0 )
570  res *= mult;
571  _exp /= 2;
572  if( _exp > 0 )
573  mult *= mult;
574  }
575  return res;
576  }
577 
578  template<typename P>
579  bool FactorizedPolynomial<P>::sqrt( FactorizedPolynomial<P>& _result ) const
580  {
581  CoeffType resultCoeff;
582  if( !carl::sqrt_exact( mCoefficient, resultCoeff ) )
583  return false;
584  if( !existsFactorization( *this ) )
585  {
586  _result = FactorizedPolynomial<P>( resultCoeff );
587  return true;
588  }
589  if( factorizedTrivially() )
590  {
591  PolyType sqrtOfPolynomial;
592  if( polynomial().sqrt( sqrtOfPolynomial ) )
593  {
594  _result = FactorizedPolynomial<P>( std::move( sqrtOfPolynomial ), mpCache );
595  content().setNewFactors( _result, 1, _result, 1 );
596  rehash();
597  _result *= resultCoeff;
598  assert( computePolynomial( _result ).pow( 2 ) == computePolynomial( *this ) );
599  return true;
600  }
601  return false;
602  }
603  Factorization<P> resultFactorization;
604  for( const auto& factor : content().factorization() )
605  {
606  if( factor.second % 2 == 0 )
607  {
608  resultFactorization.insert( std::pair<FactorizedPolynomial<P>, carl::exponent>( factor.first, factor.second/2 ) );
609  }
610  else
611  {
612  FactorizedPolynomial<P> sqrtOfFactor;
613  if( !factor.first.sqrt( sqrtOfFactor ) )
614  return false;
615  resultFactorization.insert( std::pair<FactorizedPolynomial<P>, carl::exponent>( sqrtOfFactor, factor.second ) );
616  }
617  }
618  _result = FactorizedPolynomial<P>( std::move( resultFactorization ), resultCoeff, mpCache );
619  assert( computePolynomial( _result ).pow( 2 ) == computePolynomial( *this ) );
620  return true;
621  }
622 
623  template<typename P>
624  DivisionResult<FactorizedPolynomial<P>> FactorizedPolynomial<P>::divideBy( const FactorizedPolynomial<P>& _divisor ) const
625  {
626  static_assert(is_field_type<CoeffType>::value, "Division only defined for field coefficients");
627  if( _divisor.is_one() || this->is_zero() )
628  return DivisionResult<FactorizedPolynomial<P>>( *this, FactorizedPolynomial<P>() );
629  if( this->is_constant() && _divisor.is_constant() )
630  return DivisionResult<FactorizedPolynomial<P>>( FactorizedPolynomial<P>( this->coefficient()/_divisor.coefficient() ), FactorizedPolynomial<P>() );
631  else if( _divisor.is_constant() )
632  return DivisionResult<FactorizedPolynomial<P>>( (*this)/_divisor.coefficient(), FactorizedPolynomial<P>() );
633  else if( this->is_constant() )
634  return DivisionResult<FactorizedPolynomial<P>>( FactorizedPolynomial<P>(), *this );
635  std::pair<FactorizedPolynomial<P>,FactorizedPolynomial<P>> ret = lazyDiv( *this, _divisor );
636  // TODO: maybe apply gcd, which might make this operation more expensive but a long-term investment
637  DivisionResult<P> dr = ret.first.polynomial().divideBy( ret.second.polynomial() );
638  // TODO: Could we calculate the quotient and remainder directly on the factorizations?
639  FactorizedPolynomial<P> q = dr.quotient.is_constant() ? FactorizedPolynomial<P>( dr.quotient.constant_part() ) : FactorizedPolynomial<P>( dr.quotient, mpCache );
640  FactorizedPolynomial<P> r = dr.remainder.is_constant() ? FactorizedPolynomial<P>( dr.remainder.constant_part() ) : FactorizedPolynomial<P>( dr.remainder, mpCache );
641  return DivisionResult<FactorizedPolynomial<P>>( q, r );
642  }
643 
644  template<typename P>
645  template<typename C, EnableIf<is_field_type<C>>>
646  FactorizedPolynomial<P> FactorizedPolynomial<P>::divideBy( const CoeffType& _divisor ) const
647  {
648  FactorizedPolynomial<P> result( *this );
649  result /= _divisor;
650  return result;
651  }
652 
653  template<typename P>
654  template<typename C, EnableIf<is_field_type<C>>>
655  bool FactorizedPolynomial<P>::divideBy( const FactorizedPolynomial<P>& _divisor, FactorizedPolynomial<P>& _quotient ) const
656  {
657  static_assert(is_field_type<C>::value, "Division only defined for field coefficients");
658  DivisionResult<FactorizedPolynomial<P>> dr = this->divideBy( _divisor );
659  if( dr.remainder.is_zero() )
660  {
661  _quotient = dr.quotient;
662  return true;
663  }
664  return false;
665  }
666 
667  template<typename P>
668  bool operator==( const FactorizedPolynomial<P>& _lhs, const FactorizedPolynomial<P>& _fpolyB )
669  {
670  if( _lhs.pCache() == nullptr && _fpolyB.pCache() == nullptr )
671  {
672  return _lhs.coefficient() == _fpolyB.coefficient();
673  }
674  else if( _lhs.pCache() != nullptr && _fpolyB.pCache() != nullptr )
675  {
676  if ( _lhs.coefficient() == _fpolyB.coefficient() )
677  return _lhs.content() == _fpolyB.content();
678  }
679  return false;
680  }
681 
682  template <typename P>
683  bool operator==( const FactorizedPolynomial<P>& _lhs, const typename FactorizedPolynomial<P>::CoeffType& _rhs )
684  {
685  return !existsFactorization( _lhs ) && _lhs.coefficient() == _rhs;
686  }
687 
688  template<typename P>
689  bool operator<( const FactorizedPolynomial<P>& _lhs, const typename P::CoeffType& _rhs )
690  {
691  return !existsFactorization( _lhs ) && _lhs.coefficient() < _rhs;
692  }
693 
694  template<typename P>
695  bool operator<( const FactorizedPolynomial<P>& _lhs, const FactorizedPolynomial<P>& _rhs )
696  {
697  ASSERT_CACHE_EQUAL( _lhs.pCache(), _rhs.pCache() );
698 
699  if( _lhs.pCache() == nullptr && _rhs.pCache() == nullptr )
700  {
701  return _lhs.coefficient() < _rhs.coefficient();
702  }
703  else if( _lhs.pCache() != nullptr && _rhs.pCache() != nullptr )
704  {
705  if( _lhs.content() == _rhs.content() )
706  return _lhs.coefficient() < _rhs.coefficient();
707  return _lhs.content() < _rhs.content();
708  }
709  return _lhs.pCache() == nullptr;
710  }
711 
712  template<typename P>
713  bool operator!=( const FactorizedPolynomial<P>& _lhs, const FactorizedPolynomial<P>& _rhs )
714  {
715  ASSERT_CACHE_EQUAL( _lhs.pCache(), _rhs.pCache() );
716  if( _lhs.pCache() == nullptr && _rhs.pCache() == nullptr )
717  {
718  return _lhs.coefficient() != _rhs.coefficient();
719  }
720  else if( _lhs.pCache() != nullptr && _rhs.pCache() != nullptr )
721  {
722  return !(_lhs.content() == _rhs.content());
723  }
724  return true;
725  }
726 
727  template<typename P>
728  P computePolynomial( const FactorizedPolynomial<P>& _fpoly )
729  {
730  if( _fpoly.pCache() == nullptr )
731  return std::move( P( _fpoly.coefficient() ) );
732  P result = computePolynomial( _fpoly.content() );
733  result *= _fpoly.coefficient();
734  return result;
735  }
736 
737  template<typename P>
738  Coeff<P> distributeCoefficients( Factorization<P>& _factorization )
739  {
740  Coeff<P> result(1);
741  for ( auto factor = _factorization.begin(); factor != _factorization.end(); factor++ )
742  {
743  result *= carl::pow( factor->first.coefficient(), factor->second );
744  factor->first.mCoefficient = 1;
745  }
746  return result;
747  }
748 
749  template<typename P>
750  FactorizedPolynomial<P> FactorizedPolynomial<P>::operator-() const
751  {
752  FactorizedPolynomial<P> result( *this );
753  result.mCoefficient *= Coeff<P>( -1 );
754  assert( (-computePolynomial( *this )) == computePolynomial( result ) );
755  return result;
756  }
757 
758  template<typename P>
759  FactorizedPolynomial<P> operator+( const FactorizedPolynomial<P>& _lhs, const FactorizedPolynomial<P>& _rhs )
760  {
761  ASSERT_CACHE_EQUAL( _lhs.pCache(), _rhs.pCache() );
762  _lhs.strengthenActivity();
763  _rhs.strengthenActivity();
764 
765  // Handle cases where one or both are constant
766  if( !existsFactorization( _lhs ) )
767  {
768  if( existsFactorization( _rhs ) )
769  {
770  FactorizedPolynomial<P> result( _rhs.polynomial() * _rhs.coefficient() + _lhs.coefficient(), _rhs.pCache() );
771  assert( computePolynomial( _lhs ) + computePolynomial( _rhs ) == computePolynomial( result ) );
772  return result;
773  }
774  else
775  {
776  FactorizedPolynomial<P> result( _lhs.coefficient() + _rhs.coefficient() );
777  assert( computePolynomial( _lhs ) + computePolynomial( _rhs ) == computePolynomial( result ) );
778  return result;
779  }
780  }
781  else if( !existsFactorization( _rhs ) )
782  {
783  FactorizedPolynomial<P> result( _lhs.polynomial() * _lhs.coefficient() + _rhs.coefficient(), _lhs.pCache() );
784  assert( computePolynomial( _lhs ) + computePolynomial( _rhs ) == computePolynomial( result ) );
785  return result;
786  }
787 
788  Coeff<P> coefficientCommon = Coeff<P>(carl::gcd( carl::get_num( _lhs.coefficient() ), carl::get_num( _rhs.coefficient() ) )) /
789  Coeff<P>(carl::lcm( carl::get_denom( _lhs.coefficient() ), carl::get_denom( _rhs.coefficient() ) ));
790  Coeff<P> coefficientRestA = _lhs.coefficient() / coefficientCommon;
791  Coeff<P> coefficientRestB = _rhs.coefficient() / coefficientCommon;
792 
793  if (carl::is_zero(coefficientCommon))
794  {
795  FactorizedPolynomial<P> result;
796  assert( computePolynomial( _lhs ) + computePolynomial( _rhs ) == computePolynomial( result ) );
797  return result;
798  }
799 
800  Factorization<P> factorizationRestA, factorizationRestB;
801  assert( existsFactorization( _lhs ) );
802  const Factorization<P>& factorizationA = _lhs.factorization();
803  assert( existsFactorization( _rhs ) );
804  const Factorization<P>& factorizationB = _rhs.factorization();
805 
806  //Compute common divisor as factor of result
807  Factorization<P> resultFactorization = commonDivisor( factorizationA, factorizationB, factorizationRestA, factorizationRestB );
808 
809  //Compute remaining sum
810  P sum;
811  if( resultFactorization.empty() )
812  {
813  sum = _lhs.polynomial() * coefficientRestA;
814  sum += _rhs.polynomial() * coefficientRestB;
815  }
816  else
817  {
818  sum = computePolynomial( factorizationRestA ) * coefficientRestA;
819  sum += computePolynomial( factorizationRestB ) * coefficientRestB;
820  }
821  if ( carl::is_zero(sum) )
822  {
823  FactorizedPolynomial<P> result;
824  assert( computePolynomial( _lhs ) + computePolynomial( _rhs ) == computePolynomial( result ) );
825  return result;
826  }
827  else
828  {
829  if ( sum.is_constant() )
830  coefficientCommon *= sum.constant_part();
831  else
832  {
833  FactorizedPolynomial<P> fpolySum( sum, _lhs.pCache() );
834  coefficientCommon *= fpolySum.coefficient();
835  fpolySum.mCoefficient = Coeff<P>(1);
836  resultFactorization.insert( std::pair<FactorizedPolynomial<P>, carl::exponent>( fpolySum, 1 ) );
837  }
838  }
839  FactorizedPolynomial<P> result( std::move( resultFactorization ), coefficientCommon, FactorizedPolynomial<P>::chooseCache( _lhs.pCache(), _rhs.pCache() ) );
840  assert( computePolynomial( _lhs ) + computePolynomial( _rhs ) == computePolynomial( result ) );
841  return result;
842  }
843 
844  template<typename P>
845  FactorizedPolynomial<P> operator+( const FactorizedPolynomial<P>& _fpoly, const typename FactorizedPolynomial<P>::CoeffType& _coef )
846  {
847  FactorizedPolynomial<P> result( _fpoly );
848  return std::move( result += _coef );
849  }
850 
851  template<typename P>
852  FactorizedPolynomial<P>& FactorizedPolynomial<P>::operator+=( const CoeffType& _coef )
853  {
854  FactorizedPolynomial<P> result = *this + FactorizedPolynomial<P>( _coef );
855  return *this = result;
856  }
857 
858  template<typename P>
859  FactorizedPolynomial<P>& FactorizedPolynomial<P>::operator+=( const FactorizedPolynomial<P>& _fpoly )
860  {
861  FactorizedPolynomial<P> result = *this + _fpoly;
862  return *this = result;
863  }
864 
865  template<typename P>
866  FactorizedPolynomial<P> operator-( const FactorizedPolynomial<P>& _fpoly, const typename FactorizedPolynomial<P>::CoeffType& _coef )
867  {
868  FactorizedPolynomial<P> result( _fpoly );
869  return std::move( result -= _coef );
870  }
871 
872  template<typename P>
873  FactorizedPolynomial<P>& FactorizedPolynomial<P>::operator-=( const CoeffType& _coef )
874  {
875  FactorizedPolynomial<P> result = *this + FactorizedPolynomial<P>( CoeffType( -1 ) * _coef );
876  assert( computePolynomial( *this ) - _coef == computePolynomial( result ) );
877  return *this = result;
878  }
879 
880  template<typename P>
881  FactorizedPolynomial<P>& FactorizedPolynomial<P>::operator-=( const FactorizedPolynomial<P>& _fpoly )
882  {
883  FactorizedPolynomial<P> result = *this - _fpoly;
884  return *this = result;
885  }
886 
887  template<typename P>
888  FactorizedPolynomial<P> operator-( const FactorizedPolynomial<P>& _lhs, const FactorizedPolynomial<P>& _rhs )
889  {
890  ASSERT_CACHE_EQUAL( _lhs.pCache(), _rhs.pCache() );
891  const Coeff<P>& coefficient = -_rhs.coefficient();
892  if ( existsFactorization( _rhs ) )
893  {
894  FactorizedPolynomial<P> result = _lhs + FactorizedPolynomial<P>( std::move( Factorization<P>( _rhs.factorization() ) ), coefficient, _rhs.pCache() );
895  assert( computePolynomial( _lhs ) - computePolynomial( _rhs ) == computePolynomial( result ) );
896  return result;
897  }
898  else
899  {
900  FactorizedPolynomial<P> result = _lhs + FactorizedPolynomial<P>( coefficient );
901  assert( computePolynomial( _lhs ) - computePolynomial( _rhs ) == computePolynomial( result ) );
902  return result;
903  }
904  }
905 
906  template<typename P>
907  FactorizedPolynomial<P> operator*( const FactorizedPolynomial<P>& _fpoly, const typename FactorizedPolynomial<P>::CoeffType& _coeff )
908  {
909  if( carl::is_zero( _coeff ) )
910  {
911  FactorizedPolynomial<P> result;
912  assert( computePolynomial( _fpoly ) * _coeff == computePolynomial( result ) );
913  return result;
914  }
915  FactorizedPolynomial<P> result( _fpoly );
916  result.mCoefficient *= _coeff;
917  assert( computePolynomial( _fpoly ) * _coeff == computePolynomial( result ) );
918  return result;
919  }
920 
921  template<typename P>
922  FactorizedPolynomial<P> operator*( const FactorizedPolynomial<P>& _lhs, const FactorizedPolynomial<P>& _rhs )
923  {
924  ASSERT_CACHE_EQUAL( _lhs.pCache(), _rhs.pCache() );
925  _lhs.strengthenActivity();
926  _rhs.strengthenActivity();
927 
928  if( _lhs.is_zero() || _rhs.is_zero() )
929  {
930  FactorizedPolynomial<P> result;
931  assert( computePolynomial( _lhs ) * computePolynomial( _rhs ) == computePolynomial( result ) );
932  return result;
933  }
934 
935  if( !existsFactorization( _rhs ) )
936  {
937  FactorizedPolynomial<P> result( _lhs );
938  result.mCoefficient *= _rhs.mCoefficient;
939  assert( computePolynomial( _lhs ) * computePolynomial( _rhs ) == computePolynomial( result ) );
940  return result;
941  }
942  else if( !existsFactorization( _lhs ) )
943  {
944  FactorizedPolynomial<P> result( _rhs );
945  result.mCoefficient *= _lhs.mCoefficient;
946  assert( computePolynomial( _lhs ) * computePolynomial( _rhs ) == computePolynomial( result ) );
947  return result;
948  }
949 
950  Factorization<P> resultFactorization;
951  const Factorization<P>& factorizationA = _lhs.factorization();
952  const Factorization<P>& factorizationB = _rhs.factorization();
953  auto factorA = factorizationA.begin();
954  auto factorB = factorizationB.begin();
955  while( factorA != factorizationA.end() && factorB != factorizationB.end() )
956  {
957  if( factorA->first == factorB->first )
958  {
959  resultFactorization.insert( resultFactorization.end(), std::pair<FactorizedPolynomial<P>, carl::exponent>(factorA->first, factorA->second + factorB->second ) );
960  factorA++;
961  factorB++;
962  }
963  else if( factorA->first < factorB->first )
964  {
965  resultFactorization.insert( resultFactorization.end(), *factorA );
966  factorA++;
967  }
968  else
969  {
970  resultFactorization.insert( resultFactorization.end(), *factorB );
971  factorB++;
972  }
973  }
974  while ( factorA != factorizationA.end() )
975  {
976  resultFactorization.insert( resultFactorization.end(), *factorA );
977  factorA++;
978  }
979  while ( factorB != factorizationB.end() )
980  {
981  resultFactorization.insert( resultFactorization.end(), *factorB );
982  factorB++;
983  }
984 
985  Coeff<P> coefficientResult = _lhs.coefficient() * _rhs.coefficient();
986  //TODO needed?
987  coefficientResult *= distributeCoefficients( resultFactorization );
988  FactorizedPolynomial<P> result( std::move( resultFactorization ), coefficientResult, FactorizedPolynomial<P>::chooseCache( _lhs.pCache(), _rhs.pCache() ) );
989  assert( computePolynomial( _lhs ) * computePolynomial( _rhs ) == computePolynomial( result ) );
990  return result;
991  }
992 
993  template<typename P>
994  FactorizedPolynomial<P>& FactorizedPolynomial<P>::operator*=( const CoeffType& _coef )
995  {
996  if( carl::is_zero( _coef ) && mpCache != nullptr )
997  {
998  mpCache->dereg( mCacheRef );
999  mCacheRef = CACHE::NO_REF;
1000  mpCache = nullptr;
1001  }
1002  mCoefficient *= _coef;
1003  ASSERT_CACHE_REF_LEGAL( (*this) );
1004  return *this;
1005  }
1006 
1007  template<typename P>
1008  FactorizedPolynomial<P>& FactorizedPolynomial<P>::operator*=( const FactorizedPolynomial<P>& _fpoly )
1009  {
1010  FactorizedPolynomial<P> result = *this * _fpoly;
1011  return *this = result;
1012  }
1013 
1014  template<typename P>
1015  FactorizedPolynomial<P>& FactorizedPolynomial<P>::operator/=( const CoeffType& _coef )
1016  {
1017  assert( !carl::is_zero( _coef ) );
1018  FactorizedPolynomial<P> tmp = *this;
1019  this->mCoefficient /= _coef;
1020  assert(computePolynomial(tmp) * (static_cast<CoeffType>(1)/_coef) == computePolynomial(*this));
1021  return *this;
1022  }
1023 
1024  template<typename P>
1025  FactorizedPolynomial<P>& FactorizedPolynomial<P>::operator/=( const FactorizedPolynomial<P>& _fpoly )
1026  {
1027  assert( !_fpoly.is_zero() );
1028  FactorizedPolynomial<P> result = this->quotient( _fpoly );
1029  assert( computePolynomial( *this ) / computePolynomial( _fpoly ) == computePolynomial( result ) );
1030  return *this = result;
1031  }
1032 
1033  template<typename P>
1034  FactorizedPolynomial<P> FactorizedPolynomial<P>::quotient(const FactorizedPolynomial<P>& _fdivisor) const
1035  {
1036  assert( !_fdivisor.is_zero() );
1037  if( is_zero() )
1038  {
1039  FactorizedPolynomial<P> result;
1040  assert( carl::quotient(computePolynomial( *this ), computePolynomial( _fdivisor ) ) == computePolynomial( result ) );
1041  return result;
1042  }
1043  FactorizedPolynomial<P> result = lazyDiv( *this, _fdivisor ).first;
1044  assert( carl::quotient(computePolynomial( *this ), computePolynomial( _fdivisor ) ) == computePolynomial( result ) );
1045  return result;
1046  }
1047 
1048  template<typename P>
1049  std::string FactorizedPolynomial<P>::toString( bool _infix, bool /* _friendlyVarNames */ ) const
1050  {
1051  if( existsFactorization( *this ) )
1052  {
1053  std::stringstream result;
1054  if( _infix )
1055  {
1056  if( mCoefficient != Coeff<P>( 1 ) )
1057  {
1058  result << mCoefficient << " * (";
1059  }
1060  result << content();
1061  if( mCoefficient != Coeff<P>( 1 ) )
1062  result << ")";
1063  }
1064  else
1065  {
1066  bool withCoeff = mCoefficient != Coeff<P>( 1 );
1067  if( withCoeff )
1068  {
1069  result << "(* " << mCoefficient << " ";
1070  }
1071  result << content();
1072  if( withCoeff )
1073  result << ")";
1074  }
1075  return result.str();
1076  }
1077  std::stringstream s;
1078  s << mCoefficient;
1079  return s.str();
1080  }
1081 
1082  template<typename P>
1083  FactorizedPolynomial<P> quotient( const FactorizedPolynomial<P>& _fpolyA, const FactorizedPolynomial<P>& _fpolyB )
1084  {
1085  assert( !_fpolyB.is_zero() );
1086  return _fpolyA.quotient( _fpolyB );
1087  }
1088 
1089  template<typename P>
1090  std::pair<FactorizedPolynomial<P>,FactorizedPolynomial<P>> lazyDiv( const FactorizedPolynomial<P>& _fpolyA, const FactorizedPolynomial<P>& _fpolyB )
1091  {
1092  assert( !_fpolyB.is_zero() );
1093  if( _fpolyA.is_zero() )
1094  return std::make_pair( FactorizedPolynomial<P>(), FactorizedPolynomial<P>( Coeff<P>( 1 ) ) );
1095  ASSERT_CACHE_EQUAL( _fpolyA.pCache(), _fpolyB.pCache() );
1096  _fpolyA.strengthenActivity();
1097  _fpolyB.strengthenActivity();
1098  // Handle cases where one or both are constant
1099  if( !existsFactorization( _fpolyB ) )
1100  {
1101  FactorizedPolynomial<P> fPolyASimplified( _fpolyA );
1102  assert( !carl::is_zero(_fpolyB.mCoefficient) );
1103  fPolyASimplified.mCoefficient /= _fpolyB.mCoefficient;
1104  auto result = std::make_pair( fPolyASimplified, FactorizedPolynomial<P>( Coeff<P>( 1 ) ) );
1105  assert( computePolynomial( result.first ) * computePolynomial( _fpolyB ) == computePolynomial( result.second ) * computePolynomial( _fpolyA ) );
1106  return result;
1107  }
1108  else if( !existsFactorization( _fpolyA ) )
1109  {
1110  FactorizedPolynomial<P> fPolyASimplified( _fpolyA );
1111  assert( !carl::is_zero(_fpolyB.mCoefficient) );
1112  fPolyASimplified.mCoefficient /= _fpolyB.mCoefficient;
1113  FactorizedPolynomial<P> fPolyBSimplified( _fpolyB );
1114  assert( !carl::is_zero(_fpolyA.mCoefficient) );
1115  fPolyBSimplified.mCoefficient = Coeff<P>( 1 );
1116  auto result = std::make_pair( fPolyASimplified, fPolyBSimplified );
1117  assert( computePolynomial( result.first ) * computePolynomial( _fpolyB ) == computePolynomial( result.second ) * computePolynomial( _fpolyA ) );
1118  return result;
1119  }
1120  Factorization<P> resultFactorizationA;
1121  Factorization<P> resultFactorizationB;
1122  const Factorization<P>& factorizationA = _fpolyA.factorization();
1123  const Factorization<P>& factorizationB = _fpolyB.factorization();
1124  auto factorA = factorizationA.begin();
1125  auto factorB = factorizationB.begin();
1126  while( factorA != factorizationA.end() && factorB != factorizationB.end() )
1127  {
1128  if( factorA->first == factorB->first )
1129  {
1130  if( factorA->second > factorB->second )
1131  resultFactorizationA.insert( resultFactorizationA.end(), std::pair<FactorizedPolynomial<P>, carl::exponent>(factorA->first, factorA->second - factorB->second ) );
1132  else if( factorA->second < factorB->second )
1133  resultFactorizationB.insert( resultFactorizationB.end(), std::pair<FactorizedPolynomial<P>, carl::exponent>(factorA->first, factorB->second - factorA->second ) );
1134  factorA++;
1135  factorB++;
1136  }
1137  else if( factorA->first < factorB->first )
1138  {
1139  resultFactorizationA.insert( resultFactorizationA.end(), *factorA );
1140  factorA++;
1141  }
1142  else
1143  {
1144  resultFactorizationB.insert( resultFactorizationB.end(), *factorB );
1145  factorB++;
1146  }
1147  }
1148  while ( factorA != factorizationA.end() )
1149  {
1150  resultFactorizationA.insert( resultFactorizationA.end(), *factorA );
1151  factorA++;
1152  }
1153  while ( factorB != factorizationB.end() )
1154  {
1155  resultFactorizationB.insert( resultFactorizationB.end(), *factorB );
1156  factorB++;
1157  }
1158  Coeff<P> coefficientResultA = _fpolyA.coefficient() / _fpolyB.coefficient();
1159  auto cache = FactorizedPolynomial<P>::chooseCache( _fpolyA.pCache(), _fpolyB.pCache() );
1160  FactorizedPolynomial<P> resultA( std::move( resultFactorizationA ), coefficientResultA, cache );
1161  FactorizedPolynomial<P> resultB( std::move( resultFactorizationB ), Coeff<P>( 1 ), cache );
1162  assert( computePolynomial( resultA ) * computePolynomial( _fpolyB ) == computePolynomial( resultB ) * computePolynomial( _fpolyA ) );
1163  return std::make_pair( resultA, resultB );
1164  }
1165 
1166  template<typename P>
1167  FactorizedPolynomial<P> lcm( const FactorizedPolynomial<P>& _fpolyA, const FactorizedPolynomial<P>& _fpolyB )
1168  {
1169  assert( !_fpolyA.is_zero() && !_fpolyB.is_zero() );
1170  ASSERT_CACHE_EQUAL( _fpolyA.pCache(), _fpolyB.pCache() );
1171  _fpolyA.strengthenActivity();
1172  _fpolyB.strengthenActivity();
1173  bool rehashFPolyA = false;
1174  bool rehashFPolyB = false;
1175  Coeff<P> coefficientLCM = Coeff<P>(carl::lcm( carl::get_num(_fpolyA.coefficient()), carl::get_num(_fpolyB.coefficient()) )) /
1176  Coeff<P>(carl::gcd( carl::get_denom(_fpolyA.coefficient()), carl::get_denom(_fpolyB.coefficient()) ));
1177  // Handle cases where one or both are constant
1178  if( !existsFactorization( _fpolyB ) )
1179  {
1180  FactorizedPolynomial<P> result( _fpolyA );
1181  result.mCoefficient = coefficientLCM;
1182  assert( carl::is_zero(carl::remainder(computePolynomial( result ), computePolynomial( _fpolyA ) )) );
1183  assert( carl::is_zero(carl::remainder(computePolynomial( result ), computePolynomial( _fpolyB ) )) );
1184  return result;
1185  }
1186  else if( !existsFactorization( _fpolyA ) )
1187  {
1188  FactorizedPolynomial<P> result( _fpolyB );
1189  result.mCoefficient = coefficientLCM;
1190  assert( carl::is_zero(carl::remainder(computePolynomial( result ), computePolynomial( _fpolyA ) )) );
1191  assert( carl::is_zero(carl::remainder(computePolynomial( result ), computePolynomial( _fpolyB ) )) );
1192  return result;
1193  }
1194  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", "Compute LCM of " << _fpolyA << " and " << _fpolyB );
1195 
1196  //Both polynomials are not constant
1197  Factorization<P> restAFactorization, restBFactorization;
1198  Coeff<P> c( 0 );
1199  gcd( _fpolyA.content(), _fpolyB.content(), restAFactorization, restBFactorization, c, rehashFPolyA, rehashFPolyB );
1200  if( c != Coeff<P>( 0 ) )
1201  coefficientLCM *= c;
1202 
1203  if( rehashFPolyA )
1204  _fpolyA.rehash();
1205  if( rehashFPolyB )
1206  _fpolyB.rehash();
1207 
1208  //Compute lcm as A*restB
1209  Factorization<P> lcmFactorization;
1210  if( existsFactorization( _fpolyA ) )
1211  lcmFactorization.insert( _fpolyA.factorization().begin(), _fpolyA.factorization().end() );
1212  for ( auto factor = restBFactorization.begin(); factor != restBFactorization.end(); factor++ )
1213  {
1214  lcmFactorization.insert( *factor );
1215  }
1216 
1217  coefficientLCM *= distributeCoefficients( lcmFactorization );
1218  FactorizedPolynomial<P> result( std::move( lcmFactorization ), coefficientLCM, _fpolyA.pCache() );
1219  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", "LCM of " << _fpolyA << " and " << _fpolyB << ": " << result);
1220  assert( carl::is_zero(carl::remainder(computePolynomial( result ), computePolynomial( _fpolyA ) )) );
1221  assert( carl::is_zero(carl::remainder(computePolynomial( result ), computePolynomial( _fpolyB ) )) );
1222  return result;
1223  }
1224 
1225  template<typename P>
1226  FactorizedPolynomial<P> commonMultiple( const FactorizedPolynomial<P>& _fpolyA, const FactorizedPolynomial<P>& _fpolyB )
1227  {
1228  assert( !_fpolyA.is_zero() && !_fpolyB.is_zero() );
1229  ASSERT_CACHE_EQUAL( _fpolyA.pCache(), _fpolyB.pCache() );
1230  _fpolyA.strengthenActivity();
1231  _fpolyB.strengthenActivity();
1232 
1233  Coeff<P> coefficientLCM = Coeff<P>(carl::lcm( carl::get_num(_fpolyA.coefficient()), carl::get_num(_fpolyB.coefficient()) )) /
1234  Coeff<P>(carl::gcd( carl::get_denom(_fpolyA.coefficient()), carl::get_denom(_fpolyB.coefficient()) ));
1235  // Handle cases where one or both are constant
1236  if( !existsFactorization( _fpolyB ) )
1237  {
1238  FactorizedPolynomial<P> result( _fpolyA );
1239  result.mCoefficient = coefficientLCM;
1240  assert( carl::is_zero(carl::remainder(computePolynomial( result ), computePolynomial( _fpolyA ) )) );
1241  assert( carl::is_zero(carl::remainder(computePolynomial( result ), computePolynomial( _fpolyB ) )) );
1242  return result;
1243  }
1244  else if( !existsFactorization( _fpolyA ) )
1245  {
1246  FactorizedPolynomial<P> result( _fpolyB );
1247  result.mCoefficient = coefficientLCM;
1248  assert( carl::is_zero(carl::remainder(computePolynomial( result ), computePolynomial( _fpolyA ) )) );
1249  assert( carl::is_zero(carl::remainder(computePolynomial( result ), computePolynomial( _fpolyB ) )) );
1250  return result;
1251  }
1252 
1253  Factorization<P> cmFactorization;
1254  const Factorization<P>& factorizationA = _fpolyA.factorization();
1255  const Factorization<P>& factorizationB = _fpolyB.factorization();
1256  auto factorA = factorizationA.begin();
1257  auto factorB = factorizationB.begin();
1258  while( factorA != factorizationA.end() && factorB != factorizationB.end() )
1259  {
1260  if( factorA->first == factorB->first )
1261  {
1262  cmFactorization.insert( cmFactorization.end(), factorA->second > factorB->second ? *factorA : *factorB );
1263  factorA++;
1264  factorB++;
1265  }
1266  else if( factorA->first < factorB->first )
1267  {
1268  cmFactorization.insert( cmFactorization.end(), *factorA );
1269  factorA++;
1270  }
1271  else
1272  {
1273  cmFactorization.insert( cmFactorization.end(), *factorB );
1274  factorB++;
1275  }
1276  }
1277  while ( factorA != factorizationA.end() )
1278  {
1279  cmFactorization.insert( cmFactorization.end(), *factorA );
1280  factorA++;
1281  }
1282  while ( factorB != factorizationB.end() )
1283  {
1284  cmFactorization.insert( cmFactorization.end(), *factorB );
1285  factorB++;
1286  }
1287  FactorizedPolynomial<P> result( std::move( cmFactorization ), coefficientLCM, FactorizedPolynomial<P>::chooseCache( _fpolyA.pCache(), _fpolyB.pCache() ) );
1288  assert( carl::is_zero(carl::remainder(computePolynomial( result ), computePolynomial( _fpolyA ) )) );
1289  assert( carl::is_zero(carl::remainder(computePolynomial( result ), computePolynomial( _fpolyB ) )) );
1290  return result;
1291  }
1292 
1293  template<typename P>
1294  FactorizedPolynomial<P> commonDivisor( const FactorizedPolynomial<P>& _fpolyA, const FactorizedPolynomial<P>& _fpolyB )
1295  {
1296  assert( !_fpolyA.is_zero() && !_fpolyB.is_zero() );
1297  ASSERT_CACHE_EQUAL( _fpolyA.pCache(), _fpolyB.pCache() );
1298  _fpolyA.strengthenActivity();
1299  _fpolyB.strengthenActivity();
1300  Coeff<P> coefficientCommon = Coeff<P>(carl::gcd( carl::get_num( _fpolyA.coefficient() ), carl::get_num( _fpolyB.coefficient() ) )) /
1301  Coeff<P>(carl::lcm( carl::get_denom( _fpolyA.coefficient() ), carl::get_denom( _fpolyB.coefficient() ) ));
1302  // Handle cases where one or both are constant
1303  if( !existsFactorization( _fpolyA ) || !existsFactorization( _fpolyB ) )
1304  {
1305  FactorizedPolynomial<P> result( coefficientCommon );
1306  assert( carl::is_zero(carl::remainder(computePolynomial( _fpolyA ), computePolynomial( result ) )) );
1307  assert( carl::is_zero(carl::remainder(computePolynomial( _fpolyB ), computePolynomial( result ) )) );
1308  return result;
1309  }
1310 
1311  Factorization<P> cdFactorization;
1312  const Factorization<P>& factorizationA = _fpolyA.factorization();
1313  const Factorization<P>& factorizationB = _fpolyB.factorization();
1314  auto factorA = factorizationA.begin();
1315  auto factorB = factorizationB.begin();
1316  while( factorA != factorizationA.end() && factorB != factorizationB.end() )
1317  {
1318  if( factorA->first == factorB->first )
1319  {
1320  cdFactorization.insert( cdFactorization.end(), factorA->second < factorB->second ? *factorA : *factorB );
1321  factorA++;
1322  factorB++;
1323  }
1324  else if( factorA->first < factorB->first )
1325  factorA++;
1326  else
1327  factorB++;
1328  }
1329  FactorizedPolynomial<P> result( std::move( cdFactorization ), coefficientCommon, FactorizedPolynomial<P>::chooseCache( _fpolyA.pCache(), _fpolyB.pCache() ) );
1330  assert( carl::is_zero(carl::remainder(computePolynomial( _fpolyA ), computePolynomial( result ) )) );
1331  assert( carl::is_zero(carl::remainder(computePolynomial( _fpolyB ), computePolynomial( result ) )) );
1332  return result;
1333  }
1334 
1335  template<typename P>
1336  Factorization<P> commonDivisor( const Factorization<P>& _fFactorizationA, const Factorization<P>& _fFactorizationB, Factorization<P>& _fFactorizationRestA, Factorization<P>& _fFactorizationRestB)
1337  {
1338  assert( !_fFactorizationA.empty() && !_fFactorizationB.empty() );
1339  Factorization<P> resultFactorization;
1340  _fFactorizationRestA.clear();
1341  _fFactorizationRestB.clear();
1342  auto factorA = _fFactorizationA.begin();
1343  auto factorB = _fFactorizationB.begin();
1344  while( factorA != _fFactorizationA.end() && factorB != _fFactorizationB.end() )
1345  {
1346  if( factorA->first == factorB->first )
1347  {
1348  if( factorA->second < factorB->second )
1349  {
1350  resultFactorization.insert( resultFactorization.end(), *factorA );
1351  _fFactorizationRestB.insert( _fFactorizationRestB.end(), std::make_pair( factorB->first, factorB->second - factorA->second ) );
1352  }
1353  else if( factorA->second > factorB->second )
1354  {
1355  resultFactorization.insert( resultFactorization.end(), *factorB );
1356  _fFactorizationRestA.insert( _fFactorizationRestA.end(), std::make_pair( factorA->first, factorA->second - factorB->second ) );
1357  }
1358  else
1359  {
1360  resultFactorization.insert( resultFactorization.end(), *factorA );
1361  }
1362  factorA++;
1363  factorB++;
1364  }
1365  else if( factorA->first < factorB->first )
1366  {
1367  _fFactorizationRestA.insert( _fFactorizationRestA.end(), *factorA );
1368  factorA++;
1369  }
1370  else
1371  {
1372  _fFactorizationRestB.insert( _fFactorizationRestB.end(), *factorB );
1373  factorB++;
1374  }
1375  }
1376  while ( factorA != _fFactorizationA.end() )
1377  {
1378  _fFactorizationRestA.insert( _fFactorizationRestA.end(), *factorA );
1379  factorA++;
1380  }
1381  while ( factorB != _fFactorizationB.end() )
1382  {
1383  _fFactorizationRestB.insert( _fFactorizationRestB.end(), *factorB );
1384  factorB++;
1385  }
1386  assert( computePolynomial( _fFactorizationA ) == computePolynomial( resultFactorization ) * computePolynomial( _fFactorizationRestA ) );
1387  assert( computePolynomial( _fFactorizationB ) == computePolynomial( resultFactorization ) * computePolynomial( _fFactorizationRestB ) );
1388  return resultFactorization;
1389  }
1390 
1391  template<typename P>
1392  FactorizedPolynomial<P> gcd( const FactorizedPolynomial<P>& _fpolyA, const FactorizedPolynomial<P>& _fpolyB )
1393  {
1394  assert( !_fpolyA.is_zero() && !_fpolyB.is_zero() );
1395  FactorizedPolynomial<P> restA, restB;
1396  return gcd( _fpolyA, _fpolyB, restA, restB );
1397  }
1398 
1399  template<typename P>
1400  FactorizedPolynomial<P> gcd( const FactorizedPolynomial<P>& _fpolyA, const FactorizedPolynomial<P>& _fpolyB, FactorizedPolynomial<P>& _fpolyRestA, FactorizedPolynomial<P>& _fpolyRestB )
1401  {
1402  assert( !_fpolyA.is_zero() && !_fpolyB.is_zero() );
1403  _fpolyA.strengthenActivity();
1404  _fpolyB.strengthenActivity();
1405 
1406  Coeff<P> coefficientCommon = Coeff<P>(carl::gcd( carl::get_num( _fpolyA.coefficient() ), carl::get_num( _fpolyB.coefficient() ) )) /
1407  Coeff<P>(carl::lcm( carl::get_denom( _fpolyA.coefficient() ), carl::get_denom( _fpolyB.coefficient() ) ));
1408  Coeff<P> coefficientRestA = _fpolyA.coefficient() / coefficientCommon;
1409  Coeff<P> coefficientRestB = _fpolyB.coefficient() / coefficientCommon;
1410 
1411  //Handle cases where one or both are constant
1412  if ( !existsFactorization( _fpolyA ) )
1413  {
1414  if ( existsFactorization( _fpolyB ) )
1415  {
1416  _fpolyRestA = FactorizedPolynomial<P>( coefficientRestA );
1417  _fpolyRestB = FactorizedPolynomial<P>( std::move( Factorization<P>( _fpolyB.factorization() ) ), coefficientRestB, _fpolyB.pCache() );
1418  FactorizedPolynomial<P> result( coefficientCommon );
1419  assert( computePolynomial( _fpolyA ) == computePolynomial( result ) * computePolynomial( _fpolyRestA ) );
1420  assert( computePolynomial( _fpolyB ) == computePolynomial( result ) * computePolynomial( _fpolyRestB ) );
1421  return result;
1422  }
1423  else
1424  {
1425  _fpolyRestA = FactorizedPolynomial<P>( coefficientRestA );
1426  _fpolyRestB = FactorizedPolynomial<P>( coefficientRestB );
1427  FactorizedPolynomial<P> result( coefficientCommon );
1428  assert( computePolynomial( _fpolyA ) == computePolynomial( result ) * computePolynomial( _fpolyRestA ) );
1429  assert( computePolynomial( _fpolyB ) == computePolynomial( result ) * computePolynomial( _fpolyRestB ) );
1430  return result;
1431  }
1432  }
1433  else if ( !existsFactorization( _fpolyB ) )
1434  {
1435  _fpolyRestA = FactorizedPolynomial<P>( std::move( Factorization<P>( _fpolyA.factorization() ) ), coefficientRestA, _fpolyA.pCache());
1436  _fpolyRestB = FactorizedPolynomial<P>( coefficientRestB );
1437  FactorizedPolynomial<P> result( coefficientCommon );
1438  assert( computePolynomial( _fpolyA ) == computePolynomial( result ) * computePolynomial( _fpolyRestA ) );
1439  assert( computePolynomial( _fpolyB ) == computePolynomial( result ) * computePolynomial( _fpolyRestB ) );
1440  return result;
1441  }
1442 
1443  //Both polynomials are not constant
1444  Factorization<P> restAFactorization, restBFactorization;
1445  Coeff<P> c( 0 );
1446  bool rehashFPolyA = false;
1447  bool rehashFPolyB = false;
1448  Factorization<P> gcdFactorization = gcd( _fpolyA.content(), _fpolyB.content(), restAFactorization, restBFactorization, c, rehashFPolyA, rehashFPolyB );
1449 
1450  if( c != Coeff<P>( 0 ) )
1451  coefficientCommon *= c;
1452  if( rehashFPolyA )
1453  _fpolyA.rehash();
1454  if( rehashFPolyB )
1455  _fpolyB.rehash();
1456 
1457  coefficientRestA *= distributeCoefficients( restAFactorization );
1458  coefficientRestB *= distributeCoefficients( restBFactorization );
1459  coefficientCommon *= distributeCoefficients( gcdFactorization );
1460  _fpolyRestA = FactorizedPolynomial<P>( std::move( restAFactorization ), coefficientRestA, _fpolyA.pCache());
1461  _fpolyRestB = FactorizedPolynomial<P>( std::move( restBFactorization ), coefficientRestB, _fpolyA.pCache());
1462 
1463  FactorizedPolynomial<P> result( std::move( gcdFactorization ), coefficientCommon, _fpolyA.pCache() );
1464  assert( computePolynomial( _fpolyA ) == computePolynomial( result ) * computePolynomial( _fpolyRestA ) );
1465  assert( computePolynomial( _fpolyB ) == computePolynomial( result ) * computePolynomial( _fpolyRestB ) );
1466  return result;
1467  }
1468 
1469  template<typename P>
1470  Factors<FactorizedPolynomial<P>> factor( const FactorizedPolynomial<P>& _fpoly )
1471  {
1472  return factor( _fpoly.content(), _fpoly.coefficient() );
1473  }
1474 
1475  template <typename P>
1476  std::ostream& operator<<( std::ostream& _out, const FactorizedPolynomial<P>& _fpoly )
1477  {
1478  _out << _fpoly.toString();
1479  return _out;
1480  }
1481 } // namespace carl