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 () {
85  num /= den.leadingCoefficient();
86  den /= den.leadingCoefficient();
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  **/
98  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> () {
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  **/
115  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> (const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b) : den(b.den), num(b.num), PROFILING(b.PROFILING),
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  **/
150  ~UnivariateRationalFunction<UnivariatePolynomialOverField,Field> () {}
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 
204  inline bool isMatlabOutput() {
205  if (outputFormatting == "MATLAB_OUTPUT")
206  return true;
207  else
208  return false;
209  }
210 
211  inline void setMatlabOutput() {
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 
265  inline UnivariateRationalFunction<UnivariatePolynomialOverField,Field> operator+ (const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b) const {
266  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> r(*this);
267  return (r += b);
268  }
269  inline UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& operator+= (const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b) {
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 
292  inline UnivariateRationalFunction<UnivariatePolynomialOverField,Field> operator- (const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b) const {
293  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> r(*this);
294  return (r -= b);
295  }
296  inline UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& operator-= (const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b) {
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 
319  inline UnivariateRationalFunction<UnivariatePolynomialOverField,Field> operator- () const {
320  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> r (-num, den);
321  return r;
322  }
323 
324  /**
325  * Overload operator ^
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 {
331  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> res;
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) {
338  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> x (*this);
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  /**
356  * Overload operator ^=
357  * replace xor operation by exponentiation
358  *
359  * @param e: The exponentiation, e > 0
360  **/
361  inline UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& operator^= (long long int e) {
362  *this = *this ^ e;
363  return *this;
364  }
365 
366  inline UnivariateRationalFunction<UnivariatePolynomialOverField,Field> inverse() const {
367  if (num.isZero()) {
368  std::cout << "BPAS error: division by zero from UnivariateRationalFunction<UnivariatePolynomialOverField,Field> inverse()" << std::endl;
369  exit(1);
370  }
371  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> r(den, num);
372  r.normalize();
373  return r;
374  }
375 
376  inline UnivariateRationalFunction<UnivariatePolynomialOverField,Field> operator* (const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b) const {
377  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> r(*this);
378  return (r *= b);
379  }
380  inline UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& operator*= (const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b) {
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 
399  inline UnivariateRationalFunction<UnivariatePolynomialOverField,Field> operator/ (const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b) const {
400  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> r(*this);
401  return (r /= b);
402  }
403  inline UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& operator/= (const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b) {
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  }
413  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> e (b.den, b.num);
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;
423  num /= den.leadingCoefficient();
424  den /= den.leadingCoefficient();
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 
455  inline UnivariateRationalFunction<UnivariatePolynomialOverField,Field> unitCanonical(UnivariateRationalFunction<UnivariatePolynomialOverField,Field>* u = NULL, UnivariateRationalFunction<UnivariatePolynomialOverField,Field>* v = NULL) const {
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);
461  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> ret;
462  ret.num = temp;
463  ret.den = dc;
464  if (u != NULL || v!= NULL) {
465  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> temp2;
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  /**
479  * Overload operator =
480  *
481  * @param b: A rational function
482  **/
483  inline UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& operator= (const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b) {
484  if (this != &b) {
485  num = b.num;
486  den = b.den;
487  characteristic = b.characteristic;
488  }
489  return *this;
490  }
491 
492  ExpressionTree convertToExpressionTree() const {
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;
513  Field lc = tempDen.leadingCoefficient();
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  */
534  UnivariateRationalFunction<UnivariatePolynomialOverField, Field> gcd(const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b) const {
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  */
543  inline Factors<UnivariateRationalFunction> squareFree() const {
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 
551  UnivariateRationalFunction<UnivariatePolynomialOverField, Field> euclideanSize() const {
552  std::cerr << "UnivariateRationalFunction::euclideanSize NOT YET IMPLEMENTED" << std::endl;
553  //TODO
554  return *this;
555 
556  }
557 
558  UnivariateRationalFunction<UnivariatePolynomialOverField, Field> euclideanDivision(const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b,
559  UnivariateRationalFunction<UnivariatePolynomialOverField,Field>* q = NULL) const {
560  std::cerr << "UnivariateRationalFunction::euclideanDivision NOT YET IMPLEMENTED" << std::endl;
561  //TODO
562  return *this;
563 
564  }
565 
566  UnivariateRationalFunction<UnivariatePolynomialOverField, Field> quotient(const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b) const {
567  return (*this / b);
568  }
569 
570  UnivariateRationalFunction<UnivariatePolynomialOverField, Field> remainder(const UnivariateRationalFunction<UnivariatePolynomialOverField,Field>& b) const {
571  return UnivariateRationalFunction<UnivariatePolynomialOverField, Field>(0,1);
572  }
573 
574  UnivariateRationalFunction<UnivariatePolynomialOverField, Field> extendedEuclidean(const UnivariateRationalFunction<UnivariatePolynomialOverField, Field>& b,
575  UnivariateRationalFunction<UnivariatePolynomialOverField, Field>* s = NULL,
576  UnivariateRationalFunction<UnivariatePolynomialOverField, Field>* t = NULL) const {
577  std::cerr << "UnivariateRationalFunction::extendedEuclidean NOT YET IMPLEMENTED" << std::endl;
578  //TODO
579  return *this;
580  }
581 
582  inline UnivariateRationalFunction<UnivariatePolynomialOverField, Field> operator%(const UnivariateRationalFunction<UnivariatePolynomialOverField, Field>& b) const {
583  UnivariateRationalFunction<UnivariatePolynomialOverField, Field> ret = this->remainder(b);
584  return ret;
585  }
586 
587  inline UnivariateRationalFunction<UnivariatePolynomialOverField, Field>& operator%=(const UnivariateRationalFunction<UnivariatePolynomialOverField, Field>& b) {
588  *this = this->remainder(b);
589  }
590 
591 
592  void hermiteReduce(std::vector< UnivariateRationalFunction<UnivariatePolynomialOverField,Field> > *g, UnivariateRationalFunction<UnivariatePolynomialOverField,Field> *h) {
593  std::vector<UnivariatePolynomialOverField> gg,hh;
594  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> temp;
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;
660  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> temp;
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;
698  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> temp;
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;
744  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> temp;
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;
789  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> temp;
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;
834  UnivariateRationalFunction<UnivariatePolynomialOverField,Field> temp;
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::remainder
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > remainder(const UnivariateRationalFunction< UnivariatePolynomialOverField, Field > &b) const
Get the remainder of *this and b.
Definition: urationalfunction.h:570
UnivariateRationalFunction::unitCanonical
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
SparseUnivariatePolynomial
A univariate polynomial over an arbitrary BPASRing represented sparsely.
Definition: BigPrimeField.hpp:21
UnivariateRationalFunction::canonicalize
void canonicalize()
Canonicalize this fraction, reducing terms as needed.
Definition: urationalfunction.h:418
UnivariateRationalFunction::euclideanSize
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > euclideanSize() const
Get the euclidean size of *this.
Definition: urationalfunction.h:551
BPASRationalFunction
An abstract class defining the interface of a rational function.
Definition: polynomial.h:186
ExpressionTree
An ExpressionTree encompasses various forms of data that can be expressed generically as a binary tre...
Definition: ExpressionTree.hpp:17
UnivariateRationalFunction::numerator
UnivariatePolynomialOverField numerator() const
Get the fraction&#39;s numerator.
Definition: urationalfunction.h:243
UnivariateRationalFunction::print
void print(std::ostream &ostream) const
Overload stream operator <<.
Definition: urationalfunction.h:504
UnivariateRationalFunction::convertToExpressionTree
ExpressionTree convertToExpressionTree() const
Convert this to an expression tree.
Definition: urationalfunction.h:492
UnivariateRationalFunction::operator%=
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > & operator%=(const UnivariateRationalFunction< UnivariatePolynomialOverField, Field > &b)
Assign *this to be the remainder of *this and b.
Definition: urationalfunction.h:587
UnivariateRationalFunction::operator%
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > operator%(const UnivariateRationalFunction< UnivariatePolynomialOverField, Field > &b) const
Get the remainder of *this and b;.
Definition: urationalfunction.h:582
UnivariateRationalFunction::zero
void zero()
Make *this ring element zero.
Definition: urationalfunction.h:433
Factors
A simple data structure for encapsulating a collection of Factor elements.
Definition: Factors.hpp:95
UnivariateRationalFunction::gcd
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > gcd(const UnivariateRationalFunction< UnivariatePolynomialOverField, Field > &b) const
BPASGCDDomain, BPASEuclideanDomain, BPASField virtual methods.
Definition: urationalfunction.h:534
UnivariateRationalFunction::denominator
UnivariatePolynomialOverField denominator() const
Get the fraction&#39;s denominator.
Definition: urationalfunction.h:247
UnivariateRationalFunction::inverse
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > inverse() const
Get the inverse of *this.
Definition: urationalfunction.h:366
UnivariateRationalFunction::quotient
UnivariateRationalFunction< UnivariatePolynomialOverField, Field > quotient(const UnivariateRationalFunction< UnivariatePolynomialOverField, Field > &b) const
Get the quotient of *this and b.
Definition: urationalfunction.h:566
UnivariateRationalFunction::one
void one()
Make *this ring element one.
Definition: urationalfunction.h:440
UnivariateRationalFunction::extendedEuclidean
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
UnivariateRationalFunction::isZero
bool isZero() const
Determine if *this ring element is zero, that is the additive identity.
Definition: urationalfunction.h:430
UnivariateRationalFunction
A univariate rational function templated by a unvariate polynomial over a field.
Definition: urationalfunction.h:72
Symbol
An encapsulation of a mathematical symbol.
Definition: Symbol.hpp:23
UnivariateRationalFunction::isOne
bool isOne() const
Determine if *this ring element is one, that is the multiplication identity.
Definition: urationalfunction.h:437
UnivariateRationalFunction::euclideanDivision
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
UnivariateRationalFunction::squareFree
Factors< UnivariateRationalFunction > squareFree() const
Compute squarefree factorization of *this.
Definition: urationalfunction.h:543
Generated on Thu Jul 11 2019 16:05:19 for Basic Polynomial Algebra Subprograms (BPAS) by   doxygen 1.8.13