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 #ifdef NMZ_MIC_OFFLOAD
25 #pragma offload_attribute(push, target(mic))
26 #endif
27
28 #include <cassert>
29 #include <iostream>
30 #include <sstream>
31 #include <map>
32 #include <algorithm>
33
34 #include "libnormaliz/general.h"
35 #include "libnormaliz/HilbertSeries.h"
36 #include "libnormaliz/vector_operations.h"
37 #include "libnormaliz/list_and_map_operations.h"
38 #include "libnormaliz/integer.h"
39 // #include "libnormaliz/convert.h"
40
41 #include "libnormaliz/matrix.h"
42
43 #ifdef NMZ_FLINT
44 #include "flint/flint.h"
45 #include "flint/fmpz_poly.h"
46 #endif
47
48 //---------------------------------------------------------------------------
49
50 namespace libnormaliz {
51 using std::cout;
52 using std::endl;
53 using std::flush;
54 using std::istringstream;
55 using std::ostringstream;
56 using std::pair;
57
58 //---------------------------------------------------------------------------
59
60 template <typename Integer>
permutations(const size_t & a,const size_t & b)61 Integer permutations(const size_t& a, const size_t& b) {
62 unsigned long i;
63 Integer P = 1;
64 for (i = a + 1; i <= b; i++) {
65 P *= i;
66 }
67 return P;
68 }
69
70 #ifdef NMZ_FLINT
flint_poly(fmpz_poly_t flp,const vector<mpz_class> & nmzp)71 void flint_poly(fmpz_poly_t flp, const vector<mpz_class>& nmzp) {
72 slong n = (slong)nmzp.size();
73 fmpz_poly_fit_length(flp, n);
74 for (size_t i = 0; i < nmzp.size(); ++i) {
75 fmpz_poly_set_coeff_mpz(flp, (slong)i, nmzp[i].get_mpz_t());
76 }
77 }
78
nmz_poly(vector<mpz_class> & nmzp,const fmpz_poly_t flp)79 void nmz_poly(vector<mpz_class>& nmzp, const fmpz_poly_t flp) {
80 size_t n = (size_t)fmpz_poly_length(flp);
81 nmzp.resize(n);
82 mpz_t c;
83 mpz_init(c);
84 for (size_t i = 0; i < nmzp.size(); ++i) {
85 fmpz_poly_get_coeff_mpz(c, flp, i);
86 nmzp[i] = mpz_class(c);
87 }
88 mpz_clear(c);
89 }
90 #endif
91
92 template <typename Integer>
poly_mult(const vector<Integer> & a,const vector<Integer> & b)93 vector<Integer> poly_mult(const vector<Integer>& a, const vector<Integer>& b) {
94 size_t a_size = a.size();
95 size_t b_size = b.size();
96
97 if (a_size * b_size > 1000 && a_size > 10 && b_size > 10) {
98 // omp_set_nested(1);
99 return karatsubamult(a, b);
100 // omp_set_nested(0);
101 }
102
103 vector<Integer> p(a_size + b_size - 1);
104 size_t i, j;
105 for (i = 0; i < a_size; ++i) {
106 if (a[i] == 0)
107 continue;
108 for (j = 0; j < b_size; ++j) {
109 if (b[j] == 0)
110 continue;
111 p[i + j] += a[i] * b[j];
112 }
113 }
114 return p;
115 }
116
117 #ifdef NMZ_FLINT
118 template <>
poly_mult(const vector<mpz_class> & a,const vector<mpz_class> & b)119 vector<mpz_class> poly_mult(const vector<mpz_class>& a, const vector<mpz_class>& b) {
120 size_t a_size = a.size();
121 size_t b_size = b.size();
122
123 vector<mpz_class> p(a_size + b_size - 1);
124 fmpz_poly_t flp1, flp2;
125 fmpz_poly_init(flp1);
126 fmpz_poly_init(flp2);
127
128 flint_poly(flp1, a);
129 flint_poly(flp2, b);
130 fmpz_poly_mul(flp1, flp1, flp2);
131 nmz_poly(p, flp1);
132
133 fmpz_poly_clear(flp1);
134 fmpz_poly_clear(flp2);
135
136 return p;
137 }
138 #endif
139
140 // division with remainder, a = q * b + r, deg(r) < deg(b), needs |leadcoef(b)| = 1
141 template <typename Integer>
poly_div(vector<Integer> & q,vector<Integer> & r,const vector<Integer> & a,const vector<Integer> & b)142 void poly_div(vector<Integer>& q, vector<Integer>& r, const vector<Integer>& a, const vector<Integer>& b) {
143 assert(b.back() != 0); // no unneeded zeros
144 assert(b.back() == 1 || b.back() == -1); // then division is always possible
145 r = a;
146 remove_zeros(r);
147 size_t b_size = b.size();
148 int degdiff = r.size() - b_size; // degree differenz
149 if (r.size() < b_size) {
150 q = vector<Integer>();
151 }
152 else {
153 q = vector<Integer>(degdiff + 1);
154 }
155 Integer divisor;
156 size_t i = 0;
157
158 while (r.size() >= b_size) {
159 divisor = r.back() / b.back();
160 q[degdiff] = divisor;
161 // r -= divisor * t^degdiff * b
162 for (i = 0; i < b_size; ++i) {
163 r[i + degdiff] -= divisor * b[i];
164 }
165 remove_zeros(r);
166 degdiff = r.size() - b_size;
167 }
168
169 return;
170 }
171
172 #ifdef NMZ_FLINT
173 template <>
poly_div(vector<mpz_class> & q,vector<mpz_class> & r,const vector<mpz_class> & a,const vector<mpz_class> & b)174 void poly_div(vector<mpz_class>& q, vector<mpz_class>& r, const vector<mpz_class>& a, const vector<mpz_class>& b) {
175 assert(b.back() != 0); // no unneeded zeros
176 assert(b.back() == 1 || b.back() == -1); // then division is always possible
177
178 fmpz_poly_t flpa, flpb, flpq, flpr;
179 fmpz_poly_init(flpa);
180 fmpz_poly_init(flpb);
181 fmpz_poly_init(flpq);
182 fmpz_poly_init(flpr);
183
184 flint_poly(flpa, a);
185 flint_poly(flpb, b);
186
187 fmpz_poly_divrem(flpq, flpr, flpa, flpb);
188 nmz_poly(q, flpq);
189 nmz_poly(r, flpr);
190
191 fmpz_poly_clear(flpa);
192 fmpz_poly_clear(flpb);
193 fmpz_poly_clear(flpq);
194 fmpz_poly_clear(flpr);
195
196 return;
197 }
198 #endif
199
200 template <typename Integer>
cyclotomicPoly(long n)201 vector<Integer> cyclotomicPoly(long n) {
202 // the static variable is initialized only once and then stored
203 static map<long, vector<Integer> > CyclotomicPoly = map<long, vector<Integer> >();
204 if (CyclotomicPoly.count(n) == 0) { // it was not computed so far
205 vector<Integer> poly, q, r;
206 for (long i = 1; i <= n; ++i) {
207 // compute needed and uncomputed factors
208 if (n % i == 0 && CyclotomicPoly.count(i) == 0) {
209 // compute the i-th poly by dividing X^i-1 by the
210 // d-th cycl.poly. with d divides i
211 poly = vector<Integer>(i + 1);
212 poly[0] = -1;
213 poly[i] = 1; // X^i - 1
214 for (long d = 1; d < i; ++d) { // <= i/2 should be ok
215 if (i % d == 0) {
216 poly_div(q, r, poly, CyclotomicPoly[d]);
217 assert(r.empty());
218 poly = q;
219 }
220 }
221 CyclotomicPoly[i] = poly;
222 // cout << i << "-th cycl. pol.: " << CyclotomicPoly[i];
223 }
224 }
225 }
226 assert(CyclotomicPoly.count(n) > 0);
227 return CyclotomicPoly[n];
228 }
229
230 #ifdef NMZ_FLINT
231 template <>
cyclotomicPoly(long n)232 vector<mpz_class> cyclotomicPoly(long n) {
233 // the static variable is initialized only once and then stored
234 static map<long, vector<mpz_class> > CyclotomicPoly = map<long, vector<mpz_class> >();
235 if (CyclotomicPoly.count(n) == 0) { // it was not computed so far
236 vector<mpz_class> poly;
237 fmpz_poly_t cyc;
238 fmpz_poly_init(cyc);
239 fmpz_poly_cyclotomic(cyc, (ulong)n);
240 nmz_poly(poly, cyc);
241 CyclotomicPoly[n] = poly;
242 fmpz_poly_clear(cyc);
243 // cout << i << "-th cycl. pol.: " << CyclotomicPoly[i];
244 }
245 assert(CyclotomicPoly.count(n) > 0);
246 return CyclotomicPoly[n];
247 }
248 #endif
249
lcm_of_keys(const map<long,denom_t> & m)250 long lcm_of_keys(const map<long, denom_t>& m) {
251 long l = 1;
252 for (const auto& it : m) {
253 if (it.second != 0)
254 l = lcm(l, it.first);
255 }
256 return l;
257 }
258
259 // compute the hsop numerator by multiplying the HS with a denominator
260 // of the form product of (1-t^i)
compute_hsop_num() const261 void HilbertSeries::compute_hsop_num() const {
262 // get the denominator as a polynomial by mutliplying the (1-t^i) terms
263 vector<mpz_class> hsop_denom_poly = vector<mpz_class>(1, 1);
264 long factor;
265 for (auto& it : hsop_denom) {
266 factor = it.first;
267 denom_t& denom_i = it.second;
268 poly_mult_to(hsop_denom_poly, factor, denom_i);
269 }
270 // cout << "new denominator as polynomial: " << hsop_denom_poly << endl;
271 vector<mpz_class> quot, remainder, cyclo_poly;
272 // first divide the new denom by the cyclo polynomials
273 for (auto& it : cyclo_denom) {
274 for (long i = 0; i < it.second; i++) {
275 cyclo_poly = cyclotomicPoly<mpz_class>(it.first);
276 // cout << "the cyclotomic polynomial is " << cyclo_poly << endl;
277 // TODO: easier polynomial division possible?
278 poly_div(quot, remainder, hsop_denom_poly, cyclo_poly);
279 // cout << "the quotient is " << quot << endl;
280 hsop_denom_poly = quot;
281 assert(remainder.size() == 0);
282 }
283 }
284 // multiply with the old numerator
285 hsop_num = poly_mult(hsop_denom_poly, cyclo_num);
286 }
287
288 //---------------------------------------------------------------------------
289
initialize()290 void HilbertSeries::initialize() {
291 is_simplified = false;
292 shift = 0;
293 verbose = false;
294 nr_coeff_quasipol = -1; // all coefficients
295 expansion_degree = -1;
296 period_bounded = true;
297 }
298
299 // Constructor, creates 0/1
HilbertSeries()300 HilbertSeries::HilbertSeries() {
301 num = vector<mpz_class>(1, 0);
302 // denom just default constructed
303 initialize();
304 }
305
306 // Constructor, creates num/denom, see class description for format
HilbertSeries(const vector<num_t> & numerator,const vector<denom_t> & gen_degrees)307 HilbertSeries::HilbertSeries(const vector<num_t>& numerator, const vector<denom_t>& gen_degrees) {
308 num = vector<mpz_class>(1, 0);
309 add(numerator, gen_degrees);
310 initialize();
311 }
312
313 // Constructor, creates num/denom, see class description for format
HilbertSeries(const vector<mpz_class> & numerator,const map<long,denom_t> & denominator)314 HilbertSeries::HilbertSeries(const vector<mpz_class>& numerator, const map<long, denom_t>& denominator) {
315 num = numerator;
316 denom = denominator;
317 initialize();
318 }
319
320 /*
321 // Constructor, string as created by to_string_rep
322 HilbertSeries::HilbertSeries(const string& str) {
323 from_string_rep(str);
324 initialize();
325 }
326 */
327
reset()328 void HilbertSeries::reset() {
329 num.clear();
330 num.push_back(0);
331 denom.clear();
332 denom_classes.clear();
333 shift = 0;
334 is_simplified = false;
335 }
336
set_nr_coeff_quasipol(long nr_coeff)337 void HilbertSeries::set_nr_coeff_quasipol(long nr_coeff) {
338 nr_coeff_quasipol = nr_coeff;
339 }
340
get_nr_coeff_quasipol() const341 long HilbertSeries::get_nr_coeff_quasipol() const {
342 return nr_coeff_quasipol;
343 }
344
set_period_bounded(bool on_off) const345 void HilbertSeries::set_period_bounded(bool on_off) const { // period_bounded is mutable
346 period_bounded = on_off;
347 }
348
get_period_bounded() const349 bool HilbertSeries::get_period_bounded() const {
350 return period_bounded;
351 }
352 // add another HilbertSeries to this
add(const vector<num_t> & num,const vector<denom_t> & gen_degrees)353 void HilbertSeries::add(const vector<num_t>& num, const vector<denom_t>& gen_degrees) {
354 vector<denom_t> sorted_gd(gen_degrees);
355 sort(sorted_gd.begin(), sorted_gd.end());
356 if (gen_degrees.size() > 0)
357 assert(sorted_gd[0] > 0); // TODO InputException?
358 poly_add_to(denom_classes[sorted_gd], num);
359 if (denom_classes.size() > DENOM_CLASSES_BOUND)
360 collectData();
361 is_simplified = false;
362 }
363
364 // add another HilbertSeries to this
operator +=(const HilbertSeries & other)365 HilbertSeries& HilbertSeries::operator+=(const HilbertSeries& other) {
366 // add denom_classes
367 for (auto& denom_class : other.denom_classes) {
368 poly_add_to(denom_classes[denom_class.first], denom_class.second);
369 }
370 // add accumulated data
371 vector<mpz_class> num_copy(other.num);
372 performAdd(num_copy, other.denom);
373 return (*this);
374 }
375
performAdd(const vector<num_t> & numerator,const vector<denom_t> & gen_degrees) const376 void HilbertSeries::performAdd(const vector<num_t>& numerator, const vector<denom_t>& gen_degrees) const {
377 map<long, denom_t> other_denom;
378 size_t i, s = gen_degrees.size();
379 for (i = 0; i < s; ++i) {
380 assert(gen_degrees[i] > 0);
381 other_denom[gen_degrees[i]]++;
382 }
383 // convert numerator to mpz
384 vector<mpz_class> other_num(numerator.size());
385 convert(other_num, numerator);
386 performAdd(other_num, other_denom);
387 }
388
389 // modifies other_num!!
performAdd(vector<mpz_class> & other_num,const map<long,denom_t> & oth_denom) const390 void HilbertSeries::performAdd(vector<mpz_class>& other_num, const map<long, denom_t>& oth_denom) const {
391 map<long, denom_t> other_denom(oth_denom); // TODO redesign, dont change other_denom
392 // adjust denominators
393 denom_t diff;
394 for (auto& it : denom) { // augment other
395 denom_t& ref = other_denom[it.first];
396 diff = it.second - ref;
397 if (diff > 0) {
398 ref += diff;
399 poly_mult_to(other_num, it.first, diff);
400 }
401 }
402 for (auto& it : other_denom) { // augment this
403 denom_t& ref = denom[it.first];
404 diff = it.second - ref;
405 if (diff > 0) {
406 ref += diff;
407 poly_mult_to(num, it.first, diff);
408 }
409 }
410 assert(denom == other_denom);
411
412 // now just add the numerators
413 poly_add_to(num, other_num);
414 remove_zeros(num);
415 is_simplified = false;
416 }
417
collectData() const418 void HilbertSeries::collectData() const {
419 if (denom_classes.empty())
420 return;
421 if (verbose)
422 verboseOutput() << "Adding " << denom_classes.size() << " denominator classes..." << flush;
423 for (auto& denom_class : denom_classes) {
424 INTERRUPT_COMPUTATION_BY_EXCEPTION
425
426 performAdd(denom_class.second, denom_class.first);
427 }
428 denom_classes.clear();
429 if (verbose)
430 verboseOutput() << " done." << endl;
431 }
432
433 // simplify, see class description
simplify() const434 void HilbertSeries::simplify() const {
435 if (is_simplified)
436 return;
437 collectData();
438 /* if (verbose) {
439 verboseOutput() << "Hilbert series before simplification: "<< endl << *this;
440 }*/
441 vector<mpz_class> q, r, poly; // polynomials
442 // In denom_cyclo we collect cyclotomic polynomials in the denominator.
443 // During this method the Hilbert series is given by num/(denom*cdenom)
444 // where denom | cdenom are exponent vectors of (1-t^i) | i-th cyclotminc poly.
445 map<long, denom_t> cdenom;
446
447 map<long, denom_t> save_denom = denom;
448 vector<mpz_class> save_num = num;
449 map<long, denom_t>::reverse_iterator rit;
450 long i;
451 for (rit = denom.rbegin(); rit != denom.rend(); ++rit) {
452 // check if we can divide the numerator by (1-t^i)
453 i = rit->first;
454 denom_t& denom_i = rit->second;
455 poly = coeff_vector<mpz_class>(i);
456 while (denom_i > 0) {
457 poly_div(q, r, num, poly);
458 if (r.size() == 0) { // numerator is divisable by poly
459 num = q;
460 denom_i--;
461 }
462 else {
463 break;
464 }
465 }
466 if (denom_i == 0)
467 continue;
468
469 // decompose (1-t^i) into cyclotomic polynomial
470 for (long d = 1; d <= i / 2; ++d) {
471 if (i % d == 0)
472 cdenom[d] += denom_i;
473 }
474 cdenom[i] += denom_i;
475 // the product of the cyclo. is t^i-1 = -(1-t^i)
476 if (denom_i % 2 == 1)
477 v_scalar_multiplication(num, mpz_class(-1));
478 } // end for
479 denom.clear();
480
481 auto it = cdenom.begin();
482 while (it != cdenom.end()) {
483 // check if we can divide the numerator by i-th cyclotomic polynomial
484
485 INTERRUPT_COMPUTATION_BY_EXCEPTION
486
487 i = it->first;
488 denom_t& cyclo_i = it->second;
489 poly = cyclotomicPoly<mpz_class>(i);
490 while (cyclo_i > 0) {
491 poly_div(q, r, num, poly);
492 if (r.empty()) { // numerator is divisable by poly
493 num = q;
494 cyclo_i--;
495 }
496 else {
497 break;
498 }
499 }
500
501 if (cyclo_i == 0) {
502 cdenom.erase(it++);
503 }
504 else {
505 ++it;
506 }
507 }
508 // done with canceling
509 // save this representation
510 cyclo_num = num;
511 cyclo_denom = cdenom;
512
513 // now collect the cyclotomic polynomials in (1-t^i) factors
514 it = cdenom.find(1);
515 if (it != cdenom.end())
516 dim = it->second;
517 else
518 dim = 0;
519 period = lcm_of_keys(cdenom);
520 if (period_bounded && period > 10 * PERIOD_BOUND) {
521 if (verbose) {
522 errorOutput() << "WARNING: Period is too big, the representation of the Hilbert series may have more than dimension "
523 "many factors in the denominator!"
524 << endl;
525 }
526 denom = save_denom;
527 num = save_num;
528 }
529 else {
530 while (true) {
531 // create a (1-t^k) factor in the denominator out of all cyclotomic poly.
532
533 INTERRUPT_COMPUTATION_BY_EXCEPTION
534
535 long k = 1;
536 bool empty = true;
537 vector<mpz_class> existing_factor(1, 1); // collects the existing cyclotomic gactors in the denom
538 for (auto& it : cdenom) { // with multiplicvity 1
539 if (it.second > 0) {
540 empty = false;
541 k = libnormaliz::lcm(k, it.first);
542 existing_factor = poly_mult(existing_factor, cyclotomicPoly<mpz_class>(it.first));
543 it.second--;
544 }
545 }
546 if (empty)
547 break;
548 denom[k]++;
549 vector<mpz_class> new_factor = coeff_vector<mpz_class>(k);
550 vector<mpz_class> quotient, dummy;
551 poly_div(quotient, dummy, new_factor, existing_factor);
552 assert(dummy.empty()); // assert remainder r is 0
553 num = poly_mult(num, quotient);
554 }
555 }
556
557 /* if (verbose) {
558 verboseOutput() << "Simplified Hilbert series: " << endl << *this;
559 }*/
560 if (!hsop_denom.empty()) {
561 compute_hsop_num();
562 }
563 else {
564 if (denom.empty()) { // this takes care of the exceptional case in wgich the series
565 hsop_num = num; // is a polynomial
566 }
567 }
568 is_simplified = true;
569 computeDegreeAsRationalFunction();
570 quasi_poly.clear();
571 }
572
computeDegreeAsRationalFunction() const573 void HilbertSeries::computeDegreeAsRationalFunction() const {
574 simplify();
575 long num_deg = num.size() - 1 + shift;
576 long denom_deg = 0;
577 for (auto& it : denom) {
578 denom_deg += it.first * it.second;
579 }
580 degree = num_deg - denom_deg;
581 }
582
getDegreeAsRationalFunction() const583 long HilbertSeries::getDegreeAsRationalFunction() const {
584 simplify();
585 return degree;
586 }
587
getPeriod() const588 long HilbertSeries::getPeriod() const {
589 simplify();
590 return period;
591 }
592
isHilbertQuasiPolynomialComputed() const593 bool HilbertSeries::isHilbertQuasiPolynomialComputed() const {
594 return is_simplified && !quasi_poly.empty();
595 }
596
resetHilbertQuasiPolynomial()597 void HilbertSeries::resetHilbertQuasiPolynomial() {
598 quasi_poly.clear();
599 }
600
getHilbertQuasiPolynomial() const601 const vector<vector<mpz_class> >& HilbertSeries::getHilbertQuasiPolynomial() const {
602 computeHilbertQuasiPolynomial();
603 if (quasi_poly.empty())
604 throw NotComputableException("HilbertQuasiPolynomial");
605 return quasi_poly;
606 }
607
getHilbertQuasiPolynomialDenom() const608 mpz_class HilbertSeries::getHilbertQuasiPolynomialDenom() const {
609 computeHilbertQuasiPolynomial();
610 if (quasi_poly.empty())
611 throw NotComputableException("HilbertQuasiPolynomial");
612 return quasi_denom;
613 }
614
computeHilbertQuasiPolynomial() const615 void HilbertSeries::computeHilbertQuasiPolynomial() const {
616 if (isHilbertQuasiPolynomialComputed() || nr_coeff_quasipol == 0)
617 return;
618 simplify();
619
620 vector<long> denom_vec = to_vector(denom);
621 if (nr_coeff_quasipol > (long)denom_vec.size()) {
622 if (verbose)
623 verboseOutput() << "Number of coeff of quasipol too large. Reset to deault value." << endl;
624 nr_coeff_quasipol = -1;
625 }
626
627 if (period_bounded && period > PERIOD_BOUND) {
628 if (verbose) {
629 errorOutput() << "WARNING: We skip the computation of the Hilbert-quasi-polynomial because the period " << period
630 << " is too big!" << endl;
631 errorOutput() << "Rerun with NO_PERIOD_BOUND" << endl;
632 }
633 return;
634 }
635 if (verbose && period > 1) {
636 verboseOutput() << "Computing Hilbert quasipolynomial of period " << period << " ..." << flush;
637 }
638 long i, j;
639 // period und dim encode the denominator
640 // now adjust the numerator
641 long num_size = num.size();
642 vector<mpz_class> norm_num(num_size); // normalized numerator
643 for (i = 0; i < num_size; ++i) {
644 norm_num[i] = num[i];
645 }
646 map<long, denom_t>::reverse_iterator rit;
647 long d;
648 vector<mpz_class> r;
649 for (rit = denom.rbegin(); rit != denom.rend(); ++rit) {
650 INTERRUPT_COMPUTATION_BY_EXCEPTION
651
652 d = rit->first;
653 // nothing to do if it already has the correct t-power
654 if (d != period) {
655 // norm_num *= (1-t^p / 1-t^d)^denom[d]
656 // first by multiply: norm_num *= (1-t^p)^denom[d]
657 poly_mult_to(norm_num, period, rit->second);
658
659 // then divide: norm_num /= (1-t^d)^denom[d]
660 for (i = 0; i < rit->second; ++i) {
661 poly_div(norm_num, r, norm_num, coeff_vector<mpz_class>(d));
662 assert(r.empty()); // assert remainder r is 0
663 }
664 }
665 }
666
667 // determine the common period of the coefficients that will be computed and printed
668 long reduced_period;
669 if (nr_coeff_quasipol >= 0) {
670 reduced_period = 1;
671 for (long j = 0; j < nr_coeff_quasipol; ++j)
672 reduced_period = lcm(reduced_period, denom_vec[j]);
673 }
674 else
675 reduced_period = period;
676
677 // cut numerator into period many pieces and apply standard method
678 // we make only reduced_period many components
679 quasi_poly = vector<vector<mpz_class> >(reduced_period);
680 long nn_size = norm_num.size();
681 for (j = 0; j < reduced_period; ++j) {
682 quasi_poly[j].reserve(dim);
683 }
684 for (i = 0; i < nn_size; ++i) {
685 if (i % period < reduced_period)
686 quasi_poly[i % period].push_back(norm_num[i]);
687 }
688
689 #pragma omp parallel for
690 for (j = 0; j < reduced_period; ++j) {
691 INTERRUPT_COMPUTATION_BY_EXCEPTION
692
693 quasi_poly[j] = compute_polynomial(quasi_poly[j], dim);
694 }
695
696 // substitute t by t/period:
697 // dividing by period^dim and multipling the coeff with powers of period
698 mpz_class pp = 1;
699 for (i = dim - 2; i >= 0; --i) {
700 pp *= period; // p^i ok, it is p^(dim-1-i)
701 for (j = 0; j < reduced_period; ++j) {
702 quasi_poly[j][i] *= pp;
703 }
704 } // at the end pp=p^dim-1
705 // the common denominator for all coefficients, dim! * pp
706 quasi_denom = permutations<mpz_class>(1, dim) * pp;
707 // substitute t by t-j
708 for (j = 0; j < reduced_period; ++j) {
709 // X |--> X - (j + shift)
710 linear_substitution<mpz_class>(quasi_poly[j], j + shift); // replaces quasi_poly[j]
711 }
712 // divide by gcd //TODO operate directly on vector
713 Matrix<mpz_class> QP(quasi_poly);
714 mpz_class g = QP.matrix_gcd();
715 g = libnormaliz::gcd(g, quasi_denom);
716 quasi_denom /= g;
717 QP.scalar_division(g);
718 // we use a normed shift, so that the cylcic shift % period always yields a non-negative integer
719 long normed_shift = -shift;
720 while (normed_shift < 0)
721 normed_shift += reduced_period;
722 for (j = 0; j < reduced_period; ++j) {
723 quasi_poly[j] = QP[(j + normed_shift) % reduced_period]; // QP[ (j - shift) % p ]
724 }
725
726 long delete_coeff = 0;
727 if (nr_coeff_quasipol >= 0)
728 delete_coeff = (long)quasi_poly[0].size() - nr_coeff_quasipol;
729
730 for (auto& i : quasi_poly) // delete coefficients that have not been computed completely
731 for (long j = 0; j < delete_coeff; ++j)
732 i[j] = 0;
733
734 if (verbose && period > 1) {
735 verboseOutput() << " done." << endl;
736 }
737 }
738
739 // expands the series to degree to_degree
compute_expansion() const740 void HilbertSeries::compute_expansion() const {
741 expansion.clear();
742 vector<mpz_class> denom_expansion = expand_denom();
743 expansion = poly_mult(num, denom_expansion);
744 if ((long)expansion.size() > expansion_degree + 1)
745 expansion.resize(expansion_degree + 1);
746 }
747
getExpansion() const748 vector<mpz_class> HilbertSeries::getExpansion() const {
749 compute_expansion();
750 return expansion;
751 }
752
get_expansion_degree() const753 long HilbertSeries::get_expansion_degree() const {
754 return expansion_degree;
755 }
756
set_expansion_degree(long degree)757 void HilbertSeries::set_expansion_degree(long degree) {
758 expansion_degree = degree;
759 }
760
expand_denom() const761 vector<mpz_class> HilbertSeries::expand_denom() const {
762 vector<long> denom_vec = to_vector(denom);
763 vector<mpz_class> result(1, 1); // the constant 1
764 for (long i : denom_vec) {
765 vector<mpz_class> this_factor = expand_inverse(i, expansion_degree);
766 result = poly_mult(result, this_factor);
767 if ((long)result.size() > expansion_degree + 1)
768 result.resize(expansion_degree + 1);
769 }
770 return result;
771 }
772
773 // computes the series expansion of 1/(1-t^e)
expand_inverse(size_t exponent,long to_degree)774 vector<mpz_class> expand_inverse(size_t exponent, long to_degree) {
775 vector<mpz_class> expansion(to_degree + 1, 0);
776 for (long i = 0; i <= to_degree; i += exponent)
777 expansion[i] = 1;
778 return expansion;
779 }
780
781 // returns the numerator, repr. as vector of coefficients, the h-vector
getNum() const782 const vector<mpz_class>& HilbertSeries::getNum() const {
783 simplify();
784 return num;
785 }
786 // returns the denominator, repr. as a map of the exponents of (1-t^i)^e
getDenom() const787 const map<long, denom_t>& HilbertSeries::getDenom() const {
788 simplify();
789 return denom;
790 }
791
792 // returns the numerator, repr. as vector of coefficients
getCyclotomicNum() const793 const vector<mpz_class>& HilbertSeries::getCyclotomicNum() const {
794 simplify();
795 return cyclo_num;
796 }
797 // returns the denominator, repr. as a map of the exponents of (1-t^i)^e
getCyclotomicDenom() const798 const map<long, denom_t>& HilbertSeries::getCyclotomicDenom() const {
799 simplify();
800 return cyclo_denom;
801 }
802
getHSOPDenom() const803 const map<long, denom_t>& HilbertSeries::getHSOPDenom() const {
804 simplify();
805 return hsop_denom;
806 }
807
getHSOPNum() const808 const vector<mpz_class>& HilbertSeries::getHSOPNum() const {
809 simplify();
810 assert(v_non_negative(hsop_num));
811 return hsop_num;
812 }
813
814 // shift
setShift(long s)815 void HilbertSeries::setShift(long s) {
816 if (shift != s) {
817 is_simplified = false;
818 // remove quasi-poly //TODO could also be adjusted
819 quasi_poly.clear();
820 quasi_denom = 1;
821 shift = s;
822 }
823 }
824
setHSOPDenom(vector<denom_t> new_denom)825 void HilbertSeries::setHSOPDenom(vector<denom_t> new_denom) {
826 hsop_denom = count_in_map<long, denom_t>(new_denom);
827 }
828
setHSOPDenom(map<long,denom_t> new_denom)829 void HilbertSeries::setHSOPDenom(map<long, denom_t> new_denom) {
830 hsop_denom = new_denom;
831 }
832
getShift() const833 long HilbertSeries::getShift() const {
834 return shift;
835 }
836
adjustShift()837 void HilbertSeries::adjustShift() {
838 collectData();
839 size_t adj = 0; // adjust shift by
840 while (adj < num.size() && num[adj] == 0)
841 adj++;
842 if (adj > 0) {
843 shift += adj;
844 num.erase(num.begin(), num.begin() + adj);
845 if (cyclo_num.size() != 0) {
846 assert(cyclo_num.size() >= adj);
847 cyclo_num.erase(cyclo_num.begin(), cyclo_num.begin() + adj);
848 }
849 }
850 }
851
852 /*
853 // methods for textual transfer of a Hilbert Series
854 string HilbertSeries::to_string_rep() const {
855 collectData();
856 ostringstream s;
857
858 s << num.size() << " ";
859 s << num;
860 vector<denom_t> denom_vector(to_vector(denom));
861 s << denom_vector.size() << " ";
862 s << denom_vector;
863 return s.str();
864 }
865
866 void HilbertSeries::from_string_rep(const string& input) {
867 istringstream s(input);
868 long i, size;
869
870 s >> size;
871 num.resize(size);
872 for (i = 0; i < size; ++i) {
873 s >> num[i];
874 }
875
876 vector<denom_t> denom_vector;
877 s >> size;
878 denom_vector.resize(size);
879 for (i = 0; i < size; ++i) {
880 s >> denom_vector[i];
881 }
882
883 denom = count_in_map<long, denom_t>(denom_vector);
884 is_simplified = false;
885 }
886
887 // writes in a human readable format
888 ostream& operator<<(ostream& out, const HilbertSeries& HS) {
889 HS.collectData();
890 out << "(";
891 // i == 0
892 if (HS.num.size() > 0)
893 out << " " << HS.num[0];
894 if (HS.shift != 0)
895 out << "*t^" << HS.shift;
896 for (size_t i = 1; i < HS.num.size(); ++i) {
897 if (HS.num[i] == 1)
898 out << " +t^" << i + HS.shift;
899 else if (HS.num[i] == -1)
900 out << " -t^" << i + HS.shift;
901 else if (HS.num[i] > 0)
902 out << " +" << HS.num[i] << "*t^" << i + HS.shift;
903 else if (HS.num[i] < 0)
904 out << " -" << -HS.num[i] << "*t^" << i + HS.shift;
905 }
906 out << " ) / (";
907 if (HS.denom.empty()) {
908 out << " 1";
909 }
910 for (const auto& it : HS.denom) {
911 if (it.second != 0)
912 out << " (1-t^" << it.first << ")^" << it.second;
913 }
914 out << " )" << std::endl;
915 return out;
916 }
917 */
918
919 //---------------------------------------------------------------------------
920 // polynomial operations, for polynomials repr. as vector of coefficients
921 //---------------------------------------------------------------------------
922
923 // returns the coefficient vector of 1-t^i
924 template <typename Integer>
coeff_vector(size_t i)925 vector<Integer> coeff_vector(size_t i) {
926 vector<Integer> p(i + 1, 0);
927 p[0] = 1;
928 p[i] = -1;
929 return p;
930 }
931
932 template <typename Integer>
remove_zeros(vector<Integer> & a)933 void remove_zeros(vector<Integer>& a) {
934 size_t i = a.size();
935 while (i > 0 && a[i - 1] == 0)
936 --i;
937
938 if (i < a.size()) {
939 a.resize(i);
940 }
941 }
942
943 // a += b (also possible to define the += op for vector)
944 template <typename Integer>
poly_add_to(vector<Integer> & a,const vector<Integer> & b)945 void poly_add_to(vector<Integer>& a, const vector<Integer>& b) {
946 size_t b_size = b.size();
947 if (a.size() < b_size) {
948 a.resize(b_size);
949 }
950 for (size_t i = 0; i < b_size; ++i) {
951 a[i] += b[i];
952 }
953 remove_zeros(a);
954 }
955
956 // a += b*t^m
957 template <typename Integer>
poly_add_to_tm(vector<Integer> & a,const vector<Integer> & b,long m)958 void poly_add_to_tm(vector<Integer>& a, const vector<Integer>& b, long m) {
959 size_t b_size = b.size();
960 size_t b_m = b_size + m;
961 if (a.size() < b_m) {
962 a.resize(b_m);
963 }
964 for (size_t i = 0; i < b_size; ++i) {
965 a[i + m] += b[i];
966 }
967 remove_zeros(a);
968 }
969
970 // a -= b (also possible to define the -= op for vector)
971 template <typename Integer>
poly_sub_to(vector<Integer> & a,const vector<Integer> & b)972 void poly_sub_to(vector<Integer>& a, const vector<Integer>& b) {
973 size_t b_size = b.size();
974 if (a.size() < b_size) {
975 a.resize(b_size);
976 }
977 for (size_t i = 0; i < b_size; ++i) {
978 a[i] -= b[i];
979 }
980 remove_zeros(a);
981 }
982
983 // a *= t^m
984 template <typename Integer>
poly_mult_by_tm(vector<Integer> & a,long m)985 void poly_mult_by_tm(vector<Integer>& a, long m) {
986 long a_ori_size = a.size();
987 a.resize(a_ori_size + m);
988 for (long i = a_ori_size - 1; i >= 0; --i)
989 a[i + m] = a[i];
990 for (long i = 0; i < m; ++i)
991 a[i] = 0;
992 }
993
994 // a * b
995
996 /* template<typename Integer>
997 vector<Integer> old_poly_mult(const vector<Integer>& a, const vector<Integer>& b) {
998 size_t a_size = a.size();
999 size_t b_size = b.size();
1000
1001
1002 vector<Integer> p( a_size + b_size - 1 );
1003 size_t i,j;
1004 for (i=0; i<a_size; ++i) {
1005 if (a[i] == 0) continue;
1006 for (j=0; j<b_size; ++j) {
1007 if (b[j] == 0) continue;
1008 p[i+j] += a[i]*b[j];
1009 }
1010 }
1011 return p;
1012 }*/
1013
1014 template <typename Integer>
karatsubamult(const vector<Integer> & a,const vector<Integer> & b)1015 vector<Integer> karatsubamult(const vector<Integer>& a, const vector<Integer>& b) {
1016 size_t a_size = a.size();
1017 size_t b_size = b.size();
1018 if (a_size * b_size <= 1000 || a_size <= 10 || b_size <= 10) {
1019 return poly_mult(a, b);
1020 }
1021
1022 size_t m = (a_size + 1) / 2;
1023 if (2 * m < (b_size + 1)) {
1024 m = (b_size + 1) / 2;
1025 }
1026
1027 vector<Integer> f0(m), f1(m), g0(m), g1(m);
1028 for (size_t i = 0; i < m && i < a_size; ++i)
1029 f0[i] = a[i];
1030 for (size_t i = m; i < a_size; ++i)
1031 f1[i - m] = a[i];
1032 for (size_t i = 0; i < m && i < b_size; ++i)
1033 g0[i] = b[i];
1034 for (size_t i = m; i < b_size; ++i)
1035 g1[i - m] = b[i];
1036 remove_zeros(f0);
1037 remove_zeros(f1);
1038 remove_zeros(g0);
1039 remove_zeros(g1);
1040
1041 vector<Integer> sf = f0;
1042 vector<Integer> sg = g0;
1043
1044 vector<Integer> mix;
1045 vector<Integer> h00;
1046 vector<Integer> h11;
1047
1048 #pragma omp parallel // num_threads(3)
1049 {
1050 #pragma omp single nowait
1051 {
1052 h00 = karatsubamult(f0, g0); // h00 = f0 * g0
1053 }
1054
1055 #pragma omp single nowait
1056 {
1057 h11 = karatsubamult(f1, g1); // h11 = f1 * g1
1058 }
1059
1060 #pragma omp single nowait
1061 {
1062 poly_add_to(sf, f1); // f0+f1
1063 poly_add_to(sg, g1); // g0 + g1
1064 mix = karatsubamult(sf, sg); // (f0 + f1)*(g0 + g1)
1065 }
1066
1067 } // parallel
1068
1069 f0.clear();
1070 g0.clear();
1071 f1.clear();
1072 g1.clear();
1073
1074 poly_sub_to(mix, h00); // mix = mix - f0*g0
1075 poly_sub_to(mix, h11); // mix = mix - f1*g1
1076
1077 poly_add_to_tm(h00, mix, m);
1078 poly_add_to_tm(h00, h11, 2 * m);
1079
1080 return h00;
1081 }
1082
1083 // a *= (1-t^d)^e
1084 template <typename Integer>
poly_mult_to(vector<Integer> & a,long d,long e)1085 void poly_mult_to(vector<Integer>& a, long d, long e) {
1086 assert(d > 0);
1087 assert(e >= 0);
1088 long i;
1089 a.reserve(a.size() + d * e);
1090 while (e > 0) {
1091 a.resize(a.size() + d);
1092 for (i = a.size() - 1; i >= d; --i) {
1093 a[i] -= a[i - d];
1094 }
1095 e--;
1096 }
1097 }
1098
1099 //---------------------------------------------------------------------------
1100 // computing the Hilbert polynomial from h-vector
1101 //---------------------------------------------------------------------------
1102
1103 // The algorithm follows "Cohen-Macaulay rings", 4.1.5 and 4.1.9.
1104 // The E_vector is the vector of higher multiplicities.
1105 // It is assumed that (d-1)! is used as a common denominator in the calling routine.
1106
1107 template <typename Integer>
compute_e_vector(vector<Integer> Q,int dim)1108 vector<Integer> compute_e_vector(vector<Integer> Q, int dim) {
1109 size_t j;
1110 int i;
1111 vector<Integer> E_Vector(dim, 0);
1112 // cout << "QQQ " << Q;
1113 // Q.resize(dim+1);
1114 int bound = Q.size();
1115 if (bound > dim)
1116 bound = dim;
1117 for (i = 0; i < bound; i++) {
1118 for (j = 0; j < Q.size() - i; j++) {
1119 E_Vector[i] += Q[j];
1120 }
1121 E_Vector[i] /= permutations<Integer>(1, i);
1122 for (j = 1; j < Q.size() - i; j++) {
1123 Q[j - 1] = static_cast<unsigned long>(j) * Q[j];
1124 }
1125 }
1126 return E_Vector;
1127 }
1128
1129 //---------------------------------------------------------------------------
1130
1131 template <typename Integer>
compute_polynomial(vector<Integer> h_vector,int dim)1132 vector<Integer> compute_polynomial(vector<Integer> h_vector, int dim) {
1133 // handle dimension 0
1134 if (dim == 0)
1135 return vector<Integer>(dim);
1136
1137 vector<Integer> Hilbert_Polynomial = vector<Integer>(dim);
1138 int i, j;
1139
1140 Integer mult_factor;
1141 vector<Integer> E_Vector = compute_e_vector(h_vector, dim);
1142 vector<Integer> C(dim, 0);
1143 C[0] = 1;
1144 for (i = 0; i < dim; i++) {
1145 mult_factor = permutations<Integer>(i, dim);
1146 if (((dim - 1 - i) % 2) == 0) {
1147 for (j = 0; j < dim; j++) {
1148 Hilbert_Polynomial[j] += mult_factor * E_Vector[dim - 1 - i] * C[j];
1149 }
1150 }
1151 else {
1152 for (j = 0; j < dim; j++) {
1153 Hilbert_Polynomial[j] -= mult_factor * E_Vector[dim - 1 - i] * C[j];
1154 }
1155 }
1156 for (j = dim - 1; 0 < j; j--) {
1157 C[j] = (unsigned long)(i + 1) * C[j] + C[j - 1];
1158 }
1159 C[0] = permutations<Integer>(1, i + 1);
1160 }
1161
1162 return Hilbert_Polynomial;
1163 }
1164
1165 //---------------------------------------------------------------------------
1166
1167 // substitutes t by (t-a), overwrites the polynomial!
1168 template <typename Integer>
linear_substitution(vector<Integer> & poly,const Integer & a)1169 void linear_substitution(vector<Integer>& poly, const Integer& a) {
1170 long deg = poly.size() - 1;
1171 // Iterated division by (t+a)
1172 for (long step = 0; step < deg; ++step) {
1173 for (long i = deg - 1; i >= step; --i) {
1174 poly[i] -= a * poly[i + 1];
1175 }
1176 // the remainders are the coefficients of the transformed polynomial
1177 }
1178 }
1179
1180 //---------------------------------------------------------------------------
IntegrationData()1181 IntegrationData::IntegrationData() {
1182 }
1183
set_nr_coeff_quasipol(long nr_coeff)1184 void IntegrationData::set_nr_coeff_quasipol(long nr_coeff) {
1185 weighted_Ehrhart_series.first.set_nr_coeff_quasipol(nr_coeff);
1186 }
1187
set_expansion_degree(long degree)1188 void IntegrationData::set_expansion_degree(long degree) {
1189 weighted_Ehrhart_series.first.set_expansion_degree(degree);
1190 }
1191
getPolynomial() const1192 string IntegrationData::getPolynomial() const {
1193 return polynomial;
1194 }
1195
getDegreeOfPolynomial() const1196 long IntegrationData::getDegreeOfPolynomial() const {
1197 return degree_of_polynomial;
1198 }
1199
setDegreeOfPolynomial(const long d)1200 void IntegrationData::setDegreeOfPolynomial(const long d) {
1201 degree_of_polynomial = d;
1202 }
1203
IntegrationData(const string & poly)1204 IntegrationData::IntegrationData(const string& poly) {
1205 polynomial = poly;
1206 polynomial_is_homogeneous = false; // to be on the safe side
1207 }
1208
isWeightedEhrhartQuasiPolynomialComputed() const1209 bool IntegrationData::isWeightedEhrhartQuasiPolynomialComputed() const {
1210 return weighted_Ehrhart_series.first.isHilbertQuasiPolynomialComputed();
1211 }
1212
getWeightedEhrhartQuasiPolynomial() const1213 const vector<vector<mpz_class> >& IntegrationData::getWeightedEhrhartQuasiPolynomial() const {
1214 return weighted_Ehrhart_series.first.getHilbertQuasiPolynomial();
1215 }
1216
resetHilbertQuasiPolynomial()1217 void IntegrationData::resetHilbertQuasiPolynomial() {
1218 weighted_Ehrhart_series.first.resetHilbertQuasiPolynomial();
1219 }
1220
getExpansion() const1221 vector<mpz_class> IntegrationData::getExpansion() const {
1222 return weighted_Ehrhart_series.first.getExpansion();
1223 }
1224
computeWeightedEhrhartQuasiPolynomial()1225 void IntegrationData::computeWeightedEhrhartQuasiPolynomial() {
1226 weighted_Ehrhart_series.first.computeHilbertQuasiPolynomial();
1227 }
1228
getWeightedEhrhartQuasiPolynomialDenom() const1229 mpz_class IntegrationData::getWeightedEhrhartQuasiPolynomialDenom() const {
1230 return weighted_Ehrhart_series.first.getHilbertQuasiPolynomialDenom() * weighted_Ehrhart_series.second;
1231 }
1232
1233 // the following 4 functions are nit used in Normaliz, bur provided for interfaces
getNum_ZZ() const1234 const vector<mpz_class>& IntegrationData::getNum_ZZ() const {
1235 return weighted_Ehrhart_series.first.getNum();
1236 }
1237
getDenom() const1238 const map<long, denom_t>& IntegrationData::getDenom() const {
1239 return weighted_Ehrhart_series.first.getDenom();
1240 }
1241
getCyclotomicNum_ZZ() const1242 const vector<mpz_class>& IntegrationData::getCyclotomicNum_ZZ() const {
1243 return weighted_Ehrhart_series.first.getCyclotomicNum();
1244 }
1245
getCyclotomicDenom() const1246 const map<long, denom_t>& IntegrationData::getCyclotomicDenom() const {
1247 return weighted_Ehrhart_series.first.getCyclotomicDenom();
1248 }
1249
getWeightedEhrhartSeries() const1250 const pair<HilbertSeries, mpz_class>& IntegrationData::getWeightedEhrhartSeries() const {
1251 return weighted_Ehrhart_series;
1252 }
1253
getIntegral() const1254 mpq_class IntegrationData::getIntegral() const {
1255 return integral;
1256 }
1257
getEuclideanIntegral() const1258 nmz_float IntegrationData::getEuclideanIntegral() const {
1259 return euclidean_integral;
1260 }
1261
getNumeratorCommonDenom() const1262 mpz_class IntegrationData::getNumeratorCommonDenom() const {
1263 return weighted_Ehrhart_series.second;
1264 }
1265
getVirtualMultiplicity() const1266 mpq_class IntegrationData::getVirtualMultiplicity() const {
1267 return virtual_multiplicity;
1268 }
1269
setIntegral(const mpq_class I)1270 void IntegrationData::setIntegral(const mpq_class I) {
1271 integral = I;
1272 }
1273
setEuclideanIntegral(const nmz_float I)1274 void IntegrationData::setEuclideanIntegral(const nmz_float I) {
1275 euclidean_integral = I;
1276 }
1277
setVirtualMultiplicity(const mpq_class I)1278 void IntegrationData::setVirtualMultiplicity(const mpq_class I) {
1279 virtual_multiplicity = I;
1280 }
1281
setWeightedEhrhartSeries(const pair<HilbertSeries,mpz_class> & E)1282 void IntegrationData::setWeightedEhrhartSeries(const pair<HilbertSeries, mpz_class>& E) {
1283 weighted_Ehrhart_series = E;
1284 weighted_Ehrhart_series.first.adjustShift();
1285 }
1286
setHomogeneity(const bool hom)1287 void IntegrationData::setHomogeneity(const bool hom) {
1288 polynomial_is_homogeneous = hom;
1289 }
1290
isPolynomialHomogeneous() const1291 bool IntegrationData::isPolynomialHomogeneous() const {
1292 return polynomial_is_homogeneous;
1293 }
1294
1295 } // end namespace libnormaliz
1296
1297 #ifdef NMZ_MIC_OFFLOAD
1298 #pragma offload_attribute(pop)
1299 #endif
1300