1 /*
2 * testLibraryExamples.cpp
3 *
4 * Created on: Jul 12, 2011
5 * Author: bedutra
6 * This files contains examples on how to use latte functions in a
7 * library setting.
8 *
9 * We will always be working with the 3-polytope that is a
10 * square pyramid sitting on top of a cube for the examples.
11 */
12
13 #include <iostream>
14 #include <fstream>
15 #include <sstream>
16 #include <cstdlib>
17 #include <cassert>
18
19 #include "valuation.h"
20 #include "count.h"
21 #include "latte_cddlib.h"
22
23
24 //global option on how to build the polytopes.
25 int USE_FILES = 0;
26 bool PRINT_POLYHEDRA = false ;
27
28
29 /*****************************************************************************/
30 /********** start of function headers. Nothing interesting here **************/
31 /*****************************************************************************/
32
33
34 using namespace std;
35
36 typedef enum
37 {
38 vrepToCones, vrepToVertices, hrepToCones, hrepToVertices
39 } ReadProblemType;
40
41 //write h or v files
42 void building_h_file();
43 void building_v_file();
44
45 //Reading the polyhedron
46 Polyhedron * building_polyhedron(const ReadProblemType readProblemType,
47 BarvinokParameters * parms);
48 Polyhedron * building_polyhedron_with_files(const ReadProblemType readProblemType,
49 BarvinokParameters * parms);
50 Polyhedron * building_polyhedron_without_files(const ReadProblemType readProblemType,
51 BarvinokParameters * parms);
52
53
54 //finding volumes and integrals.
55 void finding_integral_polynomial_coneDecomp();
56 void finding_integral_linearForm_triangulation();
57 void finding_volume_coneDecomp();
58 void finding_volume_triangulation();
59
60 //counting
61 void finding_lattice_count();
62 void finding_ehrhart_polynomial();
63 void finding_ehrhart_taylor();
64
65 //helper functions to make the example polytope and integrand.
66 dd_MatrixPtr getHrepMatrix();
67 dd_MatrixPtr getVrepMatrix();
68 linFormSum getLinearForms();
69 monomialSum getPolynomial();
70
71 //A B C D E F G H I J K L M N O P Q R S T U V W X Y Z.
72
73
74 string FILE_BASE;
75
main(int argc,char * argv[])76 int main(int argc, char *argv[])
77 {
78 /* //feel free to comment this back in if you want.
79 cout << "Hello, this file is designed to help developers use LattE as\n"
80 << " a library. I recommend commenting out the examples you do\n"
81 << " not care to see and un-commenting the print statements\n"
82 << " to the examples you want to see\n" << endl;
83
84 cout << "Press enter. Lots of things will be printed to the screen\n"
85 << " depending on how many examples you have compiled." << endl;
86 cin.get(); //stop, so user can read the message above.
87
88 */
89
90 //global variable. Used to make temp. files.
91 FILE_BASE = argv[0];
92
93
94 /* //I should delete this block
95 int file, num;
96 for(cin >> file >> num; file != -1; cin >> file >> num)
97 {
98 BarvinokParameters *params = new BarvinokParameters;
99 USE_FILES = file;
100 switch(num)
101 {
102 case 1:
103 cout << "doing h to cones" << endl;
104 building_polyhedron(hrepToCones, params);
105 break;
106 case 2:
107 cout << "doing h to v" << endl;
108 building_polyhedron(hrepToVertices, params);
109 break;
110 case 3:
111 cout << "doing v to cones" << endl;
112 building_polyhedron(vrepToCones, params);
113 break;
114 case 4:
115 cout << "doing v to v" << endl;
116 building_polyhedron(vrepToVertices, params);
117 break;
118 }
119 delete params;
120 }
121 */
122
123
124
125 //Reading in a polytope. The "easy" way.
126 //Note that these "building" functions are leaking memory because I
127 // am not deleting the returned polytopes.
128 BarvinokParameters *params = new BarvinokParameters;
129 building_polyhedron(hrepToCones, params);
130 building_polyhedron(hrepToVertices, params);
131 building_polyhedron(vrepToCones, params);
132 building_polyhedron(vrepToVertices, params);
133
134
135
136
137 //volume examples
138 finding_volume_coneDecomp();
139 finding_volume_triangulation();
140
141 //integral examples.
142 finding_integral_polynomial_coneDecomp();
143 finding_integral_linearForm_triangulation();
144
145 //counting examples. Note that the rational functions are
146 // written to the .rat file.
147 finding_lattice_count();
148 finding_ehrhart_polynomial();
149 finding_ehrhart_taylor();
150
151 cout << "ALL TESTS PASSED" << endl;
152 return 0;
153 }//main()
154
155
building_h_file()156 void building_h_file()
157 {
158 string hFile = FILE_BASE + ".hrep";
159 dd_MatrixPtr M;
160
161 M = getHrepMatrix();
162 fstream hOut;
163 hOut.open(hFile.c_str(), fstream::out);
164 hOut << M->rowsize << " " << M->colsize << endl;
165 for (int i = 0; i < M->rowsize; ++i)
166 {
167 for (int j = 0; j < M->colsize; ++j)
168 {
169 hOut << M->matrix[i][j] << " ";
170 }
171 hOut << endl;
172 }
173 hOut.close();
174
175 dd_FreeMatrix(M);
176
177 }///building_h_file
178
building_v_file()179 void building_v_file()
180 {
181 string vFile = FILE_BASE + ".vrep";
182 dd_MatrixPtr M;
183
184 //make the v file
185 M = getVrepMatrix();
186 fstream vOut;
187 vOut.open(vFile.c_str(), fstream::out);
188 vOut << M->rowsize << " " << M->colsize << endl;
189 for (int i = 0; i < M->rowsize; ++i)
190 {
191 for (int j = 0; j < M->colsize; ++j)
192 vOut << M->matrix[i][j] << " ";
193 vOut << endl;
194 }
195 vOut.close();
196
197 dd_FreeMatrix(M);
198
199 }//building_v_file()
200
201 /**
202 * The easy way to have latte construct a polyhedron from your v or h
203 * representation is to write it to a file.
204 *
205 *
206 */
building_polyhedron(const ReadProblemType readProblemType,BarvinokParameters * params)207 Polyhedron * building_polyhedron(const ReadProblemType readProblemType,
208 BarvinokParameters * params)
209 {
210
211 if( USE_FILES == 1 )
212 return building_polyhedron_with_files(readProblemType, params);
213 else
214 return building_polyhedron_without_files(readProblemType, params);
215
216 }
217
building_polyhedron_with_files(const ReadProblemType readProblemType,BarvinokParameters * params)218 Polyhedron * building_polyhedron_with_files(const ReadProblemType readProblemType,
219 BarvinokParameters * params)
220 {
221 string hFile = FILE_BASE + ".hrep";
222 string vFile = FILE_BASE + ".vrep";
223
224 ReadPolyhedronData readPolyData;
225 Polyhedron *Poly;
226
227 if (readProblemType == hrepToCones || readProblemType == hrepToVertices)
228 {
229 building_h_file();
230 } else
231 {
232 building_v_file();
233 }
234
235
236 //start from the h file, find the vertex representation.
237 if (readProblemType == hrepToVertices)
238 {
239
240 strcpy(readPolyData.dualApproach, "yes"); //compute the vertices
241
242 //note that at this point, you can set additional options.
243 //readPolyData.parse_option("--cdd"); //if the input is in a cdd-style file.
244 //readPolyData.parse_option("--compute-vertex-cones=cdd"); //find the tangent cones with cdd
245 //readPolyData.parse_option("--compute-vertex-cones=4ti2");// or 4ti2
246 //readPolyData.parse_option("--redundancy-check=none"); //check h-rep for redundancies or not
247 //readPolyData.parse_option("--redundancy-check=cddlib");
248 readPolyData.parse_option("--redundancy-check=full-cddlib");
249 //readPolyData.parse_option("--vrep"); //if the input file is a latte style v-rep file.
250
251 readPolyData.parse_option(hFile.c_str()); //set the file name last.
252
253
254 //params->max_determinant = 1; //if you want unimodular cones.
255
256 #ifdef HAVE_FORTYTWO_LIB
257 parse_standard_triangulation_option("--triangulation=4ti2", params);
258 parse_standard_dualization_option("--dualization=4ti2", params);
259 #endif
260
261 //ok that's all the set up we need. Read in the polyhedron.
262 Poly = readPolyData.read_polyhedron(params);
263
264 params->Number_of_Variables = Poly->numOfVars;
265 params->File_Name = (char *) readPolyData.filename.c_str();
266
267 //hey, the next print statement shows we don't have vertices. We have the dual h-rep.
268 if ( PRINT_POLYHEDRA )
269 Poly->printPolyhedron();
270 // homogenized: true: cones represent the homogenization of the polyhedron.
271 // false:cones represent the supporting cones of all vertices of the polyhedron.
272 // dualized: if the cone is dualized.
273 // NOTE: In the Polyhedron class, vertices are saved with the 'leanding 1' at the end (v, 1)
274 // For hyperplanes, they are saved as 'ax - b < 0'. This is opposite for how they are
275 // saved in the v- and h-rep files.
276
277 //to fix this, dualize the cone again.
278 dualizeCones(Poly->cones, Poly->numOfVars, params);
279 Poly->dualized = false;
280 if ( PRINT_POLYHEDRA )
281 Poly->printPolyhedron();
282
283 }
284 //start from the h file, find the tangent cones
285 else if (readProblemType == hrepToCones)
286 {
287 strcpy(readPolyData.dualApproach, "no"); //compute the cones
288
289 readPolyData.parse_option("--redundancy-check=full-cddlib");
290 readPolyData.parse_option(hFile.c_str()); //set the file name last.
291
292 #ifdef HAVE_FORTYTWO_LIB
293 parse_standard_triangulation_option("--triangulation=4ti2", params);
294 parse_standard_dualization_option("--dualization=4ti2", params);
295 #endif
296
297 //ok that's all the set up we need. Read in the polyhedron.
298 Poly = readPolyData.read_polyhedron(params);
299
300 params->Number_of_Variables = Poly->numOfVars;
301 params->File_Name = (char *) readPolyData.filename.c_str();
302
303 //hey, the next print statement shows we don't have the rays of the tangent cones,
304 // we have the facets of the tangent cones.
305 if ( PRINT_POLYHEDRA )
306 Poly->printPolyhedron();
307 //to fix this, dualize the cone again to compute the rays.
308 //but then the cone is dualized. So dualize back (it just swaps the ray/facet lists).
309 dualizeCones(Poly->cones, Poly->numOfVars, params);
310 dualizeCones(Poly->cones, Poly->numOfVars, params);
311 if ( PRINT_POLYHEDRA )
312 Poly->printPolyhedron();
313 }
314 //start from the v-rep file, find the vertices.
315 else if (readProblemType == vrepToVertices)
316 {
317 strcpy(readPolyData.dualApproach, "yes"); //compute the vertices
318
319 readPolyData.parse_option(vFile.c_str()); //set the file name last.
320 readPolyData.parse_option("--vrep"); //the input file is a latte style v-rep file.
321
322 #ifdef HAVE_FORTYTWO_LIB
323 parse_standard_triangulation_option("--triangulation=4ti2", params);
324 parse_standard_dualization_option("--dualization=4ti2", params);
325 #endif
326
327 //ok that's all the set up we need. Read in the polyhedron.
328 Poly = readPolyData.read_polyhedron(params);
329
330 params->Number_of_Variables = Poly->numOfVars;
331 params->File_Name = (char *) readPolyData.filename.c_str();
332
333 //the polyhedron is in the form you think it is!
334 if ( PRINT_POLYHEDRA )
335 Poly->printPolyhedron();
336
337 } else if (readProblemType == vrepToCones)
338 {
339
340 strcpy(readPolyData.dualApproach, "no"); //compute the cones
341
342 readPolyData.parse_option(vFile.c_str()); //set the file name last.
343 readPolyData.parse_option("--vrep"); //the input file is a latte style v-rep file.
344
345 #ifdef HAVE_FORTYTWO_LIB
346 parse_standard_triangulation_option("--triangulation=4ti2", params);
347 parse_standard_dualization_option("--dualization=4ti2", params);
348 #endif
349
350 //ok that's all the set up we need. Read in the polyhedron.
351 Poly = readPolyData.read_polyhedron(params);
352
353 params->Number_of_Variables = Poly->numOfVars;
354 params->File_Name = (char *) readPolyData.filename.c_str();
355
356 //the polyhedron is in the form you think it is!
357 if ( PRINT_POLYHEDRA )
358 Poly->printPolyhedron();
359
360 } else
361 {
362 cerr << "unknown type" << endl;
363 exit(1);
364 }
365
366 freeListVector(readPolyData.templistVec);
367 freeListVector(readPolyData.matrix);
368 return Poly;
369 }//building_polyhedron()
370
371
building_polyhedron_without_files(const ReadProblemType readProblemType,BarvinokParameters * params)372 Polyhedron * building_polyhedron_without_files(const ReadProblemType readProblemType,
373 BarvinokParameters * params)
374 {
375 string hFile = FILE_BASE + ".hrep";
376 string vFile = FILE_BASE + ".vrep";
377
378 ReadPolyhedronData readPolyData;
379 Polyhedron *Poly;
380
381 //start from the h file, find the vertex representation.
382 if (readProblemType == hrepToVertices)
383 {
384 dd_MatrixPtr M;
385 M = getHrepMatrix();
386 //note that at this point, you can set additional options.
387 //readPolyData.parse_option("--compute-vertex-cones=cdd"); //find the tangent cones with cdd
388 //readPolyData.parse_option("--compute-vertex-cones=4ti2");// or 4ti2
389 //readPolyData.parse_option("--redundancy-check=none"); //check h-rep for redundancies or not
390 //readPolyData.parse_option("--redundancy-check=cddlib");
391 readPolyData.parse_option("--redundancy-check=full-cddlib");
392
393 #ifdef HAVE_FORTYTWO_LIB
394 parse_standard_triangulation_option("--triangulation=4ti2", params);
395 parse_standard_dualization_option("--dualization=4ti2", params);
396 #endif
397
398 //ok that's all the set up we need. Read in the polyhedron.
399 Poly = readPolyData.read_polyhedron(M, params, ReadPolyhedronData::computeVertices);
400
401 if ( PRINT_POLYHEDRA )
402 Poly->printPolyhedron();
403 }
404 //start from the h file, find the tangent cones
405 else if (readProblemType == hrepToCones)
406 {
407 dd_MatrixPtr M;
408 M = getHrepMatrix();
409 readPolyData.parse_option("--redundancy-check=full-cddlib");
410
411 #ifdef HAVE_FORTYTWO_LIB
412 parse_standard_triangulation_option("--triangulation=4ti2", params);
413 parse_standard_dualization_option("--dualization=4ti2", params);
414 #endif
415
416 //ok that's all the set up we need. Read in the polyhedron.
417 Poly = readPolyData.read_polyhedron(M, params, ReadPolyhedronData::computePrimalCones);
418
419 if ( PRINT_POLYHEDRA )
420 Poly->printPolyhedron();
421 }
422 //start from the v-rep file, find the vertices.
423 else if (readProblemType == vrepToVertices)
424 {
425 dd_MatrixPtr M;
426 M = getVrepMatrix();
427 readPolyData.parse_option("--redundancy-check=full-cddlib");
428
429 #ifdef HAVE_FORTYTWO_LIB
430 parse_standard_triangulation_option("--triangulation=4ti2", params);
431 parse_standard_dualization_option("--dualization=4ti2", params);
432 #endif
433
434 //ok that's all the set up we need. Read in the polyhedron.
435 Poly = readPolyData.read_polyhedron(M, params, ReadPolyhedronData::computeVertices);
436
437 if ( PRINT_POLYHEDRA )
438 Poly->printPolyhedron();
439 } else if (readProblemType == vrepToCones)
440 {
441 dd_MatrixPtr M;
442 M = getVrepMatrix();
443 readPolyData.parse_option("--redundancy-check=full-cddlib");
444
445 #ifdef HAVE_FORTYTWO_LIB
446 parse_standard_triangulation_option("--triangulation=4ti2", params);
447 parse_standard_dualization_option("--dualization=4ti2", params);
448 #endif
449
450 //ok that's all the set up we need. Read in the polyhedron.
451 Poly = readPolyData.read_polyhedron(M, params, ReadPolyhedronData::computePrimalCones);
452
453 if ( PRINT_POLYHEDRA )
454 Poly->printPolyhedron();
455 } else
456 {
457 cerr << "unknown type" << endl;
458 exit(1);
459 }
460
461 freeListVector(readPolyData.templistVec);
462 freeListVector(readPolyData.matrix);
463
464
465 return Poly;
466 }//building_polyhedron()
467
468
469
470
471
472 /**
473 * Integrate a polynomial over a polytope using the cone decomposition method.
474 *
475 * Integrating via the triangulation method is very similar. I will not give
476 * an explicit example using the triangulation method on polynomials.
477 *
478 */
finding_integral_polynomial_coneDecomp()479 void finding_integral_polynomial_coneDecomp()
480 {
481 RationalNTL ans;
482 BarvinokParameters *params = new BarvinokParameters;
483 Polyhedron * poly = building_polyhedron(hrepToCones, params);
484 PolytopeValuation polytopeValuation(poly, *params);
485 monomialSum polynomial = getPolynomial();
486
487 //the cone-decomposition method needs to start from tangent-cones
488 ans = polytopeValuation.findIntegral(polynomial,
489 PolytopeValuation::integratePolynomialAsLinearFormCone);
490 //If a dilation is needed, the polyhedron will now be dilated.
491 //If you wanted to use the triangulation method, you just need to preplace
492 // LawrenceIntegration with TriangulationIntegration. Also, you could
493 // ask for a v-rep by replacing hrepToCones with vrepToVertices because
494 // the triangulation method can start from both polyhedron descriptions.
495
496 cout << "Integral using the cone decomposition method is \n"
497 << "Rational: " << ans << '\n' << "Real : " << ans.to_RR()
498 << endl;
499
500 assert( ans == RationalNTL("734124064/2079"));
501
502 delete poly;
503 delete params;
504 }
505
506 /**
507 * Computes the integral of powers of linear forms using the triangulation
508 * method.
509 */
finding_integral_linearForm_triangulation()510 void finding_integral_linearForm_triangulation()
511 {
512 RationalNTL ans;
513 BarvinokParameters * params = new BarvinokParameters;
514 Polyhedron * poly = building_polyhedron(hrepToVertices, params);
515 PolytopeValuation polytopeValuation(poly, *params);
516 linFormSum lForms = getLinearForms();
517
518 ans = polytopeValuation.findIntegral(lForms,
519 PolytopeValuation::integrateLinearFormTriangulation);
520
521 cout << "Integral using the triangulation method is \n" << "Rational: "
522 << ans << '\n' << "Real : " << ans.to_RR() << endl;
523
524 //note the string, number too big for an int.
525 assert( ans == RationalNTL("2448974416/15"));
526
527 delete poly;
528 delete params;
529 }//finding_integral_polynomial_coneDecomp
530
531
532 /**
533 * Computes the volume using the cone decomposition method.
534 */
finding_volume_coneDecomp()535 void finding_volume_coneDecomp()
536 {
537 RationalNTL ans;
538 BarvinokParameters *params = new BarvinokParameters;
539
540 Polyhedron * poly = building_polyhedron(vrepToCones, params);
541
542 PolytopeValuation polytopeValuation(poly, *params);
543
544
545 //the cone-decomposition method needs to start from tangent-cones
546 ans = polytopeValuation.findVolume(PolytopeValuation::volumeCone);
547 //if a dilation is needed, the polyhedron will now be dilated.
548
549 cout << "Volume using the cone decomposition method is \n" << "Rational: "
550 << ans << '\n' << "Real : " << ans.to_RR() << endl;
551
552 assert( ans == RationalNTL(208,3));
553
554 delete poly;
555 delete params;
556 }//finding_volume_coneDecomp
557
558
559 /**
560 * Computes the volume using the triangulation method.
561 */
finding_volume_triangulation()562 void finding_volume_triangulation()
563 {
564 RationalNTL ans;
565 BarvinokParameters *params = new BarvinokParameters;
566 Polyhedron * poly = building_polyhedron(vrepToVertices, params);
567 PolytopeValuation polytopeValuation(poly, *params);
568
569 //note, the triangulation method could start from the tangent-cones or a v-rep.
570 ans = polytopeValuation.findVolume(PolytopeValuation::volumeTriangulation);
571 //the original polytope is not lost if we started from a v-rep.
572
573 cout << "Volume using the triangulation method is \n" << "Rational: "
574 << ans << '\n' << "Real : " << ans.to_RR() << endl;
575
576 assert( ans == RationalNTL(208,3));
577
578 delete poly;
579 delete params;
580 }//finding_volume_triangulation
581
582
583 /**
584 * Computes the number of lattice points in the polytope.
585 * Unlike in the other examples, we cannot start from the Polhedron, we have
586 * to start from the command line!
587 *
588 * The improvement to using this method to find the count instead of simply
589 * calling system with the LattE count command is that the output is directly
590 * returned to you here. If you where to call a system command, you would have
591 * to parse the output which could case errors if the LattE output every changes.
592 */
finding_lattice_count()593 void finding_lattice_count()
594 {
595 string fileName = FILE_BASE + ".hrep";
596 building_h_file();
597
598 char * argv[] = { "./dummyName", "--count-lattice-points", "" };
599 argv[2] = (char *) fileName.c_str();
600
601 CountAnswerContainer ans;
602 ans = mainCountDriver(3, argv);
603
604 cout << "The number of lattice points is " << ans.numLaticePoints << endl;
605
606 assert( ans.numLaticePoints == to_ZZ(126));
607 }//finding_lattice_count
608
609 /**
610 * Finds the ehrhart polynomial
611 */
finding_ehrhart_polynomial()612 void finding_ehrhart_polynomial()
613 {
614 string fileName = FILE_BASE + ".hrep";
615 building_h_file();
616
617 char * argv[] = { "./dummyName", "--ehrhart-polynomial", "" };
618 argv[2] = (char *) fileName.c_str();
619
620 CountAnswerContainer ans;
621 ans = mainCountDriver(3, argv);
622
623 const char * ehrhartPoly[4] = { "1", "35/3", "44", "208/3" };
624
625 for (int i = 0; i < 4; ++i)
626 assert( ans.ehrhart_coefficients[i] == mpq_class(ehrhartPoly[i]));
627
628 }//finding_ehrhart_polynomial
629
630
631 /*
632 * Finds the first 5 terms in the taylor expansion.
633 */
finding_ehrhart_taylor()634 void finding_ehrhart_taylor()
635 {
636 string fileName = FILE_BASE + ".hrep";
637 building_h_file();
638
639 char * argv[] = { "./dummyName", "--ehrhart-taylor=5", "" };
640 argv[2] = (char *) fileName.c_str();
641
642 CountAnswerContainer ans;
643 ans = mainCountDriver(3, argv);
644
645 const char * taylorPoly[6] = { "1", "126", "755", "2304", "5189t", "9826" };
646
647 for (int i = 0; i < 6; ++i)
648 assert( ans.seriesExpansion[i] == to_ZZ(taylorPoly[i]));
649
650 }//finding_ehrhart_taylor
651
652 /**
653 * Simple function to set up a dd_matrix of facets.
654 */
getHrepMatrix()655 dd_MatrixPtr getHrepMatrix()
656 {
657 int h[][9] = { { 0, 1, 0, 0 }, { 0, 0, 1, 0 }, { 0, 0, 0, 1 }, { 4, -1, 0,
658 0 }, { 4, 0, -1, 0 }, { 32, 4, 0, -8 }, { 48, 0, -4, -8 }, { 48,
659 -4, 0, -8 }, { 32, 0, 4, -8 } }; //latte-style h rep (0 <= b - Ax)
660
661
662 dd_MatrixPtr matrix =
663 dd_CreateMatrix(9 /*num rows*/, 4 /*num of vars + 1*/);
664
665 matrix->representation = dd_Inequality; //h-rep.
666
667 for (int i = 0; i < 9; ++i)
668 for (int j = 0; j < 4; ++j)
669 dd_set(matrix->matrix[i][j], mpq_class(h[i][j]).get_mpq_t());
670
671 return matrix;
672 }//getHrepMatrix
673
674 /**
675 * Create a v-rep matrix for the polytope with vertices
676 * 0 0 0
677 * 4 0 0
678 * 0 4 0
679 * 4 4 0
680 * 0 0 4
681 * 4 0 4
682 * 0 4 4
683 * 4 4 4
684 * (the above is a 3-cube)
685 * and
686 * (2,2,5)
687 * (the above is a square pyramid on top of a cube).
688 */
getVrepMatrix()689 dd_MatrixPtr getVrepMatrix()
690 {
691 int v[][4] = { { 1, 0, 0 }, { 1, 4, 0 }, { 1, 0, 4 }, { 1, 4, 4 } }; //I'm to lazy to fill up the matrix one element at a time, so save a part of the vertices here.
692
693 dd_MatrixPtr matrix =
694 dd_CreateMatrix(9 /*num vertices*/, 4 /*num of vars + 1*/);
695
696 matrix->representation = dd_Generator; //v-rep.
697
698 for (int z = 0; z <= 1; ++z)
699 {
700 for (int i = 0; i < 4; ++i)
701 {
702 for (int j = 0; j < 3; ++j)
703 dd_set(matrix->matrix[4 * z + i][j],
704 mpq_class(v[i][j]).get_mpq_t());
705 dd_set(matrix->matrix[4 * z + i][3], mpq_class(z * 4).get_mpq_t());
706 }
707 }
708 dd_set(matrix->matrix[8][0], mpq_class(1).get_mpq_t());
709 dd_set(matrix->matrix[8][1], mpq_class(2).get_mpq_t());
710 dd_set(matrix->matrix[8][2], mpq_class(2).get_mpq_t());
711 dd_set(matrix->matrix[8][3], mpq_class(5).get_mpq_t());
712
713 return matrix;
714 }//getVrepMatrix
715
716
717 /**
718 * Returns the powers of linear forms
719 * 1/2(1x +2y + 3z)^2 + (30x +20y + 9z)^3 + 1
720 */
getLinearForms()721 linFormSum getLinearForms()
722 {
723 linFormSum lform;
724 lform.varCount = 3;
725 vec_ZZ l;
726 l.SetLength(3);
727
728 FormLoadConsumer<RationalNTL> loader;
729 loader.setFormSum(lform);
730
731 //the coefficient must be adjusted by the (degree!)
732
733 l[0] = 1;
734 l[1] = 2;
735 l[2] = 3;
736 loader.ConsumeLinForm(RationalNTL(2, 2), 2, l);// add 1/2(1x +2y + 3z)^2
737 //2!=2
738
739 l[0] = 30;
740 l[1] = 20;
741 l[2] = 9;
742 loader.ConsumeLinForm(RationalNTL(6, 1), 3, l);// add (30x +20y + 9z)^3
743 //3!=6.
744
745 l[0] = 8;
746 l[1] = 11;
747 l[2] = 123456789;
748 loader.ConsumeLinForm(RationalNTL(1, 1), 0, l);// add 1 = 1*( anything)^0
749
750
751 return lform;
752 }//getLinearForms
753
754 /**
755 * Returns the polynomial
756 * 1/4 x^2y^3z^5 + 3/4 x^2y^2z^2 + 1/9
757 */
getPolynomial()758 monomialSum getPolynomial()
759 {
760 monomialSum polynomial;
761 polynomial.varCount = 3;
762 int exp[3];
763
764 MonomialLoadConsumer<RationalNTL> loader;
765 loader.setMonomialSum(polynomial);
766
767 exp[0] = 2;
768 exp[1] = 3;
769 exp[2] = 5;
770 loader.ConsumeMonomial(RationalNTL(1, 4), exp); //add 1/4 x^2y^3z^5
771
772 exp[0] = 2;
773 exp[1] = 2;
774 exp[2] = 2;
775 loader.ConsumeMonomial(RationalNTL(3, 4), exp); //add 3/4 x^2y^2z^2
776
777 exp[0] = 0;
778 exp[1] = 0;
779 exp[2] = 0;
780 loader.ConsumeMonomial(RationalNTL(1, 9), exp); //add 1/9
781
782 return polynomial;
783 }//getPolynomial
784