1 /* barvinok.cpp -- Barvinok's decomposition of a cone.
2
3 Copyright 2002, 2003 Ruriko Yoshida
4 Copyright 2006, 2007 Matthias Koeppe
5
6 This file is part of LattE.
7
8 LattE is free software; you can redistribute it and/or modify it
9 under the terms of the version 2 of the GNU General Public License
10 as published by the Free Software Foundation.
11
12 LattE is distributed in the hope that it will be useful, but
13 WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with LattE; if not, write to the Free Software Foundation,
19 Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
20 */
21
22 #include <list>
23 #include <vector>
24
25 #include <fstream>
26 #include <cstdlib>
27 #include <cstring>
28 #include <cassert>
29 #include <math.h>
30 #include <time.h>
31
32 #include "ComputeOmega.h"
33 #include "barvinok.h"
34 #include "../ramon.h"
35 #include "../RudyResNTL.h"
36 #include "rational.h"
37 #include "convert.h"
38 #include "dual.h"
39 #include "config.h"
40 #include "Irrational.h"
41 #include "triangulation/triangulate.h"
42 #ifdef HAVE_EXPERIMENTS
43 #include "barvinok/SubspaceAvoidingDecomposition.h"
44 #endif
45
46 #undef SHOWDETS
47 #undef SHOWCONES
48 //#define SHOWDETS 1
49 //#define SHOWCONES
50 #ifdef SHOWCONES
51 #include "print.h"
52 #endif
53
BarvinokParameters()54 BarvinokParameters::BarvinokParameters() :
55 substitution(PolynomialSubstitution),
56 decomposition(DualDecomposition),
57 triangulation(
58 #ifdef HAVE_FORTYTWO_LIB
59 RegularTriangulationWith4ti2
60 #else
61 RegularTriangulationWithCdd
62 #endif
63 ),
64 triangulation_max_height(10000),
65 triangulation_bias(-1),
66 nonsimplicial_subdivision(false),
67 triangulation_special_cone(NULL),
68 triangulation_prescribed_height_data(NULL),
69 triangulation_assume_fulldim(true),
70 dualization(
71 #ifdef HAVE_FORTYTWO_LIB
72 DualizationWith4ti2
73 #else
74 DualizationWithCdd
75 #endif
76 ),
77 shortvector(LatteLLL),
78 smithform(
79 #ifdef HAVE_LIDIA
80 LidiaSmithForm
81 #else
82 IlioSmithForm
83 #endif
84 ),
85 max_determinant(0),
86 File_Name(NULL),
87 Number_of_Variables(0),
88 Flags(0),
89 Cone_Index(0),
90 total_time("Total time", true),
91 read_time("Time for reading and preprocessing"),
92 vertices_time("Time for computing vertices and supporting cones"),
93 irrationalize_time("Time for irrationalizing general cones"),
94 dualize_time("Time for dualizing general cones"),
95 triangulate_time("Time for triangulating cones into simplicial cones"),
96 decompose_time("Time for Barvinok decomposition and residue calculation"),
97 num_triangulations_with_trivial_heights(0),
98 num_triangulations_with_dependent_heights(0),
99 num_triangulations(0)
100 {}
101
~BarvinokParameters()102 BarvinokParameters::~BarvinokParameters()
103 {}
104
print_statistics(ostream & s)105 void BarvinokParameters::print_statistics(ostream &s)
106 {
107 s << read_time << vertices_time << irrationalize_time << dualize_time
108 << triangulate_time << decompose_time << total_time;
109 }
110
deliver_number_of_lattice_points(const ZZ & number)111 void BarvinokParameters::deliver_number_of_lattice_points(const ZZ &number)
112 {
113 cerr << endl << "**** The number of lattice points is: "; cerr.flush();
114 cout << number; cout.flush();
115 cerr << " ****" << endl; cerr.flush();
116 cout << endl;
117 ofstream out("numOfLatticePoints");
118 out << number << endl;
119 }
120
print_statistics(ostream & s)121 void Single_Cone_Parameters::print_statistics(ostream &s)
122 {
123 BarvinokParameters::print_statistics(s);
124 s << "Total number of simplicial cones: " <<
125 Total_Simplicial_Cones << endl;
126 if (max_determinant != 0) {
127 s << "Total number of "
128 << (max_determinant == 1 ? "unimodular" : "low-index")
129 << " cones: " << Total_Uni_Cones << endl;
130 }
131 s << "Maximum depth of the decomposition tree: " << Max_Depth << endl;
132 }
133
SetConsumer(ConeConsumer * a_consumer)134 void DelegatingSingleConeParameters::SetConsumer(ConeConsumer *a_consumer)
135 {
136 consumer = a_consumer;
137 }
138
ConsumeCone(listCone * cone)139 int DelegatingSingleConeParameters::ConsumeCone(listCone *cone)
140 {
141 assert(consumer != NULL);
142 return consumer->ConsumeCone(cone);
143 }
144
145 /* Note: We are dealing with the "Row space" of the
146 input matrix due to NTL. */
147
148 /**********************************************************************/
CheckOmega(const mat_ZZ & U,vec_ZZ & Z)149 vec_ZZ CheckOmega( const mat_ZZ & U, vec_ZZ & Z){
150
151 int m;
152 m = U.NumCols();
153 mat_ZZ Dummy;
154 Dummy.SetDims(m + 1, m);
155 Dummy[0] = Z;
156 ZZ d;
157
158 for(int i = 0; i < m; i++)
159 Dummy[i + 1] = U[i];
160
161 mat_ZZ dummy;
162 image(d, Dummy, dummy);
163
164 int flag = 1, number = 0;
165
166 for(int i = 0; i <= m; i++)
167 if(dummy[0][i] >= 0) number++;
168 if(number == (m + 1)) flag = 0;
169
170 if(flag != 0){
171 number = 0;
172 for(int i = 0; i <= m; i++)
173 if(dummy[0][i] <= 0) number++;
174 if(number == (m + 1)) flag = 0;
175 }
176 if(flag == 0){
177 Z = - Z;
178 }
179 Dummy.kill();
180 dummy.kill();
181 return Z;
182
183 }
184 /**********************************************************************/
185
MatrixGCD(mat_ZZ & B,long & m)186 void MatrixGCD(mat_ZZ & B, long & m){
187 vec_ZZ gcds;
188 gcds.SetLength(m);
189 for(int i = 1; i <= m; i++)
190 for(int j = 1; j <= m; j++)
191 if(B(i, j) != 0)
192 gcds[i-1] = GCD(gcds[i-1], B(i, j));
193 for(int i = 1; i <= m; i++)
194 for(int j = 1; j <= m; j++)
195 if(B(i, j) != 0)
196 B(i, j) = B(i, j) / gcds[i-1];
197
198 }
199 /**********************************************************************/
200
201 /* Likewise barvinok_Single, but the cone is given as a listCone;
202 the function consumes the cone. */
203 int
204 barvinok_DFS(listCone *cone, Single_Cone_Parameters *Parameters);
205
206 int
barvinok_Single(mat_ZZ B,Single_Cone_Parameters * Parameters,const Vertex * vertex)207 barvinok_Single(mat_ZZ B, Single_Cone_Parameters *Parameters,
208 const Vertex *vertex)
209 {
210 //cerr << "barvinok_Single Called." << endl;;
211
212 long m, n;
213 m = B.NumRows();
214 n = B.NumCols();
215
216 if (m != n) {
217 cerr << "Input must be square (have " << m << " rows, "
218 << n << " cols). " << endl;
219 exit(2);
220 }
221
222 ZZ D = determinant(B);
223
224 if( D == 0)
225 {
226 cerr << "Input must be linearly independent. " << endl;
227 exit(3);
228 }
229
230 Parameters->Total_Simplicial_Cones++;
231
232 /* The following routine is to get the minimal
233 integral generators for the cone. */
234
235 MatrixGCD(B, m);
236
237 listCone *dummy = createListCone();
238 dummy->coefficient = 1;
239 dummy->determinant = D;
240 dummy->vertex = new Vertex(*vertex);
241 dummy->rays = transformArrayBigVectorToListVector(B, m, n);
242
243 switch (Parameters->decomposition) {
244 case BarvinokParameters::DualDecomposition:
245 // Keep the dual cones during Barvinok decomposition
246 computeDetAndFacetsOfSimplicialCone(dummy, n);
247 break;
248 case BarvinokParameters::IrrationalPrimalDecomposition:
249 // Do Barvinok decomposition on the primal cones.
250 dualizeCone(dummy, Parameters->Number_of_Variables, Parameters);
251 irrationalizeCone(dummy, Parameters->Number_of_Variables);
252 break;
253 case BarvinokParameters::IrrationalAllPrimalDecomposition:
254 // We already have primal cones; do Barvinok decomposition
255 // on them.
256 computeDetAndFacetsOfSimplicialCone(dummy, n);
257 break;
258 default:
259 cerr << "Unknown BarvinokParameters::decomposition" << endl;
260 abort();
261 }
262
263 int result;
264 result = barvinok_DFS(dummy, Parameters);
265
266 return result;
267 }
268
269 /* Let GENERATOR and MAT be the same matrix, with determinant DET.
270 Copy the vector Z into each row of the matrix (we are dealing with
271 the row space) and compute the determinant of the resulting matrix;
272 store the determinants in DETS[]. When any determinant is
273 larger-or-equal than DET in absolute value, return false;
274 otherwise return true.
275
276 If NONDECREASE_OK is false (the default), stop the computation
277 immediately when the result is known to be false.
278 */
279 static bool
computeAndCheckDeterminants(const mat_ZZ & generator,const ZZ & Det,const vec_ZZ & Z,int m,mat_ZZ & mat,vec_ZZ & Dets,bool nondecrease_ok=false)280 computeAndCheckDeterminants(const mat_ZZ &generator, const ZZ &Det,
281 const vec_ZZ &Z, int m,
282 mat_ZZ &mat, vec_ZZ &Dets, bool nondecrease_ok = false)
283 {
284 bool decrease = true;
285 ZZ absDet = abs(Det);
286 for (int i = 1; i <= m; i++) {
287 /* Copy in the row */
288 for(int j = 1; j <= m; j++)
289 mat(i, j) = Z(j);
290 /* Compute and store the determinant. */
291 determinant(Dets[i - 1], mat);
292 /* Restore the original row */
293 for(int j = 1; j <= m; j++)
294 mat(i, j) = generator(i, j);
295 if (abs(Dets[i - 1]) >= absDet) {
296 decrease = false;
297 if (!nondecrease_ok) return false;
298 }
299 }
300 return decrease;
301 }
302
303 /* Decompose the cone spanned by GENERATOR (which has determinant DET)
304 according to Barvinok's theory into M (the dimension) many CONES
305 and store their determinants in DETS.
306
307 Entries with Det[i] == 0 have Cones[i] == NULL (we don't generate
308 lower-dimensional cones).
309 */
310 bool
barvinokStep(const listCone * Cone,vector<listCone * > & Cones,vec_ZZ & Dets,int m,Single_Cone_Parameters * Parameters)311 barvinokStep(const listCone *Cone,
312 vector <listCone *> &Cones, vec_ZZ &Dets,
313 int m, Single_Cone_Parameters *Parameters)
314 {
315 #ifdef SHOWCONES
316 cerr << "############ Decomposing: ############" << endl;
317 printCone((listCone*)Cone, m);
318 cerr << "************** Results: **************" << endl;
319 #endif
320 mat_ZZ generator = createConeDecMatrix(Cone, m, m);
321 mat_ZZ dual = createFacetMatrix(Cone, m, m);
322 /* ComputeOmega(const mat_ZZ &, long& ) computes
323 an integral vector in the parallelogram. */
324 mat_ZZ mat;
325 vec_ZZ Z;
326 switch (Parameters->shortvector) {
327 case BarvinokParameters::LatteLLL: {
328 Z = ComputeOmega(generator, dual, m, 0, 0);
329 Z = CheckOmega(generator, Z);
330
331 mat = generator;
332 // FIXME: The determinants actually do not need to be computed;
333 // they are given by the vector L(index) in ComputeOmega...
334 // --mkoeppe, Fri Nov 24 17:01:37 MET 2006
335 bool success
336 = computeAndCheckDeterminants(generator, Cone->determinant, Z,
337 m, mat, Dets);
338 if (!success) {
339 cerr << "Second loop... " << endl;
340 Z = ComputeOmega(generator, dual, m, 2, 2);
341 Z = CheckOmega(generator, Z);
342 success = computeAndCheckDeterminants(generator, Cone->determinant, Z,
343 m, mat, Dets);
344 assert(success);
345 }
346 break;
347 }
348 case BarvinokParameters::SubspaceAvoidingLLL: {
349 #ifdef HAVE_EXPERIMENTS
350 Z = ComputeShortVectorAvoidingSubspace(generator, dual);
351 Z = CheckOmega(generator, Z);
352 mat = generator;
353 bool decrease
354 = computeAndCheckDeterminants(generator, Cone->determinant, Z,
355 m, mat, Dets, true);
356 if (!decrease)
357 return false;
358 #else
359 cerr << "SubspaceAvoidingLLL not compiled in, sorry." << endl;
360 exit(1);
361 #endif
362 break;
363 }
364 default:
365 assert(0);
366 }
367
368 for(int i = 0; i < m; i++) {
369 if (Dets[i] == 0)
370 Cones[i] = NULL;
371 else {
372 Cones[i] = createListCone();
373 {
374 /* Create the rays: */
375 /* Copy in the row */
376 for(int j = 1; j <= m; j++)
377 mat(i+1, j) = Z(j);
378 Cones[i]->rays
379 = transformArrayBigVectorToListVector(mat, m, m);
380 /* Restore the original row */
381 for(int j = 1; j <= m; j++)
382 mat(i+1, j) = generator(i+1, j);
383 }
384 Cones[i]->determinant = Dets[i];
385 {
386 /* Compute the sign: */
387 long signDet = sign(Cone->determinant);
388 long signDeti = sign(Dets[i]);
389 Cones[i]->coefficient = Cone->coefficient * signDet * signDeti;
390 }
391 Cones[i]->vertex = new Vertex(*Cone->vertex);
392 computeDetAndFacetsOfSimplicialCone(Cones[i], m);
393 #ifdef SHOWCONES
394 printCone(Cones[i], m);
395 #endif
396 }
397 }
398 return true;
399 }
400
401 static int
deliver_cone(listCone * C,Single_Cone_Parameters * Parameters)402 deliver_cone(listCone *C, Single_Cone_Parameters *Parameters)
403 {
404 Parameters->Total_Uni_Cones += 1;
405 if ( Parameters->Total_Uni_Cones % 1000 == 0)
406 cerr << Parameters->Total_Uni_Cones
407 << (Parameters->max_determinant == 0
408 ? " simplicial cones done."
409 : (Parameters->max_determinant == 1
410 ? " unimodular cones done."
411 : " low-index cones done."))
412 << endl;
413 switch (Parameters->decomposition) {
414 case BarvinokParameters::DualDecomposition:
415 dualizeCone(C, Parameters->Number_of_Variables, Parameters);
416 return Parameters->ConsumeCone(C);
417 case BarvinokParameters::IrrationalPrimalDecomposition:
418 case BarvinokParameters::IrrationalAllPrimalDecomposition:
419 return Parameters->ConsumeCone(C);
420 default:
421 cerr << "Unknown BarvinokParameters::decomposition" << endl;
422 abort();
423 }
424 }
425
426 static ZZ
criterion_abs_det(listCone * C,Single_Cone_Parameters * Parameters)427 criterion_abs_det(listCone *C, Single_Cone_Parameters *Parameters)
428 {
429 switch (Parameters->decomposition) {
430 case BarvinokParameters::DualDecomposition:
431 return abs(C->dual_determinant);
432 case BarvinokParameters::IrrationalPrimalDecomposition:
433 case BarvinokParameters::IrrationalAllPrimalDecomposition:
434 return abs(C->determinant);
435 default:
436 cerr << "Unknown BarvinokParameters::decomposition" << endl;
437 abort();
438 }
439 }
440
barvinok_DFS(listCone * C,Single_Cone_Parameters * Parameters)441 int barvinok_DFS(listCone *C, Single_Cone_Parameters *Parameters)
442 {
443 if (Parameters->Current_Depth > Parameters->Max_Depth)
444 Parameters->Max_Depth = Parameters->Current_Depth;
445
446 ZZ absDet = criterion_abs_det(C, Parameters);
447
448 if (absDet == 0) {
449 cerr << "barvinok_DFS: Det = 0." << endl;
450 return 1;
451 }
452 else {
453 switch (Parameters->decomposition) {
454 case BarvinokParameters::DualDecomposition:
455 break;
456 case BarvinokParameters::IrrationalPrimalDecomposition:
457 case BarvinokParameters::IrrationalAllPrimalDecomposition:
458 checkConeIrrational(C, Parameters->Number_of_Variables);
459 break;
460 default:
461 cerr << "Unknown BarvinokParameters::decomposition";
462 abort();
463 }
464 if (Parameters->max_determinant == 0
465 || absDet <= Parameters->max_determinant)
466 return deliver_cone(C, Parameters);
467 }
468
469 //cerr << "barvinok_DFS: non-uni cone." << endl;
470
471 int result = 1;
472 long m = Parameters->Number_of_Variables;
473
474 vec_ZZ Dets;
475 Dets.SetLength(m);
476 vector<listCone *> cones1(m);
477
478 bool success = barvinokStep(C, cones1, Dets, m, Parameters);
479 if (!success) {
480 cerr << "Unable to decompose cone with index " << absDet;
481 if (absDet <= 200000) { // "Emergency" max-index
482 cerr << ", enumerating it." << endl;
483 return deliver_cone(C, Parameters);
484 }
485 else {
486 cerr << ", giving up." << endl;
487 exit(1);
488 }
489 }
490
491 ZZ max;
492 max = -1;
493
494 #ifdef SHOWDETS
495 cerr << "Level " << Parameters->Current_Depth << ": Index " << absDet << " -> ";
496 #endif
497 for(int i = 0; i < m; i++)
498 {
499 Dets[i] = abs(Dets[i]);
500 if(Dets[i] > max)
501 max = Dets[i];
502
503 if (Dets[i] > 0) {
504 Parameters->Current_Simplicial_Cones_Total ++;
505
506 #ifdef SHOWDETS
507 cerr << criterion_abs_det(cones1[i], Parameters) << ", ";
508 #endif
509 switch (Parameters->decomposition) {
510 case BarvinokParameters::DualDecomposition:
511 break;
512 case BarvinokParameters::IrrationalPrimalDecomposition:
513 case BarvinokParameters::IrrationalAllPrimalDecomposition:
514 checkConeIrrational(cones1[i], Parameters->Number_of_Variables);
515 break;
516 default:
517 cerr << "Unknown BarvinokParameters::decomposition";
518 abort();
519 }
520 }
521 }
522 #ifdef SHOWDETS
523 cerr << endl;
524 #endif
525
526 int current;
527 ZZ min;
528
529 if (Parameters->Current_Simplicial_Cones_Total > Parameters->Max_Simplicial_Cones_Total)
530 Parameters->Max_Simplicial_Cones_Total = Parameters->Current_Simplicial_Cones_Total;
531
532 Parameters->Current_Depth++;
533
534 do {
535 min = max + 1;
536 current = -1;
537 for(int j = 0; j < m; j++) {
538 if(Dets[j] < min && Dets[j] != 0)
539 {
540 current = j;
541 min = Dets[j];
542 }
543 }
544 if (current >= 0) {
545 Dets[current] = 0; // mark done
546 if(barvinok_DFS(cones1[current], Parameters) == -1)
547 result = -1;
548 Parameters->Current_Simplicial_Cones_Total--;
549 }
550 } while (current >= 0 && result == 1);
551 Parameters->Current_Depth--;
552 freeCone(C);
553 return result;
554 }
555
556 /*
557 The first step is to triangulate a cone into simplicial cones.
558 Then, by using Barvinok's decomposition, we decompose each
559 simplicial cone into unimodular cones.
560 */
561
562 int
barvinokDecomposition_Single(listCone * cone,Single_Cone_Parameters * Parameters)563 barvinokDecomposition_Single (listCone *cone,
564 Single_Cone_Parameters *Parameters)
565 {
566 int status = 1;
567 listCone *triang = triangulateCone(cone, Parameters->Number_of_Variables, Parameters);
568 Parameters->decompose_time.start();
569 listCone *t;
570 for (t = triang; t!=NULL; t=t->rest) {
571 int num_rays = lengthListVector(t->rays);
572 assert(num_rays == Parameters->Number_of_Variables);
573 mat_ZZ B = createConeDecMatrix(t, num_rays, Parameters->Number_of_Variables);
574 if ((status = barvinok_Single(B, Parameters, t->vertex)) == -1)
575 goto BAILOUT;
576 }
577 BAILOUT:
578 Parameters->decompose_time.stop();
579 freeListCone(triang);
580 return status;
581 }
582
583