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 //---------------------------------------------------------------------------
25 #include <vector>
26
27 #include "libnormaliz/reduction.h"
28 #include "libnormaliz/vector_operations.h"
29
30 namespace libnormaliz {
31 using namespace std;
32
33 //---------------------------------------------------------------------------
34
35 template <typename Integer>
Candidate(const vector<Integer> & v,const vector<Integer> & val,long sd)36 Candidate<Integer>::Candidate(const vector<Integer>& v, const vector<Integer>& val, long sd)
37 : cand(v), values(val), sort_deg(sd), reducible(true), original_generator(false) {
38 }
39
40 //---------------------------------------------------------------------------
41
42 template <typename Integer>
Candidate(const vector<Integer> & v,const Full_Cone<Integer> & C)43 Candidate<Integer>::Candidate(const vector<Integer>& v, const Full_Cone<Integer>& C) : cand(v) {
44 compute_values_deg(C);
45 reducible = true;
46 original_generator = false;
47 }
48
49 //---------------------------------------------------------------------------
50
51 template <typename Integer>
Candidate(const vector<Integer> & v,size_t max_size)52 Candidate<Integer>::Candidate(const vector<Integer>& v, size_t max_size) {
53 cand = v;
54 values.resize(max_size, 0);
55 sort_deg = 0;
56 reducible = true;
57 original_generator = false;
58 }
59
60 //---------------------------------------------------------------------------
61
62 template <typename Integer>
Candidate(size_t cand_size,size_t val_size)63 Candidate<Integer>::Candidate(size_t cand_size, size_t val_size) {
64 // cand=v;
65 values.resize(val_size, 0);
66 cand.resize(cand_size, 0);
67 sort_deg = 0;
68 reducible = true;
69 original_generator = false;
70 }
71
72 //---------------------------------------------------------------------------
73
74 template <typename Integer>
sum(const Candidate<Integer> & C,const Candidate<Integer> & D)75 Candidate<Integer> sum(const Candidate<Integer>& C, const Candidate<Integer>& D) {
76 Candidate<Integer> the_sum = C;
77 the_sum.cand = v_add(the_sum.cand, D.cand);
78 the_sum.values = v_add(the_sum.values, D.values);
79 the_sum.sort_deg += D.sort_deg;
80 the_sum.original_generator = false;
81 the_sum.reducible = true;
82 return the_sum;
83 }
84
85 template Candidate<long> sum(const Candidate<long>&, const Candidate<long>&);
86 template Candidate<long long> sum(const Candidate<long long>&, const Candidate<long long>&);
87 template Candidate<mpz_class> sum(const Candidate<mpz_class>&, const Candidate<mpz_class>&);
88 #ifdef ENFNORMALIZ
89 template Candidate<renf_elem_class> sum(const Candidate<renf_elem_class>&, const Candidate<renf_elem_class>&);
90 #endif
91
92 //---------------------------------------------------------------------------
93
94 template <typename Integer>
compute_values_deg(const Full_Cone<Integer> & C)95 void Candidate<Integer>::compute_values_deg(const Full_Cone<Integer>& C) {
96 C.Support_Hyperplanes.MxV(values, cand);
97 convert(sort_deg, v_scalar_product(cand, C.Sorting));
98 if (C.do_module_gens_intcl || C.hilbert_basis_rec_cone_known) // necessary to make all monoid generators subtractible
99 sort_deg *= 2;
100 }
101
102 template <>
compute_values_deg(const Full_Cone<renf_elem_class> & C)103 void Candidate<renf_elem_class>::compute_values_deg(const Full_Cone<renf_elem_class>& C) {
104 assert(false);
105 }
106
107 //---------------------------------------------------------------------------
108
109 template <typename Integer>
CandidateList()110 CandidateList<Integer>::CandidateList() : tmp_candidate(0, 0) {
111 verbose = false;
112 dual = false;
113 last_hyp = 0;
114 }
115
116 //---------------------------------------------------------------------------
117
118 template <typename Integer>
CandidateList(bool dual_algorithm)119 CandidateList<Integer>::CandidateList(bool dual_algorithm) : tmp_candidate(0, 0) {
120 verbose = false;
121 dual = dual_algorithm;
122 last_hyp = 0;
123 }
124
125 // size_t NrCompVect=0;
126 // size_t NrCompVal=0;
127
128 //---------------------------------------------------------------------------
129
130 //---------------------------------------------------------------------------
131
132 template <typename Integer>
divide_sortdeg_by2()133 void CandidateList<Integer>::divide_sortdeg_by2() {
134 for (auto& r : Candidates)
135 r.sort_deg /= 2;
136 }
137
138 template <typename Integer>
is_reducible(const vector<Integer> & values,const long sort_deg) const139 bool CandidateList<Integer>::is_reducible(const vector<Integer>& values, const long sort_deg) const {
140 long sd;
141 /* if(dual)
142 sd=sort_deg;
143 else */
144 sd = sort_deg / 2;
145 size_t kk = 0;
146 for (const auto& r : Candidates) {
147 /* #pragma omp atomic
148 NrCompVect++;
149 #pragma omp atomic
150 NrCompVal++; */
151 if (sd < r.sort_deg) {
152 return (false);
153 }
154 /* #pragma omp atomic
155 NrCompVal++;*/
156 size_t i = 0;
157 if (values[kk] < r.values[kk])
158 continue;
159 for (; i < values.size(); ++i) {
160 /* #pragma omp atomic
161 NrCompVal++; */
162 if (values[i] < r.values[i]) {
163 kk = i;
164 break;
165 }
166 }
167 if (i == values.size()) {
168 return (true);
169 }
170 }
171 return (false);
172 }
173
174 //---------------------------------------------------------------------------
175
176 template <typename Integer>
is_reducible(Candidate<Integer> & c) const177 bool CandidateList<Integer>::is_reducible(Candidate<Integer>& c) const {
178 /*if(dual && c.in_HB)
179 c.reducible=false;
180 else */
181 c.reducible = is_reducible(c.values, c.sort_deg);
182 return (c.reducible);
183 }
184
185 //---------------------------------------------------------------------------
186
187 /* // not used at present
188 template <typename Integer>
189 bool CandidateList<Integer>::is_reducible(vector<Integer> v, Candidate<Integer>& cand, const Full_Cone<Integer>& C) const {
190 cand = Candidate<Integer>(v, C);
191 return ((*this).is_reducible(cand));
192 }
193 */
194
195 //---------------------------------------------------------------------------
196
197 // Fourth version with parallelization and tables
198 template <typename Integer>
reduce_by(CandidateList<Integer> & Reducers)199 void CandidateList<Integer>::reduce_by(CandidateList<Integer>& Reducers) {
200 size_t cpos, csize = Candidates.size();
201
202 CandidateTable<Integer> ReducerTable(Reducers);
203
204 bool skip_remaining = false;
205 ;
206 std::exception_ptr tmp_exception;
207
208 // This parallel region cannot throw a NormalizException
209 #pragma omp parallel private(cpos) firstprivate(ReducerTable)
210 {
211 auto c = Candidates.begin();
212 cpos = 0;
213
214 #pragma omp for schedule(dynamic)
215 for (size_t k = 0; k < csize; ++k) {
216 for (; k > cpos; ++cpos, ++c)
217 ;
218 for (; k < cpos; --cpos, --c)
219 ;
220
221 if (skip_remaining)
222 continue;
223 try {
224 INTERRUPT_COMPUTATION_BY_EXCEPTION
225
226 ReducerTable.is_reducible(*c);
227
228 } catch (const std::exception&) {
229 tmp_exception = std::current_exception();
230 skip_remaining = true;
231 #pragma omp flush(skip_remaining)
232 }
233 }
234
235 } // end parallel
236
237 if (!(tmp_exception == 0))
238 std::rethrow_exception(tmp_exception);
239
240 // erase reducibles
241 for (auto c = Candidates.begin(); c != Candidates.end();) {
242 if ((*c).reducible)
243 c = Candidates.erase(c);
244 else // continue
245 ++c;
246 }
247 }
248
249 //---------------------------------------------------------------------------
250
251 template <typename Integer>
auto_reduce()252 void CandidateList<Integer>::auto_reduce() {
253 if (empty())
254 return;
255
256 sort_by_deg();
257 auto_reduce_sorted();
258 }
259
260 //---------------------------------------------------------------------------
261
262 template <typename Integer>
auto_reduce_sorted()263 void CandidateList<Integer>::auto_reduce_sorted() {
264 // uses generations defined by degrees
265
266 if (empty())
267 return;
268
269 CandidateList<Integer> Irreducibles(dual), CurrentReducers(dual);
270 Integer irred_degree;
271 size_t cs = Candidates.size();
272 if (verbose && cs > 1000) {
273 verboseOutput() << "auto-reduce " << cs << " candidates, degrees <= ";
274 }
275
276 while (!Candidates.empty()) {
277 irred_degree = Candidates.begin()->sort_deg * 2 - 1;
278 if (verbose && cs > 1000) {
279 verboseOutput() << irred_degree << " " << flush;
280 }
281 auto c = Candidates.begin();
282 while (c != Candidates.end() && c->sort_deg <= irred_degree) {
283 ++c; // find location for splicing
284 }
285 CurrentReducers.Candidates.splice(CurrentReducers.Candidates.begin(), Candidates, Candidates.begin(), c);
286 reduce_by(CurrentReducers);
287 Irreducibles.Candidates.splice(Irreducibles.Candidates.end(), CurrentReducers.Candidates);
288 }
289 if (verbose && cs > 1000) {
290 verboseOutput() << endl;
291 }
292 Candidates.splice(Candidates.begin(), Irreducibles.Candidates);
293 }
294
295 #ifdef ENFNORMALIZ
296 template <>
auto_reduce_sorted()297 void CandidateList<renf_elem_class>::auto_reduce_sorted() {
298 assert(false);
299 }
300 #endif
301 //---------------------------------------------------------------------------
302
303 template <typename Integer>
reduce_by_and_insert(Candidate<Integer> & cand,const CandidateList<Integer> & Reducers)304 bool CandidateList<Integer>::reduce_by_and_insert(Candidate<Integer>& cand, const CandidateList<Integer>& Reducers) {
305 bool irred = !Reducers.is_reducible(cand);
306 if (irred)
307 Candidates.push_back(cand);
308 return irred;
309 }
310
311 //---------------------------------------------------------------------------
312
313 template <typename Integer>
reduce_by_and_insert(const vector<Integer> & v,Full_Cone<Integer> & C,CandidateList<Integer> & Reducers)314 bool CandidateList<Integer>::reduce_by_and_insert(const vector<Integer>& v,
315 Full_Cone<Integer>& C,
316 CandidateList<Integer>& Reducers) {
317 Candidate<Integer> cand(v, C);
318 return reduce_by_and_insert(cand, Reducers);
319 }
320
321 //---------------------------------------------------------------------------
322
323 template <typename Integer>
unique_vectors()324 void CandidateList<Integer>::unique_vectors() {
325 assert(dual);
326
327 if (empty())
328 return;
329
330 // sort_by_val();
331
332 auto h_start = Candidates.begin();
333
334 h_start++;
335 for (auto h = h_start; h != Candidates.end();) {
336 auto prev = h;
337 prev--;
338 if (h->values == prev->values) // since cone may not be pointed in the dual , vectors
339 h = Candidates.erase(h); // must be made unique modulo the unit group
340 else // values gives standard embedding
341 ++h;
342 }
343 }
344
345 //---------------------------------------------------------------------------
346
347 template <typename Integer>
deg_compare(const Candidate<Integer> & a,const Candidate<Integer> & b)348 bool deg_compare(const Candidate<Integer>& a, const Candidate<Integer>& b) {
349 return (a.sort_deg < b.sort_deg);
350 }
351
352 //---------------------------------------------------------------------------
353
354 template <typename Integer>
val_compare(const Candidate<Integer> & a,const Candidate<Integer> & b)355 bool val_compare(const Candidate<Integer>& a, const Candidate<Integer>& b) {
356 if (a.sort_deg < b.sort_deg)
357 return (true);
358 if (a.sort_deg == b.sort_deg) {
359 if (a.values < b.values)
360 return true;
361 if (a.values == b.values)
362 return a.mother < b.mother;
363 }
364 return false;
365 }
366
367 //---------------------------------------------------------------------------
368
369 template <typename Integer>
sort_by_deg()370 void CandidateList<Integer>::sort_by_deg() {
371 Candidates.sort(deg_compare<Integer>);
372 }
373
374 //---------------------------------------------------------------------------
375
376 template <typename Integer>
sort_by_val()377 void CandidateList<Integer>::sort_by_val() {
378 Candidates.sort(val_compare<Integer>);
379 }
380
381 //---------------------------------------------------------------------------
382
383 template <typename Integer>
clear()384 void CandidateList<Integer>::clear() {
385 Candidates.clear();
386 }
387
388 //---------------------------------------------------------------------------
389
390 template <typename Integer>
size()391 size_t CandidateList<Integer>::size() {
392 return Candidates.size();
393 }
394
395 //---------------------------------------------------------------------------
396
397 template <typename Integer>
empty()398 bool CandidateList<Integer>::empty() {
399 return Candidates.empty();
400 }
401
402 /*
403 template<typename Integer>
404 void CandidateList<Integer>::search(){
405
406 vector<Integer> TESTV(6);
407 TESTV[0]=2;
408 TESTV[1]=6;
409 TESTV[2]=0;
410 TESTV[3]=15;
411 TESTV[4]=18;
412 TESTV[5]=2;
413
414 for(const auto& h : Candidates){
415 if(h.cand==TESTV){
416 assert(h.cand!=TESTV);
417 }
418
419 }
420
421 }
422 // */
423
424 //---------------------------------------------------------------------------
425
426 template <typename Integer>
merge(CandidateList<Integer> & NewCand)427 void CandidateList<Integer>::merge(CandidateList<Integer>& NewCand) {
428 Candidates.merge(NewCand.Candidates, deg_compare<Integer>);
429 }
430
431 //---------------------------------------------------------------------------
432
433 template <typename Integer>
merge_by_val(CandidateList<Integer> & NewCand)434 void CandidateList<Integer>::merge_by_val(CandidateList<Integer>& NewCand) {
435 list<Candidate<Integer>*> dummy;
436 merge_by_val_inner(NewCand, false, dummy);
437 }
438
439 //---------------------------------------------------------------------------
440
441 template <typename Integer>
merge_by_val(CandidateList<Integer> & NewCand,list<Candidate<Integer> * > & New_Elements)442 void CandidateList<Integer>::merge_by_val(CandidateList<Integer>& NewCand, list<Candidate<Integer>*>& New_Elements) {
443 CandidateList<Integer> dummy;
444 merge_by_val_inner(NewCand, true, New_Elements);
445 }
446
447 //---------------------------------------------------------------------------
448
449 template <typename Integer>
merge_by_val_inner(CandidateList<Integer> & NewCand,bool collect_new_elements,list<Candidate<Integer> * > & New_Elements)450 void CandidateList<Integer>::merge_by_val_inner(CandidateList<Integer>& NewCand,
451 bool collect_new_elements,
452 list<Candidate<Integer>*>& New_Elements) {
453 CandidateList<Integer> Coll;
454 Coll.dual = dual;
455 Coll.last_hyp = last_hyp;
456
457 while (!empty() || !NewCand.empty()) {
458 if (NewCand.empty()) {
459 Coll.Candidates.splice(Coll.Candidates.begin(), Candidates);
460 break;
461 }
462
463 if (empty()) {
464 typename list<Candidate<Integer> >::reverse_iterator h;
465 if (collect_new_elements) {
466 for (h = NewCand.Candidates.rbegin(); h != NewCand.Candidates.rend(); ++h)
467 New_Elements.push_front(&(*h));
468 }
469 Coll.Candidates.splice(Coll.Candidates.begin(), NewCand.Candidates);
470 break;
471 }
472
473 if (NewCand.Candidates.back().values == Candidates.back().values) { // if equal, new is erased
474 if (NewCand.Candidates.back().mother < Candidates.back().mother)
475 Candidates.back().mother = NewCand.Candidates.back().mother;
476 NewCand.Candidates.pop_back();
477 continue;
478 }
479
480 if (val_compare<Integer>(Candidates.back(), NewCand.Candidates.back())) { // old is smaller, new must be inserteed
481 if (collect_new_elements) {
482 New_Elements.push_front(&(NewCand.Candidates.back()));
483 }
484 Coll.Candidates.splice(Coll.Candidates.begin(), NewCand.Candidates, --NewCand.Candidates.end());
485 continue;
486 }
487
488 Coll.Candidates.splice(Coll.Candidates.begin(), Candidates, --Candidates.end()); // the remaining case
489 }
490
491 splice(Coll); // Coll moved to this
492 }
493
494 //---------------------------------------------------------------------------
495
496 template <typename Integer>
push_back(const Candidate<Integer> & cand)497 void CandidateList<Integer>::push_back(const Candidate<Integer>& cand) {
498 // cout << cand;
499 Candidates.push_back(cand);
500 }
501
502 //---------------------------------------------------------------------------
503
504 template <typename Integer>
extract(list<vector<Integer>> & V_List)505 void CandidateList<Integer>::extract(list<vector<Integer> >& V_List) {
506 for (const auto& c : Candidates)
507 V_List.push_back(c.cand);
508 }
509
510 //---------------------------------------------------------------------------
511
512 template <typename Integer>
splice(CandidateList<Integer> & NewCand)513 void CandidateList<Integer>::splice(CandidateList<Integer>& NewCand) {
514 Candidates.splice(Candidates.begin(), NewCand.Candidates);
515 }
516
517 //---------------------------------------------------------------------------
518
519 template <typename Integer>
CandidateTable(CandidateList<Integer> & CandList)520 CandidateTable<Integer>::CandidateTable(CandidateList<Integer>& CandList) {
521 for (auto& c : CandList.Candidates)
522 ValPointers.push_back(make_pair(c.sort_deg, &(c.values)));
523 dual = CandList.dual;
524 last_hyp = CandList.last_hyp;
525 }
526
527 //---------------------------------------------------------------------------
528
529 template <typename Integer>
CandidateTable(bool dual,size_t last_hyp)530 CandidateTable<Integer>::CandidateTable(bool dual, size_t last_hyp) {
531 this->dual = dual;
532 this->last_hyp = last_hyp;
533 }
534
535 //---------------------------------------------------------------------------
536
537 template <typename Integer>
is_reducible(Candidate<Integer> & c)538 bool CandidateTable<Integer>::is_reducible(Candidate<Integer>& c) {
539 c.reducible = is_reducible(c.values, c.sort_deg);
540 return (c.reducible);
541 }
542
543 //---------------------------------------------------------------------------
544
545 template <typename Integer>
is_reducible(const vector<Integer> & values,const long sort_deg)546 bool CandidateTable<Integer>::is_reducible(const vector<Integer>& values, const long sort_deg) {
547 long sd;
548 /* if(dual)
549 sd=sort_deg;
550 else */
551 sd = sort_deg / 2;
552 size_t kk = 0;
553 for (auto r = ValPointers.begin(); r != ValPointers.end(); ++r) {
554 /* #pragma omp atomic
555 NrCompVect++;
556 #pragma omp atomic
557 NrCompVal++;*/
558 if (sd < (long)r->first) {
559 return (false);
560 }
561 /* #pragma omp atomic
562 NrCompVal++;*/
563 size_t i = 0;
564 if (values[kk] < (*(r->second))[kk])
565 continue;
566 for (; i < values.size(); ++i) {
567 /* #pragma omp atomic
568 NrCompVal++; */
569 if (values[i] < (*(r->second))[i]) {
570 kk = i;
571 break;
572 }
573 }
574 if (i == values.size()) {
575 ValPointers.splice(ValPointers.begin(), ValPointers, r);
576 return (true);
577 }
578 }
579 return (false);
580 }
581
582 //---------------------------------------------------------------------------
583
584 template <typename Integer>
is_reducible_unordered(Candidate<Integer> & c)585 bool CandidateTable<Integer>::is_reducible_unordered(Candidate<Integer>& c) {
586 c.reducible = is_reducible_unordered(c.values, c.sort_deg);
587 return (c.reducible);
588 }
589
590 //---------------------------------------------------------------------------
591
592 template <typename Integer>
is_reducible_unordered(const vector<Integer> & values,const long sort_deg)593 bool CandidateTable<Integer>::is_reducible_unordered(const vector<Integer>& values, const long sort_deg) {
594 long sd;
595 if (dual)
596 sd = sort_deg;
597 else
598 sd = sort_deg / 2;
599 size_t kk = 0;
600 for (auto r = ValPointers.begin(); r != ValPointers.end(); ++r) {
601 if (sd <= (long)r->first) {
602 continue; // in the ordered version we can say: return(false);
603 }
604 // #pragma omp atomic
605 // redcounter++;
606 vector<Integer>* reducer = r->second;
607 if (values[last_hyp] < (*(r->second))[last_hyp])
608 continue;
609 size_t i = 0;
610 if (values[kk] < (*(r->second))[kk])
611 continue;
612 for (; i < last_hyp; ++i)
613 if (values[i] < (*reducer)[i]) {
614 kk = i;
615 break;
616 }
617 if (i == last_hyp) {
618 ValPointers.splice(ValPointers.begin(), ValPointers, r);
619 return (true);
620 }
621 }
622 return (false);
623 }
624
625 #ifndef NMZ_MIC_OFFLOAD // offload with long is not supported
626 template class CandidateList<long>;
627 template class CandidateTable<long>;
628 template class Candidate<long>;
629 #endif
630
631 template class CandidateList<long long>;
632 template class CandidateTable<long long>;
633 template class Candidate<long long>;
634
635 template class CandidateList<mpz_class>;
636 template class CandidateTable<mpz_class>;
637 template class Candidate<mpz_class>;
638
639 #ifdef ENFNORMALIZ
640 template class CandidateList<renf_elem_class>;
641 template class CandidateTable<renf_elem_class>;
642 template class Candidate<renf_elem_class>;
643 #endif
644
645 size_t redcounter = 0;
646
647 } // namespace libnormaliz
648