 |
GiNaCRA
0.6.4
|
|
GiNaCRA
0.6.4
|
00001 /* 00002 * GiNaCRA - GiNaC Real Algebra package 00003 * Copyright (C) 2010-2012 Ulrich Loup, Joachim Redies, Sebastian Junges 00004 * 00005 * This file is part of GiNaCRA. 00006 * 00007 * GiNaCRA is free software: you can redistribute it and/or modify 00008 * it under the terms of the GNU General Public License as published by 00009 * the Free Software Foundation, either version 3 of the License, or 00010 * (at your option) any later version. 00011 * 00012 * GiNaCRA is distributed in the hope that it will be useful, 00013 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 * GNU General Public License for more details. 00016 * 00017 * You should have received a copy of the GNU General Public License 00018 * along with GiNaCRA. If not, see <http://www.gnu.org/licenses/>. 00019 * 00020 */ 00021 00022 00023 #ifndef GINACRA_MULTIVARIATEPOLYNOMIALFACTORY_H 00024 #define GINACRA_MULTIVARIATEPOLYNOMIALFACTORY_H 00025 00026 // #define GINACRA_MULTIVARIATEPOLYNOMIAL_CHECK_STAIRCASEBORDER // check whether a Groebner basis defines a finite set of monomials under the staircase (as it is the case for special GBs) 00027 // #define GINACRA_MULTIVARIATEMONOMIAL_CHECK_SPECIALGROEBNERBASIS // check for non-zero-dimensional system 00028 00040 #include "MultivariatePolynomial.h" 00041 #include "MultivariateTerm.h" 00042 #include "MultivariateMonomial.h" 00043 #include "MultivariateCoefficient.h" 00044 #include "operators.h" 00045 #include <algorithm> 00046 00047 namespace GiNaCRA 00048 { 00050 // Auxiliary structs // 00052 00055 template<class monomialOrdering> 00056 struct symbol_is_greatereq_tdegree: 00057 public binary_function<symbol, symbol, bool> 00058 { 00059 MultivariatePolynomial<monomialOrdering> polynomial; 00060 00061 symbol_is_greatereq_tdegree( const MultivariatePolynomial<monomialOrdering>& polynomial ): 00062 polynomial( polynomial ) 00063 {} 00064 00070 bool operator ()( const symbol& a, const symbol& b ) const 00071 { 00072 return polynomial.tdegree( a ) >= polynomial.tdegree( b ); 00073 } 00074 }; 00075 00076 class MultivariatePolynomialFactory 00077 { 00078 private: 00079 00080 public: 00081 00083 // MultivariatePolynomial // 00085 00090 template<class monomialOrdering> 00091 static void checkForSpecialGroebnerBasis( const list<MultivariatePolynomial<monomialOrdering> >& basis ) throw ( invalid_argument ) 00092 { 00093 if( basis.size() <= 0 ) 00094 throw invalid_argument( "MultivariatePolynomial::monomialsUnderTheStaircase: Given Groebner basis is empty." ); 00095 00096 typename list<MultivariatePolynomial<monomialOrdering> >::const_iterator i = basis.begin(); 00097 vector<symbol> variables = i->Variables(); 00098 // first iteration: gather degrees and check local properties 00099 list<int> degrees = list<int>(); 00100 int d = 0; 00101 00102 for( ; i != basis.end(); ++i ) 00103 { 00104 if( variables.size() != i->Variables().size() ) 00105 throw invalid_argument( "MultivariatePolynomial::checkForSpecialGroebnerBasis: Not all elements of the given special Groebner basis have the same number of variables." ); 00106 00107 MultivariateTerm<monomialOrdering> t = MultivariateTerm<monomialOrdering>( *i ); 00108 d = t.tdegree(); 00109 00110 if( d <= (*i - t).tdegree() ) 00111 throw invalid_argument( "MultivariatePolynomial::checkForSpecialGroebnerBasis: Not all leading monomials of the given special Groebner basis have a greater total degree than the remainders." ); 00112 00113 degrees.push_back( d ); 00114 } 00115 00116 if( variables.size() != degrees.size() ) 00117 throw invalid_argument( "MultivariatePolynomial::checkForSpecialGroebnerBasis: There are not as many elements of the given special Groebner basis the number of variables." ); 00118 00119 // second iteration: check global properties involving all degrees 00120 i = basis.begin(); 00121 00122 for( ; i != basis.end(); ++i ) 00123 { 00124 MultivariatePolynomial<monomialOrdering> r = *i - MultivariateMonomial<monomialOrdering>( *i ); 00125 list<int>::const_iterator j = degrees.begin(); 00126 d = *j; // previous degree 00127 ++j; 00128 00129 for( ; j != degrees.end(); ++j ) 00130 { 00131 if( d < *j ) 00132 throw invalid_argument( "MultivariatePolynomial::checkForSpecialGroebnerBasis: Given special Groebner basis does not have decreasing degrees of its leading monomials." ); 00133 } 00134 00135 j = degrees.begin(); 00136 00137 for( vector<symbol>::const_iterator v = variables.begin(); v != variables.end(); ++v ) 00138 { 00139 if( *v != i->Variable() ) 00140 { 00141 if( i->degree( *v ) >= *j ) 00142 throw invalid_argument( "MultivariatePolynomial::checkForSpecialGroebnerBasis: Not all leading terms have a greater degree than all remainders." ); 00143 } 00144 00145 ++j; 00146 } 00147 } 00148 } 00149 00157 template<class monomialOrdering> 00158 static const MultivariatePolynomial<monomialOrdering> generalPositionTemplate( const vector<symbol>& variables, 00159 unsigned degree, 00160 numeric index ) 00161 { 00162 list<symbol> l = sortVariables( variables ); 00163 ex p = numeric( 1 ); 00164 00165 for( vector<symbol>::const_iterator i = l.begin(); i != l.end(); ++i ) 00166 { 00167 p += index * pow( *i, degree ); 00168 index *= index; // index^(i+1) 00169 } 00170 00171 return MultivariatePolynomial<monomialOrdering>( p, l ); // , pow(*l.begin(), degree) 00172 } 00173 00182 template<class monomialOrdering> 00183 static const MultivariatePolynomial<monomialOrdering> deformation( const MultivariatePolynomial<monomialOrdering>& polynomial, 00184 const list<int>& degrees, 00185 const ex& radius, 00186 const infinitesimal_symbol& infinitesimal ) 00187 throw ( invalid_argument ) 00188 { 00189 const vector<symbol> variables = polynomial.Variables(); 00190 // construct deformation polynomial 00191 list<int>::const_iterator degree = degrees.begin(); 00192 vector<symbol>::const_iterator variable = variables.begin(); 00193 ex newPolynomial; 00194 00195 for( ; degree != degrees.end(); ++degree, ++variable ) 00196 newPolynomial += pow( *variable, *degree ); 00197 00198 variable = variables.begin(); 00199 ++variable; 00200 00201 for( ; variable != variables.end(); ++variable ) 00202 newPolynomial += pow( *variable, 2 ); 00203 00204 newPolynomial *= radius; 00205 newPolynomial -= numeric( 2 * variables.size() - 1 ); 00206 return MultivariatePolynomial<monomialOrdering>( infinitesimal * newPolynomial 00207 + (1 - infinitesimal) * static_cast<ex>(polynomial), variables ); 00208 } 00209 00219 template<class monomialOrdering> 00220 static const list<MultivariatePolynomial<monomialOrdering> > specialGroebnerBasis( 00221 const MultivariatePolynomial<monomialOrdering>& polynomial, 00222 list<int> degrees, 00223 const ex& radius, 00224 const infinitesimal_symbol& infinitesimal ) 00225 throw ( invalid_argument ) 00226 { 00227 MultivariatePolynomial<monomialOrdering> deformation = MultivariatePolynomialFactory::deformation<monomialOrdering>( polynomial, 00228 degrees, radius, 00229 infinitesimal ); 00230 list<symbol> variables = list<symbol>( polynomial.Variables() ); // Note: variables->size() == degrees.size() 00231 int k = degrees.size(); 00232 // fill the list deformationDiffs with the partial derivatives of deformation w.r.t. all variables except the first 00233 lst deformationDiffs = lst(); 00234 vector<symbol>::const_iterator v = variables.erase( variables.begin() ); // first symbol not needed 00235 00236 for( ; v != variables.end(); ++v ) 00237 deformationDiffs.append( deformation.diff( v )); 00238 00239 // 00240 // Compute first deformation (w.r.t. first variable) as a deformation enriched by an appropriately large coefficient 00241 // and reduced modulo each partial derivative of the other deformations. 00242 // Note: 00243 // * the degree sequence should be descending but the last value should be larger than 2 because this polynomial whould 00244 // then be derived to a linear polynomial which is not wanted 00245 // * The reduction is performed in order to remove the higher-degree terms in other variables of the first deformation 00246 // polynomial. This requires a reduction modulo each deformation derivative and w.r.t. their respective leading monomials. 00247 00248 // compute coefficient 00249 numeric n = 2 * numeric( degrees.size() ) - 3; 00250 lst degreesEx = lst(); 00251 list<int>::const_iterator d = degrees.begin(); 00252 00253 for( ; d != degrees.end(); ++d ) 00254 degreesEx.append( *d ); 00255 00256 ex firstDegree = *degreesEx.begin(); 00257 degreesEx.remove_first(); 00258 deformation *= pow( prod( degreesEx ), 2 ) * pow( infinitesimal, n ) * pow( radius, n * firstDegree ); // reuse deformation as the first element of the Groebner basis 00259 MultivariatePolynomial<monomialOrdering> firstDeformation = MultivariatePolynomial<monomialOrdering>( deformation, 00260 polynomial.Variables() ); 00261 // compute reductions 00262 v = variables.begin(); 00263 d = degrees.erase( degrees.begin() ); 00264 lst::const_iterator def = deformationDiffs.begin(); 00265 00266 for( ; def != deformationDiffs.end(); ++def ) 00267 { 00268 ex m = pow( *v, *d ); 00269 MultivariatePolynomial<monomialOrdering> g( *def, polynomial.Variables() ); 00270 firstDeformation = firstDeformation.reduction( g, m ); // omitted (cf. p. 491 in ISBN 0-387-94090-1 and ISBN-13: 978-3642069642): .reduction(g, m); 00271 ++v; 00272 ++d; 00273 } 00274 00275 // fill the final output list 00276 list<MultivariatePolynomial<monomialOrdering> > basis = list<MultivariatePolynomial<monomialOrdering> >(); 00277 basis.push_back( firstDeformation ); 00278 00279 for( lst::const_iterator def = deformationDiffs.begin(); def != deformationDiffs.end(); ++def ) 00280 basis.push_back( MultivariatePolynomial<monomialOrdering>( *def, polynomial.Variables() )); 00281 00282 return basis; 00283 } 00284 00289 template<class monomialOrdering> 00290 static const list<int> synthesizeSpecialGroebnerBasisDegrees( const MultivariatePolynomial<monomialOrdering>& polynomial ) 00291 { 00292 list<int> degrees = list<int>(); 00293 vector<symbol> variables = vector<symbol>( polynomial.Variables() ); // Note: variables->size() == degrees.size() 00294 variables.sort( symbol_is_greatereq_tdegree<monomialOrdering>( polynomial )); // sort ascendingly by degree (=> first variable corresponds to max. total degree) 00295 00296 for( vector<symbol>::const_iterator v = variables.begin(); v != variables.end(); ++v ) 00297 { 00298 int d = max( polynomial.tdegree( *v ), 3 ); // do not go below a degree value of 4 (next even of 3) because these terms should be seperated from the squares in the deformation polynomial 00299 degrees.push_back( d % 2 == 0 ? d + 2 : d + 1 ); // next even number strictly greater than d 00300 } 00301 00302 return degrees; 00303 } 00304 00314 template<class monomialOrdering> 00315 static const list<MultivariatePolynomial<monomialOrdering> > synthesizeSpecialGroebnerBasis( 00316 const MultivariatePolynomial<monomialOrdering>& polynomial, 00317 const ex& radius, 00318 const infinitesimal_symbol& infinitesimal ) 00319 throw ( invalid_argument ) 00320 { 00321 return specialGroebnerBasis( polynomial, synthesizeSpecialGroebnerBasisDegrees( polynomial ), radius, infinitesimal ); 00322 } 00323 00325 // MultivariateMonomial // 00327 00332 template<class monomialOrdering> 00333 static const list<MultivariateMonomial<monomialOrdering> > cornersOfTheStaircase( 00334 const list<MultivariatePolynomial<monomialOrdering> >& gb ) 00335 throw ( invalid_argument ) 00336 { 00337 if( gb.size() <= 0 ) 00338 throw invalid_argument( "cornersOfTheStaircase: Given Groebner basis is empty." ); 00339 00340 list<MultivariateMonomial<monomialOrdering> > lmons = list<MultivariateMonomial<monomialOrdering> >(); 00341 00342 for( typename list<MultivariatePolynomial<monomialOrdering> >::const_iterator i = gb.begin(); i != gb.end(); ++i ) 00343 lmons.push_back( MultivariateMonomial<monomialOrdering>( *i )); // constructs MultivariateMonomial directly from MultivariateMonomial.Lmon 00344 00345 return lmons; 00346 } 00347 00352 template<class monomialOrdering> 00353 static const list<MultivariateMonomial<monomialOrdering> > monomialsUnderTheStaircase( 00354 const list<MultivariatePolynomial<monomialOrdering> >& gb ) 00355 throw ( invalid_argument ) 00356 { 00357 if( gb.size() <= 0 ) 00358 throw invalid_argument( "MultivariateMonomial::monomialsUnderTheStaircase: Given Groebner basis is empty." ); 00359 00360 // collect the degrees and corresponding main symbols for each main variable 00361 list<int> degrees = list<int>(); 00362 list<symbol> vars = vector<symbol>(); 00363 00364 for( typename list<MultivariatePolynomial<monomialOrdering> >::const_iterator i = gb.begin(); i != gb.end(); ++i ) 00365 { 00366 vars.push_back( i->Variable() ); 00367 degrees.push_back( i->lmon().degree( i->Variable() )); 00368 } 00369 00370 #ifdef GINACRA_MULTIVARIATEMONOMIAL_CHECK_SPECIALGROEBNERBASIS 00371 // check input to avoid non-termination (in case of non-zero-dimensional system) 00372 MultivariatePolynomialFactory::checkForSpecialGroebnerBasis<monomialOrdering>( gb ); 00373 #endif 00374 // construct monomials under the staircase 00375 list<MultivariateMonomial<monomialOrdering> > monomials = list<MultivariateMonomial<monomialOrdering> >(); 00376 monomials.push_back( MultivariateMonomial<monomialOrdering>( ex( 1 ), vars )); 00377 list<int>::const_iterator degree = degrees.begin(); 00378 00379 for( vector<symbol>::const_iterator var = vars.begin(); var != vars.end(); ++var ) 00380 { 00381 #ifdef GINACRA_MULTIVARIATEMONOMIAL_DEBUG 00382 cout << "MultivariateMonomial::monomialsUnderTheStaircase: (" << *var << ", " << *degree << ")" << endl; 00383 #endif 00384 searchMonomialsUnderTheStaircase( var, vars.end(), degree, MultivariateMonomial<monomialOrdering>( ex( 1 ), vars ), monomials ); 00385 ++degree; 00386 } 00387 00388 return monomials; 00389 } 00390 00397 template<class monomialOrdering> 00398 static const list<MultivariateTerm<monomialOrdering> > termsUnderTheStaircase( 00399 const list<MultivariateMonomial<monomialOrdering> >& monomials, 00400 const list<MultivariatePolynomial<monomialOrdering> >& gb, 00401 ex& coeff = *(new ex( 1 ))) 00402 throw ( invalid_argument ) 00403 { 00404 if( monomials.size() <= 0 || gb.size() <= 0 ) 00405 throw invalid_argument( "termsUnderTheStaircase: Given list of monomials under the staircase or the given special Groebner basis is empty." ); 00406 00407 list<MultivariateTerm<monomialOrdering> > terms = list<MultivariateTerm<monomialOrdering> >(); 00408 // fetch common leading coefficient of the given basis 00409 coeff = coefficientLcm( gb ); 00410 00411 // multiply appropriate powers of common coefficient to the monomials under the staircase 00412 for( typename list<MultivariateMonomial<monomialOrdering> >::const_iterator i = monomials.begin(); i != monomials.end(); ++i ) 00413 terms.push_back( MultivariateTerm<monomialOrdering>( pow( coeff, i->tdegree() ) * static_cast<ex>(*i), i->Variables() )); // since *i is a monomial, tdegree should determine the sum of exponents contained in *i 00414 00415 return terms; 00416 } 00417 00422 template<class monomialOrdering> 00423 static const list<MultivariateMonomial<monomialOrdering> > borderOfTheStaircase( 00424 const list<MultivariateMonomial<monomialOrdering> >& corners ) 00425 throw ( invalid_argument ) 00426 { 00427 if( corners.size() <= 0 ) 00428 throw invalid_argument( "borderOfTheStaircase: Given list of corners of the staircase is empty." ); 00429 00430 list<MultivariateMonomial<monomialOrdering> > monomials = list<MultivariateMonomial<monomialOrdering> >(); 00431 00432 for( typename list<MultivariateMonomial<monomialOrdering> >::const_iterator corners_current = corners.begin(); 00433 corners_current != corners.end(); ++corners_current ) 00434 searchBorderOfTheStaircase( corners_current, corners.begin(), corners.end(), *corners_current, monomials ); 00435 00436 return monomials; 00437 } 00438 00444 template<class monomialOrdering> 00445 static const list<MultivariateTerm<monomialOrdering> > borderTermsOfTheStaircase( 00446 const list<MultivariateMonomial<monomialOrdering> >& border, 00447 const list<MultivariatePolynomial<monomialOrdering> >& gb ) 00448 throw ( invalid_argument ) 00449 { 00450 if( border.size() <= 0 || gb.size() <= 0 ) 00451 throw invalid_argument( "borderTermsOfTheStaircase: Given border of the staircase or the given special Groebner basis is empty." ); 00452 00453 list<MultivariateTerm<monomialOrdering> > terms = list<MultivariateTerm<monomialOrdering> >(); 00454 // fetch common leading coefficient of the given basis 00455 const ex coeff = coefficientLcm( gb ); 00456 00457 // multiply appropriate powers of common coefficient to the monomials under the staircase 00458 for( typename list<MultivariateMonomial<monomialOrdering> >::const_iterator i = border.begin(); i != border.end(); ++i ) 00459 terms.push_back( MultivariateTerm<monomialOrdering>( pow( coeff, i->tdegree() ) * static_cast<ex>(*i), i->Variables() )); // since *i is a monomial, tdegree should determine the sum of exponents contained in *i 00460 00461 return terms; 00462 } 00463 00464 // Auxiliary methods 00465 00470 template<class monomialOrdering> 00471 static MultivariatePolynomial<monomialOrdering> coefficientLcm( const list<MultivariatePolynomial<monomialOrdering> >& gb ) 00472 { 00473 lst lcoeffs = lst(); 00474 for( typename list<MultivariatePolynomial<monomialOrdering> >::const_iterator i = gb.begin(); i != gb.end(); ++i ) 00475 lcoeffs.append( i->lcoeff() ); 00476 return MultivariatePolynomial<monomialOrdering>( lcm( lcoeffs ), gb.front().Variables() ); 00477 } 00478 00487 template<class monomialOrdering> 00488 static void searchMonomialsUnderTheStaircase( vector<symbol>::const_iterator variables_current, 00489 const vector<symbol>::const_iterator& variables_end, 00490 list<int>::const_iterator degrees_current, 00491 const MultivariateMonomial<monomialOrdering>& monomial, 00492 list<MultivariateMonomial<monomialOrdering> >& monomials ) 00493 throw ( invalid_argument ) 00494 { 00495 #ifdef GINACRA_MULTIVARIATEMONOMIAL_DEBUG 00496 cout << "MultivariateMonomial::searchMonomialsUnderTheStaircase: degree " << *degrees_current << " symbol " << *variables_current 00497 << endl; 00498 #endif 00499 00500 // terminate in case no new symbol is available for multiplication 00501 if( variables_current == variables_end ) 00502 return; 00503 00504 // multiply with current variable 00505 const MultivariateMonomial<monomialOrdering> m = monomial.mul( MultivariateMonomial<monomialOrdering>( *variables_current, 00506 monomial.Variables() )); 00507 00508 // terminate in case the degree bound for this variable is met 00509 if( m.degree( ex_to<symbol>( *variables_current )) >= *degrees_current ) 00510 return; 00511 00512 // add new monomial 00513 monomials.push_back( m ); 00514 00515 // try to multiply the next variables to the monomial 00516 for( ; variables_current != variables_end; ++variables_current, ++degrees_current ) 00517 searchMonomialsUnderTheStaircase( variables_current, variables_end, degrees_current, m, monomials ); 00518 } 00519 00528 template<class monomialOrdering> 00529 static void searchBorderOfTheStaircase( typename list<MultivariateMonomial<monomialOrdering> >::const_iterator corners_current, 00530 const typename list<MultivariateMonomial<monomialOrdering> >::const_iterator& corners_begin, 00531 const typename list<MultivariateMonomial<monomialOrdering> >::const_iterator& corners_end, 00532 const MultivariateMonomial<monomialOrdering>& monomial, 00533 list<MultivariateMonomial<monomialOrdering> >& monomials ) 00534 throw ( invalid_argument ) 00535 { 00536 // append current products of corners 00537 monomials.push_back( monomial ); 00538 // append all missing combinations of the current monomial (i.e. corner product) multiplied with other variables 00539 list<symbol> variables = vector<symbol>(); 00540 list<int> degrees = list<int>(); 00541 00542 for( typename list<MultivariateMonomial<monomialOrdering> >::const_iterator c = corners_begin; c != corners_end; ++c ) 00543 { 00544 if( monomial != *c ) 00545 { 00546 // corner is current monomial 00547 const_iterator p = monomial.begin(); 00548 00549 for( ; p != monomial.end(); ++p ) 00550 { 00551 if( *p == *c ) // corner occurs in current monomial 00552 break; 00553 } 00554 00555 if( p == monomial.end() ) 00556 { 00557 // corner does not occur in current monomial, so add corresponding variable and max. degree 00558 variables.push_back( c->Variable() ); 00559 degrees.push_back( c->degree() ); 00560 } 00561 } 00562 } 00563 00564 list<int>::const_iterator deg = degrees.begin(); 00565 00566 for( vector<symbol>::const_iterator var = variables.begin(); var != variables.end(); ++deg, ++var ) 00567 searchMonomialsUnderTheStaircase( var, variables.end(), deg, monomial, monomials ); 00568 00569 // proceed with the next corner combinations 00570 while( true ) 00571 { 00572 ++corners_current; // choose next corner 00573 00574 if( corners_current == corners_end ) 00575 return; 00576 00577 searchBorderOfTheStaircase( corners_current, corners_begin, corners_end, monomial.mul( *corners_current ), monomials ); 00578 } 00579 } 00580 00581 }; 00582 // class MultivariatePolynomialFactory 00583 00584 } // namespace GiNaC 00585 00586 #endif // GINACRA_MULTIVARIATEPOLYNOMIALFACTORY_H
1.8.0