1 #include "field_rationalfunctions.h"
2
3 #include <memory>
4 #include <assert.h>
5 #include <cstdio> /* Always include cstdio before gmp.h.*/
6 #include <gmp.h>
7 #include <sstream>
8 #include <iostream>
9
10 #include "termorder.h"
11 #include "division.h"
12 #include "buchberger.h"
13 #include "printer.h"
14
15 #include "log.h"
16
17 int FieldElementRationalFunctionsLiving;
18 /**
19 * This field is the field of _univariate_ rational functions.
20 */
21 class FieldElementRationalFunction : public FieldElementImplementation
22 {
23 public:
24 Polynomial p,q;
FieldElementRationalFunction(FieldImplementation & a)25 FieldElementRationalFunction(FieldImplementation &a):
26 FieldElementImplementation(a),
27 p(((FieldRationalFunctionsImplementation*)&a)->getPolynomialRing()),
28 q(Term(((FieldRationalFunctionsImplementation*)&a)->getPolynomialRing().getField().zHomomorphism(1),Monomial(((FieldRationalFunctionsImplementation*)&a)->getPolynomialRing(),IntegerVector(1))))
29 {
30 FieldElementRationalFunctionsLiving++;
31 }
FieldElementRationalFunction(FieldImplementation & a,int n_)32 FieldElementRationalFunction(FieldImplementation &a,int n_):
33 FieldElementImplementation(a),
34 p(Term(((FieldRationalFunctionsImplementation*)&a)->getPolynomialRing().getField().zHomomorphism(n_),Monomial(((FieldRationalFunctionsImplementation*)&a)->getPolynomialRing(),IntegerVector(1)))),
35 q(Term(((FieldRationalFunctionsImplementation*)&a)->getPolynomialRing().getField().zHomomorphism(1),Monomial(((FieldRationalFunctionsImplementation*)&a)->getPolynomialRing(),IntegerVector(1))))
36 {
37 if(n_==0)p=Polynomial(((FieldRationalFunctionsImplementation*)&a)->getPolynomialRing());
38 /* fprintf(stderr,"constructing\n");
39 AsciiPrinter(Stderr).printPolynomial(p);
40 AsciiPrinter(Stderr).printPolynomial(q);
41 fprintf(stderr,"\n");*/
42 FieldElementRationalFunctionsLiving++;
43 }
normalize()44 void normalize()
45 {
46 PolynomialSet g(p.getRing());
47 g.push_back(p);
48 g.push_back(q);
49 LexicographicTermOrder T;
50 buchberger(&g,T);
51 minimize(&g);
52 assert(g.size()==1);
53 Polynomial r=g.front();
54 PolynomialSet Q(p.getRing());
55 PolynomialSet P(p.getRing());
56 division(p,g,T,&P);
57 division(q,g,T,&Q);
58 p=*P.begin();
59 q=*Q.begin();
60 assert(!q.isZero());
61 // FieldElement a=q.terms.begin()->second;
62 // TermMap temp=q.terms.end();
63 FieldElement a=q.terms.rbegin()->second;
64 p*=a;
65 q*=a.inverse();
66 // assert(0);//make q's leading coefficient 1
67 }
FieldElementRationalFunction(FieldImplementation & a,Polynomial const & p_,Polynomial const & q_)68 FieldElementRationalFunction(FieldImplementation &a, Polynomial const &p_, Polynomial const &q_):
69 FieldElementImplementation(a),
70 p(p_),
71 q(q_)
72 {
73
74 FieldElementRationalFunctionsLiving++;
75 }
76 /* FieldElementRationalFunction(FieldImplementation &a, mpq_t *n_):FieldElementImplementation(a)
77 {
78 FieldElementRationalsLiving++;
79 mpq_init(value);
80 mpq_set(value,*n_);
81 }*/
~FieldElementRationalFunction()82 virtual ~FieldElementRationalFunction()
83 {
84 FieldElementRationalsLiving--;
85 // mpq_clear(value);
86 }
operator =(const FieldElementRationalFunction & a)87 FieldElementRationalFunction& operator=(const FieldElementRationalFunction& a)
88 {
89 assert(0);
90 /* const FieldElementRational *A=(const FieldElementRational*)&a;
91 if (this != A) {
92 mpq_clear(value);
93 mpz_init_set(mpq_numref(value), mpq_numref(a.value));
94 mpz_init_set(mpq_denref(value), mpq_denref(a.value));
95 }
96 */
97 return *this;
98 }
99
operator *=(const FieldElementImplementation & a)100 void operator*=(const FieldElementImplementation &a)
101 {
102 const FieldElementRationalFunction *A=(const FieldElementRationalFunction*)&a;
103 assert(A);
104
105 // TimerScope ts(&rationalTimer);
106 p*=A->p;
107 q*=A->q;
108 normalize();
109 // mpq_mul(value,value,A->value);
110 }
operator +=(const FieldElementImplementation & a)111 void operator+=(const FieldElementImplementation &a)
112 {
113 const FieldElementRationalFunction *A=(const FieldElementRationalFunction*)&a;
114 assert(A);
115
116 p=p*A->q+A->p*q;
117 q=A->q*q;
118 normalize();
119 }
madd(const FieldElementImplementation & a,const FieldElementImplementation & b)120 void madd(const FieldElementImplementation &a,const FieldElementImplementation &b)
121 {
122 const FieldElementRationalFunction *A=(const FieldElementRationalFunction*)&a;
123 const FieldElementRationalFunction *B=(const FieldElementRationalFunction*)&b;
124 assert(A);
125 assert(B);
126
127 p=p*(A->q*B->q)+(A->p*B->p)*q;
128 q=A->q*B->q*q;
129 normalize();
130 }
131
132 FieldElementRationalFunction *one() const;
isZero() const133 bool isZero() const
134 {
135 return p.isZero();
136 }
137
sum(const FieldElementImplementation & b) const138 FieldElementRationalFunction *sum(const FieldElementImplementation &b)const
139 {
140 const FieldElementRationalFunction *B=(const FieldElementRationalFunction*)&b;
141
142 FieldElementRationalFunction *r= new FieldElementRationalFunction(*getField(),p*B->q+B->p*q,B->q*q);
143
144 return r;
145 }
difference(const FieldElementImplementation & b) const146 FieldElementRationalFunction *difference(const FieldElementImplementation &b)const
147 {
148 const FieldElementRationalFunction *B=(const FieldElementRationalFunction*)&b;
149 FieldElementRationalFunction *r= new FieldElementRationalFunction(*getField(),p*B->q-B->p*q,B->q*q);
150 return r;
151 }
negation() const152 FieldElementRationalFunction *negation()const
153 {
154 FieldElementRationalFunction *r= new FieldElementRationalFunction(*getField(),p-p-p,q);
155
156 return r;
157 }
inverse() const158 FieldElementImplementation *inverse()const
159 {
160 if(isZero())
161 {
162 AsciiPrinter P(Stderr);
163 P.printString("Error inverting FieldElement: ");
164 // P.printFieldElement(*this);
165 P.printString("\n");
166 assert(0);
167 }
168
169 FieldElementRationalFunction *r= new FieldElementRationalFunction(*getField(),q,p);
170
171 return r;
172 }
173
sign() const174 int sign()const
175 {
176 if(isZero())return 0;
177 return p.terms.rbegin()->second.sign();
178 }
179
LaTeXTranslator(const string & s)180 static string LaTeXTranslator(const string &s)
181 {
182 int startIndex=0;
183 string sign;
184 if(s[0]=='-')
185 {
186 sign=string("-");
187 startIndex=1;
188 }
189 int slashIndex=-1;
190 for(int i=startIndex;i<s.length();i++)if(s[i]=='/')slashIndex=i;
191 if(slashIndex==-1)
192 return string(s);
193
194 return sign+string("{").append(s,startIndex,slashIndex-startIndex)+string("\\over ").append(s,slashIndex+1,s.length()-slashIndex-1)+string("}");
195 }
196
toString(bool writeIfOne=true,bool alwaysWriteSign=false,bool latexMode=false) const197 std::string toString(bool writeIfOne=true, bool alwaysWriteSign=false, bool latexMode=false /*, bool mathMode=true*/) const
198 {
199 stringstream s;
200
201 // s << "(" <<p.toString(latexMode) << "/" << q.toString(latexMode) << ")";
202 s << "(";
203 s<<p.toString(latexMode);
204 s << "/";
205 s << q.toString(latexMode);
206 s << ")";
207 /*
208 s<<"(";
209 FieldElementRational *tempOne=one();
210 FieldElementRational *temp=difference(*tempOne);
211 bool isOne=temp->isZero();
212 //fprintf(Stderr,"DELETE\n");
213 delete temp;
214 temp=sum(*tempOne);
215 bool isMinusOne=temp->isZero();
216 //fprintf(Stderr,"DELETE\n");
217 delete temp;
218 //fprintf(Stderr,"DELETE\n");
219 delete tempOne;
220
221 if(!writeIfOne && isOne)
222 {
223 if(alwaysWriteSign)return std::string("+");
224 return std::string("");
225 }
226 if(!writeIfOne && isMinusOne)
227 return std::string("-");
228 static char s[1290*1000];
229 // mpq_get_str(s,10,value); //// CHECK BUFFER SIZE!!!!
230 // Changed to make code gmp 3.1.1 compatible
231 mpz_get_str(s,10,mpq_numref(value)); //CHECK BUFFER SIZE!!!!
232 string S(s);
233 if(mpz_cmp_ui(mpq_denref(value),1)!=0)
234 {
235 mpz_get_str(s,10,mpq_denref(value)); //CHECK BUFFER SIZE!!!!
236 S=S+string("/")+string(s);
237 }
238
239 if(latexMode)S=LaTeXTranslator(S);
240
241 if(alwaysWriteSign && mpq_sgn(value)!=-1)
242 return std::string("+")+S;
243 return S;
244 */
245 return s.str();
246 }
247
copy() const248 FieldElementRationalFunction *copy()const
249 {
250 FieldElementRationalFunction *r= new FieldElementRationalFunction(*getField());
251 // fprintf(Stderr,"NEW\n");
252 /*
253 mpq_clear(r->value);
254 mpz_init_set(mpq_numref(r->value), mpq_numref(value));
255 mpz_init_set(mpq_denref(r->value), mpq_denref(value));
256 */
257 r->p=p;
258 r->q=q;
259
260 return r;
261 }
262 };
263
getPolynomialRing() const264 PolynomialRing FieldRationalFunctionsImplementation::getPolynomialRing()const
265 {
266 return thePolynomialRing;
267 }
268
isRationals() const269 bool FieldRationalFunctionsImplementation::isRationals()const
270 {
271 return false;
272 }
273
getCharacteristic() const274 int FieldRationalFunctionsImplementation::getCharacteristic()const
275 {
276 return this->coefficientField.getCharacteristic();
277 }
278
279
FieldRationalFunctionsImplementation(Field const & f_,string const & parameterName_)280 FieldRationalFunctionsImplementation::FieldRationalFunctionsImplementation(Field const &f_, string const ¶meterName_):
281 coefficientField(f_),
282 parameterName(parameterName_),
283 thePolynomialRing(f_,1)
284 {
285 vector<string> l;
286 l.push_back(parameterName_);
287 thePolynomialRing= PolynomialRing(f_,l);
288 }
289
toString() const290 std::string FieldRationalFunctionsImplementation::toString()const
291 {
292 stringstream s;
293 s<< coefficientField.toString() << "(" << parameterName << ")";
294 return s.str();
295 }
296
zHomomorphismImplementation(int n)297 FieldElementImplementation *FieldRationalFunctionsImplementation::zHomomorphismImplementation(int n)
298 {//indsaet igen
299 /* if(n==0)
300 {
301 static FieldElementImplementation *p;
302 if(p==0)p=new FieldElementRational(*this,0);
303 p->refCount++;
304 return p;
305 }
306 else
307 if(n==1)
308 {
309 static FieldElementImplementation *p;
310 if(p==0)p=new FieldElementRational(*this,1);
311 p->refCount++;
312 return p;
313 }
314 else
315 if(n==-1)
316 {
317 static FieldElementImplementation *p;
318 if(p==0)p=new FieldElementRational(*this,-1);
319 p->refCount++;
320 return p;
321 }
322 */
323 FieldElementImplementation *ret=new FieldElementRationalFunction(*this,n);
324 // fprintf(Stderr,"NEW\n");
325 // ret->refCount++;
326 return ret;
327 }
328
zHomomorphism(int n)329 FieldElement FieldRationalFunctionsImplementation::zHomomorphism(int n)
330 {
331 // fprintf(Stderr,"NEW\n");
332 // return FieldElement(new FieldElementRational(*this,n));
333 return FieldElement(zHomomorphismImplementation(n));
334 }
335
name()336 const char *FieldRationalFunctionsImplementation::name()
337 {
338 return "Rational functions";
339 }
340
one() const341 FieldElementRationalFunction *FieldElementRationalFunction::one() const
342 {
343 // fprintf(Stderr,"NEW\n");
344 return new FieldElementRationalFunction(*getField(),1);
345 }
346
347
348
FieldRationalFunctions(Field const & coefficientField,string const & parameterName)349 FieldRationalFunctions::FieldRationalFunctions(Field const &coefficientField, string const ¶meterName):
350 Field(new FieldRationalFunctionsImplementation(coefficientField,parameterName))
351 {
352 }
353
354
exponent(int power)355 FieldElement FieldRationalFunctions::exponent(int power)
356 {
357 #if USESHAREDPTR
358 FieldRationalFunctionsImplementation *imp=dynamic_cast<FieldRationalFunctionsImplementation*>(implementingObject.get());
359 #else
360 FieldRationalFunctionsImplementation *imp=dynamic_cast<FieldRationalFunctionsImplementation*>(implementingObject);
361 #endif
362 IntegerVector v(1);
363 v[0]=power;
364 Polynomial p=Term(imp->getPolynomialRing().getField().zHomomorphism(1),Monomial(imp->getPolynomialRing(),v));
365 Polynomial q=Term(imp->getPolynomialRing().getField().zHomomorphism(1),Monomial(imp->getPolynomialRing(),IntegerVector(1)));
366
367 return new FieldElementRationalFunction(*imp, p, q);
368 }
369
370
fromCoefficientField(FieldElement const & c)371 FieldElement FieldRationalFunctions::fromCoefficientField(FieldElement const &c)
372 {
373 #if USESHAREDPTR
374 FieldRationalFunctionsImplementation *imp=dynamic_cast<FieldRationalFunctionsImplementation*>(implementingObject.get());
375 #else
376 FieldRationalFunctionsImplementation *imp=dynamic_cast<FieldRationalFunctionsImplementation*>(implementingObject);
377 #endif
378 IntegerVector v(1);
379 Polynomial p=Term(c,Monomial(imp->getPolynomialRing(),v));
380 Polynomial q=Term(imp->getPolynomialRing().getField().zHomomorphism(1),Monomial(imp->getPolynomialRing(),IntegerVector(1)));
381
382 return new FieldElementRationalFunction(*imp, p, q);
383 }
384
385
386
substitute(FieldElement const & e,FieldElement const & tvalue) const387 FieldElement FieldRationalFunctions::substitute(FieldElement const &e, FieldElement const &tvalue)const
388 {
389 #if USESHAREDPTR
390 FieldElementRationalFunction *eimp=dynamic_cast<FieldElementRationalFunction*>(e.implementingObject.get());
391 #else
392 FieldElementRationalFunction *eimp=dynamic_cast<FieldElementRationalFunction*>(e.implementingObject);
393 #endif
394 FieldElement p=eimp->p.evaluate(tvalue);
395 FieldElement q=eimp->q.evaluate(tvalue);
396
397 assert(!q.isZero());
398
399 return p*q.inverse();
400 }
401
402 #include "field_rationals.h"
403 #include <iostream>
testRationalFunctionField()404 void testRationalFunctionField()
405 {
406 FieldRationalFunctions F(Q,string("t"));
407 AsciiPrinter(Stderr).printField(F);
408
409
410 FieldElement a=F.zHomomorphism(2);
411
412 FieldElement b=F.exponent(3);
413
414 a=b*a;
415 #if USESHAREDPTR
416 AsciiPrinter(Stderr).printPolynomial(((dynamic_cast<FieldElementRationalFunction*>((a.implementingObject.get())))->p));
417 #else
418 AsciiPrinter(Stderr).printPolynomial(((dynamic_cast<FieldElementRationalFunction*>((a.implementingObject)))->p));
419 #endif
420 /* // cerr<< ((FieldElementRationalFunction*)&(a.implementingObject))-> p.toString(false);
421 cerr<< (dynamic_cast<FieldElementRationalFunction*>(a.implementingObject))-> p.toString(false);
422 */
423
424
425 cerr<<a.implementingObject->toString();
426 }
427