Basic Polynomial Algebra Subprograms (BPAS)  v. 1.652
urationalfunction.h
1 #ifndef _URATIONALFUNCTION_H_
2 #define _URATIONALFUNCTION_H_
3
4 #include <complex>
5 #include <mps/mpc.h>
6 #include "../polynomial.h"
7 #include "rationalfunction_euclideanmethods.h"
8 #include "rationalfunction_symbolicintegration.h"
9 #include "multiprecision_rootfinding.h"
10 #include "rationalfunction_integrationpostprocessing.h"
11 #include "rationalfunction_symbolicnumericintegration.h"
12 #include "rationalfunction_integrationprinting.h"
13 #include "rationalfunction_complexrationalnumberordering.h"
14
15 struct IntegralTerm
16 {
17  int Sindex;
18  int tindex;
19
20  bool operator()(const IntegralTerm a,const IntegralTerm b) const {
21  if (a.Sindex > b.Sindex)
22  return 1;
23  else if (a.Sindex < b.Sindex)
24  return 0;
25  else {
26  if (a.tindex > b.tindex)
27  return 1;
28  else
29  return 0;
30  }
31  }
32
33  bool operator<(const IntegralTerm b) const {
34  if (this->Sindex < b.Sindex)
35  return 1;
36  else if (this->Sindex > b.Sindex)
37  return 0;
38  else {
39  if (this->tindex < b.tindex)
40  return 1;
41  else
42  return 0;
43  }
44  }
45
46  inline friend std::ostream& operator<< (std::ostream &out, IntegralTerm b) {
47  out << "[" << b.Sindex << "," << b.tindex << "]";
48  return out;
49  }
50 };
51
52 struct complexMPF
53 {
54  mpc_t c;
55 };
56
57 typedef std::map<int, IntegralTerm> RootIntegralTermMap;
58 typedef RootIntegralTermMap::const_iterator RITMIter;
59 typedef std::multimap<IntegralTerm, int> IntegralTermRootMap;
60 typedef IntegralTermRootMap::const_iterator ITRMIter;
61 typedef std::multimap<ComplexRationalNumber, int, CompareByRealThenReverseImaginary> residueRootIndexMap;
62 typedef residueRootIndexMap::const_iterator RRIMIter;
63 typedef std::multimap<int, int> TermRootMap;
64 typedef TermRootMap::const_iterator TRMIter;
65
66 /**
67  * A univariate rational function templated by a unvariate polynomial over a field.
68  * The univariate polynomial and the coefficient BPASField must be passed separately
69  * and explicitly.
70  */
71 template <class UnivariatePolynomialOverField, class Field>
72 class UnivariateRationalFunction : public BPASRationalFunction<UnivariatePolynomialOverField, UnivariateRationalFunction<UnivariatePolynomialOverField,Field>>,
73  private Derived_from<Field,BPASField<Field>> {
74  private:
75  UnivariatePolynomialOverField den;
76  UnivariatePolynomialOverField num;
77
78  bool PROFILING;
79  bool ERROR_ANALYSIS;
80  bool PFD_LOG_PART;
81  bool floatingPointPrinting;
82  std::string outputFormatting;
83
84  inline void normalize () {
87  }
88  public:
89  mpz_class characteristic;
90  static bool isPrimeField;
91  static bool isSmallPrimeField;
92  static bool isComplexField;
93  /**
94  * Construct the zero univariate rational function
95  *
96  * @param
97  **/
99  den.one();
100  num.zero();
101  PROFILING = false;
102  ERROR_ANALYSIS = false;
103  PFD_LOG_PART = false;
104  floatingPointPrinting = false;
105  outputFormatting = "MAPLE_OUTPUT";
106  UnivariatePolynomialOverField e;
107  characteristic = e.characteristic;
108  };
109
110  /**
111  * Copy constructor
112  *
113  * @param b: A rational function
114  **/
116 ERROR_ANALYSIS(b.ERROR_ANALYSIS), PFD_LOG_PART(b.PFD_LOG_PART), floatingPointPrinting(b.floatingPointPrinting), outputFormatting(b.outputFormatting) {
117  characteristic = b.characteristic;
118  }
119
120  /**
121  *
122  * @param a: the numerator
123  * @param b: the denominator
124  **/
125  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> (UnivariatePolynomialOverField a, UnivariatePolynomialOverField b) {
126  if (a.variable() != b.variable()) {
127  std::cout << "BPAS error: numerator and denominator must have the same variable." << std::endl;
128  exit(1);
129  }
130  if (b.isZero()) {
131  std::cout << "BPAS error: denominator is zero from UnivariateRationalFunction<UnivariatePolynomialOverField,Field>" << std::endl;
132  exit(1);
133  }
134  num = a;
135  den = b;
136  PROFILING = false;
137  ERROR_ANALYSIS = false;
138  PFD_LOG_PART = false;
139  floatingPointPrinting = false;
140  UnivariatePolynomialOverField e;
141  characteristic = e.characteristic;
142  outputFormatting = "MAPLE_OUTPUT";
143  }
144
145  /**
146  * Destroy the rational function
147  *
148  * @param
149  **/
151
152  inline void setVariableName(Symbol name) {
153  num.setVariableName(name);
154  den.setVariableName(name);
155  }
156
157  inline Symbol variable() {
158  return num.variable();
159  }
160
161  inline bool isProfiling() {
162  return PROFILING;
163  }
164
165  inline void setProfiling(bool a) {
166  PROFILING = a;
167  }
168
169  inline bool isAnalyzingError() {
170  return ERROR_ANALYSIS;
171  }
172
173  inline void setAnalyzingError(bool a) {
174  ERROR_ANALYSIS = a;
175  }
176
177  inline bool isPFDLogPart() {
178  return PFD_LOG_PART;
179  }
180
181  inline void setPFDLogPart(bool a) {
182  PFD_LOG_PART = a;
183  }
184
185  inline bool isFloatingPointPrinting() {
186  return floatingPointPrinting;
187  }
188
189  inline void setFloatingPointPrinting(bool a) {
190  floatingPointPrinting = a;
191  }
192
193  inline bool isMapleOutput() {
194  if (outputFormatting == "MAPLE_OUTPUT")
195  return true;
196  else
197  return false;
198  }
199
200  inline void setMapleOutput() {
201  outputFormatting = "MAPLE_OUTPUT";
202  }
203
205  if (outputFormatting == "MATLAB_OUTPUT")
206  return true;
207  else
208  return false;
209  }
210
212  outputFormatting = "MATLAB_OUTPUT";
213  }
214
215  inline void setNumerator(const UnivariatePolynomialOverField& b) {
216  if (num.variable() != b.variable()) {
217  std::cout << "BPAS error: numerator and denominator must have the same variable." << std::endl;
218  exit(1);
219  }
220  num = b;
221  canonicalize();
222  }
223
224  inline void setDenominator(const UnivariatePolynomialOverField& b) {
225  if (num.variable() != b.variable()) {
226  std::cout << "BPAS error: numerator and denominator must have the same variable." << std::endl;
227  exit(1);
228  }
229  den = b;
230  canonicalize();
231  }
232
233  inline void set(const UnivariatePolynomialOverField& a, const UnivariatePolynomialOverField& b) {
234  if (a.variable() != b.variable()) {
235  std::cout << "BPAS error: numerator and denominator must have the same variable." << std::endl;
236  exit(1);
237  }
238  num = a;
239  den = b;
240  canonicalize();
241  }
242
243  inline UnivariatePolynomialOverField numerator() const {
244  return num;
245  }
246
247  inline UnivariatePolynomialOverField denominator() const {
248  return den;
249  }
250
251  inline Field evaluate(const Field& c) {
252  Field output;
253  output = num.evaluate(c);
254  output /= den.evaluate(c);
255  return output;
256  }
257
258  inline bool operator!= (const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b) const {
259  return (!(num == b.num) || !(den == b.den));
260  }
261  inline bool operator== (const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b) const {
262  return ((num == b.num) && (den == b.den));
263  }
264
267  return (r += b);
268  }
270  if (num.variable()!= b.num.variable()) {
271  std::cout << "BPAS: error, trying to add between RationalFunction[" << num.variable() << "] and RationalFunction[" << b.num.variable() << "]." << std::endl;
272  exit(1);
273  }
274  UnivariatePolynomialOverField g, r1(den), r2(b.den);
275  g = den.gcd(b.den);
276  r1 /= g;
277  r2 /= g;
278  den = r1 * r2;
279  r1 *= b.num;
280  r2 *= num;
281  r1 += r2; // r1 = ar2 + cr1;
282  r2 = r1.gcd(g);
283  r1 /= r2; // r1 = e;
284  g /= r2; // g = g';
285  den *= g;
286  num = r1;
287  normalize();
288
289  return *this;
290  }
291
294  return (r -= b);
295  }
297  if (num.variable()!= b.num.variable()) {
298  std::cout << "BPAS: error, trying to subtract between RationalFunction[" << num.variable() << "] and RationalFunction[" << b.num.variable() << "]." << std::endl;
299  exit(1);
300  }
301  UnivariatePolynomialOverField g, r1(den), r2(b.den);
302  g = den.gcd(b.den);
303  r1 /= g;
304  r2 /= g;
305  den = r1 * r2;
306  r1 *= -b.num;
307  r2 *= num;
308  r1 += r2; // r1 = ar2 + cr1;
309  r2 = r1.gcd(g);
310  r1 /= r2; // r1 = e;
311  g /= r2; // g = g';
312  den *= g;
313  num = r1;
314  normalize();
315
316  return *this;
317  }
318
321  return r;
322  }
323
324  /**
326  * replace xor operation by exponentiation
327  *
328  * @param e: The exponentiation, e > 0
329  **/
330  inline UnivariateRationalFunction<UnivariatePolynomialOverField,Field> operator^ (long long int e) const {
332  res.setVariableName(num.variable());
333  if (isZero() || isOne() || e == 1)
334  res = *this;
335  else if (e == 2)
336  res = *this * *this;
337  else if (e > 2) {
339  res.one();
340
341  while (e) {
342  if (e % 2) { res *= x; }
343  x = x * x;
344  e >>= 1;
345  }
346  }
347  else if (!e)
348  res.one();
349  else {
350  res = *this;
351  }
352  return res;
353  }
354
355  /**
357  * replace xor operation by exponentiation
358  *
359  * @param e: The exponentiation, e > 0
360  **/
362  *this = *this ^ e;
363  return *this;
364  }
365
367  if (num.isZero()) {
368  std::cout << "BPAS error: division by zero from UnivariateRationalFunction<UnivariatePolynomialOverField,Field> inverse()" << std::endl;
369  exit(1);
370  }
372  r.normalize();
373  return r;
374  }
375
378  return (r *= b);
379  }
381  if (num.variable()!= b.num.variable()) {
382  std::cout << "BPAS: error, trying to multiply between RationalFunction[" << num.variable() << "] and RationalFunction[" << b.num.variable() << "]." << std::endl;
383  exit(1);
384  }
385
386  UnivariatePolynomialOverField g1, g2, r;
387  g1 = num.gcd(b.den);
388  g2 = den.gcd(b.num);
389  num /= g1;
390  r = b.num / g2;
391  num *= r;
392  den /= g2;
393  r = b.den / g1;
394  den *= r;
395  normalize();
396  return *this;
397  }
398
401  return (r /= b);
402  }
404  if (b.isZero()) {
405  std::cout << "BPAS error: division by zero from UnivariateRationalFunction<UnivariatePolynomialOverField,Field> operator/=" << std::endl;
406  exit(1);
407  }
408
409  if (num.variable()!= b.num.variable()) {
410  std::cout << "BPAS: error, trying to divide between RationalFunction[" << num.variable() << "] and RationalFunction[" << b.num.variable() << "]." << std::endl;
411  exit(1);
412  }
414  *this *= e;
415  return *this;
416  }
417
418  inline void canonicalize() {
419  UnivariatePolynomialOverField temp;
420  temp = num.gcd(den);
421  num /= temp;
422  den /= temp;
425  //UnivariateRationalFunction<UnivariatePolynomialOverField,Field> ret;
426  //ret = this->unitCanonical();
427  //*this = ret;
428  }
429
430  inline bool isZero() const {
431  return num.isZero();
432  }
433  inline void zero() {
434  num.zero();
435  den.one();
436  }
437  inline bool isOne() const {
438  return num.isOne() && den.isOne();
439  }
440  inline void one() {
441  num.one();
442  den.one();
443  }
444  inline bool isNegativeOne() const {
445  return (num.isNegativeOne() && den.isOne()) || (num.isOne() && den.isNegativeOne());
446  }
447  inline void negativeOne() {
448  num.negativeOne();
449  den.one();
450  }
451  inline int isConstant() const {
452  return num.isConstant() && den.isConstant();
453  }
454
456  UnivariatePolynomialOverField du,dc,dv,g,temp;
457  g = num.gcd(den);
458  temp = den/g;
459  dc = temp.unitCanonical(&du,&dv);
460  temp = du*(num/g);
462  ret.num = temp;
463  ret.den = dc;
464  if (u != NULL || v!= NULL) {
466  temp2.one();
467  if (u != NULL) {
468  *u = temp2;
469  }
470  if (v != NULL) {
471  *v = temp2;
472  }
473  }
474
475  return ret;
476  }
477
478  /**
480  *
481  * @param b: A rational function
482  **/
484  if (this != &b) {
485  num = b.num;
486  den = b.den;
487  characteristic = b.characteristic;
488  }
489  return *this;
490  }
491
493  std::cerr << "UnivariateRationalFunction::convertToExpressionTree NOT YET IMPLEMENTED" << std::endl;
494  //TODO
495  return ExpressionTree();
496  }
497
498  /**
499  * Overload stream operator <<
500  *
501  * @param out: Stream object
502  * @param b: The rational function
503  **/
504  void print(std::ostream& ostream) const {
505  ostream << "(" << num << ")/(" << den << ")";
506  }
507
508  /*UnivariateRationalFunction<UnivariatePolynomialOverField, Field> unitCanonical(UnivariateRationalFunction<UnivariatePolynomialOverField,Field>* u = NULL,
509  UnivariateRationalFunction<UnivariatePolynomialOverField,Field>* v = NULL) const {
510  UnivariatePolynomialOverField temp;
511  temp = num.gcd(den);
512  UnivariatePolynomialOverField tempDen = den / temp;
514  temp *= lc;
515  if (u != NULL) {
516  //TODO
517  *u = UnivariateRationalFunction<UnivariatePolynomialOverField,Field>();
518  }
519  if (v != NULL) {
520  *v = UnivariateRationalFunction<UnivariatePolynomialOverField,Field>(temp, temp);
521  }
522
523  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> ret = *this;
524  ret.canonicalize();
525  return ret;
526  }*/
527
528
529  /** BPASGCDDomain, BPASEuclideanDomain, BPASField virtual methods **/
530
531  /**
532  * Get this GCD of *this and b.
533  */
535  std::cerr << "UnivariateRationalFunction::gcd NOT YET IMPLEMENTED" << std::endl;
536  //TODO
537  return *this;
538  }
539
540  /**
541  * Compute squarefree factorization of *this
542  */
544  std::cerr << "UnivariateRationalFunction::squareFree NOT YET IMPLEMENTED" << std::endl;
545  //TODO
546  std::vector<UnivariateRationalFunction> ret;
547  ret.push_back(*this);
548  return ret;
549  }
550
552  std::cerr << "UnivariateRationalFunction::euclideanSize NOT YET IMPLEMENTED" << std::endl;
553  //TODO
554  return *this;
555
556  }
557
560  std::cerr << "UnivariateRationalFunction::euclideanDivision NOT YET IMPLEMENTED" << std::endl;
561  //TODO
562  return *this;
563
564  }
565
567  return (*this / b);
568  }
569
572  }
573
577  std::cerr << "UnivariateRationalFunction::extendedEuclidean NOT YET IMPLEMENTED" << std::endl;
578  //TODO
579  return *this;
580  }
581
584  return ret;
585  }
586
588  *this = this->remainder(b);
589  }
590
591
593  std::vector<UnivariatePolynomialOverField> gg,hh;
595
596  temp.setVariableName(num.variable());
597  h->setVariableName(num.variable());
598  g->clear();
599  _hermiteReduce<UnivariatePolynomialOverField,Field>(num,den,&gg,&hh);
600  int i(0);
601  while (i<gg.size()) {
602  temp.set(gg.at(i),gg.at(i+1));
603  g->push_back(temp);
604  i += 2;
605  }
606  temp.set(hh.at(0),hh.at(1));
607  *h = temp;
608  }
609
610  void integrateRationalFunctionLogPart(std::vector< SparseUnivariatePolynomial<UnivariatePolynomialOverField> > *S, std::vector<UnivariatePolynomialOverField> *U) {
611  U->clear();
612  S->clear();
613
614  _integrateRationalFunctionLogPart<UnivariatePolynomialOverField,Field>(S,U,num,den,PROFILING);
615  }
616
617  void differentiate() {
618  /* Destructive rational function differentiation */
619  // input a/d
620  // output (a'*d-a*d')/d^2 = a'*e-a*f/d*e; e=d/g, f=d'/g, g=gcd(d,d')
621
622  UnivariatePolynomialOverField D(den); // D = d
623  UnivariatePolynomialOverField dD(den);
624  UnivariatePolynomialOverField temp;
625  //std::cout << "here?" << std::endl;
626  dD.differentiate(1); // dD = d'
627  //std::cout << "dD = " << dD << std::endl;
628  temp = D.gcd(dD); // temp = g
629  //std::cout << "temp = " << temp << std::endl;
630  D /= temp; // D = e
631  //std::cout << "D = " << D << std::endl;
632  dD /= temp; // dD = f
633  //std::cout << "dD = " << dD << std::endl;
634  temp = -num; // temp = -a
635  //std::cout << "temp = " << temp << std::endl;
636  temp *= dD; // temp = -a*f
637  //std::cout << "temp = " << temp << std::endl;
638  dD = num; // dD = a
639  //std::cout << "dD = " << dD << std::endl;
640  dD.differentiate(1); // dD = a'
641  //std::cout << "dD = " << dD << std::endl;
642  dD *= D; // dD = a'*e
643  //std::cout << "dD = " << dD << std::endl;
644  temp += dD; // temp = a'*e-a*f
645  //std::cout << "temp = " << temp << std::endl;
646  D *= den; // D = d*e
647  //std::cout << "D = " << D << std::endl;
648
649  //std::cout << "here?" << std::endl;
650  num = temp;
651  den = D;
652  canonicalize();
653  }
654
655  void integrate(UnivariatePolynomialOverField *P, std::vector< UnivariateRationalFunction<UnivariatePolynomialOverField,Field> > *g, std::vector<UnivariatePolynomialOverField> *U, std::vector< SparseUnivariatePolynomial<UnivariatePolynomialOverField> > *S) {
656  g->clear();
657  U->clear();
658  S->clear();
659  std::vector<UnivariatePolynomialOverField> G;
661  temp.setVariableName(num.variable());
662
663  // Profiling variables
664  unsigned long long start;
665  float elapsed = 0;
666
667  if (PROFILING){
668  std::cout << "integrate" << std::endl;
669  std::cout << "--------------------------------------" << std::endl;
670  startTimer(&start);
671  }
672
673  _integrateRationalFunction<UnivariatePolynomialOverField,Field>(num,den,P,&G,U,S,PROFILING);
674
675  if (PROFILING){
676  stopTimer(&start,&elapsed);
677  std::cout << "--------------------------------------" << std::endl;
678  std::cout << "integrate runtime: " << elapsed << " s" << std::endl;
679  }
680
681  int i(0);
682  while (i<G.size()) {
683  temp.set(G.at(i),G.at(i+1));
684  g->push_back(temp);
685  i += 2;
686  }
687
688  }
689
690  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) {
691  P->zero();
692  g->clear();
693  lg->clear();
694  Lg->clear();
695  atn->clear();
696  Atn->clear();
697  std::vector<UnivariatePolynomialOverField> G;
699  temp.setVariableName(num.variable());
700
701  std::cout << "[realSymbolicNumericIntegrate (snInt): Symbolic-Numeric Integration with BPAS and MPSolve]" << std::endl;
702  std::cout << "[Integration method: Hermite reduction, LRT integration]" << std::endl;
703  std::cout << "Starting..." << std::endl;
704
705  // Profiling variables
706  unsigned long long start;
707  float elapsed = 0;
708
709  if (PROFILING){
710  std::cout << "--------------------------------------" << std::endl;
711  startTimer(&start);
712  }
713
714  _realSNIntegrate<UnivariatePolynomialOverField,Field>(num,den,P,&G,lg,Lg,atn,Atn,prec,PROFILING,PFD_LOG_PART,ERROR_ANALYSIS);
715
716  if (PROFILING){
717  stopTimer(&start,&elapsed);
718  std::cout << "--------------------------------------" << std::endl;
719  std::cout << "realSymbolicNumericIntegrate runtime: " << elapsed << " s" << std::endl;
720  std::fstream fs;
721  fs.open("perftiming.txt",std::fstream::in | std::fstream::out | std::fstream::app);
722  fs << elapsed << std::endl;
723  fs.close();
724  }
725
726  int i(0);
727  while (i<G.size()) {
728  temp.set(G.at(i),G.at(i+1));
729  g->push_back(temp);
730  i += 2;
731  }
732
733  }
734
735  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> *Atn1, std::vector<UnivariatePolynomialOverField> *Atn2, int prec) {
736  P->zero();
737  g->clear();
738  lg->clear();
739  Lg->clear();
740  atn->clear();
741  Atn1->clear();
742  Atn2->clear();
743  std::vector<UnivariatePolynomialOverField> G;
745  temp.setVariableName(num.variable());
746
747  std::cout << "[realSymbolicNumericIntegrate (snInt): Symbolic-Numeric Integration with BPAS and MPSolve]" << std::endl;
748  std::cout << "[Integration method: Hermite reduction, LRT integration]" << std::endl;
749  std::cout << "Starting..." << std::endl;
750
751  // Profiling variables
752  unsigned long long start;
753  float elapsed = 0;
754
755  if (PROFILING){
756  std::cout << "--------------------------------------" << std::endl;
757  startTimer(&start);
758  }
759
760  _realSNIntegrate<UnivariatePolynomialOverField,Field>(num,den,P,&G,lg,Lg,atn,Atn1,Atn2,prec,PROFILING,PFD_LOG_PART,ERROR_ANALYSIS);
761
762  if (PROFILING){
763  stopTimer(&start,&elapsed);
764  std::cout << "--------------------------------------" << std::endl;
765  std::cout << "realSymbolicNumericIntegrate runtime: " << elapsed << " s" << std::endl;
766  std::fstream fs;
767  fs.open("perftiming.txt",std::fstream::in | std::fstream::out | std::fstream::app);
768  fs << elapsed << std::endl;
769  fs.close();
770  }
771
772  int i(0);
773  while (i<G.size()) {
774  temp.set(G.at(i),G.at(i+1));
775  g->push_back(temp);
776  i += 2;
777  }
778
779  }
780
781  void realSymbolicNumericIntegratePFD(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) {
782  P->zero();
783  g->clear();
784  lg->clear();
785  Lg->clear();
786  atn->clear();
787  Atn->clear();
788  std::vector<UnivariatePolynomialOverField> G;
790  temp.setVariableName(num.variable());
791
792  std::cout << "[realSymbolicNumericIntegratePFD (snIntPFD): Symbolic-Numeric Integration with BPAS and MPSolve]" << std::endl;
793  std::cout << "[Integration method: Hermite reduction, PFD integration]" << std::endl;
794  std::cout << "Starting..." << std::endl;
795
796  // Profiling variables
797  unsigned long long start;
798  float elapsed = 0;
799
800  if (PROFILING){
801  std::cout << "--------------------------------------" << std::endl;
802  startTimer(&start);
803  }
804
805  _realSNIntegratePFD<UnivariatePolynomialOverField,Field>(num,den,P,&G,lg,Lg,atn,Atn,prec,PROFILING,PFD_LOG_PART,ERROR_ANALYSIS);
806
807  if (PROFILING){
808  stopTimer(&start,&elapsed);
809  std::cout << "--------------------------------------" << std::endl;
810  std::cout << "realSymbolicNumericIntegratePFD runtime: " << elapsed << " s" << std::endl;
811  std::fstream fs;
812  fs.open("perftiming.txt",std::fstream::in | std::fstream::out | std::fstream::app);
813  fs << elapsed << std::endl;
814  fs.close();
815  }
816
817  int i(0);
818  while (i<G.size()) {
819  temp.set(G.at(i),G.at(i+1));
820  g->push_back(temp);
821  i += 2;
822  }
823
824  }
825
826  void realSymbolicNumericIntegrateSimplePFD(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) {
827  P->zero();
828  g->clear();
829  lg->clear();
830  Lg->clear();
831  atn->clear();
832  Atn->clear();
833  std::vector<UnivariatePolynomialOverField> G;
835  temp.setVariableName(num.variable());
836
837  std::cout << "[realSymbolicNumericIntegratePFD (snIntPFD): Symbolic-Numeric Integration with BPAS and MPSolve]" << std::endl;
838  std::cout << "[Integration method: Hermite reduction, PFD integration]" << std::endl;
839  std::cout << "Starting..." << std::endl;
840
841  // Profiling variables
842  unsigned long long start;
843  float elapsed = 0;
844
845  if (PROFILING){
846  std::cout << "--------------------------------------" << std::endl;
847  startTimer(&start);
848  }
849
850  _realSNIntegrateSimplePFD<UnivariatePolynomialOverField,Field>(num,den,P,&G,lg,Lg,atn,Atn,prec,PROFILING,PFD_LOG_PART,ERROR_ANALYSIS);
851
852  if (PROFILING){
853  stopTimer(&start,&elapsed);
854  std::cout << "--------------------------------------" << std::endl;
855  std::cout << "realSymbolicNumericIntegratePFD runtime: " << elapsed << " s" << std::endl;
856  std::fstream fs;
857  fs.open("perftiming.txt",std::fstream::in | std::fstream::out | std::fstream::app);
858  fs << elapsed << std::endl;
859  fs.close();
860  }
861
862  int i(0);
863  while (i<G.size()) {
864  temp.set(G.at(i),G.at(i+1));
865  g->push_back(temp);
866  i += 2;
867  }
868
869  }
870
871  /*void realSymbolicNumericIntegratePFD(UnivariatePolynomialOverField *P, std::vector< UnivariateRationalFunction<UnivariatePolynomialOverField,Field> > *g, std::vector<Field> *lg, std::vector<UnivariatePolynomialOverField> *Lg, std::vector<Field> *atn, std::vector<UnivariatePolynomialOverField> *Atn1, std::vector<UnivariatePolynomialOverField> *Atn2, int prec) {
872  P->zero();
873  g->clear();
874  lg->clear();
875  Lg->clear();
876  atn->clear();
877  Atn1->clear();
878  Atn2->clear();
879  std::vector<UnivariatePolynomialOverField> G;
880  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> temp;
881  temp.setVariableName(num.variable());
882
883  std::cout << "[realSymbolicNumericIntegratePFD (snIntPFD): Symbolic-Numeric Integration with BPAS and MPSolve]" << std::endl;
884  std::cout << "[Integration method: Hermite reduction, PFD integration]" << std::endl;
885  std::cout << "Starting..." << std::endl;
886
887  // Profiling variables
888  unsigned long long start;
889  float elapsed = 0;
890
891  if (PROFILING){
892  std::cout << "--------------------------------------" << std::endl;
893  startTimer(&start);
894  }
895
896  _realSNIntegratePFD<UnivariatePolynomialOverField,Field>(num,den,P,&G,lg,Lg,atn,Atn1,Atn2,prec,PROFILING,PFD_LOG_PART,ERROR_ANALYSIS);
897
898  if (PROFILING){
899  stopTimer(&start,&elapsed);
900  std::cout << "--------------------------------------" << std::endl;
901  std::cout << "realSymbolicNumericIntegratePFD runtime: " << elapsed << " s" << std::endl;
902  }
903
904  int i(0);
905  while (i<G.size()) {
906  temp.set(G.at(i),G.at(i+1));
907  g->push_back(temp);
908  i += 2;
909  }
910
911  }*/
912
913  void printIntegral(UnivariatePolynomialOverField &P, std::vector< UnivariateRationalFunction<UnivariatePolynomialOverField,Field> > &g, std::vector<UnivariatePolynomialOverField> &U, std::vector< SparseUnivariatePolynomial<UnivariatePolynomialOverField> > &S){
914  std::vector<UnivariatePolynomialOverField> G;
915  for (int i=0; i<g.size(); i++) {
916  G.push_back(g.at(i).num);
917  G.push_back(g.at(i).den);
918  }
919  _printFormalIntegral<UnivariatePolynomialOverField,Field>(num,den,P,G,U,S, false, floatingPointPrinting, false);
920  }
921
922  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){
923  std::vector<UnivariatePolynomialOverField> G;
924  for (int i=0; i<g.size(); i++) {
925  G.push_back(g.at(i).num);
926  G.push_back(g.at(i).den);
927  }
928  std::vector<UnivariatePolynomialOverField> empty;
929  _printIntegral<UnivariatePolynomialOverField,Field>(num,den,P,G,lg,Lg,atn,Atn,empty, false, floatingPointPrinting, false, outputFormatting);
930  }
931
932  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> &Atn1, std::vector<UnivariatePolynomialOverField> &Atn2){
933  std::vector<UnivariatePolynomialOverField> G;
934  for (int i=0; i<g.size(); i++) {
935  G.push_back(g.at(i).num);
936  G.push_back(g.at(i).den);
937  }
938  _printIntegral<UnivariatePolynomialOverField,Field>(num,den,P,G,lg,Lg,atn,Atn1,Atn2, false, floatingPointPrinting, false, outputFormatting);
939  }
940
941  void realSymbolicNumericIntegrate(int prec) {
942  UnivariatePolynomialOverField P;
943  std::vector< UnivariateRationalFunction<UnivariatePolynomialOverField,Field> > g;
944  std::vector<Field> lg, atn;
945  std::vector<UnivariatePolynomialOverField> Lg, Atn1, Atn2;
946
947  realSymbolicNumericIntegrate(&P, &g, &lg, &Lg, &atn, &Atn1, &Atn2, prec);
948  printIntegral(P, g, lg, Lg, atn, Atn1, Atn2);
949  }
950
951  void integrate() {
952  UnivariatePolynomialOverField P;
953  std::vector< UnivariateRationalFunction<UnivariatePolynomialOverField,Field> > g;
954  std::vector<UnivariatePolynomialOverField> U;
955  std::vector< SparseUnivariatePolynomial<UnivariatePolynomialOverField> > S;
956
957  integrate(&P, &g, &U, &S);
958  printIntegral(P, g, U, S);
959  }
960 };
961 #endif
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > remainder(const UnivariateRationalFunction< UnivariatePolynomialOverField, Field > &b) const
Get the remainder of *this and b.
Definition: urationalfunction.h:570
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > unitCanonical(UnivariateRationalFunction< UnivariatePolynomialOverField, Field > *u=NULL, UnivariateRationalFunction< UnivariatePolynomialOverField, Field > *v=NULL) const
Obtain the unit normal (a.k.a canonical associate) of an element.
Definition: urationalfunction.h:455
A univariate polynomial over an arbitrary BPASRing represented sparsely.
Definition: BigPrimeField.hpp:21
void canonicalize()
Canonicalize this fraction, reducing terms as needed.
Definition: urationalfunction.h:418
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > euclideanSize() const
Get the euclidean size of *this.
Definition: urationalfunction.h:551
An abstract class defining the interface of a rational function.
Definition: polynomial.h:186
An ExpressionTree encompasses various forms of data that can be expressed generically as a binary tre...
Definition: ExpressionTree.hpp:17
UnivariatePolynomialOverField numerator() const
Get the fraction&#39;s numerator.
Definition: urationalfunction.h:243
void print(std::ostream &ostream) const
Definition: urationalfunction.h:504
ExpressionTree convertToExpressionTree() const
Convert this to an expression tree.
Definition: urationalfunction.h:492
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > & operator%=(const UnivariateRationalFunction< UnivariatePolynomialOverField, Field > &b)
Assign *this to be the remainder of *this and b.
Definition: urationalfunction.h:587
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > operator%(const UnivariateRationalFunction< UnivariatePolynomialOverField, Field > &b) const
Get the remainder of *this and b;.
Definition: urationalfunction.h:582
void zero()
Make *this ring element zero.
Definition: urationalfunction.h:433
A simple data structure for encapsulating a collection of Factor elements.
Definition: Factors.hpp:95
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > gcd(const UnivariateRationalFunction< UnivariatePolynomialOverField, Field > &b) const
BPASGCDDomain, BPASEuclideanDomain, BPASField virtual methods.
Definition: urationalfunction.h:534
UnivariatePolynomialOverField denominator() const
Get the fraction&#39;s denominator.
Definition: urationalfunction.h:247
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > inverse() const
Get the inverse of *this.
Definition: urationalfunction.h:366
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > quotient(const UnivariateRationalFunction< UnivariatePolynomialOverField, Field > &b) const
Get the quotient of *this and b.
Definition: urationalfunction.h:566
void one()
Make *this ring element one.
Definition: urationalfunction.h:440
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > extendedEuclidean(const UnivariateRationalFunction< UnivariatePolynomialOverField, Field > &b, UnivariateRationalFunction< UnivariatePolynomialOverField, Field > *s=NULL, UnivariateRationalFunction< UnivariatePolynomialOverField, Field > *t=NULL) const
Perofrm the extended euclidean division on *this and b.
Definition: urationalfunction.h:574
bool isZero() const
Determine if *this ring element is zero, that is the additive identity.
Definition: urationalfunction.h:430
A univariate rational function templated by a unvariate polynomial over a field.
Definition: urationalfunction.h:72
An encapsulation of a mathematical symbol.
Definition: Symbol.hpp:23
bool isOne() const
Determine if *this ring element is one, that is the multiplication identity.
Definition: urationalfunction.h:437
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > euclideanDivision(const UnivariateRationalFunction< UnivariatePolynomialOverField, Field > &b, UnivariateRationalFunction< UnivariatePolynomialOverField, Field > *q=NULL) const
Perform the eucldiean division of *this and b.
Definition: urationalfunction.h:558
Factors< UnivariateRationalFunction > squareFree() const
Compute squarefree factorization of *this.
Definition: urationalfunction.h:543