carl  25.04
Computer ARithmetic Library
PolynomialFactorizationPair.tpp
Go to the documentation of this file.
1 /*
2  * File: PolynomialFactorizationPair.tpp
3  * Author: Florian Corzilius
4  *
5  * Created on September 5, 2014, 3:57 PM
6  */
7 
8 #pragma once
9 
10 #include "PolynomialFactorizationPair.h"
11 
12 #include "FactorizedPolynomial.h"
13 #include <carl-logging/carl-logging.h>
14 #include <carl-arith/poly/umvpoly/functions/Definiteness.h>
15 #include <carl-arith/poly/umvpoly/functions/GCD.h>
16 
17 namespace carl
18 {
19  template <typename P>
20  std::string factorizationToString( const Factorization<P>& _factorization, bool _infix, bool _friendlyVarNames )
21  {
22  if( _factorization.empty() )
23  {
24  return "1";
25  }
26  else
27  {
28  std::string result;
29  if( _infix )
30  {
31  for( auto polyExpPair = _factorization.begin(); polyExpPair != _factorization.end(); ++polyExpPair )
32  {
33  if( polyExpPair != _factorization.begin() )
34  result += " * ";
35  assert( polyExpPair->second > 0 );
36  result += "(" + polyExpPair->first.toString( true, _friendlyVarNames ) + ")";
37  if( polyExpPair->second > 1 )
38  {
39  std::stringstream s;
40  s << polyExpPair->second;
41  result += "^" + s.str();
42  }
43  }
44  }
45  else
46  {
47  if( _factorization.size() == 1 && _factorization.begin()->second == 1 )
48  {
49  return _factorization.begin()->first.toString( false, _friendlyVarNames );
50  }
51  result += "(*";
52  for( auto polyExpPair = _factorization.begin(); polyExpPair != _factorization.end(); ++polyExpPair )
53  {
54  assert( polyExpPair->second > 0 );
55  for( size_t i = 0; i < polyExpPair->second; ++i )
56  result += " " + polyExpPair->first.toString( false, _friendlyVarNames );
57  }
58  result += ")";
59  }
60  return result;
61  }
62  }
63 
64  template <typename P>
65  std::ostream& operator<<( std::ostream& _out, const Factorization<P>& _factorization )
66  {
67  _out << factorizationToString( _factorization );
68  return _out;
69  }
70 
71  template<typename P>
72  bool factorizationsEqual( const Factorization<P>& _factorizationA, const Factorization<P>& _factorizationB )
73  {
74  auto iterA = _factorizationA.begin();
75  auto iterB = _factorizationB.begin();
76  while( iterA != _factorizationA.end() && iterB != _factorizationB.end() )
77  {
78  if( iterA->second != iterB->second || !(iterA->first == iterB->first) )
79  break;
80  ++iterA; ++iterB;
81  }
82  return iterA == _factorizationA.end() && iterB == _factorizationB.end();
83  }
84 
85  template<typename P>
86  PolynomialFactorizationPair<P>::PolynomialFactorizationPair( Factorization<P>&& _factorization, P* _polynomial ):
87  mHash( 0 ),
88  mMutex(),
89  mFactorization( std::move( _factorization ) ),
90  mpPolynomial( _polynomial ),
91  mIrreducible( -1 )
92  {
93  if ( mpPolynomial == nullptr )
94  {
95  if ( mFactorization.size() == 1 )
96  {
97  // No factorization -> set polynomial
98  assert( existsFactorization( mFactorization.begin()->first ) );
99  mpPolynomial = new P( carl::pow(*mFactorization.begin()->first.content().mpPolynomial, mFactorization.begin()->second ) );
100  assert( mpPolynomial != nullptr );
101  }
102  }
103  else
104  {
105  if ( is_one(*mpPolynomial) )
106  {
107  assert( mFactorization.size() == 0);
108  mpPolynomial = nullptr;
109  }
110  }
111 
112  // Check correctness
113  assert( mpPolynomial == nullptr || carl::is_one(mpPolynomial->coprime_factor()) );
114  assert( mpPolynomial == nullptr || mFactorization.empty() || assertFactorization() );
115 
116  rehash();
117  }
118 
119  template<typename P>
120  PolynomialFactorizationPair<P>::~PolynomialFactorizationPair()
121  {
122  delete mpPolynomial;
123  }
124 
125  template<typename P>
126  void PolynomialFactorizationPair<P>::rehash() const
127  {
128  std::lock_guard<std::recursive_mutex> lock( mMutex );
129  if( mpPolynomial == nullptr )
130  {
131  assert( mFactorization.empty() || mFactorization.size() > 1 );
132  mHash = 0;
133  for( auto polyExpPair = mFactorization.begin(); polyExpPair != mFactorization.end(); ++polyExpPair )
134  {
135  mHash = (mHash << 5) | (mHash >> (sizeof(size_t)*8 - 5));
136  mHash ^= std::hash<FactorizedPolynomial<P>>()( polyExpPair->first );
137  mHash = (mHash << 5) | (mHash >> (sizeof(size_t)*8 - 5));
138  mHash ^= polyExpPair->second;
139  }
140  }
141  else
142  {
143  mHash = std::hash<P>()( *mpPolynomial );
144  }
145  }
146 
147  template<typename P>
148  bool operator==( const PolynomialFactorizationPair<P>& _polyFactA, const PolynomialFactorizationPair<P>& _polyFactB )
149  {
150  if( &_polyFactA == &_polyFactB )
151  return true;
152  //TODO fix
153  //if ( _polyFactA.mHash != _polyFactB.mHash )
154  // return false;
155  std::lock_guard<std::recursive_mutex> lockA( _polyFactA.mMutex );
156  std::lock_guard<std::recursive_mutex> lockB( _polyFactB.mMutex );
157  if( _polyFactA.mpPolynomial != nullptr && _polyFactB.mpPolynomial != nullptr )
158  {
159  return *_polyFactA.mpPolynomial == *_polyFactB.mpPolynomial;
160  }
161  else
162  {
163  if( factorizationsEqual( _polyFactA.factorization(), _polyFactB.factorization() ) )
164  return true;
165  else
166  {
167  // There is no way around this );
168  if( _polyFactA.mpPolynomial == nullptr )
169  _polyFactA.mpPolynomial = new P( computePolynomial( _polyFactA.factorization() ) );
170  if( _polyFactB.mpPolynomial == nullptr )
171  _polyFactB.mpPolynomial = new P( computePolynomial( _polyFactB.factorization() ) );
172  return *_polyFactA.mpPolynomial == *_polyFactB.mpPolynomial;
173  }
174  }
175  }
176 
177  template<typename P>
178  bool operator<( const PolynomialFactorizationPair<P>& _polyFactA, const PolynomialFactorizationPair<P>& _polyFactB )
179  {
180  if( &_polyFactA == &_polyFactB )
181  return false;
182  std::lock_guard<std::recursive_mutex> lockA( _polyFactA.mMutex );
183  std::lock_guard<std::recursive_mutex> lockB( _polyFactB.mMutex );
184  if( _polyFactA.mpPolynomial != nullptr && _polyFactB.mpPolynomial != nullptr )
185  {
186  return *_polyFactA.mpPolynomial < *_polyFactB.mpPolynomial;
187  }
188  else
189  {
190  auto iterA = _polyFactA.factorization().begin();
191  auto iterB = _polyFactB.factorization().begin();
192  while( iterA != _polyFactA.factorization().end() && iterB != _polyFactB.factorization().end() )
193  {
194  if( iterA->first < iterB->first )
195  {
196  return true;
197  }
198  else if( iterA->first == iterB->first )
199  {
200  if( iterA->second < iterB->second )
201  {
202  return true;
203  }
204  else if( iterA->second > iterB->second )
205  {
206  return false;
207  }
208  }
209  else
210  {
211  return false;
212  }
213  ++iterA; ++iterB;
214  }
215  return iterA == _polyFactA.factorization().end();
216  }
217  }
218 
219  template<typename P>
220  bool canBeUpdated( const PolynomialFactorizationPair<P>& _toUpdate, const PolynomialFactorizationPair<P>& _updateWith )
221  {
222  if( &_toUpdate == &_updateWith )
223  return false;
224  std::lock_guard<std::recursive_mutex> lockA( _toUpdate.mMutex );
225  std::lock_guard<std::recursive_mutex> lockB( _updateWith.mMutex );
226  assert( _toUpdate.hash() == _updateWith.hash() && _toUpdate == _updateWith );
227  if( _toUpdate.mpPolynomial == nullptr && _updateWith.mpPolynomial != nullptr )
228  return true;
229  assert( _updateWith.mpPolynomial == nullptr || (*_toUpdate.mpPolynomial) == (*_updateWith.mpPolynomial) );
230  return !_updateWith.factorization().empty() && !factorizationsEqual( _toUpdate.factorization(), _updateWith.factorization() );
231  }
232 
233  template<typename P>
234  void update( PolynomialFactorizationPair<P>& _toUpdate, PolynomialFactorizationPair<P>& _updateWith )
235  {
236  assert( canBeUpdated( _toUpdate, _updateWith ) ); // This assertion only ensures efficient use this method.
237  assert( &_toUpdate != &_updateWith );
238  assert( _toUpdate.mpPolynomial == nullptr || _updateWith.mpPolynomial == nullptr || *_toUpdate.mpPolynomial == *_updateWith.mpPolynomial );
239  std::lock_guard<std::recursive_mutex> lockA( _toUpdate.mMutex );
240  std::lock_guard<std::recursive_mutex> lockB( _updateWith.mMutex );
241  if( _toUpdate.mpPolynomial == nullptr && _updateWith.mpPolynomial != nullptr )
242  _toUpdate.mpPolynomial = _updateWith.mpPolynomial;
243  // The factorization of the PolynomialFactorizationPair to update which can be empty, if constructed freshly by a polynomial.
244  if( _toUpdate.factorizedTrivially() && !_updateWith.factorizedTrivially() )
245  {
246  _toUpdate.mFactorization = _updateWith.mFactorization;
247  }
248  if( !factorizationsEqual( _toUpdate.factorization(), _updateWith.factorization() ) )
249  {
250  // Calculating the gcd refines both factorizations to the same factorization
251  bool refineA = false;
252  bool refineB = false;
253  Factorization<P> restA, restB;
254  typename P::CoeffType c( 0 );
255  gcd( _toUpdate, _updateWith, restA, restB, c, refineA, refineB );
256  assert( c == typename P::CoeffType( 0 ) || c == typename P::CoeffType( 1 ) );
257  }
258  }
259 
260  template<typename P>
261  P computePolynomial( const Factorization<P>& _factorization )
262  {
263  P result( 1 );
264  for (auto ft = _factorization.begin(); ft != _factorization.end(); ft++ )
265  {
266  result *= carl::pow(computePolynomial( ft->first ), ft->second );
267  }
268  return result;
269  }
270 
271  template<typename P>
272  P computePolynomial( const PolynomialFactorizationPair<P>& _pfPair )
273  {
274  if( _pfPair.mpPolynomial != nullptr )
275  return *_pfPair.mpPolynomial;
276  return computePolynomial( _pfPair.factorization() );
277  }
278 
279  template<typename P>
280  typename P::CoeffType PolynomialFactorizationPair<P>::flattenFactorization() const
281  {
282  typename P::CoeffType result( 0 );
283  if ( factorizedTrivially() )
284  {
285  return result;
286  }
287  std::lock_guard<std::recursive_mutex> lock( mMutex );
288  for( auto ft = mFactorization.begin(); ft != mFactorization.end(); )
289  {
290  if( ft->first.content().factorizedTrivially() )
291  {
292  if ( ft->first.coefficient() != 1 )
293  {
294  if( result == typename P::CoeffType( 0 ) )
295  result = typename P::CoeffType( 1 );
296  carl::exponent e = ft->second;
297  assert( e != 0 );
298  result *= carl::pow( ft->first.coefficient(), e );
299  ft->first.setCoefficient( 1 );
300  }
301  ++ft;
302  }
303  else
304  {
305  // Update factorization
306  if( result == typename P::CoeffType( 0 ) )
307  result = typename P::CoeffType( 1 );
308  const Factorization<P>& partFactorization = ft->first.factorization();
309  assert( carl::is_positive( ft->first.coefficient() ) );
310  carl::exponent e = ft->second;
311  assert( e != 0 );
312  result *= carl::pow( ft->first.coefficient(), e );
313 
314  for( auto partFactor = partFactorization.begin(); partFactor != partFactorization.end(); partFactor++ )
315  {
316  mFactorization.insert( std::pair<FactorizedPolynomial<P>, carl::exponent>( partFactor->first, partFactor->second * e ) );
317  }
318  ft = mFactorization.erase(ft);
319  }
320  }
321  assert( assertFactorization() );
322  return result;
323  }
324 
325  template<typename P>
326  bool PolynomialFactorizationPair<P>::isIrreducible() const
327  {
328  if ( mIrreducible != -1 )
329  return mIrreducible == 1;
330 
331  assert( mpPolynomial != nullptr );
332  if ( mpPolynomial->is_linear() )
333  {
334  mIrreducible = 1;
335  return true;
336  }
337  Definiteness definiteness = carl::definiteness(*mpPolynomial);
338  if ( definiteness == Definiteness::POSITIVE || definiteness == Definiteness::NEGATIVE )
339  {
340  mIrreducible = 1;
341  return true;
342  }
343  mIrreducible = 0;
344  return false;
345  }
346 
347  template<typename P>
348  void PolynomialFactorizationPair<P>::setNewFactors( const FactorizedPolynomial<P>& _fpolyA, carl::exponent exponentA, const FactorizedPolynomial<P>& _fpolyB, carl::exponent exponentB ) const
349  {
350  assert( mFactorization.size() == 1 );
351  assert( factorizedTrivially() );
352  assert( !_fpolyA.is_one() );
353  assert( !_fpolyB.is_one() );
354  //assert( carl::is_one(_fpolyA.coefficient()) );
355  //assert( carl::is_one(_fpolyB.coefficient()) );
356  assert( exponentA > 0 );
357  assert( exponentB > 0 );
358  mFactorization.clear();
359  if( _fpolyA == _fpolyB )
360  {
361  mFactorization.insert ( std::pair<FactorizedPolynomial<P>, carl::exponent>( _fpolyA, exponentA+exponentB ) );
362  }
363  else
364  {
365  mFactorization.insert ( std::pair<FactorizedPolynomial<P>, carl::exponent>( _fpolyA, exponentA ) );
366  mFactorization.insert ( std::pair<FactorizedPolynomial<P>, carl::exponent>( _fpolyB, exponentB ) );
367  }
368  assert( mpPolynomial != nullptr );
369  assert( *mpPolynomial == computePolynomial( mFactorization ) );
370  assert( assertFactorization() );
371  }
372 
373  template<typename P>
374  Factorization<P> gcd( const PolynomialFactorizationPair<P>& _pfPairA, const PolynomialFactorizationPair<P>& _pfPairB, Factorization<P>& _restA, Factorization<P>& _restB, typename P::CoeffType& _coeff, bool& _pfPairARefined, bool& _pfPairBRefined )
375  {
376  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", "****************************************************" );
377  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", "Compute GCD (internal) of " << _pfPairA << " and " << _pfPairB );
378  if( &_pfPairA == &_pfPairB )
379  return _pfPairA.factorization();
380  std::lock_guard<std::recursive_mutex> lockA( _pfPairA.mMutex );
381  std::lock_guard<std::recursive_mutex> lockB( _pfPairB.mMutex );
382  _coeff = typename P::CoeffType( 1 );
383  _pfPairARefined = false;
384  _pfPairBRefined = false;
385  _restA.clear();
386  _restB.clear();
387  // First flatten the factorizations of both polynomials
388  typename P::CoeffType cA = _pfPairA.flattenFactorization();
389  if( cA != typename P::CoeffType( 0 ) )
390  {
391  // Refinement took place -> set flag which indicates to update the first given factorized polynomial in its cache
392  _pfPairARefined = true;
393  _coeff *= cA;
394  }
395  typename P::CoeffType cB = _pfPairB.flattenFactorization();
396  if( cB != typename P::CoeffType( 0 ) )
397  {
398  // Refinement took place -> set flag which indicates to update the second given factorized polynomial in its cache
399  _pfPairBRefined = true;
400  _coeff *= cB;
401  }
402 // Factorization<P> factorizationA = _pfPairA.mFactorization;
403 // Factorization<P> factorizationB = _pfPairB.mFactorization;
404 // Factorization<P> result;
405  Factorization<P> factorizationA;
406  Factorization<P> factorizationB;
407  Factorization<P> result = commonDivisor( _pfPairA.mFactorization, _pfPairB.mFactorization, factorizationA, factorizationB );
408  while ( !factorizationA.empty() )
409  {
410  //Consider first factor in currently not checked factorization of A
411  FactorizedPolynomial<P> factorA = factorizationA.begin()->first;
412  carl::exponent exponentA = factorizationA.begin()->second;
413  factorizationA.erase( factorizationA.begin() );
414  bool breaked = false;
415  while ( !factorA.is_one() && !factorizationB.empty() )
416  {
417  // If the first factor is has a non trivial factorization -> flatten (this could be introduced during this loop even if we flatted beforehand)
418  const Factorization<P>& currentFactorizationA = factorA.factorization();
419  if( currentFactorizationA.size() != 1 )
420  {
421  factorizationA.insert( currentFactorizationA.begin(), currentFactorizationA.end() );
422  breaked = true;
423  break;
424  } else if ( currentFactorizationA.begin()->second > 1 ) {
425  // Polynomial has exponent > 1
426  exponentA *= currentFactorizationA.begin()->second;
427  factorA = currentFactorizationA.begin()->first;
428  }
429  // Take the second factor
430  FactorizedPolynomial<P> factorB = factorizationB.begin()->first;
431  carl::exponent exponentB = factorizationB.begin()->second;
432  factorizationB.erase( factorizationB.begin() );
433  // If the second factor is has a non trivial factorization -> flatten (this could be introduced during this loop even if we flatted beforehand)
434  const Factorization<P>& currentFactorizationB = factorB.factorization();
435  if( currentFactorizationB.size() != 1 )
436  {
437  factorizationB.insert( currentFactorizationB.begin(), currentFactorizationB.end() );
438  continue;
439  } else if ( currentFactorizationB.begin()->second > 1 ) {
440  // Polynomial has exponent > 1
441  exponentB *= currentFactorizationB.begin()->second;
442  factorB = currentFactorizationB.begin()->first;
443  }
444  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", "" );
445  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", "####################################################" );
446  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", "For first factorization: " );
447  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", " current factor: (" << factorA << ")^" << exponentA );
448  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", " open remainder: " << factorizationA );
449  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", " closed remainder: " << _restA );
450  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", "For second factorization: " );
451  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", " current factor: (" << factorB << ")^" << exponentB );
452  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", " open remainder: " << factorizationB );
453  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", " closed remainder: " << _restB );
454  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", "already found gcd is:" << result );
455  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", "" );
456  // Check if the factors are equal (even if we removed equal factor by commonDivisor() beforehand, common factors could be introduced during this loop)
457  if( factorA == factorB )
458  {
459  //Common factor found
460  carl::exponent exponentCommon = exponentA < exponentB ? exponentA : exponentB;
461  result.insert( std::pair<FactorizedPolynomial<P>, carl::exponent>( factorA, exponentCommon ) );
462  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", __LINE__ << ": add (" << factorA << ")^" << exponentCommon << " to gcd: " << result );
463  if (exponentA > exponentCommon)
464  {
465  factorizationA.insert( std::pair<FactorizedPolynomial<P>, carl::exponent>( factorA, exponentA-exponentCommon ) );
466  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", __LINE__ << ": add (" << factorA << ")^" << (exponentA-exponentCommon) << " to first open remainder: " << factorizationA );
467  }
468  if (exponentB > exponentCommon)
469  {
470  //Ignore FactorB as it has no remaining common factor with current FactorB
471  _restB.insert( std::pair<FactorizedPolynomial<P>, carl::exponent>( factorB, exponentB-exponentCommon ) );
472  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", __LINE__ << ": add (" << factorB << ")^" << (exponentB-exponentCommon) << " to second open remainder: " << _restB );
473  }
474  factorA = FactorizedPolynomial<P>( carl::constant_one<typename P::CoeffType>::get() );
475  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", __LINE__ << ": set current factor to: (1)^1" );
476  }
477  else
478  {
479  P polGCD, polA, polB;
480  assert( existsFactorization( factorA ) );
481  assert( existsFactorization( factorB ) );
482  if ( factorA.content().isIrreducible() && factorB.content().isIrreducible() )
483  polGCD = P( 1 );
484  else
485  {
486  //Compute GCD of factors
487  assert( factorA.content().mpPolynomial != nullptr );
488  assert( factorB.content().mpPolynomial != nullptr );
489  polA = *factorA.content().mpPolynomial;
490  polB = *factorB.content().mpPolynomial;
491  polGCD = carl::gcd( polA, polB );
492  if (carl::is_negative(polGCD.lcoeff())) {
493  polGCD = -polGCD;
494  }
495  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", __LINE__ << ": GCD of " << polA << " and " << polB << ": " << polGCD);
496  }
497 
498  if (is_one(polGCD))
499  {
500  //Ignore FactorB as it has no common factor with current FactorA
501  _restB.insert( std::pair<FactorizedPolynomial<P>, carl::exponent>( factorB, exponentB ) );
502  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", __LINE__ << ": add (" << factorB << ")^" << exponentB << " to second finished remainder: " << _restB );
503  }
504  else
505  {
506  //New common factor
507  P remainA = carl::quotient(polA, polGCD );
508  P remainB = carl::quotient(polB, polGCD );
509  carl::exponent exponentCommon = exponentA < exponentB ? exponentA : exponentB;
510  std::shared_ptr<Cache<PolynomialFactorizationPair<P>>> cache = factorA.pCache();
511  //Set new part of GCD
512  FactorizedPolynomial<P> gcdResult( polGCD, cache );
513  result.insert( std::pair<FactorizedPolynomial<P>, carl::exponent>( gcdResult, exponentCommon ) );
514  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", __LINE__ << ": add (" << gcdResult << ")^" << exponentCommon << " to gcd: " << result );
515  if (is_one(remainA))
516  {
517  if ( exponentA > exponentCommon )
518  {
519  exponentA -= exponentCommon;
520  }
521  else
522  {
523  factorA = FactorizedPolynomial<P>( carl::constant_one<typename P::CoeffType>::get() );
524  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", __LINE__ << ": set current factor to: (1)^1" );
525  }
526  }
527  else
528  {
529  //Set new factorization
530  FactorizedPolynomial<P> polRemainA( remainA, cache );
531  factorA.content().setNewFactors( gcdResult, 1, polRemainA, 1 );
532  _pfPairARefined = true;
533  factorA = polRemainA;
534  //Add remaining factorization
535  if (exponentA > exponentCommon)
536  {
537  factorizationA.insert( std::pair<FactorizedPolynomial<P>, carl::exponent>( gcdResult, exponentA-exponentCommon ) );
538  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", __LINE__ << ": add (" << gcdResult << ")^" << (exponentA-exponentCommon) << " to first open remainder: " << factorizationA );
539  }
540  }
541 
542  //Part of FactorB remains
543  if ( exponentB > exponentCommon )
544  {
545  factorizationB.insert( std::pair<FactorizedPolynomial<P>, carl::exponent>( gcdResult, exponentB-exponentCommon) );
546  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", __LINE__ << ": add (" << gcdResult << ")^" << (exponentB-exponentCommon) << " to second open remainder: " << factorizationB );
547  }
548 
549  if (!is_one(remainB))
550  {
551  //Set new factorization
552  FactorizedPolynomial<P> polRemainB( remainB, cache );
553  factorB.content().setNewFactors( gcdResult, 1, polRemainB, 1 );
554  _pfPairBRefined = true;
555  //Ignore remaining factorization as it has no common factor with FactorA anymore
556  factorizationB.insert( std::pair<FactorizedPolynomial<P>, carl::exponent>( polRemainB, exponentB) );
557  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", __LINE__ << ": add (" << polRemainB << ")^" << exponentB << " to second open remainder: " << _restB );
558  }
559  }
560  }
561  } //End of inner while
562  if( breaked )
563  continue;
564  //Insert remaining factorA into rest
565  if( !factorA.is_one() )
566  {
567  _restA.insert( std::pair<FactorizedPolynomial<P>, carl::exponent>( factorA, exponentA ) );
568  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", __LINE__ << ": add (" << factorA << ")^" << exponentA << " to first closed remainder: " << _restA );
569  }
570  //Reset factorizationB
571  factorizationB.insert( _restB.begin(), _restB.end() );
572  _restB.clear();
573  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", __LINE__ << ": shift second closed remainder to second open remainder: " << factorizationB );
574  } //End of outer while
575  _restB = factorizationB;
576  cA = _pfPairA.flattenFactorization();
577  if( cA != typename P::CoeffType( 0 ) )
578  {
579  _pfPairARefined = true;
580  _coeff *= cA;
581  }
582  cB = _pfPairB.flattenFactorization();
583  if( cB != typename P::CoeffType( 0 ) )
584  {
585  _pfPairBRefined = true;
586  _coeff *= cB;
587  }
588  assert( carl::is_one(_coeff) );
589  // Check correctness
590  assert( _pfPairA.assertFactorization() );
591  assert( _pfPairB.assertFactorization() );
592  CARL_LOG_DEBUG( "carl.core.factorizedpolynomial", "GCD (internal) of " << _pfPairA << " and " << _pfPairB << ": " << result << " with rests " << _restA << " and " << _restB );
593  assert( computePolynomial( result ) * computePolynomial( _restA ) == computePolynomial( _pfPairA ) );
594  assert( computePolynomial( result ) * computePolynomial( _restB ) == computePolynomial( _pfPairB ) );
595  assert( carl::is_one(carl::gcd( computePolynomial( _restA ), computePolynomial( _restB ))));
596  return result;
597  }
598 
599  template<typename P>
600  Factors<FactorizedPolynomial<P>> factor( const PolynomialFactorizationPair<P>& _pfPair, const typename P::CoeffType& _coeff )
601  {
602  typename P::CoeffType constantFactor = _coeff;
603  bool allFactorsIrreducible = false;
604  Factors<FactorizedPolynomial<P>> result;
605  if( !_pfPair.mFactorization.empty() )
606  {
607  _pfPair.mMutex.lock();
608  while( !allFactorsIrreducible )
609  {
610  _pfPair.flattenFactorization();
611  allFactorsIrreducible = true;
612  for( auto ft = _pfPair.mFactorization.begin(); ft != _pfPair.mFactorization.end(); )
613  {
614  assert( existsFactorization( ft->first ) );
615  if( ft->first.content().isIrreducible() )
616  {
617  ++ft;
618  }
619  else
620  {
621  assert( ft->first.factorization().size() == 1 );
622  assert( carl::is_one( ft->first.coefficient() ) );
623  carl::exponent e = ft->second;
624  assert( e != 0 );
625  Factors<typename FactorizedPolynomial<P>::PolyType> factorFactorization = carl::factor( ft->first.polynomial() );
626  Factorization<P> refinement;
627  for( const auto& pt : factorFactorization )
628  {
629  if( pt.first.is_constant() )
630  constantFactor *= pt.first.constant_part();
631  else
632  {
633  FactorizedPolynomial<P> fp(pt.first, ft->first.pCache());
634  constantFactor *= carl::pow( fp.coefficient(), pt.second );
635  refinement.insert( std::pair<FactorizedPolynomial<P>, carl::exponent>( fp.coprime_coefficients(), pt.second ) );
636  }
637  }
638  ft = _pfPair.mFactorization.erase(ft);
639 
640  for( auto partFactor = refinement.begin(); partFactor != refinement.end(); partFactor++ )
641  {
642  assert( existsFactorization( partFactor->first ) );
643  auto insertResult = _pfPair.mFactorization.insert( std::pair<FactorizedPolynomial<P>, carl::exponent>( partFactor->first, partFactor->second * e ) );
644  if( insertResult.first->first.factorization().size() > 1 ) // Note that factorization() takes care about the flattening
645  {
646  assert( insertResult.first->first.content().mIrreducible == 0 );
647  allFactorsIrreducible = false;
648  }
649  else
650  {
651  insertResult.first->first.content().mIrreducible = 1;
652  }
653  }
654  }
655  }
656  }
657  _pfPair.mMutex.unlock();
658  for( auto ft = _pfPair.mFactorization.begin(); ft != _pfPair.mFactorization.end(); ft++ )
659  {
660  result.insert( std::pair<FactorizedPolynomial<P>, carl::exponent>( ft->first, ft->second ) );
661  }
662  }
663  result.insert( std::pair<FactorizedPolynomial<P>, unsigned>( FactorizedPolynomial<P>(constantFactor), 1 ) );
664  return result;
665  }
666 
667  template <typename P>
668  std::ostream& operator<<(std::ostream& _out, const PolynomialFactorizationPair<P>& _pfPair)
669  {
670  if (_pfPair.factorizedTrivially()) {
671  return _out << _pfPair.polynomial();
672  } else {
673  return _out << factorizationToString(_pfPair.factorization());
674  }
675  }
676 
677 } // namespace carl