1 /* dec.cpp -- Calling Rudy's Decomposition Procedure
2
3 Copyright 2002, 2003 Raymond Hemmecke, Ruriko Yoshida
4 Copyright 2006 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 <cassert>
23 #include <fstream>
24 #include <stdlib.h>
25 #include <time.h>
26
27 #include "cone.h"
28 #include "print.h"
29 #include "ramon.h"
30 #include "latte_ntl_integer.h"
31 #include "RudyResNTL.h"
32 #include "PolyTree.h"
33 #include "flags.h"
34 #include "convert.h"
35 #include "dual.h"
36 #include "genFunction/piped.h"
37 #include "Residue.h"
38
39 #define MODULUS 1000000000
40 #define Exponent_Ten_Power 10
41
Collecting_Single_Cone_Parameters()42 Collecting_Single_Cone_Parameters::Collecting_Single_Cone_Parameters()
43 : Decomposed_Cones(NULL)
44 {
45 }
46
Collecting_Single_Cone_Parameters(const BarvinokParameters & params)47 Collecting_Single_Cone_Parameters::Collecting_Single_Cone_Parameters(const BarvinokParameters ¶ms)
48 : Single_Cone_Parameters(params), Decomposed_Cones(NULL)
49 {
50 }
51
ConsumeCone(listCone * cone)52 int Collecting_Single_Cone_Parameters::ConsumeCone(listCone *cone)
53 {
54 assert(cone->rest == NULL);
55 cone->rest = Decomposed_Cones;
56 Decomposed_Cones = cone;
57 return 1; // means "success, please continue"
58 }
59
60 listCone*
decomposeCones(listCone * cones,bool dualize,BarvinokParameters & param)61 decomposeCones(listCone *cones, bool dualize,
62 BarvinokParameters ¶m)
63 {
64 int numOfConesDecomposed,numOfAllCones;
65 listCone *tmp;
66
67 Collecting_Single_Cone_Parameters parameters(param);
68
69 if (dualize) {
70 dualizeCones(cones, param.Number_of_Variables, ¶m);
71 }
72
73 cerr << "Decomposing all cones.\n";
74 numOfConesDecomposed=0;
75 numOfAllCones=lengthListCone(cones);
76
77 parameters.Cone_Index = 0;
78 tmp=cones;
79 while (tmp) {
80 //parameters.Cone = tmp;
81 int result = barvinokDecomposition_Single(tmp,
82 ¶meters);
83 assert(result >= 0);
84
85 tmp=tmp->rest;
86 numOfConesDecomposed++;
87 if (numOfConesDecomposed==50*(numOfConesDecomposed/50)) {
88 cerr << numOfConesDecomposed << " / " << numOfAllCones << " done.\n";
89 }
90 parameters.Cone_Index++;
91 }
92
93 cerr << "All cones have been decomposed.\n";
94 cerr << lengthListCone(parameters.Decomposed_Cones) << " cones in total.\n";
95
96 return parameters.Decomposed_Cones;
97 }
98
99 listCone*
decomposeCones(listCone * cones,int numOfVars,unsigned int Flags,const char * File_Name,int max_determinant,bool dualize,BarvinokParameters::DecompositionType decomposition,bool debug_triangulation)100 decomposeCones(listCone *cones, int numOfVars, unsigned int Flags,
101 const char *File_Name, int max_determinant,
102 bool dualize,
103 BarvinokParameters::DecompositionType decomposition,
104 bool debug_triangulation)
105 {
106 Collecting_Single_Cone_Parameters parameters;
107 parameters.Flags = Flags;
108 parameters.Number_of_Variables = numOfVars;
109 parameters.max_determinant = max_determinant;
110 parameters.File_Name = File_Name;
111 parameters.decomposition = decomposition;
112 parameters.debug_triangulation = debug_triangulation;
113 return decomposeCones(cones, dualize, parameters);
114 }
115
116 /* ----------------------------------------------------------------- */
117
118
119
120 vec_ZZ
guess_generic_vector(int numOfVars)121 guess_generic_vector(int numOfVars)
122 {
123 vec_ZZ result;
124 result.SetLength(numOfVars);
125 int i;
126 for (i = 0; i < numOfVars; i++)
127 result[i] = (rand () % MODULUS) * ((rand() % 2) * 2 - 1);
128 return result;
129 }
130
131 void
InitializeComputation()132 Generic_Vector_Single_Cone_Parameters::InitializeComputation()
133 {
134 generic_vector = guess_generic_vector(Number_of_Variables);
135 }
136
137 void
barvinokDecomposition_List(listCone * cones,Generic_Vector_Single_Cone_Parameters & Parameters)138 barvinokDecomposition_List(listCone *cones,
139 Generic_Vector_Single_Cone_Parameters &Parameters)
140 {
141 do {
142 Parameters.InitializeComputation();
143 try {
144 listCone *cone;
145 int index;
146 for (cone = cones, index = 0; cone != NULL; cone = cone->rest, index++) {
147 int status;
148 status = barvinokDecomposition_Single(cone, &Parameters);
149 if (status < 0) {
150 static NotGenericException not_generic;
151 throw not_generic; // FIXME: Later replace this return
152 // value handling by exception.
153 }
154 if (index%1 == 0)
155 cerr << index << " vertex cones done. " << endl;
156 }
157 return;
158 }
159 catch (NotGenericException) {
160 cerr << "Generic vector chosen unsuccessfully, trying again." << endl;
161 };
162 } while (1);
163 }
164
165
166
167 void
InitializeComputation()168 Standard_Single_Cone_Parameters::InitializeComputation()
169 {
170 Generic_Vector_Single_Cone_Parameters::InitializeComputation();
171 int i;
172 for (i = 0; i <= Degree_of_Taylor_Expansion; i++)
173 Taylor_Expansion_Result[i] = 0;
174 Total_Lattice_Points = 0;
175 Total_Uni_Cones = 0;
176 Cone_Index = 0;
177 Max_Depth = 0;
178 Current_Depth = 0;
179 }
180
181 int
ConsumeCone(listCone * Cone)182 Standard_Single_Cone_Parameters::ConsumeCone(listCone *Cone)
183 {
184 //cerr << "barvinok_DFS: Calculating points in Parallelepiped" << endl;
185 if ( (Flags & DUAL_APPROACH) == 0)
186 computePointsInParallelepiped(Cone, Number_of_Variables, this);
187 //printListCone(Cone, Number_of_Variables);
188
189 if (Flags & DUAL_APPROACH)
190 {
191 //cerr << "barvinok_DFS: Calling ResidueFunction_Single_Cone" << endl;
192 return ResidueFunction_Single_Cone (Cone, this);
193 }
194 else
195 {
196 return Residue_Single_Cone (Cone, Number_of_Variables, generic_vector, &Total_Lattice_Points, &Ten_Power);
197 }
198 }
199
200 vec_ZZ
decomposeAndComputeResidue(listCone * cones,int degree,bool dualize,Standard_Single_Cone_Parameters & param)201 decomposeAndComputeResidue(listCone *cones, int degree, bool dualize,
202 Standard_Single_Cone_Parameters ¶m)
203 {
204 vec_ZZ polynomialAnswer;
205 int numOfAllCones;
206 mat_ZZ mat;
207
208 if (dualize) {
209 dualizeCones(cones, param.Number_of_Variables, ¶m);
210 }
211
212 cerr << "decomposeCones_Single: Decomposing all cones. (Memory Save on)\n";
213 numOfAllCones=lengthListCone(cones);
214 cerr << numOfAllCones << " cones total to be done!";
215
216 //Set Ten_Power to 100 billion
217 // FIXME: What is magic about this number? --mkoeppe, Sat Mar 4 21:21:45 PST 2006
218 param.Ten_Power = 1;
219 for (int i = 0; i < Exponent_Ten_Power; i++)
220 param.Ten_Power *= 10;
221
222 cerr << "decomposeCones_Single: degree = " << degree << endl;
223 param.Taylor_Expansion_Result = new ZZ [degree + 1 ];
224 polynomialAnswer.SetLength(degree+1);
225
226 param.Degree_of_Taylor_Expansion = degree;
227 param.Controller = new Node_Controller(param.Number_of_Variables + 1, degree);
228
229 cerr << "Number of cones: " << numOfAllCones << endl;
230
231 barvinokDecomposition_List(cones, param);
232
233 cerr << endl << "Total Unimodular Cones: " << param.Total_Uni_Cones << endl;
234 ofstream UniOut("numOfUnimodularCones");
235 UniOut << param.Total_Uni_Cones << endl;
236 cerr << "Maximum number of simplicial cones in memory at once: " << param.Max_Simplicial_Cones_Total << endl;
237
238
239 if ( param.Flags & DUAL_APPROACH)
240 {
241 cerr << "Memory Save Mode: Taylor Expansion:" << endl;
242 if(degree > 1){
243 for (int i = 0; i<= degree; i++)
244 {
245 ZZ number;
246 number = ( param.Taylor_Expansion_Result[i] + param.Ten_Power/2)/(param.Ten_Power) ;
247 cout << number;
248 polynomialAnswer[i] = number;
249
250 if (i != 0)
251 {
252 cout << "t^" << i;
253 }
254 cout << endl;
255 }
256 }
257 else if(degree == 1){
258 Integer numOfLatticePoints
259 = ( param.Taylor_Expansion_Result[1] + param.Ten_Power/2)/param.Ten_Power;
260 param.deliver_number_of_lattice_points(numOfLatticePoints);
261 polynomialAnswer.SetLength(2);
262 polynomialAnswer[0] = 0;
263 polynomialAnswer[1] = numOfLatticePoints;
264 }
265 }
266 else {
267 param.Total_Lattice_Points = abs(param.Total_Lattice_Points);
268 ZZ number;
269 number = (param.Total_Lattice_Points + param.Ten_Power/2)/param.Ten_Power;
270 param.deliver_number_of_lattice_points(number);
271
272 polynomialAnswer.SetLength(2);
273 polynomialAnswer[0] = 0;
274 polynomialAnswer[1] = number;
275
276 }
277
278 //cerr << lengthListCone(newCones->rest) << " cones in total.\n";
279
280 delete param.Controller;
281 delete [] param.Taylor_Expansion_Result;
282 return polynomialAnswer;
283
284 }
285
decomposeCones_Single(listCone * cones,int numOfVars,int degree,unsigned int flags,char * File_Name,int max_determinant,bool dualize,BarvinokParameters::DecompositionType decomposition)286 void decomposeCones_Single (listCone *cones, int numOfVars, int degree,
287 unsigned int flags, char *File_Name, int max_determinant,
288 bool dualize,
289 BarvinokParameters::DecompositionType decomposition)
290 {
291 Standard_Single_Cone_Parameters *Barvinok_Parameters
292 = new Standard_Single_Cone_Parameters;
293
294 Barvinok_Parameters->Flags = flags;
295 Barvinok_Parameters->Number_of_Variables = numOfVars;
296 Barvinok_Parameters->max_determinant = max_determinant;
297 Barvinok_Parameters->File_Name = File_Name;
298 Barvinok_Parameters->decomposition = decomposition;
299
300 decomposeAndComputeResidue(cones, degree, dualize, *Barvinok_Parameters);
301
302 delete Barvinok_Parameters;
303 }
304 /* ----------------------------------------------------------------- */
305
306