1 #include "PolyTrie.h"
2 #include <stdio.h>
3 #include <sstream>
4
5 #define PT_DEBUG 0
6
7 //Loads a string by parsing it as a sum of monomials
8 //monomial sum: c_{1}*(x_{1}^e_{1}...x_{varCount}^e_{varCount}) + ...
9 //nested lists: [[c_{1}, [e_{1}, e_{2}, ..., e_{varCount}]], .. ]
loadMonomials(monomialSum & monomials,const string & line)10 void loadMonomials(monomialSum &monomials, const string &line)
11 {
12 monomials.termCount = 0;
13 MonomialLoadConsumer<RationalNTL>* myLoader = new MonomialLoadConsumer<
14 RationalNTL> ();
15 myLoader->setMonomialSum(monomials);
16 parseMonomials(myLoader, line);
17 delete myLoader;
18 }
19
parseMonomials(MonomialConsumer<RationalNTL> * consumer,const string & line)20 void parseMonomials(MonomialConsumer<RationalNTL>* consumer, const string &line)
21 {
22 int varCount = 0;
23 for (int i = 0; line[i] != ']'; i++)
24 {
25 varCount += (line[i] == ',');
26 }
27 if (varCount < 1)
28 {
29 cout << "line: `" << line << "'" << endl;
30 cout << "There are " << varCount << " variables, bailing." << endl;
31 return;
32 }
33 consumer->setDimension(varCount);
34
35 int termIndex, lastPos, expIndex, flag;
36 termIndex = lastPos = flag = 0; //0 means we expect coefficient, 1 means we expect exponent vector
37
38 int *exponents = new int[varCount];
39 RationalNTL coefficient;
40
41 for (size_t i = 1; i < line.length() - 1; i++) //ignore outermost square brackets
42 {
43 if (line[i] == '[')
44 {
45 switch (flag)
46 {
47 case 0: //coefficient
48 lastPos = i + 1;
49 for (; line[i] != ','; i++)
50 ;
51 coefficient = RationalNTL(
52 line.substr(lastPos, i - lastPos).c_str());
53 flag = 1;
54 break;
55 case 1: //exponent vector
56 expIndex = 0;
57 for (i++; line[i] != ']'; i++)
58 {
59 if (line[i] != ' ')
60 {
61 lastPos = i;
62 for (; line[i] != ',' && line[i] != ']'; i++)
63 ;
64 exponents[expIndex++] = atoi(line.substr(lastPos, i
65 - lastPos).c_str());
66 }
67 }
68 consumer->ConsumeMonomial(coefficient, exponents);
69 flag = 0;
70 break;
71 default: //error
72 cout << "Flag is " << flag << ", bailing." << endl;
73 return;
74 }
75 }
76 }
77
78 delete[] exponents;
79 }
80
81 //watch the magic happen
insertMonomial(const RationalNTL & coefficient,int * exponents,monomialSum & monomials)82 void insertMonomial(const RationalNTL& coefficient, int* exponents,
83 monomialSum& monomials)
84 {
85 BurstTrie<RationalNTL, int>* curTrie;
86
87 if ( coefficient == 0)
88 return;
89
90 if (monomials.termCount == 0) //need to construct the first burst trie (sorted on the first variable) and first container
91 {
92 if (PT_DEBUG)
93 {
94 cout << "Creating trie" << endl;
95 }
96 monomials.myMonomials = new BurstTrie<RationalNTL, int> ();
97 curTrie = monomials.myMonomials;
98 } else
99 {
100 curTrie = monomials.myMonomials;
101 }
102
103 curTrie->insertTerm(coefficient, exponents, 0, monomials.varCount, -1);
104
105 monomials.termCount++;
106 }
107
108 //Prints a nested list representation of our sum of monomials
109 //monomial sum: c_{1}*(x_{1}^e_{1}...x_{varCount}^e_{varCount}) + ...
110 //nested lists: [[c_{1}, [e_{1}, e_{2}, ..., e_{varCount}]], .. ]
printMonomials(const monomialSum & myPoly)111 string printMonomials(const monomialSum &myPoly)
112 {
113 BTrieIterator<RationalNTL, int>* it =
114 new BTrieIterator<RationalNTL, int> ();
115 term<RationalNTL, int>* temp;
116 it->setTrie(myPoly.myMonomials, myPoly.varCount);
117 it->begin();
118 int i = 0;
119 stringstream output(stringstream::in | stringstream::out);
120 temp = it->nextTerm();
121 do
122 {
123 if (output.str() != "")
124 {
125 output << ", ";
126 }
127 output << "[" << temp->coef << ", [";
128 for (int j = 0; j < temp->length; j++)
129 {
130 output << temp->exps[j];
131 if (j + 1 < temp->length)
132 {
133 output << ", ";
134 }
135 }
136 output << "]]";
137 temp = it->nextTerm();
138 } while (temp);
139 delete it;
140 return "[" + output.str() + "]";
141 }
142
143 //Deallocates space and nullifies internal pointers and counters
destroyMonomials(monomialSum & myPoly)144 void destroyMonomials(monomialSum &myPoly)
145 {
146 delete myPoly.myMonomials;
147 myPoly.myMonomials = NULL;
148 myPoly.termCount = myPoly.varCount = 0;
149 }
150
151 // -------------------------------------------------------------------------
152
loadLinForms(linFormSum & forms,const string line)153 void loadLinForms(linFormSum &forms, const string line)
154 {
155 forms.termCount = 0;
156 FormLoadConsumer<RationalNTL>* myLoader =
157 new FormLoadConsumer<RationalNTL> ();
158 myLoader->setFormSum(forms);
159 parseLinForms(myLoader, line);
160 delete myLoader;
161 }
162 //Loads a string by parsing it as a sum of linear forms
163 //linear form: (c_{1} / d_{1}!)[(p_{1}*x_{1} + ... p_{varCount}*x_{varCount})^d_{1}] + ...
164 //nested list: [[c_{1}, [d_{1}, [p_{1}, p_{2}, ..., p_{varCount}]], .. ]
parseLinForms(FormSumConsumer<RationalNTL> * consumer,const string & line)165 void parseLinForms(FormSumConsumer<RationalNTL>* consumer, const string& line)
166 {
167 int termIndex = 0;
168 int lastPos = 0;
169 int varCount = 0;
170 int k;
171 int flag = 0; //0 means we expect coefficient, 1 means we expect degree, 2 means we expect coefficient vector
172
173 //cout << "parseLinForms: line = " << line.c_str() << endl;
174 for (int i = 0; line[i] != ']'; i++)
175 {
176 varCount += (line[i] == ',');
177 }
178 //varCount is now the number of commas in a linear form - there is 1 less variable;
179 varCount--;
180 if (varCount < 1)
181 {
182 cout << "line: `" << line << "'" << endl;
183 cout << "There are " << varCount << " variables, bailing." << endl;
184 return;
185 }
186 consumer->setDimension(varCount);
187
188 vec_ZZ coefs;
189 coefs.SetLength(varCount);
190 int degree;
191 RationalNTL coefficient;
192
193 for (size_t i = 1; i < line.length() - 1; i++) //ignore outermost square brackets
194 {
195 if (line[i] == '[')
196 {
197 ZZ degreeFactorial;
198 switch (flag)
199 {
200 case 0: //coefficient
201 lastPos = i + 1;
202 for (; line[i] != ','; i++)
203 ;
204 coefficient = RationalNTL(
205 (line.substr(lastPos, i - lastPos).c_str()));
206 flag = 1;
207 break;
208 case 1: //degree
209 lastPos = i + 1;
210 for (; line[i] != ','; i++)
211 ;
212 degree = atoi(line.substr(lastPos, i - lastPos).c_str());
213 flag = 2;
214 break;
215 case 2: //coefficient vector
216 k = 0;
217 for (i++; line[i] != ']'; i++)
218 {
219 if (line[i] != ' ')
220 {
221 lastPos = i;
222 for (; line[i] != ',' && line[i] != ']'; i++)
223 ;
224
225 coefs[k++] = to_ZZ(
226 line.substr(lastPos, i - lastPos).c_str());
227
228 }
229 }
230 degreeFactorial = 1;
231 for (int j = 1; j <= degree; j++)
232 {
233 degreeFactorial *= j;
234 //coefficient *= j;
235 } //in linFormSum, coefficient is assumed to be divided by the factorial of the form degree
236 coefficient *= degreeFactorial;
237 consumer->ConsumeLinForm(coefficient, degree, coefs);
238 flag = 0;
239 break;
240 default: //error
241 cout << "Flag is " << flag << ", bailing." << endl;
242 return;
243 }
244 }
245 }
246 }
247
248 // Attempts to find a linear form in formSum with same degree and coefficients as those passed in
249 // if found, the linear form coefficient in formSum is incremented by coef
250 // if not found, a new linear form term is added and formSum's termCount is incremented
251 // if formSum is empty, this sets formSum's lHead and cHead variables
insertLinForm(const RationalNTL & coef,int degree,const vec_ZZ & coeffs,linFormSum & formSum)252 void insertLinForm(const RationalNTL& coef, int degree, const vec_ZZ& coeffs,
253 linFormSum& formSum) //sort on degree first or last?
254 {
255 BurstTrie<RationalNTL, ZZ> *curTrie;
256 //cout << "inserting into linear form with " << formSum.varCount << " variables" << endl;
257
258 if (coef == 0)
259 return;
260
261 if (formSum.termCount == 0) //need to construct the first burst trie (sorted on the first variable) and first container
262 {
263 formSum.myForms = new BurstTrie<RationalNTL, ZZ> ();
264 curTrie = formSum.myForms;
265 } else
266 {
267 curTrie = formSum.myForms;
268 }
269
270 ZZ* exps = new ZZ[formSum.varCount];
271 for (int i = 0; i < formSum.varCount; i++)
272 {
273 exps[i] = coeffs[i];
274 }
275 curTrie->insertTerm(coef, exps, 0, formSum.varCount, degree);
276
277 delete[] exps;
278 formSum.termCount++;
279 }
280
281 //Prints a nested list representation of our sum of linear forms
282 //linear form: (c_{1} / d_{1}!)[(p_{1}*x_{1} + ... p_{varCount}*x_{varCount})^d_{1}] + ...
283 //nested list: [[c_{1}, [d_{1}, [p_{1}, p_{2}, ..., p_{varCount}]], .. ]
printLinForms(const linFormSum & myForm)284 string printLinForms(const linFormSum &myForm)
285 {
286 BTrieIterator<RationalNTL, ZZ>* it = new BTrieIterator<RationalNTL, ZZ> ();
287 term<RationalNTL, ZZ>* temp;
288 it->setTrie(myForm.myForms, myForm.varCount);
289 it->begin();
290
291 stringstream output(stringstream::in | stringstream::out);
292 temp = it->nextTerm();
293 do
294 {
295 if (output.str() != "")
296 {
297 output << ", ";
298 }
299 output << "[" << temp->coef << ", [" << temp->degree << ", [";
300 for (int j = 0; j < temp->length; j++)
301 {
302 output << temp->exps[j];
303 if (j + 1 < temp->length)
304 {
305 output << ", ";
306 }
307 }
308 output << "]]]";
309 temp = it->nextTerm();
310 } while (temp);
311 delete it;
312 return "[" + output.str() + "]";
313 }
314
315 //Deallocates space and nullifies internal pointers and counters
destroyLinForms(linFormSum & myPoly)316 void destroyLinForms(linFormSum &myPoly)
317 {
318 if (myPoly.myForms)
319 delete myPoly.myForms;
320 myPoly.myForms = NULL;
321 myPoly.termCount = myPoly.varCount = 0;
322 }
323
324 // -------------------------------------------------------------------------
325
loadLinFormProducts(linFormProductSum & forms,const string line)326 void loadLinFormProducts(linFormProductSum &forms, const string line)
327 {
328 forms.varCount = 0;
329 FormProductLoadConsumer<RationalNTL>* myLoader =
330 new FormProductLoadConsumer<RationalNTL> ();
331 myLoader->setFormProductSum(forms);
332 parseLinFormProducts(myLoader, line);
333 delete myLoader;
334 }
335 //Loads a string by parsing it as a sum of linear forms
336 //linear form: a * (c^b * f^e * h^g) + ...
337 //nested list: [ [a, [[b, [c]], [e, [f]], [g, [h]]]], ... ]
parseLinFormProducts(FormProductLoadConsumer<RationalNTL> * consumer,const string & line)338 void parseLinFormProducts(FormProductLoadConsumer<RationalNTL>* consumer, const string& line)
339 {
340 int termIndex = 0;
341 int lastPos = 0;
342 int varCount = 0;
343 int k;
344 int flag = 0; //0 means we expect coefficient, 1 means we expect degree, 2 means we expect coefficient vector
345
346 //cout << "parseLinForms: line = " << line.c_str() << endl;
347 for (int i = 0; line[i] != ']'; i++)
348 {
349 varCount += (line[i] == ',');
350 }
351 //varCount is now the number of commas in a linear form - there is 1 less variable;
352 varCount--;
353 if (varCount < 1)
354 {
355 cout << "line: `" << line << "'" << endl;
356 cout << "There are " << varCount << " variables, error." << endl;
357 exit(1);
358 }
359 consumer->setDimension(varCount);
360
361 vec_ZZ coefs;
362 coefs.SetLength(varCount);
363 int degree;
364 RationalNTL coefficient;
365
366 int productIndex = 0;
367 for (int i = 1; i < line.length() - 1; i++) //ignore outermost square brackets
368 {
369 if (line[i] == '[')
370 {
371 ZZ degreeFactorial;
372 switch (flag)
373 {
374 case 0: //coefficient
375 lastPos = i + 1;
376 for (; line[i] != ','; i++)
377 ;
378 coefficient = RationalNTL(
379 (line.substr(lastPos, i - lastPos).c_str()));
380 flag = 1;
381 break;
382 case 1: //start of a product of linear forms [[p, [l]], ...],
383 productIndex = consumer->initializeNewProduct();
384 //cout << "new index:" << productIndex << endl;
385 flag = 2;
386 break;
387 case 2: //start of a power of linear form [p, [l]]
388 lastPos = i + 1;
389 for (; line[i] != ','; i++)
390 ;
391 degree = atoi(line.substr(lastPos, i - lastPos).c_str());
392 flag = 3;
393 break;
394 case 3: //coefficient vector
395 k = 0;
396 for (i++; line[i] != ']'; i++)
397 {
398 if (line[i] != ' ')
399 {
400 lastPos = i;
401 for (; line[i] != ',' && line[i] != ']'; i++)
402 ;
403
404 coefs[k++] = to_ZZ(
405 line.substr(lastPos, i - lastPos).c_str());
406
407 }
408 }//end of exponent vector.
409
410 //cout << "inserting " << coefficient << " ";
411 //for(int w = 0; w < coefs.length(); ++w)
412 // cout << coefs[w] << ", ";
413 //cout << "into " << productIndex << endl;
414 consumer->ConsumeLinForm(productIndex, coefficient, degree, coefs);
415 coefficient = 1; //the other coefficients in this product are 1.
416
417 for(; line[i] != ']'; ++i); //end of linear form list: [[p, [l]] ]
418 //for(; line[i] != ']'; ++i); //end of linear form list: [[[p, [l]] ]
419 //check to see if there are any more products.
420 while(true)
421 {
422 ++i;
423 if ( line[i] == ',')
424 {
425 flag = 2; //read next power of linear form into the current product.
426 break;
427 }
428 else if ( line[i] == ']')
429 {
430 flag = 0; //new product
431 break;
432 }
433 }
434 break;
435 default: //error
436 cout << "Flag is " << flag << ", error." << endl;
437 exit(1);
438 }
439 }
440 }
441 }
442
443
444
445
446 //Deallocates space and nullifies internal pointers and counters
destroyLinFormProducts(linFormProductSum & myProd)447 void destroyLinFormProducts(linFormProductSum &myProd)
448 {
449 for(int i = 0; i < myProd.myFormProducts.size(); ++i)
450 {
451 destroyLinForms(myProd.myFormProducts[i]);
452 }
453 myProd.myFormProducts.clear();
454 }
455
456
457 /**
458 * Unlike printLinForms(), this function is more for debugging.
459 */
printLinFormProducts(const linFormProductSum & plf)460 string printLinFormProducts(const linFormProductSum &plf)
461 {
462 stringstream out;
463 for(int i = 0; i < plf.myFormProducts.size(); ++i)
464 {
465 cout << i << " started" << endl;
466 cout << printLinForms(plf[i]).c_str() << endl;
467 out << "Term " << i << " " << printLinForms(plf[i]) << "\n";
468 cout << i << " finished" << endl;
469 }
470 return out.str();
471 }
472
473
474 // -------------------------------------------------------------------------
475
476 //INPUT: monomial specified by myPoly.coefficientBlocks[mIndex / BLOCK_SIZE].data[mIndex % BLOCK_SIZE]
477 // and myPoly.exponentBlocks[mIndex / BLOCK_SIZE].data[mIndex % BLOCK_SIZE]
478 //OUTPUT: lForm now also contains the linear decomposition of this monomial
479 // note: all linear form coefficients assumed to be divided by their respective |M|!, and the form is assumed to be of power M
decompose(BTrieIterator<RationalNTL,int> * it,linFormSum & lForm)480 void decompose(BTrieIterator<RationalNTL, int>* it, linFormSum &lForm)
481 {
482 //cout << "decomposing " << lForm.varCount << " variables" << endl;
483 term<RationalNTL, int>* temp;
484 //BTrieIterator<ZZ, int>* it = new BTrieIterator<ZZ, int>();
485
486
487 //it->setTrie(myPoly.myMonomials, myPoly.varCount);
488 it->begin();
489
490 //cout << "in decompost()::\n";
491 temp = it->nextTerm();
492
493 do
494 {
495 //cout << "monomial " << temp->coef;
496 //for(int i = 0; i < temp->length; ++i)
497 // cout << temp->exps[i];
498 //cout << "\n";
499 decompose(temp, lForm);
500 temp = it->nextTerm();
501 } while (temp);
502
503 //delete myPoly.myMonomials;
504 }
505
decompose(term<RationalNTL,int> * myTerm,linFormSum & lForm)506 void decompose(term<RationalNTL, int>* myTerm, linFormSum &lForm)
507 {
508 vec_ZZ myExps;
509 myExps.SetLength(lForm.varCount);
510
511 //April 27. 2011 Brandon: I don't think this if is ever true. If I insert [[[1,[0,0]]], I get a monomial of length 2 still.
512 //This or the "string to polynomial" function is not correct. I think we can delete this if statement.
513 if (myTerm->length == 0) //constant
514 {
515 for (int j = 0; j < lForm.varCount; j++)
516 {
517 myExps[j] = 0;
518 }
519 insertLinForm(myTerm->coef, 0, myExps, lForm);
520 return;
521 }
522
523
524 ZZ formsCount = to_ZZ(myTerm->exps[0] + 1); //first exponent
525 int totalDegree = myTerm->exps[0];
526 for (int i = 1; i < myTerm->length; i++)
527 {
528 formsCount *= myTerm->exps[i] + 1;
529 totalDegree += myTerm->exps[i];
530 }
531 formsCount--;
532
533 //If this is zero, then the term is a constant [c,[0,0,0,0...0]]
534 if ( formsCount == 0)
535 {
536 for (int j = 0; j < lForm.varCount; j++)
537 {
538 myExps[j] = 0;
539 }
540 insertLinForm(myTerm->coef, 0, myExps, lForm);
541 return;
542 }//if formsCount
543
544 //cout << "At most " << formsCount << " linear forms will be required for this decomposition." << endl;
545 //cout << "Total degree is " << totalDegree << endl;
546
547 int* p = new int[lForm.varCount];
548 int* counter = new int[lForm.varCount];
549 ZZ* binomCoeffs = new ZZ[lForm.varCount]; //for calculating the product of binomial coefficients for each linear form
550
551 RationalNTL temp;
552 int g;
553 int myIndex;
554 for (int i = 0; i < lForm.varCount; i++)
555 {
556 counter[i] = 0;
557 binomCoeffs[i] = to_ZZ(1);
558 }
559 for (ZZ i = to_ZZ(1); i <= formsCount; i++)
560 {
561 //cout << "i is " << i << endl;
562 counter[0] += 1;
563 for (myIndex = 0; counter[myIndex] > myTerm->exps[myIndex]; myIndex++)
564 {
565 counter[myIndex] = 0;
566 binomCoeffs[myIndex] = to_ZZ(1);
567 counter[myIndex + 1] += 1;
568 }
569 binomCoeffs[myIndex] *= myTerm->exps[myIndex] - counter[myIndex] + 1; // [n choose k] = [n choose (k - 1) ] * (n - k + 1)/k
570 binomCoeffs[myIndex] /= counter[myIndex];
571 //cout << "counter is: " << counter << endl;
572
573 //find gcd of all elements and calculate the product of binomial coefficients for the linear form
574 g = counter[0];
575 int parity = totalDegree - counter[0];
576 p[0] = counter[0];
577 temp = binomCoeffs[0];
578 if (lForm.varCount > 1)
579 {
580 for (int k = 1; k < lForm.varCount; k++)
581 {
582 p[k] = counter[k];
583 g = GCD(g, p[k]);
584 parity -= p[k];
585 temp *= binomCoeffs[k];
586 }
587 }
588
589 //calculate coefficient
590 temp *= myTerm->coef;
591 if ((parity % 2) == 1)
592 {
593 temp *= to_ZZ(-1);
594 } // -1 ^ [|M| - (p[0] + p[1] + ... p[n])], checks for odd parity using modulo 2
595
596 if (g != 1)
597 {
598 for (int k = 0; k < lForm.varCount; k++)
599 {
600 p[k] /= g;
601 }
602 temp *= power_ZZ(g, totalDegree);
603 }
604 //cout << "coefficient is " << temp << endl;
605 for (int i = 0; i < lForm.varCount; i++)
606 {
607 myExps[i] = p[i];
608 }
609 insertLinForm(temp, totalDegree, myExps, lForm);
610 }
611 delete[] p;
612 delete[] counter;
613 delete[] binomCoeffs;
614 }
615