1 /*
2 * Normaliz
3 * Copyright (C) 2007-2021 W. Bruns, B. Ichim, Ch. Soeger, U. v. d. Ohe
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 *
17 * As an exception, when this program is distributed through (i) the App Store
18 * by Apple Inc.; (ii) the Mac App Store by Apple Inc.; or (iii) Google Play
19 * by Google Inc., then that store may impose any digital rights management,
20 * device limits and/or redistribution restrictions that are required by its
21 * terms of service.
22 */
23
24 #include "libnormaliz/project_and_lift.h"
25 #include "libnormaliz/vector_operations.h"
26 #include "libnormaliz/sublattice_representation.h"
27 #include "libnormaliz/cone.h"
28
29 namespace libnormaliz {
30 using std::vector;
31
32 //---------------------------------------------------------------------------
33
34 //---------------------------------------------------------------------------
35 // computes c1*v1-c2*v2
36 template <typename Integer>
FM_comb(Integer c1,const vector<Integer> & v1,Integer c2,const vector<Integer> & v2,bool & is_zero)37 vector<Integer> FM_comb(Integer c1, const vector<Integer>& v1, Integer c2, const vector<Integer>& v2, bool& is_zero) {
38 size_t dim = v1.size();
39 vector<Integer> new_supp(dim);
40 is_zero = false;
41 size_t k = 0;
42 for (; k < dim; ++k) {
43 new_supp[k] = c1 * v1[k] - c2 * v2[k];
44 if (!check_range(new_supp[k]))
45 break;
46 }
47 Integer g = 0;
48 if (k == dim)
49 g = v_make_prime(new_supp);
50 else { // redo in GMP if necessary
51 #pragma omp atomic
52 GMP_hyp++;
53 vector<mpz_class> mpz_neg(dim), mpz_pos(dim), mpz_sum(dim);
54 convert(mpz_neg, v1);
55 convert(mpz_pos, v2);
56 for (k = 0; k < dim; k++)
57 mpz_sum[k] = convertTo<mpz_class>(c1) * mpz_neg[k] - convertTo<mpz_class>(c2) * mpz_pos[k];
58 mpz_class GG = v_make_prime(mpz_sum);
59 convert(new_supp, mpz_sum);
60 convert(g, GG);
61 }
62 if (g == 0)
63 is_zero = true;
64 return new_supp;
65 }
66
67 template <typename IntegerPL, typename IntegerRet>
order_supps(const Matrix<IntegerPL> & Supps)68 vector<size_t> ProjectAndLift<IntegerPL, IntegerRet>::order_supps(const Matrix<IntegerPL>& Supps) {
69 assert(Supps.nr_of_rows() > 0);
70 size_t dim = Supps.nr_of_columns();
71
72 vector<pair<nmz_float, size_t> > NewPos, NewNeg, NewNeutr; // to record the order of the support haperplanes
73 for (size_t i = 0; i < Supps.nr_of_rows(); ++i) {
74 if (Supps[i][dim - 1] == 0) {
75 NewNeutr.push_back(make_pair(0.0, i));
76 continue;
77 }
78 nmz_float num, den;
79 convert(num, Supps[i][0]);
80 convert(den, Supps[i][dim - 1]);
81 nmz_float quot = num / den;
82 if (Supps[i][dim - 1] > 0)
83 NewPos.push_back(make_pair(Iabs(quot), i));
84 else
85 NewNeg.push_back(make_pair(Iabs(quot), i));
86 }
87 sort(NewPos.begin(), NewPos.end());
88 sort(NewNeg.begin(), NewNeg.end());
89 NewPos.insert(NewPos.end(), NewNeutr.begin(), NewNeutr.end());
90
91 size_t min_length = NewNeg.size();
92 if (NewPos.size() < min_length)
93 min_length = NewPos.size();
94
95 vector<size_t> Order;
96
97 for (size_t i = 0; i < min_length; ++i) {
98 Order.push_back(NewPos[i].second);
99 Order.push_back(NewNeg[i].second);
100 }
101 for (size_t i = min_length; i < NewPos.size(); ++i)
102 Order.push_back(NewPos[i].second);
103 for (size_t i = min_length; i < NewNeg.size(); ++i)
104 Order.push_back(NewNeg[i].second);
105
106 assert(Order.size() == Supps.nr_of_rows());
107
108 /* for(size_t i=0;i<Order.size();++i)
109 cout << Supps[Order[i]][dim-1] << " ";
110 cout << endl;*/
111
112 return Order;
113 }
114 //---------------------------------------------------------------------------
115 template <typename IntegerPL, typename IntegerRet>
compute_projections(size_t dim,size_t down_to,vector<dynamic_bitset> & Ind,vector<dynamic_bitset> & Pair,vector<dynamic_bitset> & ParaInPair,size_t rank,bool only_projections)116 void ProjectAndLift<IntegerPL, IntegerRet>::compute_projections(size_t dim,
117 size_t down_to,
118 vector<dynamic_bitset>& Ind,
119 vector<dynamic_bitset>& Pair,
120 vector<dynamic_bitset>& ParaInPair,
121 size_t rank,
122 bool only_projections) {
123 INTERRUPT_COMPUTATION_BY_EXCEPTION
124
125 const Matrix<IntegerPL>& Supps = AllSupps[dim];
126
127 size_t dim1 = dim - 1;
128
129 if (verbose)
130 verboseOutput() << "embdim " << dim << " inequalities " << Supps.nr_of_rows() << endl;
131
132 if (dim == down_to)
133 return;
134
135 // Supps.pretty_print(cout);
136 // cout << Ind;
137
138 // cout << "SSS" << Ind.size() << " " << Ind;
139
140 // We now augment the given cone by the last basis vector and its negative
141 // Afterwards we project modulo the subspace spanned by them
142
143 vector<key_t> Neg, Pos; // for the Fourier-Motzkin elimination of inequalities
144 Matrix<IntegerPL> SuppsProj(0, dim); // for the support hyperplanes of the projection
145 Matrix<IntegerPL> EqusProj(0, dim); // for the equations (both later minimized)
146
147 // First we make incidence vectors with the given generators
148 vector<dynamic_bitset> NewInd; // for the incidence vectors of the new hyperplanes
149 vector<dynamic_bitset> NewPair; // for the incidence vectors of the new hyperplanes
150 vector<dynamic_bitset> NewParaInPair; // for the incidence vectors of the new hyperplanes
151
152 dynamic_bitset TRUE;
153 if (!is_parallelotope) {
154 TRUE.resize(Ind[0].size());
155 TRUE.set();
156 }
157
158 vector<bool> IsEquation(Supps.nr_of_rows());
159
160 bool rank_goes_up = false; // if we add the last unit vector
161 size_t PosEquAt = 0; // we memorize the positions of pos/neg equations if rank goes up
162 size_t NegEquAt = 0;
163
164 for (size_t i = 0; i < Supps.nr_of_rows(); ++i) {
165 if (!is_parallelotope && Ind[i] == TRUE)
166 IsEquation[i] = true;
167
168 if (Supps[i][dim1] == 0) { // already independent of last coordinate
169 no_crunch = false;
170 if (IsEquation[i])
171 EqusProj.append(Supps[i]); // is equation
172 else {
173 SuppsProj.append(Supps[i]); // neutral support hyperplane
174 if (!is_parallelotope)
175 NewInd.push_back(Ind[i]);
176 else {
177 NewPair.push_back(Pair[i]);
178 NewParaInPair.push_back(ParaInPair[i]);
179 }
180 }
181 continue;
182 }
183 if (IsEquation[i])
184 rank_goes_up = true;
185 if (Supps[i][dim1] > 0) {
186 if (IsEquation[i])
187 PosEquAt = i;
188 Pos.push_back(i);
189 continue;
190 }
191 Neg.push_back(i);
192 if (IsEquation[i])
193 NegEquAt = i;
194 }
195
196 // cout << "Nach Pos/Neg " << EqusProj.nr_of_rows() << " " << Pos.size() << " " << Neg.size() << endl;
197
198 // now the elimination, matching Pos and Neg
199
200 bool skip_remaining;
201 std::exception_ptr tmp_exception;
202
203 if (rank_goes_up) {
204 assert(!is_parallelotope);
205
206 for (size_t p : Pos) { // match pos and neg equations
207 if (!IsEquation[p])
208 continue;
209 IntegerPL PosVal = Supps[p][dim1];
210 for (size_t n : Neg) {
211 if (!IsEquation[n])
212 continue;
213 IntegerPL NegVal = Supps[n][dim1];
214 bool is_zero;
215 // cout << Supps[p];
216 // cout << Supps[n];
217 vector<IntegerPL> new_equ = FM_comb(PosVal, Supps[n], NegVal, Supps[p], is_zero);
218 // cout << "zero " << is_zero << endl;
219 // cout << "=====================" << endl;
220 if (is_zero)
221 continue;
222 EqusProj.append(new_equ);
223 }
224 }
225
226 for (size_t p : Pos) { // match pos inequalities with a negative equation
227 if (IsEquation[p])
228 continue;
229 IntegerPL PosVal = Supps[p][dim1];
230 IntegerPL NegVal = Supps[NegEquAt][dim1];
231 vector<IntegerPL> new_supp(dim);
232 bool is_zero;
233 new_supp = FM_comb(PosVal, Supps[NegEquAt], NegVal, Supps[p], is_zero);
234 /* cout << Supps[NegEquAt];
235 cout << Supps[p];
236 cout << new_supp;
237 cout << "zero " << is_zero << endl;
238 cout << "+++++++++++++++++++++" << endl; */
239 if (is_zero) // cannot happen, but included for analogy
240 continue;
241 SuppsProj.append(new_supp);
242 NewInd.push_back(Ind[p]);
243 }
244
245 for (size_t n : Neg) { // match neg inequalities with a positive equation
246 if (IsEquation[n])
247 continue;
248 IntegerPL PosVal = Supps[PosEquAt][dim1];
249 IntegerPL NegVal = Supps[n][dim1];
250 vector<IntegerPL> new_supp(dim);
251 bool is_zero;
252 new_supp = FM_comb(PosVal, Supps[n], NegVal, Supps[PosEquAt], is_zero);
253 /* cout << Supps[PosEquAt];
254 cout << Supps[n];
255 cout << new_supp;
256 cout << "zero " << is_zero << endl;
257 cout << "=====================" << endl;*/
258
259 if (is_zero) // cannot happen, but included for analogy
260 continue;
261 SuppsProj.append(new_supp);
262 NewInd.push_back(Ind[n]);
263 }
264 }
265
266 // cout << "Nach RGU " << EqusProj.nr_of_rows() << " " << SuppsProj.nr_of_rows() << endl;
267
268 if (!rank_goes_up && !is_parallelotope) { // must match pos and neg hyperplanes
269
270 // cout << "Pos " << Pos.size() << " Neg " << Neg.size() << " Supps " << SuppsProj.nr_of_rows() << endl;
271
272 skip_remaining = false;
273
274 size_t min_nr_vertices = rank - 2;
275 /*if(rank>=3){
276 min_nr_vertices=1;
277 for(long i=0;i<(long) rank -3;++i)
278 min_nr_vertices*=2;
279
280 }*/
281
282 #pragma omp parallel for schedule(dynamic)
283 for (size_t i = 0; i < Pos.size(); ++i) {
284 if (skip_remaining)
285 continue;
286
287 try {
288 size_t p = Pos[i];
289 IntegerPL PosVal = Supps[p][dim1];
290 vector<key_t> PosKey;
291 for (size_t k = 0; k < Ind[i].size(); ++k)
292 if (Ind[p][k])
293 PosKey.push_back(k);
294
295 for (size_t n : Neg) {
296 INTERRUPT_COMPUTATION_BY_EXCEPTION
297
298 // // to give a facet of the extended cone
299 // match incidence vectors
300 dynamic_bitset incidence(TRUE.size());
301 size_t nr_match = 0;
302 vector<key_t> CommonKey;
303 for (unsigned int k : PosKey)
304 if (Ind[n][k]) {
305 incidence[k] = true;
306 CommonKey.push_back(k);
307 nr_match++;
308 }
309 if (rank >= 2 && nr_match < min_nr_vertices) // cannot make subfacet of augmented cone
310 continue;
311
312 bool IsSubfacet = true;
313 for (size_t k = 0; k < Supps.nr_of_rows(); ++k) {
314 if (k == p || k == n || IsEquation[k])
315 continue;
316 bool contained = true;
317 for (unsigned int j : CommonKey) {
318 if (!Ind[k][j]) {
319 contained = false;
320 break;
321 }
322 }
323 if (contained) {
324 IsSubfacet = false;
325 break;
326 }
327 }
328 if (!IsSubfacet)
329 continue;
330 //}
331
332 IntegerPL NegVal = Supps[n][dim1];
333 vector<IntegerPL> new_supp(dim);
334 bool is_zero;
335 new_supp = FM_comb(PosVal, Supps[n], NegVal, Supps[p], is_zero);
336 if (is_zero) // linear combination is 0
337 continue;
338
339 if (nr_match == TRUE.size()) { // gives an equation
340 #pragma omp critical(NEWEQ)
341 EqusProj.append(new_supp);
342 continue;
343 }
344 #pragma omp critical(NEWSUPP)
345 {
346 SuppsProj.append(new_supp);
347 NewInd.push_back(incidence);
348 }
349 }
350
351 } catch (const std::exception&) {
352 tmp_exception = std::current_exception();
353 skip_remaining = true;
354 #pragma omp flush(skip_remaining)
355 }
356 }
357
358 } // !rank_goes_up && !is_parallelotope
359
360 if (!(tmp_exception == 0))
361 std::rethrow_exception(tmp_exception);
362
363 if (!rank_goes_up && is_parallelotope) { // must match pos and neg hyperplanes
364
365 size_t codim = dim1 - 1; // the minimal codim a face of the original cone must have
366 // in order to project to a subfacet of the current one
367 size_t original_dim = Pair[0].size() + 1;
368 size_t max_number_containing_factes = original_dim - codim;
369
370 skip_remaining = false;
371
372 size_t nr_pos = Pos.size();
373 // if(nr_pos>10000)
374 // nr_pos=10000;
375 size_t nr_neg = Neg.size();
376 // if(nr_neg>10000)
377 // nr_neg=10000;
378
379 #pragma omp parallel for schedule(dynamic)
380 for (size_t i = 0; i < nr_pos; ++i) {
381 if (skip_remaining)
382 continue;
383
384 try {
385 INTERRUPT_COMPUTATION_BY_EXCEPTION
386
387 size_t p = Pos[i];
388 IntegerPL PosVal = Supps[p][dim1];
389
390 for (size_t j = 0; j < nr_neg; ++j) {
391 size_t n = Neg[j];
392 dynamic_bitset IntersectionPair(Pair[p].size());
393 size_t nr_hyp_intersection = 0;
394 bool in_parallel_hyperplanes = false;
395 bool codim_too_small = false;
396
397 for (size_t k = 0; k < Pair[p].size(); ++k) { // run over all pairs
398 if (Pair[p][k] || Pair[n][k]) {
399 nr_hyp_intersection++;
400 IntersectionPair[k] = true;
401 if (nr_hyp_intersection > max_number_containing_factes) {
402 codim_too_small = true;
403 break;
404 }
405 }
406 if (Pair[p][k] && Pair[n][k]) {
407 if (ParaInPair[p][k] != ParaInPair[n][k]) {
408 in_parallel_hyperplanes = true;
409 break;
410 }
411 }
412 }
413 if (in_parallel_hyperplanes || codim_too_small)
414 continue;
415
416 dynamic_bitset IntersectionParaInPair(Pair[p].size());
417 for (size_t k = 0; k < ParaInPair[p].size(); ++k) {
418 if (Pair[p][k])
419 IntersectionParaInPair[k] = ParaInPair[p][k];
420 else if (Pair[n][k])
421 IntersectionParaInPair[k] = ParaInPair[n][k];
422 }
423
424 // we must nevertheless use the comparison test
425 bool IsSubfacet = true;
426 if (!no_crunch) {
427 for (size_t k = 0; k < Supps.nr_of_rows(); ++k) {
428 if (k == p || k == n || IsEquation[k])
429 continue;
430 bool contained = true;
431
432 for (size_t u = 0; u < IntersectionPair.size(); ++u) {
433 if (Pair[k][u] && !IntersectionPair[u]) { // hyperplane k contains facet of Supp
434 contained = false; // not our intersection
435 continue;
436 }
437 if (Pair[k][u] && IntersectionPair[u]) {
438 if (ParaInPair[k][u] != IntersectionParaInPair[u]) { // they are contained in parallel
439 contained = false; // original facets
440 continue;
441 }
442 }
443 }
444
445 if (contained) {
446 IsSubfacet = false;
447 break;
448 }
449 }
450 }
451 if (!IsSubfacet)
452 continue;
453
454 IntegerPL NegVal = Supps[n][dim1];
455 bool dummy;
456 vector<IntegerPL> new_supp = FM_comb(PosVal, Supps[n], NegVal, Supps[p], dummy);
457 #pragma omp critical(NEWSUPP)
458 {
459 SuppsProj.append(new_supp);
460 NewPair.push_back(IntersectionPair);
461 NewParaInPair.push_back(IntersectionParaInPair);
462 }
463 }
464
465 } catch (const std::exception&) {
466 tmp_exception = std::current_exception();
467 skip_remaining = true;
468 #pragma omp flush(skip_remaining)
469 }
470 }
471
472 if (!(tmp_exception == 0))
473 std::rethrow_exception(tmp_exception);
474
475 } // !rank_goes_up && is_parallelotope
476
477 // cout << "Nach FM " << EqusProj.nr_of_rows() << " " << SuppsProj.nr_of_rows() << endl;
478
479 Ind.clear(); // no longer needed
480
481 EqusProj.resize_columns(dim1); // cut off the trailing 0
482 SuppsProj.resize_columns(dim1); // project hyperplanes
483
484 // Equations have not yet been appended to support hypwerplanes
485 EqusProj.row_echelon(); // reduce equations
486 // cout << "Nach eche " << EqusProj.nr_of_rows() << endl;
487 /* for(size_t i=0;i<EqusProj.nr_of_rows(); ++i)
488 cout << EqusProj[i]; */
489 SuppsProj.append(EqusProj); // append them as pairs of inequalities
490 EqusProj.scalar_multiplication(-1);
491 SuppsProj.append(EqusProj);
492 AllNrEqus[dim1] = EqusProj.nr_of_rows();
493 // We must add indictor vectors for the equations
494 for (size_t i = 0; i < 2 * EqusProj.nr_of_rows(); ++i)
495 NewInd.push_back(TRUE);
496
497 if (dim1 > 1 && !only_projections)
498 AllOrders[dim1] = order_supps(SuppsProj);
499 swap(AllSupps[dim1], SuppsProj);
500
501 size_t new_rank = dim1 - EqusProj.nr_of_rows();
502
503 compute_projections(dim - 1, down_to, NewInd, NewPair, NewParaInPair, new_rank);
504 }
505
506 //---------------------------------------------------------------------------
507 template <typename IntegerPL, typename IntegerRet>
fiber_interval(IntegerRet & MinInterval,IntegerRet & MaxInterval,const vector<IntegerRet> & base_point)508 bool ProjectAndLift<IntegerPL, IntegerRet>::fiber_interval(IntegerRet& MinInterval,
509 IntegerRet& MaxInterval,
510 const vector<IntegerRet>& base_point) {
511 size_t dim = base_point.size() + 1;
512 Matrix<IntegerPL>& Supps = AllSupps[dim];
513 vector<size_t>& Order = AllOrders[dim];
514
515 bool FirstMin = true, FirstMax = true;
516 vector<IntegerPL> LiftedGen;
517 convert(LiftedGen, base_point);
518 // cout << LiftedGen;
519 size_t check_supps = Supps.nr_of_rows();
520 if (check_supps > 1000 && dim < EmbDim && !no_relax)
521 check_supps = 1000;
522 for (size_t j = 0; j < check_supps; ++j) {
523 INTERRUPT_COMPUTATION_BY_EXCEPTION
524
525 IntegerPL Den = Supps[Order[j]].back();
526 if (Den == 0)
527 continue;
528 IntegerPL Num = -v_scalar_product_vectors_unequal_lungth(LiftedGen, Supps[Order[j]]);
529 // cout << "Num " << Num << endl;
530 IntegerRet Bound = 0;
531 // frac=(Num % Den !=0);
532 if (Den > 0) { // we must produce a lower bound of the interval
533 Bound = ceil_quot<IntegerRet, IntegerPL>(Num, Den);
534 if (FirstMin || Bound > MinInterval) {
535 MinInterval = Bound;
536 FirstMin = false;
537 }
538 }
539 if (Den < 0) { // we must produce an upper bound of the interval
540 Bound = floor_quot<IntegerRet, IntegerPL>(Num, Den);
541 if (FirstMax || Bound < MaxInterval) {
542 MaxInterval = Bound;
543 FirstMax = false;
544 }
545 }
546 if (!FirstMax && !FirstMin && MaxInterval < MinInterval)
547 return false; // interval empty
548 }
549 return true; // interval nonempty
550 }
551
552 ///---------------------------------------------------------------------------
553 template <typename IntegerPL, typename IntegerRet>
lift_points_to_this_dim(list<vector<IntegerRet>> & Deg1Proj)554 void ProjectAndLift<IntegerPL, IntegerRet>::lift_points_to_this_dim(list<vector<IntegerRet> >& Deg1Proj) {
555 if (Deg1Proj.empty())
556 return;
557 size_t dim = Deg1Proj.front().size() + 1;
558 size_t dim1 = dim - 1;
559
560 list<vector<IntegerRet> > Deg1Lifted; // to this dimension if < EmbDim
561 vector<list<vector<IntegerRet> > > Deg1Thread(omp_get_max_threads());
562 size_t max_nr_per_thread = 100000 / omp_get_max_threads();
563
564 vector<vector<num_t> > h_vec_pos_thread(omp_get_max_threads());
565 vector<vector<num_t> > h_vec_neg_thread(omp_get_max_threads());
566
567 size_t nr_to_lift = Deg1Proj.size();
568 NrLP[dim1] += nr_to_lift;
569
570 bool not_done = true;
571
572 while (not_done) {
573 // cout << "Durchgang dim " << dim << endl;
574
575 not_done = false;
576 bool message_printed = false;
577
578 bool skip_remaining;
579 std::exception_ptr tmp_exception;
580
581 skip_remaining = false;
582 int omp_start_level = omp_get_level();
583
584 #pragma omp parallel
585 {
586 int tn;
587 if (omp_get_level() == omp_start_level)
588 tn = 0;
589 else
590 tn = omp_get_ancestor_thread_num(omp_start_level + 1);
591
592 size_t nr_points_in_thread = 0;
593
594 size_t ppos = 0;
595 auto p = Deg1Proj.begin();
596 #pragma omp for schedule(dynamic)
597 for (size_t i = 0; i < nr_to_lift; ++i) {
598 if (skip_remaining)
599 continue;
600
601 for (; i > ppos; ++ppos, ++p)
602 ;
603 for (; i < ppos; --ppos, --p)
604 ;
605
606 if ((*p)[0] == 0) // point done
607 continue;
608
609 if (!not_done && verbose) {
610 #pragma omp critical
611 {
612 if (!message_printed)
613 verboseOutput() << "Lifting to dimension " << dim << endl;
614 message_printed = true;
615 }
616 }
617
618 not_done = true;
619
620 try {
621 IntegerRet MinInterval = 0, MaxInterval = 0; // the fiber over *p is an interval -- 0 to make gcc happy
622 fiber_interval(MinInterval, MaxInterval, *p);
623 // cout << "Min " << MinInterval << " Max " << MaxInterval << endl;
624 IntegerRet add_nr_Int = 0;
625 if (MaxInterval >= MinInterval)
626 add_nr_Int = 1 + MaxInterval - MinInterval;
627 long long add_nr = convertToLongLong(add_nr_Int);
628 if (dim == EmbDim && count_only && add_nr >= 1 && Congs.nr_of_rows() == 0 && Grading.size() == 0) {
629 #pragma omp atomic
630 TotalNrLP += add_nr;
631 }
632 else { // lift ppoint
633 for (IntegerRet k = MinInterval; k <= MaxInterval; ++k) {
634 INTERRUPT_COMPUTATION_BY_EXCEPTION
635
636 vector<IntegerRet> NewPoint(dim);
637 for (size_t j = 0; j < dim1; ++j)
638 NewPoint[j] = (*p)[j];
639 NewPoint[dim1] = k;
640 if (dim == EmbDim) {
641 if (!Congs.check_congruences(NewPoint))
642 continue;
643
644 #pragma omp atomic
645 TotalNrLP++;
646
647 if (!count_only)
648 Deg1Thread[tn].push_back(NewPoint);
649
650 if (Grading.size() > 0) {
651 long deg = convertToLong(v_scalar_product(Grading, NewPoint));
652 if (deg >= 0) {
653 if (deg >= (long)h_vec_pos_thread[tn].size())
654 h_vec_pos_thread[tn].resize(deg + 1);
655 h_vec_pos_thread[tn][deg]++;
656 }
657 else {
658 deg *= -1;
659 if (deg >= (long)h_vec_neg_thread[tn].size())
660 h_vec_neg_thread[tn].resize(deg + 1);
661 h_vec_neg_thread[tn][deg]++;
662 }
663 }
664 }
665 else
666 Deg1Thread[tn].push_back(NewPoint);
667 }
668 }
669
670 (*p)[0] = 0; // mark point as done
671 if (dim < EmbDim)
672 nr_points_in_thread += add_nr;
673 if (nr_points_in_thread > max_nr_per_thread) { // thread is full
674 skip_remaining = true;
675 #pragma omp flush(skip_remaining)
676 }
677
678 } catch (const std::exception&) {
679 tmp_exception = std::current_exception();
680 skip_remaining = true;
681 #pragma omp flush(skip_remaining)
682 }
683
684 } // lifting
685
686 } // pararllel
687
688 if (!(tmp_exception == 0))
689 std::rethrow_exception(tmp_exception);
690
691 for (size_t i = 0; i < Deg1Thread.size(); ++i)
692 Deg1Lifted.splice(Deg1Lifted.begin(), Deg1Thread[i]);
693
694 if (dim == EmbDim)
695 Deg1Points.splice(Deg1Points.end(), Deg1Lifted);
696
697 for (size_t i = 0; i < Deg1Thread.size(); ++i) {
698 if (h_vec_pos_thread[i].size() > h_vec_pos.size())
699 h_vec_pos.resize(h_vec_pos_thread[i].size());
700 for (size_t j = 0; j < h_vec_pos_thread[i].size(); ++j)
701 h_vec_pos[j] += h_vec_pos_thread[i][j];
702 h_vec_pos_thread[i].clear();
703 }
704
705 for (size_t i = 0; i < Deg1Thread.size(); ++i) {
706 if (h_vec_neg_thread[i].size() > h_vec_neg.size())
707 h_vec_neg.resize(h_vec_neg_thread[i].size());
708 for (size_t j = 0; j < h_vec_neg_thread[i].size(); ++j)
709 h_vec_neg[j] += h_vec_neg_thread[i][j];
710 h_vec_neg_thread[i].clear();
711 }
712
713 lift_points_to_this_dim(Deg1Lifted);
714 Deg1Lifted.clear();
715
716 } // not_done
717
718 return;
719
720 /* Deg1.pretty_print(cout);
721 cout << "*******************" << endl; */
722 }
723
724 ///---------------------------------------------------------------------------
725 template <typename IntegerPL, typename IntegerRet>
lift_point_recursively(vector<IntegerRet> & final_latt_point,const vector<IntegerRet> & latt_point_proj)726 void ProjectAndLift<IntegerPL, IntegerRet>::lift_point_recursively(vector<IntegerRet>& final_latt_point,
727 const vector<IntegerRet>& latt_point_proj) {
728 size_t dim1 = latt_point_proj.size();
729 size_t dim = dim1 + 1;
730 size_t final_dim = AllSupps.size() - 1;
731
732 IntegerRet MinInterval = 0, MaxInterval = 0; // the fiber over Deg1Proj[i] is an interval -- 0 to make gcc happy
733 fiber_interval(MinInterval, MaxInterval, latt_point_proj);
734 for (IntegerRet k = MinInterval; k <= MaxInterval; ++k) {
735 INTERRUPT_COMPUTATION_BY_EXCEPTION
736
737 vector<IntegerRet> NewPoint(dim);
738 for (size_t j = 0; j < dim1; ++j)
739 NewPoint[j] = latt_point_proj[j];
740 NewPoint[dim1] = k;
741 if (dim == final_dim && NewPoint != excluded_point) {
742 final_latt_point = NewPoint;
743 break;
744 }
745 if (dim < final_dim) {
746 lift_point_recursively(final_latt_point, NewPoint);
747 if (final_latt_point.size() > 0)
748 break;
749 }
750 }
751 }
752
753 ///---------------------------------------------------------------------------
754 template <typename IntegerPL, typename IntegerRet>
find_single_point()755 void ProjectAndLift<IntegerPL, IntegerRet>::find_single_point() {
756 size_t dim = AllSupps.size() - 1;
757 assert(dim >= 2);
758
759 vector<IntegerRet> start(1, GD);
760 vector<IntegerRet> final_latt_point;
761 lift_point_recursively(final_latt_point, start);
762 if (final_latt_point.size() > 0) {
763 SingleDeg1Point = final_latt_point;
764 if (verbose)
765 verboseOutput() << "Found point" << endl;
766 }
767 else {
768 if (verbose)
769 verboseOutput() << "No point found" << endl;
770 }
771 }
772
773 ///---------------------------------------------------------------------------
774 template <typename IntegerPL, typename IntegerRet>
compute_latt_points()775 void ProjectAndLift<IntegerPL, IntegerRet>::compute_latt_points() {
776 size_t dim = AllSupps.size() - 1;
777 assert(dim >= 2);
778
779 vector<IntegerRet> start(1, GD);
780 list<vector<IntegerRet> > start_list;
781 start_list.push_back(start);
782 lift_points_to_this_dim(start_list);
783 // cout << "TTTT " << TotalNrLP << endl;
784 NrLP[EmbDim] = TotalNrLP;
785 if (verbose) {
786 for (size_t i = 2; i < NrLP.size(); ++i)
787 verboseOutput() << "embdim " << i << " LatticePoints " << NrLP[i] << endl;
788 }
789 }
790
791 ///---------------------------------------------------------------------------
792 template <typename IntegerPL, typename IntegerRet>
compute_latt_points_float()793 void ProjectAndLift<IntegerPL, IntegerRet>::compute_latt_points_float() {
794 ProjectAndLift<nmz_float, IntegerRet> FloatLift(*this);
795 FloatLift.compute_latt_points();
796 Deg1Points.swap(FloatLift.Deg1Points);
797 TotalNrLP = FloatLift.TotalNrLP;
798 h_vec_pos = FloatLift.h_vec_pos;
799 h_vec_neg = FloatLift.h_vec_neg;
800 }
801
802 //---------------------------------------------------------------------------
803 template <typename IntegerPL, typename IntegerRet>
initialize(const Matrix<IntegerPL> & Supps,size_t rank)804 void ProjectAndLift<IntegerPL, IntegerRet>::initialize(const Matrix<IntegerPL>& Supps, size_t rank) {
805 EmbDim = Supps.nr_of_columns(); // our embedding dimension
806 AllSupps.resize(EmbDim + 1);
807 AllOrders.resize(EmbDim + 1);
808 AllNrEqus.resize(EmbDim + 1);
809 AllSupps[EmbDim] = Supps;
810 AllOrders[EmbDim] = order_supps(Supps);
811 StartRank = rank;
812 GD = 1; // the default choice
813 verbose = true;
814 is_parallelotope = false;
815 no_crunch = true;
816 use_LLL = false;
817 no_relax = false;
818 TotalNrLP = 0;
819 NrLP.resize(EmbDim + 1);
820
821 Congs = Matrix<IntegerRet>(0, EmbDim + 1);
822
823 LLL_Coordinates = Sublattice_Representation<IntegerRet>(EmbDim); // identity
824 }
825
826 //---------------------------------------------------------------------------
827 template <typename IntegerPL, typename IntegerRet>
ProjectAndLift()828 ProjectAndLift<IntegerPL, IntegerRet>::ProjectAndLift() {
829 }
830 //---------------------------------------------------------------------------
831 // General constructor
832 template <typename IntegerPL, typename IntegerRet>
ProjectAndLift(const Matrix<IntegerPL> & Supps,const vector<dynamic_bitset> & Ind,size_t rank)833 ProjectAndLift<IntegerPL, IntegerRet>::ProjectAndLift(const Matrix<IntegerPL>& Supps,
834 const vector<dynamic_bitset>& Ind,
835 size_t rank) {
836 initialize(Supps, rank);
837 StartInd = Ind;
838 }
839
840 //---------------------------------------------------------------------------
841 // Constructor for parallelotopes
842 template <typename IntegerPL, typename IntegerRet>
ProjectAndLift(const Matrix<IntegerPL> & Supps,const vector<dynamic_bitset> & Pair,const vector<dynamic_bitset> & ParaInPair,size_t rank)843 ProjectAndLift<IntegerPL, IntegerRet>::ProjectAndLift(const Matrix<IntegerPL>& Supps,
844 const vector<dynamic_bitset>& Pair,
845 const vector<dynamic_bitset>& ParaInPair,
846 size_t rank) {
847 initialize(Supps, rank);
848 is_parallelotope = true;
849 StartPair = Pair;
850 StartParaInPair = ParaInPair;
851 }
852
853 //---------------------------------------------------------------------------
854 template <typename IntegerPL, typename IntegerRet>
set_congruences(const Matrix<IntegerRet> & congruences)855 void ProjectAndLift<IntegerPL, IntegerRet>::set_congruences(const Matrix<IntegerRet>& congruences) {
856 Congs = congruences;
857 }
858
859 //---------------------------------------------------------------------------
860 template <typename IntegerPL, typename IntegerRet>
set_verbose(bool on_off)861 void ProjectAndLift<IntegerPL, IntegerRet>::set_verbose(bool on_off) {
862 verbose = on_off;
863 }
864
865 //---------------------------------------------------------------------------
866 template <typename IntegerPL, typename IntegerRet>
set_LLL(bool on_off)867 void ProjectAndLift<IntegerPL, IntegerRet>::set_LLL(bool on_off) {
868 use_LLL = on_off;
869 }
870
871 //---------------------------------------------------------------------------
872 template <typename IntegerPL, typename IntegerRet>
set_no_relax(bool on_off)873 void ProjectAndLift<IntegerPL, IntegerRet>::set_no_relax(bool on_off) {
874 no_relax = on_off;
875 }
876
877 //---------------------------------------------------------------------------
878 template <typename IntegerPL, typename IntegerRet>
set_grading_denom(const IntegerRet GradingDenom)879 void ProjectAndLift<IntegerPL, IntegerRet>::set_grading_denom(const IntegerRet GradingDenom) {
880 GD = GradingDenom;
881 }
882
883 //---------------------------------------------------------------------------
884 template <typename IntegerPL, typename IntegerRet>
set_grading(const vector<IntegerRet> & grad)885 void ProjectAndLift<IntegerPL, IntegerRet>::set_grading(const vector<IntegerRet>& grad) {
886 Grading = grad;
887 }
888
889 //---------------------------------------------------------------------------
890 template <typename IntegerPL, typename IntegerRet>
set_excluded_point(const vector<IntegerRet> & excl_point)891 void ProjectAndLift<IntegerPL, IntegerRet>::set_excluded_point(const vector<IntegerRet>& excl_point) {
892 excluded_point = excl_point;
893 }
894
895 //---------------------------------------------------------------------------
896 template <typename IntegerPL, typename IntegerRet>
set_vertices(const Matrix<IntegerPL> & Verts)897 void ProjectAndLift<IntegerPL, IntegerRet>::set_vertices(const Matrix<IntegerPL>& Verts) {
898 Vertices = Verts;
899 }
900
901 //---------------------------------------------------------------------------
902 template <typename IntegerPL, typename IntegerRet>
compute(bool all_points,bool lifting_float,bool do_only_count)903 void ProjectAndLift<IntegerPL, IntegerRet>::compute(bool all_points, bool lifting_float, bool do_only_count) {
904 // Project-and-lift for lattice points in a polytope.
905 // The first coordinate is homogenizing. Its value for polytope points ism set by GD so that
906 // a grading denominator 1=1 can be accomodated.
907 // We need only the support hyperplanes Supps and the facet-vertex incidence matrix Ind.
908 // Its rows correspond to facets.
909
910 #ifdef NMZ_EXTENDED_TESTS
911 if(!using_GMP<IntegerRet>() && !using_renf<IntegerRet>() && test_arith_overflow_proj_and_lift)
912 throw ArithmeticException(0);
913 #endif
914
915 assert(all_points || !lifting_float); // only all points allowed with float
916
917 assert(all_points || !do_only_count); // counting maks only sense for all points
918
919 if (use_LLL) {
920 LLL_coordinates_without_1st_col(LLL_Coordinates, AllSupps[EmbDim], Vertices, verbose);
921 // Note: LLL_Coordinates is of type IntegerRet.
922 Matrix<IntegerPL>
923 Aconv; // we cannot use to_sublattice_dual directly (not even with convert) since the integer types may not match
924 convert(Aconv, LLL_Coordinates.getEmbeddingMatrix());
925 // Aconv.transpose().pretty_print(cout);
926 AllSupps[EmbDim] = AllSupps[EmbDim].multiplication(Aconv.transpose());
927
928 if (Congs.nr_of_rows() > 0) { // must also transform congruences
929 vector<IntegerRet> Moduli(Congs.nr_of_rows());
930 for (size_t i = 0; i < Congs.nr_of_rows(); ++i)
931 Moduli[i] = Congs[i][Congs.nr_of_columns() - 1];
932 Matrix<IntegerRet> WithoutModuli(0, Congs.nr_of_columns() - 1);
933 for (size_t i = 0; i < Congs.nr_of_rows(); ++i) {
934 vector<IntegerRet> trans = Congs[i];
935 trans.resize(trans.size() - 1);
936 WithoutModuli.append(trans);
937 }
938 Congs = LLL_Coordinates.to_sublattice_dual(WithoutModuli);
939 Congs.insert_column(Congs.nr_of_columns(), Moduli);
940 }
941 if (Grading.size() > 0)
942 Grading = LLL_Coordinates.to_sublattice_dual_no_div(Grading);
943 }
944
945 count_only = do_only_count; // count_only belongs to *this
946
947 if (verbose)
948 verboseOutput() << "Projection" << endl;
949 compute_projections(EmbDim, 1, StartInd, StartPair, StartParaInPair, StartRank);
950 if (all_points) {
951 if (verbose)
952 verboseOutput() << "Lifting" << endl;
953 if (!lifting_float || (lifting_float && using_float<IntegerPL>())) {
954 compute_latt_points();
955 }
956 else {
957 compute_latt_points_float(); // with intermediate conversion to float
958 }
959 /* if(verbose)
960 verboseOutput() << "Number of lattice points " << TotalNrLP << endl;*/
961 }
962 else {
963 if (verbose)
964 verboseOutput() << "Try finding a lattice point" << endl;
965 find_single_point();
966 }
967
968 // cout << " POS " << h_vec_pos;
969 // cout << " NEG " << h_vec_neg;
970 }
971
972 //---------------------------------------------------------------------------
973 template <typename IntegerPL, typename IntegerRet>
compute_only_projection(size_t down_to)974 void ProjectAndLift<IntegerPL, IntegerRet>::compute_only_projection(size_t down_to) {
975 assert(down_to >= 1);
976 compute_projections(EmbDim, down_to, StartInd, StartPair, StartParaInPair, StartRank, true);
977 }
978
979 //---------------------------------------------------------------------------
980 template <typename IntegerPL, typename IntegerRet>
put_eg1Points_into(Matrix<IntegerRet> & LattPoints)981 void ProjectAndLift<IntegerPL, IntegerRet>::put_eg1Points_into(Matrix<IntegerRet>& LattPoints) {
982 while (!Deg1Points.empty()) {
983 if (use_LLL) {
984 /* cout << "ori " << Deg1Points.front();
985 cout << "tra " << LLL_Coordinates.from_sublattice(Deg1Points.front());
986 cout << "tra1 " << LLL_Coordinates.A.VxM(Deg1Points.front());*/
987 LattPoints.append(LLL_Coordinates.from_sublattice(Deg1Points.front()));
988 }
989 else
990 LattPoints.append(Deg1Points.front());
991 Deg1Points.pop_front();
992 }
993 }
994
995 //---------------------------------------------------------------------------
996 template <typename IntegerPL, typename IntegerRet>
put_single_point_into(vector<IntegerRet> & LattPoint)997 void ProjectAndLift<IntegerPL, IntegerRet>::put_single_point_into(vector<IntegerRet>& LattPoint) {
998 if (use_LLL)
999 LattPoint = LLL_Coordinates.from_sublattice(SingleDeg1Point);
1000 else
1001 LattPoint = SingleDeg1Point;
1002 }
1003
1004 //---------------------------------------------------------------------------
1005 template <typename IntegerPL, typename IntegerRet>
getNumberLatticePoints() const1006 size_t ProjectAndLift<IntegerPL, IntegerRet>::getNumberLatticePoints() const {
1007 return TotalNrLP;
1008 }
1009
1010 //---------------------------------------------------------------------------
1011 template <typename IntegerPL, typename IntegerRet>
get_h_vectors(vector<num_t> & pos,vector<num_t> & neg) const1012 void ProjectAndLift<IntegerPL, IntegerRet>::get_h_vectors(vector<num_t>& pos, vector<num_t>& neg) const {
1013 pos = h_vec_pos;
1014 neg = h_vec_neg;
1015 }
1016
1017 //---------------------------------------------------------------------------
1018 // For projection of cones
1019 template <typename IntegerPL, typename IntegerRet>
putSuppsAndEqus(Matrix<IntegerPL> & SuppsRet,Matrix<IntegerPL> & EqusRet,size_t in_dim)1020 void ProjectAndLift<IntegerPL, IntegerRet>::putSuppsAndEqus(Matrix<IntegerPL>& SuppsRet,
1021 Matrix<IntegerPL>& EqusRet,
1022 size_t in_dim) {
1023 assert(in_dim < EmbDim);
1024 assert(in_dim > 0);
1025
1026 EqusRet.resize(0, in_dim); // to make it well-defined
1027 size_t equs_start_in_row = AllSupps[in_dim].nr_of_rows() - 2 * AllNrEqus[in_dim];
1028 for (size_t i = equs_start_in_row; i < AllSupps[in_dim].nr_of_rows(); i += 2) // equations come in +- pairs
1029 EqusRet.append(AllSupps[in_dim][i]);
1030 AllSupps[in_dim].swap(SuppsRet);
1031 // SuppsRet.resize_colums(equs_start_in_row,in_dim);
1032 SuppsRet.resize(equs_start_in_row); // we must delete the superfluous rows because the transformation
1033 // to vector<vector> could else fail.
1034 }
1035 //---------------------------------------------------------------------------
1036
1037 template class ProjectAndLift<mpz_class, mpz_class>;
1038 template class ProjectAndLift<long, long long>;
1039 template class ProjectAndLift<mpz_class, long long>;
1040 template class ProjectAndLift<long long, long long>;
1041 template class ProjectAndLift<nmz_float, mpz_class>;
1042 template class ProjectAndLift<nmz_float, long long>;
1043 // template class ProjectAndLift<nmz_float, nmz_float>;
1044 #ifndef NMZ_MIC_OFFLOAD // offload with long is not supported
1045 template class ProjectAndLift<long, long>;
1046 template class ProjectAndLift<nmz_float, long>;
1047 #endif
1048
1049 #ifdef ENFNORMALIZ
1050 template class ProjectAndLift<renf_elem_class, mpz_class>;
1051 #endif
1052
1053 } // end namespace libnormaliz
1054