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