1 #include "field_rationalfunctions2.h"
2 
3 #define USEPOLYNOMIALGCD 1
4 /* For some reason, setting USEFACTORY to 1 breaks code. For example
5 ./gfan _groebnerfan --stdin ~/gfan/current/examples/3x3of3x4
6 no longer works */
7 
8 #include <memory>
9 #include <assert.h>
10 #include <gmp.h>
11 #include <sstream>
12 #include <iostream>
13 #include <time.h>
14 
15 #include "termorder.h"
16 #include "division.h"
17 #include "buchberger.h"
18 #include "saturation.h"
19 #include "printer.h"
20 #include "linalg.h"
21 
22 #if USEPOLYNOMIALGCD
23 #include "polynomialgcd.h"
24 #endif
25 
26 #include "log.h"
27 
28 int FieldElementRationalFunctions2Living;
29 
30 /**
31  * This field is the field of _multivariate_ rational functions.
32  */
33 class FieldElementRationalFunction2 : public FieldElementImplementation
34 {
gcd(Polynomial a,Polynomial b)35 	Polynomial gcd(Polynomial a, Polynomial b)
36 	{
37 		if(a.degree(IntegerVector::standardVector(1,0))<b.degree(IntegerVector::standardVector(1,0)))
38 			swap(a,b);
39 		while(!b.isZero())
40 		{
41 			LexicographicTermOrder T;
42 
43 			a.mark(T);
44 			b.mark(T);
45 			PolynomialSet B(b.getRing());B.push_back(b);
46 
47 			b=division(a,B,T);
48 			a=B.front();
49 		}
50 		return a;
51 	}
randomizedFactorTest(PolynomialRing const & r,Polynomial const & p,Polynomial const & q)52 	bool randomizedFactorTest(PolynomialRing const &r, Polynomial const &p, Polynomial const &q)
53 	//returns true if it is discovered that no factor exists. False return value does not mean anything.
54 	{
55 		return false;//DISABLED
56 		vector<Polynomial> l;
57 		l.push_back(p);
58 		l.push_back(q);
59 		int n=r.getNumberOfVariables();
60 		for(int i=0;i<n;i++)
61 		{
62 			int numberOfRetries=10;
63 			retry:
64 			numberOfRetries--;
65 			if(numberOfRetries<0) return false;//<----------------------------------------
66 			FieldVector val(r.getField(),n-1);
67 			for(int j=0;j<n-1;j++)val[j]=r.getField().zHomomorphism((rand()&127)+1);
68 			vector<Polynomial> l2;
69 			PolynomialRing r2(r.getField(),1);
70 
71 		for(int I=0;I<2;I++)
72 			{
73 				Polynomial temp(r2);
74 				for(TermMap::const_iterator j=l[I].terms.begin();j!=l[I].terms.end();j++)
75 				{
76 					FieldElement mp=j->second;
77 					for(int k=0;k<n-1;k++)
78 						for(int l=j->first.exponent[k+(k>=i)];l>0;l--)mp*=val[k];
79 					IntegerVector expo(1);expo[0]=j->first.exponent[i];
80 					temp+=Term(mp,Monomial(r2,expo));
81 				}
82 				if(temp.isZero())goto retry;
83 				if(temp.degree(IntegerVector::standardVector(1,0)!=l[I].degree(IntegerVector::standardVector(n,i)))||
84 						temp.degree(-IntegerVector::standardVector(1,0)!=l[I].degree(-IntegerVector::standardVector(n,i))))
85 					goto retry;
86 				l2.push_back(temp);
87 			}
88 			//gcd
89 			if(gcd(l2[0],l2[1]).numberOfTerms()!=1)
90 				return false;
91 		}
92 		return true;
93 	}
94 
95 
96  public:
97   Polynomial p,q;
FieldElementRationalFunction2(FieldImplementation & a)98   FieldElementRationalFunction2(FieldImplementation &a):
99     FieldElementImplementation(a),
100     p(((FieldRationalFunctions2Implementation*)&a)->getPolynomialRing()),
101     q(Term(((FieldRationalFunctions2Implementation*)&a)->getPolynomialRing().getField().zHomomorphism(1),Monomial(((FieldRationalFunctions2Implementation*)&a)->getPolynomialRing())))
102   {
103     FieldElementRationalFunctions2Living++;
104   }
FieldElementRationalFunction2(FieldImplementation & a,int n_)105   FieldElementRationalFunction2(FieldImplementation &a,int n_):
106     FieldElementImplementation(a),
107     p(Term(((FieldRationalFunctions2Implementation*)&a)->getPolynomialRing().getField().zHomomorphism(n_),Monomial(((FieldRationalFunctions2Implementation*)&a)->getPolynomialRing()))),
108     q(Term(((FieldRationalFunctions2Implementation*)&a)->getPolynomialRing().getField().zHomomorphism(1),Monomial(((FieldRationalFunctions2Implementation*)&a)->getPolynomialRing())))
109     {
110       if(n_==0)p=Polynomial(((FieldRationalFunctions2Implementation*)&a)->getPolynomialRing());
111       FieldElementRationalFunctions2Living++;
112     }
normalize(bool doExpensiveGCD=true)113   void normalize(bool doExpensiveGCD=true)
114   {
115 	  if(p.isZero()){q=p.getRing().one();return;}
116 	  //STEP 1: Scale leading term of denominator to 1
117 	  LexicographicTermOrder T;
118 	  q.mark(T);
119 	  FieldElement s=q.getMarked().c;
120 	  if(!s.isOne())
121 	  {
122 		  s=s.inverse();
123 		  q*=s;
124 		  p*=s;
125 	  }
126 
127 	  //STEP 2: If denominator is monomial, then everything is easy
128 	  if(q.isMonomial())
129 	  {
130 		  if(q.totalDegree()==0)return;
131 		  if(p.isMonomial())
132 		  {
133 			  p.mark(T);
134 			  Monomial s(p.getRing(),-min(p.getMarked().m.exponent,q.getMarked().m.exponent));
135 			  p*=s;
136 			  q*=s;
137 			  return;
138 		  }
139 	  }
140 
141 	  //STEP 3: We factor out monomial factors in both p and q and remember them for later
142 	  Monomial qm(q.getRing(),q.greatestCommonMonomialDivisor());
143 	  q.saturate();
144 	  Monomial pm(p.getRing(),p.greatestCommonMonomialDivisor());
145 	  p.saturate();
146 
147 	  IntegerVector A=min(qm.exponent,pm.exponent);
148 	  qm.exponent-=A;
149 	  pm.exponent-=A;
150 
151 	  //STEP 4: If one of our polynomials is a monomial, then we can find no more common factors
152 	  if(q.isMonomial()||p.isMonomial())
153 	  {
154 		  p*=pm;
155 		  q*=qm;
156 		  return;
157 	  }
158 
159 	  // STEP 5: We check special cases where p divides q or q divides p.
160 	  q.mark(T);
161 	  p.mark(T);
162 	  {
163 		  //if q divides p then we don't need to compute gcd
164 		  PolynomialSet R(q.getRing());
165 		  PolynomialSet Q(q.getRing());
166 		  Q.push_back(q);
167 		  if(division(p,Q,T,&R).isZero())
168 		  {
169 			  //pout<<p<<"/"<<q<<"="<<*R.begin()<<"\n";
170 
171 			  p=*R.begin()*pm;
172 			  q=Term(qm.getRing().getField().zHomomorphism(-1),qm);
173 			  return;
174 		  }
175 	  }
176 	  {
177 		  //if p divides q then we don't need to compute gcd
178 		  PolynomialSet R(p.getRing());
179 		  PolynomialSet P(p.getRing());
180 		  P.push_back(p);
181 		  if(division(q,P,T,&R).isZero())
182 		  {
183 			  //pout<<p<<"/"<<q<<"="<<*R.begin()<<"\n";
184 
185 			  q=*R.begin()*qm;
186 			  p=Term(qm.getRing().getField().zHomomorphism(-1),pm);
187 			  return;
188 		  }
189 	  }
190 
191 	  // STEP 6: This should never happen....
192 	  if((q-p).isZero())
193 	  {
194 		  p=Term(pm.getRing().getField().zHomomorphism(1),pm);
195 		  q=Term(qm.getRing().getField().zHomomorphism(1),qm);
196 		  return;
197 	  }
198 
199 	  // STEP 7: Try to prove that no factor exists by random substitution (currently disabled)
200 	  if(randomizedFactorTest(p.getRing(),p,q))
201 	  {
202 		  //		  pout<<"FACTORTESTWORKS\n";
203 		  p*=pm;
204 		  q*=qm;
205 		  return;
206 	  }
207 	  //	  assert(!q.isZero());
208 	  /*	  if(p.isZero()){
209 		  q=q.getRing().one();
210 		  return;
211 	  }
212 	   */
213 //	  pout<<"p"<<p<<"q"<<q<<"\n";//<<p*q<<"\n";
214 	  static int FACTOR,NOFACTOR;
215 
216 	  // STEP 8: Finally we have to do a gcd computation somehow
217 	  if(!doExpensiveGCD){p*=pm;
218 	  q*=qm;return;}
219 
220 	  #if USEPOLYNOMIALGCD
221 	  {//computes gcd
222 
223 //		  debug<<"gcd\n";
224 //debug<<"BLA\n";
225 
226 //		  long start=clock();
227 
228 		  Polynomial r=polynomialGCD(p,q);
229 //		  if(clock()-start>CLOCKS_PER_SEC)
230 //		  {
231 //			  debug<<"GCD TAKES TOO LONG:\n"<<p.getRing()<<"\n{"<<p<<",\n"<<q<<"}\n"<<r<<"\n\n";
232 //		  }
233 
234 
235 //		  cerr<<"RETURNNNN---------------------===============\n";
236 //		  debug<<p<<"\n";
237 //		  debug<<q<<"\n";
238 //		  debug<<r<<"\n";
239 		  PolynomialSet R(p.getRing());
240 		  R.push_back(r);
241 		  R.mark(T);
242 		//  PolynomialSet pp(p.getRing());
243 		  p.mark(T);
244 		  q.mark(T);
245 //		  debug<<p<<"\n";
246 //		  debug<<q<<"\n";
247 		  PolynomialSet Q(p.getRing());
248 		  PolynomialSet P(p.getRing());
249 		  division(p,R,T,&P);
250 		  division(q,R,T,&Q);
251 		  p=P.front();
252 		  q=Q.front();
253 //		  debug<<p<<"\n";
254 //		  debug<<q<<"\n";
255 		  p*=pm;
256 		  q*=qm;
257 
258 //		  debug<<"dcg\n";
259 //		  debug<<p<<"\n";
260 //		  debug<<q<<"\n";
261 //		  cerr<<"RETURNNNN\n";
262 	  }
263 	  #else
264 	  {//computes lcm
265 	  PolynomialSet a(p.getRing());
266 	  PolynomialSet b(q.getRing());
267 	  a.push_back(p);
268 	  b.push_back(q);
269 //	  cerr<<"Ideal Intersection\n";
270 	  PolynomialSet c=idealIntersection(a,b);
271 //	  debug<<c;
272 //	  cerr<<"Done Ideal Intersection\n";
273 	  if(c.size()!=1)
274 	  {
275 		  //		  pout<<a<<b<<c;
276 		  assert(c.size()==1);
277 	  }
278 	  //Polynomial
279 	  Polynomial r=*(c.begin());
280 
281 
282 	  //	  LexicographicTermOrder T;
283 	  {
284 		  PolynomialSet Q(p.getRing());
285 		  PolynomialSet P(p.getRing());
286 		  r.mark(T);
287 		  a.mark(T);
288 		  b.mark(T);
289 		  division(r,b,T,&P);
290 		  division(r,a,T,&Q);
291 		  if(p.numberOfTerms()!=P.begin()->numberOfTerms()){FACTOR++;//pout<<p<<"/"<<q<<":"<<*P.begin()<<"/"<<*Q.begin()<<"\n";
292 		  }
293 		  else NOFACTOR++;
294 		  p=*P.begin();
295 		  q=*Q.begin();
296 		  //pout<<"r="<<r<<"\n";
297 		  //pout<<"p="<<p<<"\n";
298 		  //pout<<"q="<<q<<"\n";
299 		  p*=pm;
300 		  q*=qm;
301 	  }
302 	  }
303 #endif
304 	 // if((NOFACTOR&(256*2-1))==0)pout<<"FAC"<<FACTOR<<","<<NOFACTOR<<"\n";
305   }
FieldElementRationalFunction2(FieldImplementation & a,Polynomial const & p_,Polynomial const & q_)306   FieldElementRationalFunction2(FieldImplementation &a, Polynomial const &p_, Polynomial const &q_):
307     FieldElementImplementation(a),
308     p(p_),
309     q(q_)
310     {
311 
312       FieldElementRationalFunctions2Living++;
313     }
~FieldElementRationalFunction2()314   virtual ~FieldElementRationalFunction2()
315     {
316       FieldElementRationalFunctions2Living--;
317     }
operator =(const FieldElementRationalFunction2 & a)318   FieldElementRationalFunction2& operator=(const FieldElementRationalFunction2& a)
319     {
320       assert(0);
321       return *this;
322     }
323 
operator *=(const FieldElementImplementation & a)324   void operator*=(const FieldElementImplementation &a)
325     {
326       const FieldElementRationalFunction2 *A=(const FieldElementRationalFunction2*)&a;
327       assert(A);
328 #if 1
329 /*      debug<<"----\n";
330       debug<<p<<" - "<<A->q<<"\n";
331       debug<<q<<" - "<<A->p<<"\n";
332       debug<<(q*A->p)<<" - "<<p*A->q<<"\n";
333 */
334       Polynomial fac1=polynomialGCD(p,A->q);
335       Polynomial fac2=polynomialGCD(q,A->p);
336 
337 //      Polynomial pTemp=p;
338 //      Polynomial qTemp=q;
339 
340 //    	  p*=A->p;
341 //    	  q*=A->q;
342 //    	  normalize();
343 //#if 1
344 //    	  int totalTerms=p.numberOfTerms()+q.numberOfTerms();
345 
346  //   	  p=pTemp;
347  //   	  q=qTemp;
348 #if 0
349       p*=A->p;
350       q*=A->q;
351 
352 //      debug<<"A\n";
353       bool err=false;
354       err|=!fac1.divides(p,&p);
355       err|=!fac1.divides(q,&q);
356       err|=!fac2.divides(p,&p);
357       err|=!fac2.divides(q,&q);
358 //      debug<<"B\n";
359       assert(!err);
360 #else
361       Polynomial tp=p.exactlyDividedBy(fac1);
362       Polynomial tAq=A->q.exactlyDividedBy(fac1);
363       Polynomial tq=q.exactlyDividedBy(fac2);
364       Polynomial tAp=A->p.exactlyDividedBy(fac2);
365 
366       p=tp*tAp;
367       q=tq*tAq;
368 #endif
369       normalize(false);
370 
371 /*      if(totalTerms!=p.numberOfTerms()+q.numberOfTerms())
372       {
373 //    	  debug<<pTemp<<"/"<<qTemp<<"\n";
374 //    	  debug<<A->p<<"/"<<A->q<<"\n";
375     	  p=pTemp;
376     	  q=qTemp;
377     	  p*=A->p;
378     	  q*=A->q;
379     	  normalize();
380 //    	  debug<<A->p<<"/"<<A->q<<"\n";
381 
382   //  	  assert(0);
383       }*/
384 //#endif
385 
386 #else
387       p*=A->p;
388       q*=A->q;
389       normalize();
390 #endif
391     }
operator +=(const FieldElementImplementation & a)392   void operator+=(const FieldElementImplementation &a)
393     {
394       const FieldElementRationalFunction2 *A=(const FieldElementRationalFunction2*)&a;
395       assert(A);
396 
397       p=p*A->q+A->p*q;
398       q=A->q*q;
399       normalize();
400     }
madd(const FieldElementImplementation & a,const FieldElementImplementation & b)401   void madd(const FieldElementImplementation &a,const FieldElementImplementation &b)
402     {
403       const FieldElementRationalFunction2 *A=(const FieldElementRationalFunction2*)&a;
404       const FieldElementRationalFunction2 *B=(const FieldElementRationalFunction2*)&b;
405       assert(A);
406       assert(B);
407 #if 0
408       p=p*(A->q*B->q)+(A->p*B->p)*q;
409       q=A->q*B->q*q;
410       normalize();
411 #else
412       Polynomial fac1=polynomialGCD(B->p,A->q);
413       Polynomial fac2=polynomialGCD(B->q,A->p);
414 
415       Polynomial tBp=B->p.exactlyDividedBy(fac1);
416       Polynomial tAq=A->q.exactlyDividedBy(fac1);
417       Polynomial tBq=B->q.exactlyDividedBy(fac2);
418       Polynomial tAp=A->p.exactlyDividedBy(fac2);
419 
420       p=p*(tBq*tAq)+(tBp*tAp)*q;
421       q=tBq*tAq*q;
422       normalize();
423 #endif
424     }
425   FieldElementRationalFunction2 *one() const;
isZero() const426   bool isZero() const
427     {
428       return p.isZero();
429     }
430 
sum(const FieldElementImplementation & b) const431   FieldElementRationalFunction2 *sum(const FieldElementImplementation &b)const
432     {
433       const FieldElementRationalFunction2 *B=(const FieldElementRationalFunction2*)&b;
434 
435       Polynomial fac1=polynomialGCD(q,B->q);
436       Polynomial nq=B->q.exactlyDividedBy(fac1)*q;
437       Polynomial np=p*B->q.exactlyDividedBy(fac1)+B->p*q.exactlyDividedBy(fac1);
438       // It is still possible that more common factors exist. These must appear as factors of fac1
439       FieldElementRationalFunction2 *r= new FieldElementRationalFunction2(*getField(),np,nq);
440       r->normalize();
441 //      FieldElementRationalFunction2 *r= new FieldElementRationalFunction2(*getField(),p*B->q-B->p*q,B->q*q);
442       return r;
443 #if 0
444       const FieldElementRationalFunction2 *B=(const FieldElementRationalFunction2*)&b;
445 
446 //      pout<<B->q<<"------"<<q<<"\n";
447 
448       Polynomial quotient(q.getRing());
449       if(B->q.divides(q,&quotient))
450       {
451 //          pout<<B->q<<"+-----"<<q<<"--"<<quotient<<"\n";
452           return new FieldElementRationalFunction2(*getField(),p-B->p*quotient,q);
453       }
454       else if(q.divides(B->q,&quotient))
455       {
456           return new FieldElementRationalFunction2(*getField(),B->p-p*quotient,B->q);
457       }
458 
459 
460       FieldElementRationalFunction2 *r= new FieldElementRationalFunction2(*getField(),p*B->q+B->p*q,B->q*q);
461 
462       return r;
463 #endif
464     }
difference(const FieldElementImplementation & b) const465   FieldElementRationalFunction2 *difference(const FieldElementImplementation &b)const
466     {
467       const FieldElementRationalFunction2 *B=(const FieldElementRationalFunction2*)&b;
468 
469       Polynomial fac1=polynomialGCD(q,B->q);
470       Polynomial nq=B->q.exactlyDividedBy(fac1)*q;
471       Polynomial np=p*B->q.exactlyDividedBy(fac1)-B->p*q.exactlyDividedBy(fac1);
472       // It is still possible that more common factors exist. These must appear as factors of fac1
473       FieldElementRationalFunction2 *r= new FieldElementRationalFunction2(*getField(),np,nq);
474       r->normalize();
475 //      FieldElementRationalFunction2 *r= new FieldElementRationalFunction2(*getField(),p*B->q-B->p*q,B->q*q);
476       return r;
477     }
negation() const478   FieldElementRationalFunction2 *negation()const
479     {
480       FieldElementRationalFunction2 *r= new FieldElementRationalFunction2(*getField(),p-p-p,q);
481 
482       return r;
483     }
inverse() const484   FieldElementImplementation *inverse()const
485   {
486     if(isZero())
487       {
488 	AsciiPrinter P(Stderr);
489 	P.printString("Error inverting FieldElement: ");
490 	P.printPolynomial(p);
491 	P.printString(" ");
492 	P.printPolynomial(q);
493 	//	P.printFieldElement(*this);
494 	P.printString("\n");
495 	assert(0);
496       }
497 
498     FieldElementRationalFunction2 *r= new FieldElementRationalFunction2(*getField(),q,p);
499 
500     return r;
501   }
502 
sign() const503   int sign()const
504   {
505 	  assert(0);//not an ordered field (yet)
506     if(isZero())return 0;
507     return p.terms.rbegin()->second.sign();
508   }
509 
LaTeXTranslator(const string & s)510   static string LaTeXTranslator(const string &s)
511   {
512 	  assert(0);//not supported yet
513 /*    int startIndex=0;
514     string sign;
515     if(s[0]=='-')
516       {
517 	sign=string("-");
518 	startIndex=1;
519       }
520     int slashIndex=-1;
521     for(int i=startIndex;i<s.length();i++)if(s[i]=='/')slashIndex=i;
522     if(slashIndex==-1)
523       return string(s);
524 
525     return sign+string("{").append(s,startIndex,slashIndex-startIndex)+string("\\over ").append(s,slashIndex+1,s.length()-slashIndex-1)+string("}");
526 */
527 	  }
528 
toString(bool writeIfOne=true,bool alwaysWriteSign=false,bool latexMode=false) const529   std::string toString(bool writeIfOne=true, bool alwaysWriteSign=false, bool latexMode=false /*, bool mathMode=true*/) const
530   {
531     stringstream s;
532 
533     //    s << "(" <<p.toString(latexMode) << "/" << q.toString(latexMode) << ")";
534     s << "(";
535     s<<p.toString(latexMode);
536     s << "/";
537     s << q.toString(latexMode);
538     s << ")";
539     return s.str();
540   }
541 
copy() const542   FieldElementRationalFunction2 *copy()const
543   {
544     FieldElementRationalFunction2 *r= new FieldElementRationalFunction2(*getField());
545     r->p=p;
546     r->q=q;
547 
548     return r;
549   }
550 };
551 
getPolynomialRing() const552 PolynomialRing FieldRationalFunctions2Implementation::getPolynomialRing()const
553 {
554   return thePolynomialRing;
555 }
556 
isRationals() const557 bool FieldRationalFunctions2Implementation::isRationals()const
558 {
559   return false;
560 }
561 
getCharacteristic() const562 int FieldRationalFunctions2Implementation::getCharacteristic()const
563 {
564 	return this->getPolynomialRing().getField().getCharacteristic();
565 }
566 
567 
FieldRationalFunctions2Implementation(PolynomialRing const & r)568 FieldRationalFunctions2Implementation::FieldRationalFunctions2Implementation(PolynomialRing const &r):
569   thePolynomialRing(r)
570 {
571 }
572 
toString() const573 std::string FieldRationalFunctions2Implementation::toString()const
574 {
575   stringstream s;
576   s<< thePolynomialRing.getField().toString() << "("<<thePolynomialRing.toStringVariableNames()<< ")";
577   return s.str();
578 }
579 
zHomomorphismImplementation(int n)580 FieldElementImplementation *FieldRationalFunctions2Implementation::zHomomorphismImplementation(int n)
581 {
582   FieldElementImplementation *ret=new FieldElementRationalFunction2(*this,n);
583   return ret;
584 }
585 
zHomomorphism(int n)586 FieldElement FieldRationalFunctions2Implementation::zHomomorphism(int n)
587 {
588   return FieldElement(zHomomorphismImplementation(n));
589 }
590 
name()591 const char *FieldRationalFunctions2Implementation::name()
592 {
593   return "Rational functions in n variables";
594 }
595 
one() const596 FieldElementRationalFunction2 *FieldElementRationalFunction2::one() const
597 {
598   return new FieldElementRationalFunction2(*getField(),1);
599 }
600 
601 
602 
FieldRationalFunctions2(PolynomialRing const & r)603 FieldRationalFunctions2::FieldRationalFunctions2(PolynomialRing const &r):
604   Field(new FieldRationalFunctions2Implementation(r))
605 {
606 }
607 
608 
polynomialToFraction(Polynomial const & p)609 FieldElement FieldRationalFunctions2::polynomialToFraction(Polynomial const &p)
610 {
611 #if USESHAREDPTR
612 	FieldRationalFunctions2Implementation *imp=dynamic_cast<FieldRationalFunctions2Implementation*>(implementingObject.get());
613 #else
614 	FieldRationalFunctions2Implementation *imp=dynamic_cast<FieldRationalFunctions2Implementation*>(implementingObject);
615 #endif
616   Polynomial q=Term(imp->getPolynomialRing().getField().zHomomorphism(1),Monomial(imp->getPolynomialRing()));
617 
618   return new FieldElementRationalFunction2(*imp, p, q);
619 }
620 
621 /*****************************************************
622  * Conversion functions
623  *****************************************************/
makeVariablesParameters(PolynomialRing const & r,int numberOfParameters)624 PolynomialRing makeVariablesParameters(PolynomialRing const &r, int numberOfParameters)
625 {
626 	assert(numberOfParameters>=0);
627 	assert(numberOfParameters<=r.getNumberOfVariables());
628 	vector<string> names(numberOfParameters);
629 	for(int i=0;i<numberOfParameters;i++)names[i]=r.getVariableName(i);
630 	PolynomialRing temp(r.getField(),names);
631 	vector<string> names2(r.getNumberOfVariables()-numberOfParameters);
632 	for(int i=0;i<names2.size();i++)names2[i]=r.getVariableName(i+numberOfParameters);
633 
634 	return PolynomialRing(FieldRationalFunctions2(temp),names2);
635 }
636 
makeVariablesParameters(PolynomialRing const & genericRing,Polynomial const & p)637 Polynomial makeVariablesParameters(PolynomialRing const &genericRing, Polynomial const &p)
638 {
639 	Polynomial ret(genericRing);
640 #if USESHAREDPTR
641 	FieldRationalFunctions2Implementation const *coefficientField=dynamic_cast<FieldRationalFunctions2Implementation const*>(genericRing.getField().implementingObject.get());
642 #else
643 	FieldRationalFunctions2Implementation const *coefficientField=dynamic_cast<FieldRationalFunctions2Implementation const*>(genericRing.getField().implementingObject);
644 #endif
645 	FieldRationalFunctions2 &DANGER=(FieldRationalFunctions2&)genericRing.getField();
646 	PolynomialRing coefRing=coefficientField->getPolynomialRing();
647 
648 	for(TermMap::const_iterator i=p.terms.begin();i!=p.terms.end();i++)
649 	{
650 		FieldElement c=i->second;
651 		IntegerVector v=i->first.exponent;
652 		IntegerVector coefficientExponent=v.subvector(0,p.getRing().getNumberOfVariables()-genericRing.getNumberOfVariables());
653 		IntegerVector monomialExponent=v.subvector(p.getRing().getNumberOfVariables()-genericRing.getNumberOfVariables(),v.size());
654 		FieldElement c2=DANGER.polynomialToFraction(Term( c,Monomial(coefRing, coefficientExponent)));//does the numerator not belong to a field?
655 		ret+=Polynomial(Term(c2,Monomial(genericRing,monomialExponent)));
656 	}
657 	return ret;
658 }
659 
makeVariablesParameters(PolynomialRing const & genericRing,PolynomialSet const & p)660 PolynomialSet makeVariablesParameters(PolynomialRing const &genericRing, PolynomialSet const &p)
661 {
662 	PolynomialSet ret(genericRing);
663 	for(PolynomialSet::const_iterator i=p.begin();i!=p.end();i++)
664 		ret.push_back(makeVariablesParameters(genericRing,*i));
665 
666 	return ret;
667 }
668