All Classes Functions Friends
urationalfunction.h
1 #ifndef _URATIONALFUNCTION_H_
2 #define _URATIONALFUNCTION_H_
3 
4 #include "../polynomial.h"
5 #include "rationalfunction_euclideanmethods.h"
6 #include "rationalfunction_symbolicintegration.h"
7 #include "multiprecision_rootfinding.h"
8 #include "rationalfunction_integrationpostprocessing.h"
9 #include "rationalfunction_symbolicnumericintegration.h"
10 #include "rationalfunction_integrationprinting.h"
11 #ifndef SERIAL
12  #include <cilktools/cilkview.h>
13 #else
14  #include <time.h>
15 #endif
16 
17 template <class UnivariatePolynomialOverField, class Field>
19  private:
20  UnivariatePolynomialOverField den;
21  UnivariatePolynomialOverField num;
22 
23  bool PROFILING;
24  bool floatingPointPrinting;
25  std::string outputFormatting;
26 
27  inline void normalize () {
28  num /= den.leadingCoefficient();
29  den /= den.leadingCoefficient();
30  }
31  public:
32  int characteristic;
33  static bool isPrimeField;
34  static bool isComplexField;
35  /**
36  * Construct the zero univariate rational function
37  *
38  * @param
39  **/
41  den.one();
42  num.zero();
43  PROFILING = false;
44  floatingPointPrinting = false;
45  outputFormatting = "MAPLE_OUTPUT";
46  UnivariatePolynomialOverField e;
47  characteristic = e.characteristic;
48  };
49 
50  /**
51  * Copy constructor
52  *
53  * @param b: A rational function
54  **/
55  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> (const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b) : den(b.den), num(b.num), PROFILING(b.PROFILING), floatingPointPrinting(b.floatingPointPrinting), outputFormatting(b.outputFormatting) {
56  characteristic = b.characteristic;
57  }
58 
59  /**
60  *
61  * @param a: the numerator
62  * @param b: the denominator
63  **/
64  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> (UnivariatePolynomialOverField a, UnivariatePolynomialOverField b) {
65  if (a.variable() != b.variable()) {
66  std::cout << "BPAS error: numerator and denominator must have the same variable." << std::endl;
67  exit(1);
68  }
69  if (b.isZero()) {
70  std::cout << "BPAS error: denominator is zero from UnivariateRationalFunction<UnivariatePolynomialOverField,Field>" << std::endl;
71  exit(1);
72  }
73  num = a;
74  den = b;
75  PROFILING = false;
76  floatingPointPrinting = false;
77  UnivariatePolynomialOverField e;
78  characteristic = e.characteristic;
79  outputFormatting = "MAPLE_OUTPUT";
80  }
81 
82  /**
83  * Destroy the rational function
84  *
85  * @param
86  **/
88 
89  inline void setVariableName(std::string name) {
90  num.setVariableName(name);
91  den.setVariableName(name);
92  }
93 
94  inline std::string variable() {
95  return num.variable();
96  }
97 
98  inline bool isProfiling() {
99  return PROFILING;
100  }
101 
102  inline void setProfiling(bool a) {
103  #ifndef SERIAL
104  PROFILING = a;
105  #endif
106  }
107 
108  inline bool isFloatingPointPrinting() {
109  return floatingPointPrinting;
110  }
111 
112  inline void setFloatingPointPrinting(bool a) {
113  floatingPointPrinting = a;
114  }
115 
116  inline bool isMapleOutput() {
117  if (outputFormatting == "MAPLE_OUTPUT")
118  return true;
119  else
120  return false;
121  }
122 
123  inline void setMapleOutput() {
124  outputFormatting = "MAPLE_OUTPUT";
125  }
126 
127  inline bool isMatlabOutput() {
128  if (outputFormatting == "MATLAB_OUTPUT")
129  return true;
130  else
131  return false;
132  }
133 
134  inline void setMatlabOutput() {
135  outputFormatting = "MATLAB_OUTPUT";
136  }
137 
138  inline void setNumerator(UnivariatePolynomialOverField& b) {
139  if (num.variable() != b.variable()) {
140  std::cout << "BPAS error: numerator and denominator must have the same variable." << std::endl;
141  exit(1);
142  }
143  num = b;
144  canonicalize();
145  }
146 
147  inline void setDenominator(UnivariatePolynomialOverField& b) {
148  if (num.variable() != b.variable()) {
149  std::cout << "BPAS error: numerator and denominator must have the same variable." << std::endl;
150  exit(1);
151  }
152  den = b;
153  canonicalize();
154  }
155 
156  inline void set(UnivariatePolynomialOverField& a, UnivariatePolynomialOverField& b) {
157  if (a.variable() != b.variable()) {
158  std::cout << "BPAS error: numerator and denominator must have the same variable." << std::endl;
159  exit(1);
160  }
161  num = a;
162  den = b;
163  canonicalize();
164  }
165 
166  inline UnivariatePolynomialOverField numerator() {
167  return num;
168  }
169 
170  inline UnivariatePolynomialOverField denominator() {
171  return den;
172  }
173 
175  return (!(num == b.num) || !(den == b.den));
176  }
178  return ((num == b.num) && (den == b.den));
179  }
180 
183  return (r += b);
184  }
186  if (num.variable()!= b.num.variable()) {
187  std::cout << "BPAS: error, trying to add between RationalFunction[" << num.variable() << "] and RationalFunction[" << b.num.variable() << "]." << std::endl;
188  exit(1);
189  }
190  UnivariatePolynomialOverField g, r1(den), r2(b.den);
191  g = den.gcd(b.den);
192  r1 /= g;
193  r2 /= g;
194  den = r1 * r2;
195  r1 *= b.num;
196  r2 *= num;
197  r1 += r2; // r1 = ar2 + cr1;
198  r2 = r1.gcd(g);
199  r1 /= r2; // r1 = e;
200  g /= r2; // g = g';
201  den *= g;
202  num = r1;
203  normalize();
204 
205  return *this;
206  }
207 
210  return (r -= b);
211  }
213  if (num.variable()!= b.num.variable()) {
214  std::cout << "BPAS: error, trying to subtract between RationalFunction[" << num.variable() << "] and RationalFunction[" << b.num.variable() << "]." << std::endl;
215  exit(1);
216  }
217  UnivariatePolynomialOverField g, r1(den), r2(b.den);
218  g = den.gcd(b.den);
219  r1 /= g;
220  r2 /= g;
221  den = r1 * r2;
222  r1 *= -b.num;
223  r2 *= num;
224  r1 += r2; // r1 = ar2 + cr1;
225  r2 = r1.gcd(g);
226  r1 /= r2; // r1 = e;
227  g /= r2; // g = g';
228  den *= g;
229  num = r1;
230  normalize();
231 
232  return *this;
233  }
234 
237  return r;
238  }
239 
240  /**
241  * Overload operator ^
242  * replace xor operation by exponentiation
243  *
244  * @param e: The exponentiation, e > 0
245  **/
248  res.setVariableName(num.variable());
249  if (isZero() || isOne() || e == 1)
250  res = *this;
251  else if (e == 2)
252  res = *this * *this;
253  else if (e > 2) {
255  res.one();
256 
257  while (e) {
258  if (e % 2) { res *= x; }
259  x = x * x;
260  e >>= 1;
261  }
262  }
263  else if (!e)
264  res.one();
265  else {
266  res = *this;
267  }
268  return res;
269  }
270  /**
271  * Overload operator ^=
272  * replace xor operation by exponentiation
273  *
274  * @param e: The exponentiation, e > 0
275  **/
277  *this = *this ^ e;
278  return *this;
279  }
280 
282  if (num.isZero()) {
283  std::cout << "BPAS error: division by zero from UnivariateRationalFunction<UnivariatePolynomialOverField,Field> inverse()" << std::endl;
284  exit(1);
285  }
287  r.normalize();
288  return r;
289  }
290 
293  return (r *= b);
294  }
296  if (num.variable()!= b.num.variable()) {
297  std::cout << "BPAS: error, trying to multiply between RationalFunction[" << num.variable() << "] and RationalFunction[" << b.num.variable() << "]." << std::endl;
298  exit(1);
299  }
300 
301  UnivariatePolynomialOverField g1, g2, r;
302  g1 = num.gcd(b.den);
303  g2 = den.gcd(b.num);
304  num /= g1;
305  r = b.num / g2;
306  num *= r;
307  den /= g2;
308  r = b.den / g1;
309  den *= r;
310  normalize();
311  return *this;
312  }
313 
316  return (r /= b);
317  }
319  if (b.isZero()) {
320  std::cout << "BPAS error: division by zero from UnivariateRationalFunction<UnivariatePolynomialOverField,Field> operator/=" << std::endl;
321  exit(1);
322  }
323 
324  if (num.variable()!= b.num.variable()) {
325  std::cout << "BPAS: error, trying to divide between RationalFunction[" << num.variable() << "] and RationalFunction[" << b.num.variable() << "]." << std::endl;
326  exit(1);
327  }
329  *this *= e;
330  return *this;
331  }
332 
333  inline void canonicalize() {
334  UnivariatePolynomialOverField temp;
335  temp = num.gcd(den);
336  num /= temp;
337  den /= temp;
338  num /= den.leadingCoefficient();
339  den /= den.leadingCoefficient();
340  }
341 
342  inline bool isZero() {
343  return num.isZero();
344  }
345  inline void zero() {
346  num.zero();
347  den.one();
348  }
349  inline bool isOne() {
350  return num.isOne() && den.isOne();
351  }
352  inline void one() {
353  num.one();
354  den.one();
355  }
356  inline bool isNegativeOne() {
357  return (num.isNegativeOne() && den.isOne()) || (num.isOne() && den.isNegativeOne());
358  }
359  inline void negativeOne() {
360  num.negativeOne();
361  den.one();
362  }
363  inline int isConstant() {
364  return num.isConstant() && den.isConstant();
365  }
366 
367  /**
368  * Overload operator =
369  *
370  * @param b: A rational function
371  **/
373  if (this != &b) {
374  num = b.num;
375  den = b.den;
376  characteristic = b.characteristic;
377  }
378  return *this;
379  }
380 
381  /**
382  * Overload stream operator <<
383  *
384  * @param out: Stream object
385  * @param b: The rational function
386  **/
387  inline friend std::ostream& operator<< (std::ostream &out, UnivariateRationalFunction<UnivariatePolynomialOverField,Field> b) {
388  out << "(" << b.num << ")/(" << b.den << ")";
389 
390  return out;
391  }
392 
394  std::vector<UnivariatePolynomialOverField> gg,hh;
396 
397  temp.setVariableName(num.variable());
398  h->setVariableName(num.variable());
399  g->clear();
400  _hermiteReduce<UnivariatePolynomialOverField,Field>(num,den,&gg,&hh);
401  int i(0);
402  while (i<gg.size()) {
403  temp.set(gg.at(i),gg.at(i+1));
404  g->push_back(temp);
405  i += 2;
406  }
407  temp.set(hh.at(0),hh.at(1));
408  *h = temp;
409  }
410 
411  void integrateRationalLogPart(std::vector< SparseUnivariatePolynomial<UnivariatePolynomialOverField> > *S, std::vector<UnivariatePolynomialOverField> *U) {
412  U->clear();
413  S->clear();
414 
415  _integrateRationalLogPart<UnivariatePolynomialOverField,Field>(S,U,num,den,PROFILING);
416  }
417 
418  void differentiate() {
419  /* Destructive rational function differentiation */
420  // input a/d
421  // output (a'*d-a*d')/d^2 = a'*e-a*f/d*e; e=d/g, f=d'/g, g=gcd(d,d')
422 
423  UnivariatePolynomialOverField D(den); // D = d
424  UnivariatePolynomialOverField dD(den);
425  UnivariatePolynomialOverField temp;
426  //std::cout << "here?" << std::endl;
427  dD.differentiate(1); // dD = d'
428  //std::cout << "dD = " << dD << std::endl;
429  temp = D.gcd(dD); // temp = g
430  //std::cout << "temp = " << temp << std::endl;
431  D /= temp; // D = e
432  //std::cout << "D = " << D << std::endl;
433  dD /= temp; // dD = f
434  //std::cout << "dD = " << dD << std::endl;
435  temp = -num; // temp = -a
436  //std::cout << "temp = " << temp << std::endl;
437  temp *= dD; // temp = -a*f
438  //std::cout << "temp = " << temp << std::endl;
439  dD = num; // dD = a
440  //std::cout << "dD = " << dD << std::endl;
441  dD.differentiate(1); // dD = a'
442  //std::cout << "dD = " << dD << std::endl;
443  dD *= D; // dD = a'*e
444  //std::cout << "dD = " << dD << std::endl;
445  temp += dD; // temp = a'*e-a*f
446  //std::cout << "temp = " << temp << std::endl;
447  D *= den; // D = d*e
448  //std::cout << "D = " << D << std::endl;
449 
450  //std::cout << "here?" << std::endl;
451  num = temp;
452  den = D;
453  canonicalize();
454  }
455 
456  void integrate(UnivariatePolynomialOverField *P, std::vector< UnivariateRationalFunction<UnivariatePolynomialOverField,Field> > *g, std::vector<UnivariatePolynomialOverField> *U, std::vector< SparseUnivariatePolynomial<UnivariatePolynomialOverField> > *S) {
457  g->clear();
458  U->clear();
459  S->clear();
460  std::vector<UnivariatePolynomialOverField> G;
462  temp.setVariableName(num.variable());
463 
464  #ifndef SERIAL
465  // Profiling variables
466  unsigned long long start;
467  unsigned long long end;
468  float elapsed = 0;
469  #else
470  // Profiling variables
471  clock_t start;
472  clock_t end;
473  float elapsed = 0;
474  #endif
475 
476  if (PROFILING){
477  std::cout << "integrate" << std::endl;
478  std::cout << "--------------------------------------" << std::endl;
479  #ifndef SERIAL
480  start = __cilkview_getticks();
481  #else
482  start = clock();
483  #endif
484  }
485 
486  _integrateRationalFunction<UnivariatePolynomialOverField,Field>(num,den,P,&G,U,S,PROFILING);
487 
488  if (PROFILING){
489  #ifndef SERIAL
490  end = __cilkview_getticks();
491  elapsed = (end - start) / 1000.f;
492  #else
493  end = clock();
494  elapsed = (float)(end - start) / CLOCKS_PER_SEC;
495  #endif
496  std::cout << "--------------------------------------" << std::endl;
497  std::cout << "integrate runtime: " << elapsed << " s" << std::endl;
498  }
499 
500  int i(0);
501  while (i<G.size()) {
502  temp.set(G.at(i),G.at(i+1));
503  g->push_back(temp);
504  i += 2;
505  }
506 
507  }
508 
509  void realSymbolicNumericIntegrate(UnivariatePolynomialOverField *P, std::vector< UnivariateRationalFunction<UnivariatePolynomialOverField,Field> > *g, std::vector<Field> *lg, std::vector<UnivariatePolynomialOverField> *Lg, std::vector<Field> *atn, std::vector<UnivariatePolynomialOverField> *Atn, int prec) {
510  P->zero();
511  g->clear();
512  lg->clear();
513  Lg->clear();
514  atn->clear();
515  Atn->clear();
516  std::vector<UnivariatePolynomialOverField> G;
518  temp.setVariableName(num.variable());
519 
520  std::cout << "[realSymbolicNumericIntegrate (snInt): Symbolic-Numeric Integration with BPAS and MPSolve]" << std::endl;
521  std::cout << "[Integration method: Hermite reduction, LRT integration]" << std::endl;
522  std::cout << "Starting..." << std::endl;
523 
524  #ifndef SERIAL
525  // Profiling variables
526  unsigned long long start;
527  unsigned long long end;
528  float elapsed = 0;
529  #else
530  // Profiling variables
531  clock_t start;
532  clock_t end;
533  float elapsed = 0;
534  #endif
535 
536  if (PROFILING){
537  std::cout << "--------------------------------------" << std::endl;
538  #ifndef SERIAL
539  start = __cilkview_getticks();
540  #else
541  start = clock();
542  #endif
543  }
544 
545  _realSNIntegrate<UnivariatePolynomialOverField,Field>(num,den,P,&G,lg,Lg,atn,Atn,prec,PROFILING);
546 
547  if (PROFILING){
548  #ifndef SERIAL
549  end = __cilkview_getticks();
550  elapsed = (end - start) / 1000.f;
551  #else
552  end = clock();
553  elapsed = (float)(end - start) / CLOCKS_PER_SEC;
554  #endif
555  std::cout << "--------------------------------------" << std::endl;
556  std::cout << "realSymbolicNumericIntegrate runtime: " << elapsed << " s" << std::endl;
557  }
558 
559  int i(0);
560  while (i<G.size()) {
561  temp.set(G.at(i),G.at(i+1));
562  g->push_back(temp);
563  i += 2;
564  }
565 
566  }
567 
568  void printIntegral(UnivariatePolynomialOverField &P, std::vector< UnivariateRationalFunction<UnivariatePolynomialOverField,Field> > &g, std::vector<UnivariatePolynomialOverField> &U, std::vector< SparseUnivariatePolynomial<UnivariatePolynomialOverField> > &S){
569  std::vector<UnivariatePolynomialOverField> G;
570  for (int i=0; i<g.size(); i++) {
571  G.push_back(g.at(i).num);
572  G.push_back(g.at(i).den);
573  }
574  _printFormalIntegral<UnivariatePolynomialOverField,Field>(num,den,P,G,U,S, false, floatingPointPrinting, false);
575  }
576 
577  void printIntegral(UnivariatePolynomialOverField &P, std::vector< UnivariateRationalFunction<UnivariatePolynomialOverField,Field> > &g, std::vector<Field> &lg, std::vector<UnivariatePolynomialOverField> &Lg, std::vector<Field> &atn, std::vector<UnivariatePolynomialOverField> &Atn){
578  std::vector<UnivariatePolynomialOverField> G;
579  for (int i=0; i<g.size(); i++) {
580  G.push_back(g.at(i).num);
581  G.push_back(g.at(i).den);
582  }
583  std::vector<UnivariatePolynomialOverField> empty;
584  _printIntegral<UnivariatePolynomialOverField,Field>(num,den,P,G,lg,Lg,atn,Atn,empty, false, floatingPointPrinting, false, outputFormatting);
585  }
586 };
587 #endif
588 /* This file is part of the BPAS library http://www.bpaslib.org
589 
590  BPAS is free software: you can redistribute it and/or modify
591  it under the terms of the GNU General Public License as published by
592  the Free Software Foundation, either version 3 of the License, or
593  (at your option) any later version.
594 
595  BPAS is distributed in the hope that it will be useful,
596  but WITHOUT ANY WARRANTY; without even the implied warranty of
597  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
598  GNU General Public License for more details.
599 
600  You should have received a copy of the GNU General Public License
601  along with BPAS. If not, see <http://www.gnu.org/licenses/>.
602 
603  Copyright: Changbo Chen <changbo.chen@hotmail.com>
604  Farnam Mansouri <mansouri.farnam@gmail.com>
605  Marc Moreno Maza <moreno@csd.uwo.ca>
606  Ning Xie <nxie6@csd.uwo.ca>
607  Yuzhen Xie <yuzhenxie@yahoo.ca>
608 
609 */
610 
611 
Definition: upolynomial.h:10
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > & operator^=(int e)
Definition: urationalfunction.h:276
Definition: polynomial.h:112
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > & operator=(UnivariateRationalFunction< UnivariatePolynomialOverField, Field > b)
Definition: urationalfunction.h:372
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > operator^(int e)
Definition: urationalfunction.h:246
Definition: urationalfunction.h:18