All Classes Functions Friends
mpolynomial.h
1 #ifndef _MPOLYNOMIAL_H_
2 #define _MPOLYNOMIAL_H_
3 
4 /**
5  * Data Structure for multivariate polynomial
6  * stored in a sparse case
7  **/
8 
9 #include "../polynomial.h"
10 
11 template <class Ring>
13  private:
14  int var; // Number of variables
15  std::string* names; // Variable names
16  // names[0] = "1" or "9" (user specified)
17  // names[1] < names[2] < ...
18  std::vector< MultivariateTerm<Ring> > terms; // Each term of the polynomial
19 
20  /* Return the first non-zero variate in the ascending order */
21  inline int firstNonZero(int* d) {
22  for (int i = 0; i < var; i++) {
23  if (d[i])
24  return i;
25  }
26  return -1;
27  }
28  /* Return the first non-equal variate in the descending order */
29  inline int checkDegs(int* d, int* b) {
30  for (int i = var-1; i > -1; --i) {
31  if (d[i] != b[i])
32  return i;
33  }
34  return -1;
35  }
36  inline bool isOrderedRing(SparseMultivariatePolynomial<Ring>& b, std::vector<int>& xs) {
37  if (names[0] != b.names[0]) { return 0; }
38  if (names[0] == "1") {
39  int k = 1;
40  for (; k <= var && k <= b.var; ++k) {
41  xs.push_back(k);
42  xs.push_back(k);
43  }
44  for (int i = k; i <= var; ++i) {
45  xs.push_back(i);
46  xs.push_back(0);
47  }
48  for (int i = k; i <= b.var; ++i) {
49  xs.push_back(0);
50  xs.push_back(i);
51  }
52  return 1;
53  }
54  if (!var) {
55  for (int i = 1; i <= b.var; ++i) {
56  xs.push_back(0);
57  xs.push_back(i);
58  }
59  return 1;
60  }
61  if (!b.var) {
62  for (int i = 1; i <= var; ++i) {
63  xs.push_back(i);
64  xs.push_back(0);
65  }
66  return 1;
67  }
68 
69  bool isFound = 0;
70  int* pos = new int[var];
71  for (int i = 1; i <= var; ++i) {
72  isFound = 0;
73  for (int j = 1; j <= b.var; ++j) {
74  if (names[i] == b.names[j]) {
75  pos[i-1] = j;
76  isFound = 1;
77  break;
78  }
79  }
80  if (!isFound) { pos[i-1] = 0; }
81  }
82 
83  isFound = 0;
84  int ak = 1, bk = pos[0];
85  if (pos[0]) {
86  for(int j = 1; j < pos[0]; ++j) {
87  xs.push_back(0);
88  xs.push_back(j);
89  }
90  xs.push_back(1);
91  xs.push_back(pos[0]);
92  }
93  for (int i = 1; i < var; ++i) {
94  if (pos[i] > bk) {
95  while (pos[i] - bk > 1 && (i > ak || !bk)) {
96  if (bk) { ak++; }
97  bk++;
98  if (names[ak] < b.names[bk]) {
99  xs.push_back(ak);
100  xs.push_back(0);
101  xs.push_back(0);
102  xs.push_back(bk);
103  }
104  else if (names[ak] > b.names[bk]) {
105  xs.push_back(0);
106  xs.push_back(bk);
107  xs.push_back(ak);
108  xs.push_back(0);
109  }
110  }
111  for(int j = bk+1; j < pos[i]; ++j) {
112  xs.push_back(0);
113  xs.push_back(j);
114  }
115  for(int j = (!bk)? ak : ak+1; j <= i; ++j) {
116  xs.push_back(j);
117  xs.push_back(0);
118  }
119  xs.push_back(i+1);
120  xs.push_back(pos[i]);
121  bk = pos[i];
122  ak = i + 1;
123  }
124  else if (pos[i] && pos[i] < bk) {
125  isFound = 1;
126  break;
127  }
128  }
129  if (!isFound) {
130  for (int i = (bk)? ak+1 : ak, j = bk+1; i <= var && j <= b.var; ++i, ++j) {
131  if (names[i] < b.names[j]) {
132  xs.push_back(i);
133  xs.push_back(0);
134  xs.push_back(0);
135  xs.push_back(j);
136  }
137  else if (names[i] > b.names[j]) {
138  xs.push_back(0);
139  xs.push_back(j);
140  xs.push_back(i);
141  xs.push_back(0);
142  }
143  ak = i, bk = j;
144  }
145  for (int i = ak+1; i <= var; ++i) {
146  xs.push_back(i);
147  xs.push_back(0);
148  }
149  for (int j = bk+1; j <= b.var; ++j) {
150  xs.push_back(0);
151  xs.push_back(j);
152  }
153  }
154 
155  delete [] pos;
156 
157  if (isFound) { xs.clear(); return 0; }
158  else { return 1; }
159  }
160 
161  /**
162  * Compare two term degrees in the same Ring
163  *
164  * Return:
165  * 1: a > b
166  * 0: a = b
167  * -1: a < b
168  **/
169  inline int compareTermDegs(MultivariateTerm<Ring>& a, MultivariateTerm<Ring>& b) {
170  for (int i = a.v-1; i > -1; --i) {
171  if (a.degs[i] < b.degs[i])
172  return -1;
173  else if (a.degs[i] > b.degs[i])
174  return 1;
175  }
176  return 0;
177  }
178  inline int compareTermDegs(MultivariateTerm<Ring>& a, MultivariateTerm<Ring>& b, std::vector<int> xs) {
179  int n = xs.size() / 2;
180  for (int i = n-1; i > -1; --i) {
181  if (xs[2*i] && xs[2*i+1]) {
182  if (a.degs[xs[2*i]-1] < b.degs[xs[2*i+1]-1])
183  return -1;
184  else if (a.degs[xs[2*i]-1] > b.degs[xs[2*i+1]-1])
185  return 1;
186  }
187  else if (!xs[2*i+1] && xs[2*i] && a.degs[xs[2*i]-1])
188  return 1;
189  else if (!xs[2*i] && xs[2*i+1] && b.degs[xs[2*i+1]-1])
190  return -1;
191  }
192  return 0;
193  }
194 
195  inline void basicOp(Ring c, bool isOp) {
196  if (!c.isZero()) {
197  if (terms.size()) {
198  bool isIt = 1;
199  for (int i = 0; i < var; ++i) {
200  if (terms[0].degs[i]) {
201  isIt = 0;
202  break;
203  }
204  }
205  if (isIt) {
206  if (isOp)
207  terms[0].coef -= c;
208  else
209  terms[0].coef += c;
210  if (terms[0].coef == 0)
211  terms.erase(terms.begin());
212  }
213  else {
214  MultivariateTerm<Ring> a;
215  if (isOp)
216  a.coef = -c;
217  else
218  a.coef = c;
219  a.v = var;
220  a.degs = new int[var];
221  for (int i = 0; i < var; ++i)
222  a.degs[i] = 0;
223  terms.insert(terms.begin(), a);
224  }
225  }
226  else {
227  MultivariateTerm<Ring> a;
228  if (isOp)
229  a.coef = -c;
230  else
231  a.coef = c;
232  a.v = var;
233  a.degs = new int[var];
234  for (int i = 0; i < var; ++i)
235  a.degs[i] = 0;
236  terms.push_back(a);
237  }
238  }
239  }
240 
241  inline void pomopo(MultivariateTerm<Ring>& t, SparseMultivariatePolynomial<Ring>& b, std::vector<int> xs) {
242  for (int i = 0; i < b.terms.size(); ++i) {
243  MultivariateTerm<Ring> a;
244  a.coef = t.coef * b.terms[i].coef;
245  a.v = var;
246  a.degs = new int[var];
247  for (int j = 0; j < var; ++j) {
248  a.degs[j] = 0;
249  if (xs[2*j])
250  a.degs[j] += t.degs[xs[2*j]-1];
251  if (xs[2*j+1])
252  a.degs[j] += b.terms[i].degs[xs[2*j+1]-1];
253  }
254  int k = 0;
255  while (k <= terms.size()) {
256  if (k == terms.size()) {
257  terms.push_back(a);
258  break;
259  }
260 
261  int s = compareTermDegs(terms[k], a);
262  if (!s) {
263  terms[k].coef += a.coef;
264  if (terms[k].coef.isZero())
265  terms.erase(terms.begin()+k);
266  break;
267  }
268  else if (s > 0) {
269  terms.insert(terms.begin()+k, a);
270  break;
271  }
272  else { k++; }
273  }
274  }
275  }
276 
277  /* Is equal to another polynomial */
278  inline bool isEqual(SparseMultivariatePolynomial<Ring>& b) {
279  if (terms.size() != b.terms.size()) { return 0; }
280 
281  std::vector<int> xs;
282  bool isOrdered = isOrderedRing(b, xs);
283  if (!isOrdered) { return 0; }
284 
285  int v = xs.size() / 2;
286  for (int i = 0; i < terms.size(); ++i) {
287  if (terms[i].coef != b.terms[i].coef)
288  return 0;
289  for (int j = 0; j < v; ++j) {
290  if (xs[2*j] && xs[2*j+1] && (terms[i].degs[xs[2*j]-1] != b.terms[i].degs[xs[2*j+1]-1]))
291  return 0;
292  else if (!xs[2*j] && xs[2*j+1] && (b.terms[i].degs[xs[2*j+1]-1] != 0))
293  return 0;
294  else if (!xs[2*j+1] && xs[2*j] && (terms[i].degs[xs[2*j]-1] != 0))
295  return 0;
296  }
297  }
298 
299  return 1;
300  }
301 
302  public:
303  int characteristic;
304  static bool isPrimeField;
305  static bool isComplexField;
306  /**
307  * Construct a multivariate polynomial
308  *
309  * @param
310  **/
312  names = new std::string[1];
313  names[0] = "1";
314  Ring e;
315  characteristic = e.characteristic;
316  }
317  /**
318  * Construct a multivariate polynomial with number of terms and variates
319  *
320  * @param v: Number of variables
321  **/
323  var = v;
324  names = new std::string[var+1];
325  names[0] = "1";
326  for (int i = 1; i <= var; ++i) {
327  std::ostringstream convert;
328  convert << var - i + 1;
329  names[i] = "_";
330  names[i] += convert.str();
331  }
332  Ring e;
333  characteristic = e.characteristic;
334  }
335  /**
336  * Construct with a variable name
337  * such that f(x) = x
338  *
339  * @param x: The variable name
340  **/
342  var = 1;
343  names = new std::string[2];
344  names[0] = "9";
345  names[1] = x;
346  MultivariateTerm<Ring> t;
347  t.coef.one();
348  t.v = 1;
349  t.degs = new int[1];
350  t.degs[0] = 1;
351  terms.push_back(t);
352  Ring e;
353  characteristic = e.characteristic;
354  }
355  /**
356  * Copy Constructor
357  *
358  * @param b: A sparse multivariate polynomial
359  **/
361  names = new std::string[var+1];
362  std::copy(b.names, b.names+var+1, names);
363  Ring e;
364  characteristic = e.characteristic;
365  }
366  /**
367  * Construct from a SUP<SMQP> polynomial
368  *
369  * @param s: The SUP<SMQP> polynomial
370  **/
372  names = new std::string[1];
373  int d = s.degree();
374  for (int k = 0; k <= d; ++k) {
376  t.var = c.var + 1;
377  t.names = new std::string[t.var+1];
378  if (c.var)
379  std::copy(c.names, c.names+t.var, t.names);
380  else
381  t.names[0] = "9";
382  t.names[t.var] = s.variable();
383  for (int i = 0; i < c.terms.size(); ++i) {
384  MultivariateTerm<Ring> a;
385  a.coef = c.terms[i].coef;
386  a.v = t.var;
387  a.degs = new int[t.var];
388  for (int j = 0; j < c.var; ++j)
389  a.degs[j] = c.terms[i].degs[j];
390  a.degs[c.var] = k;
391  t.terms.push_back(a);
392  }
393  if (k) { *this = *this + t; }
394  else { *this = t; }
395  }
396  Ring e;
397  characteristic = e.characteristic;
398  }
399  /**
400  * Destroy the polynomial
401  *
402  * @param
403  **/
405  terms.clear();
406  delete [] names;
407  }
408 
409  /**
410  * Get number of variables
411  *
412  * @param
413  **/
414  inline int numberOfVariables() {
415  return var;
416  }
417 
418  /**
419  * Get number of non-zero terms
420  *
421  * @param
422  **/
423  inline int numberOfTerms() {
424  return terms.size();
425  }
426  /**
427  * Get the degree of a variable
428  *
429  * @param x: The variable name
430  **/
431  inline int degree(std::string x) {
432  int k = 0;
433  for (int i = 1; i <= var; ++i) {
434  if (names[i] == x) {
435  k = i;
436  break;
437  }
438  }
439  if (!k) { return 0; }
440  int d = 0;
441  for (int i = terms.size()-1; i > -1; --i) {
442  if (terms[i].degs[k] > d)
443  d = terms[i].degs[k];
444  }
445  return d;
446  }
447 
448  /**
449  * Get the leading variable
450  *
451  * @param
452  **/
453  inline std::string leadingVariable() {
454  return names[var];
455  }
456  /**
457  * Get the leading term's leading variable's degree
458  *
459  * @param
460  **/
461  inline int leadingVariableDegree() {
462  int n = terms.size();
463  if (n)
464  return terms[n-1].degs[var-1];
465  else
466  return 0;
467  }
468  /**
469  * Get the leading coefficient
470  *
471  * @param
472  **/
473  inline Ring leadingCoefficient() {
474  int n = terms.size();
475  if (n)
476  return terms[n].coef;
477  else
478  return 0;
479  }
480  /**
481  * Get the leading coefficient over a variable
482  *
483  * @param x: The name of the variable
484  * @param e: The leading exponent of the variable
485  **/
487  if (e == NULL)
488  e = new int;
489  *e = 0;
490 
491  if (isConstant()) {
493  r.var = var;
494  r.names = new std::string[var+1];
495  std::copy(names, names+var+1, r.names);
496  return r;
497  }
498 
499  int k = 0;
500  for (int i = 1; i <= var; ++i) {
501  if (names[i] == x) {
502  k = i;
503  break;
504  }
505  }
506  int v = var - 1;
508  for (int i = 0; i < var; ++i) {
509  if (i < k)
510  r.names[i] = names[i];
511  else
512  r.names[i] = names[i+1];
513  }
514  if (!k) { return r; }
515  k--;
516  for (int i = 0; i < terms.size(); ++i) {
517  if (terms[i].degs[k] >= *e) {
518  if (terms[i].degs[k] > *e) {
519  *e = terms[i].degs[k];
520  r.terms.clear();
521  }
522  MultivariateTerm<Ring> t;
523  t.coef = terms[i].coef;
524  t.v = v;
525  t.degs = new int[v];
526  for (int j = 0; j < v; ++j) {
527  if (j < k)
528  t.degs[j] = terms[i].degs[j];
529  else
530  t.degs[j] = terms[i].degs[j+1];
531  }
532  r.terms.push_back(t);
533  }
534  }
535  return r;
536  }
537 
538  /**
539  * Is zero polynomial
540  *
541  * @param
542  **/
543  inline bool isZero() {
544  return !terms.size();
545  }
546  /**
547  * Zero polynomial
548  *
549  * @param
550  **/
551  inline void zero() {
552  terms.clear();
553  }
554  /**
555  * Is polynomial a constant 1
556  *
557  * @param
558  **/
559  inline bool isOne() {
560  if (terms.size() == 1) {
561  for (int i = 0; i < var; ++i) {
562  if (terms[0].degs[i] != 0)
563  return 0;
564  }
565  if (terms[0].coef.isOne())
566  return 1;
567  }
568  return 0;
569  }
570  /**
571  * Set polynomial to 1
572  *
573  * @param
574  **/
575  inline void one() {
576  terms.clear();
577  MultivariateTerm<Ring> t;
578  t.coef.one();
579  t.v = var;
580  t.degs = new int[var];
581  for (int i = 0; i < var; ++i)
582  t.degs[i] = 0;
583  terms.push_back(t);
584  }
585  /**
586  * Is polynomial a constant -1
587  *
588  * @param
589  **/
590  inline bool isNegativeOne() {
591  if (terms.size() == 1) {
592  for (int i = 0; i < var; ++i) {
593  if (terms[0].degs[i] != 0)
594  return 0;
595  }
596  if (terms[0].coef.isNegativeOne())
597  return 1;
598  }
599  return 0;
600  }
601  /**
602  * Set polynomial to -1
603  *
604  * @param
605  **/
606  inline void negativeOne() {
607  terms.clear();
608  MultivariateTerm<Ring> t;
609  t.coef.negativeOne();
610  t.v = var;
611  t.degs = new int[var];
612  for (int i = 0; i < var; ++i)
613  t.degs[i] = 0;
614  terms.push_back(t);
615  }
616  /**
617  * Is a constant
618  *
619  * @param
620  **/
621  inline int isConstant() {
622  if (!terms.size())
623  return 1;
624  else if (terms.size() == 1) {
625  for (int i = 0; i < var; ++i) {
626  if (terms[0].degs[i] != 0)
627  return 0;
628  }
629  if (terms[0].coef.isConstant() >= 0) { return 1; }
630  else { return -1; }
631  }
632  return 0;
633  }
634  /**
635  * Overload operator =
636  *
637  * @param b: A sparse multivariate polynomial
638  **/
640  if (this != &b) {
641  terms.clear();
642  delete [] names;
643 
644  var = b.numberOfVariables();
645  names = new std::string[var+1];
646  std::copy(b.names, b.names+var+1, names);
647  terms = b.terms;
648  Ring e;
649  characteristic = e.characteristic;
650  }
651  return *this;
652  }
653  /**
654  * Overload operator ==
655  *
656  * @param b: A multivariate rational polynomial
657  **/
659  return isEqual(b);
660  }
661  /**
662  * Overload operator !=
663  *
664  * @param b: A multivariate rational polynomial
665  **/
667  return !isEqual(b);
668  }
669  /**
670  * Overload operator ^
671  * replace xor operation by exponentiation
672  *
673  * @param e: The exponentiation, e > 0
674  **/
677  res.var = var;
678  res.names = new std::string[var+1];
679  std::copy(names, names+var+1, res.names);
680  res.one();
681  unsigned long int q = e / 2, r = e % 2;
682  SparseMultivariatePolynomial<Ring> p2 = *this * *this;
683  for (int i = 0; i < q; ++i)
684  res *= p2;
685  if (r) { res *= *this; }
686  return res;
687  }
688  /**
689  * Overload operator ^=
690  * replace xor operation by exponentiation
691  *
692  * @param e: The exponentiation, e > 0
693  **/
695  *this = *this ^ e;
696  return *this;
697  }
698 
699  /**
700  * Overload operator +
701  *
702  * @param b: A multivariate rational polynomial
703  **/
705  if (!terms.size()) { return b; }
706  if (!b.terms.size()) { return *this; }
707  if (isConstant()) { return (b + terms[0].coef); }
708  if (b.isConstant()) { return (*this + b.terms[0].coef); }
709  std::vector<int> xs;
710  bool isOrdered = isOrderedRing(b, xs);
711 
712  if (!isOrdered) {
713  std::cout << "BPAS: error, trying to add between Ring[";
714  for (int i = 1; i <= var; ++i) {
715  std::cout << names[i];
716  if (i < var) { std::cout << ", "; }
717  }
718  std::cout << "] and Ring[";
719  for (int i = 1; i <= b.var; ++i) {
720  std::cout << b.names[i];
721  if (i < b.var) { std::cout << ", "; }
722  }
723  std::cout << "] from SparseMultivariatePolynomial<Ring>." << std::endl;
724  exit(1);
725  }
726 
727  int v = xs.size() / 2;
729  r.names[0] = names[0];
730  if (names[0] == "9") {
731  for (int i = 1; i <= v; ++i) {
732  if (xs[2*i-2])
733  r.names[i] = names[xs[2*i-2]];
734  else
735  r.names[i] = b.names[xs[2*i-1]];
736  }
737  }
738 
739  int i = 0, j = 0;
740  while (i < terms.size() && j < b.terms.size()) {
741  MultivariateTerm<Ring> t;
742  t.v = v;
743  t.degs = new int[v];
744  int c = compareTermDegs(terms[i], b.terms[j], xs);
745  if (!c) {
746  t.coef = terms[i].coef + b.terms[j].coef;
747  if (!t.coef.isZero()) {
748  for (int k = 0; k < v; ++k) {
749  if (xs[2*k])
750  t.degs[k] = terms[i].degs[xs[2*k]-1];
751  else
752  t.degs[k] = b.terms[j].degs[xs[2*k+1]-1];
753  }
754  r.terms.push_back(t);
755  }
756  i++, j++;
757  }
758  else if (c < 0) {
759  t.coef = terms[i].coef;
760  for (int k = 0; k < v; ++k) {
761  if (xs[2*k])
762  t.degs[k] = terms[i].degs[xs[2*k]-1];
763  else
764  t.degs[k] = 0;
765  }
766  r.terms.push_back(t);
767  i++;
768  }
769  else {
770  t.coef = b.terms[j].coef;
771  for (int k = 0; k < v; ++k) {
772  if (xs[2*k+1])
773  t.degs[k] = b.terms[j].degs[xs[2*k+1]-1];
774  else
775  t.degs[k] = 0;
776  }
777  r.terms.push_back(t);
778  j++;
779  }
780  }
781  for (; i < terms.size(); ++i) {
782  MultivariateTerm<Ring> t;
783  t.v = v;
784  t.coef = terms[i].coef;
785  t.degs = new int[v];
786  for (int k = 0; k < v; ++k) {
787  if (xs[2*k])
788  t.degs[k] = terms[i].degs[xs[2*k]-1];
789  else
790  t.degs[k] = 0;
791  }
792  r.terms.push_back(t);
793  }
794  for (; j < b.terms.size(); ++j) {
795  MultivariateTerm<Ring> t;
796  t.v = v;
797  t.coef = b.terms[j].coef;
798  t.degs = new int[v];
799  for (int k = 0; k < v; ++k) {
800  if (xs[2*k+1])
801  t.degs[k] = b.terms[j].degs[xs[2*k+1]-1];
802  else
803  t.degs[k] = 0;
804  }
805  r.terms.push_back(t);
806  }
807  return r;
808  }
809  /**
810  * Overload operator +=
811  *
812  * @param b: A multivariate rational polynomial
813  **/
815  *this = *this + b;
816  return *this;
817  }
818  /**
819  * Overload operator +
820  *
821  * @param c: A constant
822  **/
825  return (r += c);
826  }
827  /**
828  * Overload operator +=
829  *
830  * @param c: A constant
831  **/
833  basicOp(c, 0);
834  return *this;
835  }
836 
838  return (p + c);
839  }
840  /**
841  * Overload operator -
842  *
843  * @param b: A multivariate rational polynomial
844  **/
846  if (!terms.size()) { return -b; }
847  if (!b.terms.size()) { return *this; }
848  if (isConstant()) { return (-b + terms[0].coef); }
849  if (b.isConstant()) { return (*this - b.terms[0].coef); }
850 
851  std::vector<int> xs;
852  bool isOrdered = isOrderedRing(b, xs);
853 
854  if (!isOrdered) {
855  std::cout << "BPAS: error, trying to subtract between Ring[";
856  for (int i = 1; i <= var; ++i) {
857  std::cout << names[i];
858  if (i < var) { std::cout << ", "; }
859  }
860  std::cout << "] and Ring[";
861  for (int i = 1; i <= b.var; ++i) {
862  std::cout << b.names[i];
863  if (i < b.var) { std::cout << ", "; }
864  }
865  std::cout << "] from SparseMultivariatePolynomial<Ring>." << std::endl;
866  exit(1);
867  }
868 
869  int v = xs.size() / 2;
871  r.names[0] = names[0];
872  if (names[0] == "9") {
873  for (int i = 1; i <= v; ++i) {
874  if (xs[2*i-2])
875  r.names[i] = names[xs[2*i-2]];
876  else
877  r.names[i] = b.names[xs[2*i-1]];
878  }
879  }
880 
881  int i = 0, j = 0;
882  while (i < terms.size() && j < b.terms.size()) {
883  MultivariateTerm<Ring> t;
884  t.v = v;
885  t.degs = new int[v];
886  int c = compareTermDegs(terms[i], b.terms[j], xs);
887  if (!c) {
888  t.coef = terms[i].coef - b.terms[j].coef;
889  if (!t.coef.isZero()) {
890  for (int k = 0; k < v; ++k) {
891  if (xs[2*k])
892  t.degs[k] = terms[i].degs[xs[2*k]-1];
893  else
894  t.degs[k] = b.terms[j].degs[xs[2*k+1]-1];
895  }
896  r.terms.push_back(t);
897  }
898  i++, j++;
899  }
900  else if (c < 0) {
901  t.coef = terms[i].coef;
902  for (int k = 0; k < v; ++k) {
903  if (xs[2*k])
904  t.degs[k] = terms[i].degs[xs[2*k]-1];
905  else
906  t.degs[k] = 0;
907  }
908  r.terms.push_back(t);
909  i++;
910  }
911  else {
912  t.coef = -b.terms[j].coef;
913  for (int k = 0; k < v; ++k) {
914  if (xs[2*k+1])
915  t.degs[k] = b.terms[j].degs[xs[2*k+1]-1];
916  else
917  t.degs[k] = 0;
918  }
919  r.terms.push_back(t);
920  j++;
921  }
922  }
923  for (; i < terms.size(); ++i) {
924  MultivariateTerm<Ring> t;
925  t.v = v;
926  t.coef = terms[i].coef;
927  t.degs = new int[v];
928  for (int k = 0; k < v; ++k) {
929  if (xs[2*k])
930  t.degs[k] = terms[i].degs[xs[2*k]-1];
931  else
932  t.degs[k] = 0;
933  }
934  r.terms.push_back(t);
935  }
936  for (; j < b.terms.size(); ++j) {
937  MultivariateTerm<Ring> t;
938  t.v = v;
939  t.coef = -b.terms[j].coef;
940  t.degs = new int[v];
941  for (int k = 0; k < v; ++k) {
942  if (xs[2*k+1])
943  t.degs[k] = b.terms[j].degs[xs[2*k+1]-1];
944  else
945  t.degs[k] = 0;
946  }
947  r.terms.push_back(t);
948  }
949  return r;
950  }
951  /**
952  * Overload operator -=
953  *
954  * @param b: A multivariate rational polynomial
955  **/
957  *this = *this - b;
958  return *this;
959  }
960  /**
961  * Overload operator -
962  *
963  * @param c: A constant
964  **/
967  return (r -= c);
968  }
969  /**
970  * Overload operator -=
971  *
972  * @param c: A constant
973  **/
975  basicOp(c, 1);
976  return *this;
977  }
979  return (-p + c);
980  }
981  /**
982  * Overload operator - (negate)
983  *
984  * @param
985  **/
988  std::copy(names, names+var+1, res.names);
989  for (int i = 0; i < terms.size(); ++i) {
990  MultivariateTerm<Ring> t;
991  t.coef = -terms[i].coef;
992  t.v = var;
993  t.degs = new int[var];
994  for (int j = 0; j < var; ++j)
995  t.degs[j] = terms[i].degs[j];
996  res.terms.push_back(t);
997  }
998  return res;
999  }
1000  /**
1001  * Overload operator *
1002  *
1003  * @param b: A multivariate rational polynomial
1004  **/
1006  if (!terms.size() || !b.terms.size()) {
1008  r.var = var;
1009  r.names = new std::string[var+1];
1010  std::copy(names, names+var+1, r.names);
1011  return r;
1012  }
1013  if (isConstant()) { return (b * terms[0].coef); }
1014  if (b.isConstant()) { return (*this * b.terms[0].coef); }
1015  std::vector<int> xs;
1016  bool isOrdered = isOrderedRing(b, xs);
1017 
1018  if (!isOrdered) {
1019  std::cout << "BPAS: error, trying to multiply between Ring[";
1020  for (int i = 1; i <= var; ++i) {
1021  std::cout << names[i];
1022  if (i < var) { std::cout << ", "; }
1023  }
1024  std::cout << "] and Ring[";
1025  for (int i = 1; i <= b.var; ++i) {
1026  std::cout << b.names[i];
1027  if (i < b.var) { std::cout << ", "; }
1028  }
1029  std::cout << "] from SparseMultivariatePolynomial<Ring>." << std::endl;
1030  exit(1);
1031  }
1032 
1033  int v = xs.size() / 2;
1035  r.names[0] = names[0];
1036  if (names[0] == "9") {
1037  for (int i = 1; i <= v; ++i) {
1038  if (xs[2*i-2])
1039  r.names[i] = names[xs[2*i-2]];
1040  else
1041  r.names[i] = b.names[xs[2*i-1]];
1042  }
1043  }
1044  for (int i = 0; i < terms.size(); ++i)
1045  r.pomopo(terms[i], b, xs);
1046  return r;
1047  }
1048  /**
1049  * Overload operator *=
1050  *
1051  * @param b: A multivariate rational polynomial
1052  **/
1054  *this = *this * b;
1055  return *this;
1056  }
1057  /**
1058  * Overload operator *
1059  *
1060  * @param e: A constant
1061  **/
1064  return (r *= e);
1065  }
1068  return (r *= e);
1069  }
1070  /**
1071  * Overload operator *=
1072  *
1073  * @param e: A constant
1074  **/
1076  if (!e.isZero() && !e.isOne()) {
1077  for (int i = 0; i < terms.size(); ++i)
1078  terms[i].coef *= e;
1079  }
1080  else if (e.isZero()) { terms.clear(); }
1081  return *this;
1082  }
1084  return (p * c);
1085  }
1086  /**
1087  * Overload operator *=
1088  *
1089  * @param e: A machine-word constant
1090  **/
1092  if (e != 0 && e != 1) {
1093  for (int i = 0; i < terms.size(); ++i)
1094  terms[i].coef *= e;
1095  }
1096  else if (e == 0) { terms.clear(); }
1097  return *this;
1098  }
1100  return (p * c);
1101  }
1102  /**
1103  * Overload operator /
1104  *
1105  * @param b: A multivariate polynomial
1106  **/
1109  return (r /= b);
1110  }
1111  /**
1112  * Overload operator /=
1113  *
1114  * @param b: A multivariate polynomial
1115  **/
1117  if (!b.terms.size()) {
1118  std::cout << "BPAS: error, dividend is zero from SparseMultivariatePolynomial<Ring>." << std::endl;
1119  exit(1);
1120  }
1121  if (b.isConstant()) { return (*this /= b.terms[0].coef); }
1122  if (isConstant()) {
1123  terms.clear();
1124  return *this;
1125  }
1126 
1127  bool isIt = 0;
1128  int s = 0;
1129  std::vector<int> xs;
1130  for (int i = 1; i <= b.var; ++i) {
1131  bool isFound = 0;
1132  for (int j = i; j <= var; ++j)
1133  if (b.names[i] == names[j]) {
1134  if (j > i + s) {
1135  for (int k = i; k < j; ++k) {
1136  xs.push_back(k);
1137  xs.push_back(0);
1138  }
1139  s = j - i;
1140  }
1141  xs.push_back(j);
1142  xs.push_back(i);
1143  isFound = 1;
1144  break;
1145  }
1146  if (!isFound) {
1147  isIt = 1;
1148  break;
1149  }
1150  }
1151  for (int i = xs.size()/2; i <= var; ++i) {
1152  xs.push_back(i);
1153  xs.push_back(0);
1154  }
1155 
1156  if (isIt) {
1157  std::cout << "BPAS: error, trying to exact divide between Ring[";
1158  for (int i = 1; i <= var; ++i) {
1159  std::cout << names[i];
1160  if (i < var)
1161  std::cout << ", ";
1162  }
1163  std::cout << "] and Ring[";
1164  for (int i = 1; i <= b.var; ++i) {
1165  std::cout << b.names[i];
1166  if (i < var)
1167  std::cout << ", ";
1168  }
1169  std::cout << "] from SparseMultivariatePolynomial<Ring>." << std::endl;
1170  exit(1);
1171  }
1172 
1174  terms.clear();
1175 
1176  MultivariateTerm<Ring> bt = b.terms[b.terms.size()-1];
1177  s = compareTermDegs(rem.terms[rem.terms.size()-1], bt, xs);
1178  while (s >= 0) {
1179  MultivariateTerm<Ring> at = rem.terms[rem.terms.size()-1];
1180  MultivariateTerm<Ring> lc, nlc;
1181  lc.coef = at.coef / bt.coef;
1182  nlc.coef = - lc.coef;
1183  lc.v = var;
1184  nlc.v = var;
1185  lc.degs = new int[var];
1186  nlc.degs = new int[var];
1187  for (int i = 0; i < var; ++i) {
1188  nlc.degs[i] = lc.degs[i] = at.degs[xs[2*i]-1] - bt.degs[xs[2*i+1]-1];
1189  if (lc.degs[i] < 0) {
1190  std::cout << "BPAS: error, not exact division from SparseMultivariatePolynomial<Ring>." << std::endl;
1191  exit(1);
1192  }
1193  }
1194  rem.pomopo(nlc, b, xs);
1195  terms.insert(terms.begin(), lc);
1196  if (rem.terms.size())
1197  s = compareTermDegs(rem.terms[rem.terms.size()-1], bt, xs);
1198  else { s = -1; }
1199  }
1200  if (!rem.isZero()) {
1201  std::cout << "BPAS: error, not exact division from SparseMultivariatePolynomial<Ring>." << std::endl;
1202  exit(1);
1203  }
1204 
1205  return *this;
1206  }
1207  /**
1208  * Overload operator /
1209  *
1210  * @param e: A constant
1211  **/
1214  return (r /= e);
1215  }
1216  /**
1217  * Overload operator /=
1218  *
1219  * @param e: A constant
1220  **/
1222  if (e.isZero()) {
1223  std::cout << "BPAS: error, dividend is zero." << std::endl;
1224  exit(1);
1225  }
1226  else if (!e.isOne()) {
1227  for (int i = 0; i < terms.size(); ++i)
1228  terms[i].coef /= e;
1229  }
1230  return *this;
1231  }
1233  if (p.isZero()) {
1234  std::cout << "BPAS: error, dividend is zero from SparseMultivariatePolynomial<Ring>." << std::endl;
1235  exit(1);
1236  }
1237 
1239  q.var = p.var;
1240  q.names = new std::string[q.var+1];
1241  std::copy(p.names, p.names+p.var+1, q.names);
1242 
1243  if (p.isConstant()) {
1244  q += c;
1245  return (q /= p.terms[0].coef);
1246  }
1247  else { return q; }
1248  }
1249 
1250  /**
1251  * Set variable names
1252  *
1253  * @param xs: Variable names
1254  **/
1255  inline void setVariableNames (std::vector<std::string> xs) {
1256  int ns = xs.size();
1257  if (ns != var) {
1258  std::cout << "BPAS: error, SMQP(" << var << "), but trying to setVariableNames with " << ns << " variables." << std::endl;
1259  exit(1);
1260  }
1261  if (names[0] == "1") {
1262  names[0] = "9";
1263  for (int i = var, j = 0; i > 0 && j < ns; --i, ++j)
1264  names[i] = xs[j];
1265  }
1266  else {
1267  bool isReorder = 0;
1268  int* pos = new int[var];
1269  for (int i = 0; i < ns; ++i) {
1270  bool isIt = 1;
1271  for (int j = 0; j < var; ++j) {
1272  if (names[j+1] == xs[ns-i-1]) {
1273  pos[i] = j;
1274  isIt = 0;
1275  break;
1276  }
1277  }
1278  if (isIt) { pos[i] = ns - i - 1; }
1279  else { isReorder = 1; }
1280  }
1281 
1282  if (isReorder) {
1284  r.var = var;
1285  r.names = new std::string[var+1];
1286  r.names[0] = "9";
1287  for (int i = var, j = 0; i > 0 && j < ns; --i, ++j)
1288  r.names[i] = xs[j];
1289 
1290  int* d = new int[var];
1291  for (int i = 0; i < terms.size(); ++i) {
1292  for (int j = 0; j < var; ++j)
1293  d[var-j-1] = terms[i].degs[pos[j]];
1294  r.setCoefficient(var, d, terms[i].coef);
1295  }
1296  delete [] d;
1297 
1298  *this = r;
1299  }
1300  else {
1301  for (int i = var, j = 0; i > 0 && j < ns; --i, ++j)
1302  names[i] = xs[j];
1303  }
1304 
1305  delete [] pos;
1306  }
1307  }
1308  /**
1309  * Get variable names
1310  *
1311  * @param
1312  **/
1313  inline std::vector<std::string> variables() {
1314  std::vector<std::string> xs;
1315  for (int i = var; i > 0; --i)
1316  xs.push_back(names[i]);
1317  return xs;
1318  }
1319  /**
1320  * Get the coefficient of one term
1321  *
1322  * @param d: The exponet of each variable
1323  * @param v: Number of variables
1324  **/
1325  inline Ring coefficient(int v, int* d) {
1326  int n = terms.size();
1327  for (int i = 0; i < n; ++i) {
1328  bool isIt = 1;
1329  for (int j = var-1, k = 0; j > -1 && k < v; --j, ++k) {
1330  if (terms[i].degs[j] != d[k]) {
1331  isIt = 0;
1332  break;
1333  }
1334  }
1335  for (int j = v-var-1; j > -1; --j) {
1336  if (d[j] != 0) {
1337  isIt = 0;
1338  break;
1339  }
1340  }
1341  if (isIt)
1342  return terms[i].coef;
1343  }
1344  return Ring();
1345  }
1346  /**
1347  * Set the coefficients of one term
1348  *
1349  * @param v: Number of variables
1350  * @param d: Its exponent of each variable
1351  * @param val: Coefficient
1352  **/
1353  inline void setCoefficient(int v, int* d, Ring val) {
1354  if (v != var) {
1355  std::cout << "BPAS: error, SMQP(" << var << "), but trying to setCoefficient with " << v << " variables." << std::endl;
1356  exit(1);
1357  }
1358  MultivariateTerm<Ring> a;
1359  a.coef = val;
1360  a.v = var;
1361  a.degs = new int[var];
1362  for (int i = 0; i < var; ++i)
1363  a.degs[i] = d[v-i-1];
1364 
1365  int i = 0;
1366  for (; i < terms.size(); ++i) {
1367  int k = compareTermDegs(a, terms[i]);
1368  if (!k) {
1369  if (val.isZero())
1370  terms.erase(terms.begin()+i);
1371  else
1372  terms[i].coef = val;
1373  break;
1374  }
1375  else if (k < 0 && !val.isZero()) {
1376  terms.insert(terms.begin()+i, a);
1377  break;
1378  }
1379  }
1380  if (i == terms.size() && !val.isZero())
1381  terms.push_back(a);
1382  }
1383 
1384  /**
1385  * Negate all the coefficients
1386  *
1387  * @param
1388  **/
1389  inline void negate() {
1390  for (int i = 0; i < terms.size(); ++i)
1391  terms[i].coef = -terms[i].coef;
1392  }
1393 
1394  /**
1395  * Convert to SUP<SMQP>
1396  *
1397  * @param x: The variable name
1398  **/
1401  r.setVariableName(x);
1402 
1403  int k = 0;
1404  for (int i = 1; i <= var; ++i)
1405  if (names[i] == x) { k = i; break; }
1406 
1407  if (k) {
1408  int d = 0;
1409  for (int i = terms.size()-1; i > -1; --i)
1410  if (terms[i].degs[k-1] > d) { d = terms[i].degs[k-1]; }
1411  int v = var - 1;
1413  for (int i = 0; i <= d; ++i) {
1414  t[i].var = v;
1415  delete [] t[i].names;
1416  t[i].names = new std::string[v+1];
1417  for (int j = 0; j < var; ++j) {
1418  if (j < k) { t[i].names[j] = names[j]; }
1419  else { t[i].names[j] = names[j+1]; }
1420  }
1421  }
1422  k--;
1423  for (int i = 0; i < terms.size(); ++i) {
1424  MultivariateTerm<Ring> a;
1425  a.v = v;
1426  a.coef = terms[i].coef;
1427  a.degs = new int[v];
1428  for (int j = 0; j < v; ++j) {
1429  if (j < k) { a.degs[j] = terms[i].degs[j]; }
1430  else { a.degs[j] = terms[i].degs[j+1]; }
1431  }
1432  t[terms[i].degs[k]].terms.push_back(a);
1433  }
1434  for (int i = 0; i <= d; ++i)
1435  r.setCoefficient(i, t[i]);
1436  delete [] t;
1437  }
1438  else { r.setCoefficient(0, *this); }
1439 
1440  return r;
1441  }
1442 
1443  /**
1444  * Overload stream operator <<
1445  *
1446  * @param out: Stream object
1447  * @param b: A multivariate rational polynoial
1448  **/
1449  inline friend std::ostream& operator<< (std::ostream& out, SparseMultivariatePolynomial<Ring> b) {
1450  int n = b.numberOfTerms();
1451  int var = b.numberOfVariables();
1452  if (!n) { out << "0"; }
1453  for (int i = 0; i < n; i++) {
1454  if (i && b.terms[i].coef.isConstant() >= 0)
1455  out << "+";
1456  else if (b.terms[i].coef.isNegativeOne())
1457  out << "-";
1458  if (!b.terms[i].coef.isOne() && !b.terms[i].coef.isNegativeOne())
1459  out << b.terms[i].coef;
1460  bool isIt = 1;
1461  int* d = b.terms[i].degs;
1462  for (int j = 0; j < b.terms[i].v; ++j) {
1463  if (d[j]) {
1464  if ((!b.terms[i].coef.isOne() && !b.terms[i].coef.isNegativeOne() && isIt) || !isIt)
1465  out << "*";
1466  out << b.names[j+1];
1467  if (d[j] > 1)
1468  out << "^" << d[j];
1469  isIt = 0;
1470  }
1471  }
1472  if (isIt && (b.terms[i].coef.isOne() || b.terms[i].coef.isNegativeOne()))
1473  out << "1";
1474  }
1475  return out;
1476  }
1477 };
1478 
1479 #endif
1480 /* This file is part of the BPAS library http://www.bpaslib.org
1481 
1482  BPAS is free software: you can redistribute it and/or modify
1483  it under the terms of the GNU General Public License as published by
1484  the Free Software Foundation, either version 3 of the License, or
1485  (at your option) any later version.
1486 
1487  BPAS is distributed in the hope that it will be useful,
1488  but WITHOUT ANY WARRANTY; without even the implied warranty of
1489  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1490  GNU General Public License for more details.
1491 
1492  You should have received a copy of the GNU General Public License
1493  along with BPAS. If not, see <http://www.gnu.org/licenses/>.
1494 
1495  Copyright: Changbo Chen <changbo.chen@hotmail.com>
1496  Farnam Mansouri <mansouri.farnam@gmail.com>
1497  Marc Moreno Maza <moreno@csd.uwo.ca>
1498  Ning Xie <nxie6@csd.uwo.ca>
1499  Yuzhen Xie <yuzhenxie@yahoo.ca>
1500 
1501 */
1502 
1503 
int degree(std::string x)
Definition: mpolynomial.h:431
Definition: upolynomial.h:10
int numberOfTerms()
Definition: mpolynomial.h:423
int isConstant()
Definition: mpolynomial.h:621
bool operator==(SparseMultivariatePolynomial< Ring > &b)
Definition: mpolynomial.h:658
bool isZero()
Definition: mpolynomial.h:543
SparseMultivariatePolynomial< Ring > operator^(int e)
Definition: mpolynomial.h:675
bool isNegativeOne()
Definition: mpolynomial.h:590
SparseUnivariatePolynomial< SparseMultivariatePolynomial< Ring > > convertToSUP(std::string x)
Definition: mpolynomial.h:1399
SparseMultivariatePolynomial< Ring > & operator*=(SparseMultivariatePolynomial< Ring > b)
Definition: mpolynomial.h:1053
int leadingVariableDegree()
Definition: mpolynomial.h:461
SparseMultivariatePolynomial< Ring > & operator/=(SparseMultivariatePolynomial< Ring > b)
Definition: mpolynomial.h:1116
void setCoefficient(int v, int *d, Ring val)
Definition: mpolynomial.h:1353
std::vector< std::string > variables()
Definition: mpolynomial.h:1313
SparseMultivariatePolynomial< Ring > & operator+=(SparseMultivariatePolynomial< Ring > b)
Definition: mpolynomial.h:814
Definition: mpolynomial.h:12
SparseMultivariatePolynomial< Ring > operator*(SparseMultivariatePolynomial< Ring > b)
Definition: mpolynomial.h:1005
void negativeOne()
Definition: mpolynomial.h:606
SparseMultivariatePolynomial< Ring > operator/(SparseMultivariatePolynomial< Ring > b)
Definition: mpolynomial.h:1107
SparseMultivariatePolynomial< Ring > & operator=(SparseMultivariatePolynomial< Ring > b)
Definition: mpolynomial.h:639
SparseMultivariatePolynomial< Ring > leadingCoefficientInVariable(std::string x, int *e=NULL)
Definition: mpolynomial.h:486
bool isOne()
Definition: mpolynomial.h:559
Definition: polynomial.h:74
void negate()
Definition: mpolynomial.h:1389
SparseMultivariatePolynomial< Ring > & operator^=(int e)
Definition: mpolynomial.h:694
void zero()
Definition: mpolynomial.h:551
void one()
Definition: mpolynomial.h:575
SparseMultivariatePolynomial< Ring > & operator-=(SparseMultivariatePolynomial< Ring > b)
Definition: mpolynomial.h:956
SparseMultivariatePolynomial< Ring > operator+(SparseMultivariatePolynomial< Ring > b)
Definition: mpolynomial.h:704
Ring leadingCoefficient()
Definition: mpolynomial.h:473
Ring coefficient(int v, int *d)
Definition: mpolynomial.h:1325
void setCoefficient(int e, Ring c)
Definition: upolynomial.h:383
void setVariableName(std::string c)
Definition: upolynomial.h:422
bool operator!=(SparseMultivariatePolynomial< Ring > &b)
Definition: mpolynomial.h:666
int numberOfVariables()
Definition: mpolynomial.h:414
void setVariableNames(std::vector< std::string > xs)
Definition: mpolynomial.h:1255
SparseMultivariatePolynomial< Ring > operator-()
Definition: mpolynomial.h:986
std::string leadingVariable()
Definition: mpolynomial.h:453