GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it

1
/*
2
 * Normaliz
3
 * Copyright (C) 2007-2014  Winfried Bruns, Bogdan Ichim, Christof Soeger
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/HilbertSeries.h"
35
#include "libnormaliz/vector_operations.h"
36
#include "libnormaliz/map_operations.h"
37
#include "libnormaliz/integer.h"
38
#include "libnormaliz/convert.h"
39
#include "libnormaliz/my_omp.h"
40

41
#include "libnormaliz/matrix.h"
42

43
//---------------------------------------------------------------------------
44

45
namespace libnormaliz {
46
using std::cout; using std::endl; using std::flush;
47
using std::istringstream; using std::ostringstream;
48
using std::pair;
49

50
long lcm_of_keys(const map<long, denom_t>& m){
51
    long l = 1;
52
    map<long, denom_t>::const_iterator it;
53
    for (it = m.begin(); it != m.end(); ++it) {
54
        if (it->second != 0)
55
            l = lcm(l,it->first);
56
    }
57
    return l;
58
}
59

60
// compute the hsop numerator by multiplying the HS with a denominator
61
// of the form product of (1-t^i)
62
void HilbertSeries::compute_hsop_num() const{
63
        // get the denominator as a polynomial by mutliplying the (1-t^i) terms
64
        vector<mpz_class> hsop_denom_poly=vector<mpz_class>(1,1);
65
        map<long,denom_t>::iterator it;
66
        long factor;
67
        for (it=hsop_denom.begin();it!=hsop_denom.end();++it){
68
            factor = it->first;
69
            denom_t& denom_i = it->second;
70
            poly_mult_to(hsop_denom_poly,factor,denom_i);
71
        }
72
        //cout << "new denominator as polynomial: " << hsop_denom_poly << endl;
73
        vector<mpz_class>  quot,remainder,cyclo_poly;
74
        //first divide the new denom by the cyclo polynomials
75
        for (auto it=cyclo_denom.begin();it!=cyclo_denom.end();++it){
76
            for(long i=0;i<it->second;i++){
77
                cyclo_poly = cyclotomicPoly<mpz_class>(it->first);
78
                //cout << "the cyclotomic polynomial is " << cyclo_poly << endl;
79
                // TODO: easier polynomial division possible?
80
                poly_div(quot,remainder,hsop_denom_poly,cyclo_poly);
81
                //cout << "the quotient is " << quot << endl;
82
                hsop_denom_poly=quot;
83
                assert(remainder.size()==0);
84
            }
85
        }
86
        // multiply with the old numerator
87
        hsop_num = poly_mult(hsop_denom_poly,cyclo_num);
88
}
89

90
//---------------------------------------------------------------------------
91

92
void HilbertSeries::initialize(){
93
    
94
    is_simplified = false;
95
    shift = 0;
96
    verbose = false;
97
    nr_coeff_quasipol=-1; // all coefficients
98
    period_bounded=true;
99
}
100

101
// Constructor, creates 0/1
102
HilbertSeries::HilbertSeries() {
103
    num = vector<mpz_class>(1,0);
104
    //denom just default constructed
105
    initialize();
106
}
107

108
// Constructor, creates num/denom, see class description for format
109
HilbertSeries::HilbertSeries(const vector<num_t>& numerator, const vector<denom_t>& gen_degrees) {
110
    num = vector<mpz_class>(1,0);
111
    add(numerator, gen_degrees);
112
    initialize();
113
}
114

115
// Constructor, creates num/denom, see class description for format
116
HilbertSeries::HilbertSeries(const vector<mpz_class>& numerator, const map<long, denom_t>& denominator) {
117
    num = numerator;
118
    denom = denominator;
119
    initialize();
120
}
121

122
// Constructor, string as created by to_string_rep
123
HilbertSeries::HilbertSeries(const string& str) {
124
    from_string_rep(str);
125
    initialize();
126
}
127

128

129
void HilbertSeries::reset() {
130
    num.clear();
131
    num.push_back(0);
132
    denom.clear();
133
    denom_classes.clear();
134
    shift = 0;
135
    is_simplified = false;
136
}
137

138
void HilbertSeries::set_nr_coeff_quasipol(long nr_coeff){
139
    nr_coeff_quasipol=nr_coeff;
140
}
141

142
long HilbertSeries::get_nr_coeff_quasipol() const{
143
    return nr_coeff_quasipol;
144
}
145

146
void HilbertSeries::set_period_bounded(bool on_off) const{ //period_bounded is mutable
147
    period_bounded=on_off;
148
}
149

150
bool HilbertSeries::get_period_bounded() const{
151
    return period_bounded;
152
}
153
// add another HilbertSeries to this
154
void HilbertSeries::add(const vector<num_t>& num, const vector<denom_t>& gen_degrees) {
155
    vector<denom_t> sorted_gd(gen_degrees);
156
    sort(sorted_gd.begin(), sorted_gd.end());
157
    if (gen_degrees.size() > 0) 
158
        assert(sorted_gd[0]>0); //TODO InputException?
159
    poly_add_to(denom_classes[sorted_gd], num);
160
    if (denom_classes.size() > DENOM_CLASSES_BOUND)
161
        collectData();
162
    is_simplified = false;
163
}
164

165

166
// add another HilbertSeries to this
167
HilbertSeries& HilbertSeries::operator+=(const HilbertSeries& other) {
168
    // add denom_classes
169
    map< vector<denom_t>, vector<num_t> >::const_iterator it;
170
    for (it = other.denom_classes.begin(); it != other.denom_classes.end(); ++it) {
171
        poly_add_to(denom_classes[it->first], it->second);
172
    }
173
    // add accumulated data
174
    vector<mpz_class> num_copy(other.num);
175
    performAdd(num_copy, other.denom);
176
    return (*this);
177
}
178

179
void HilbertSeries::performAdd(const vector<num_t>& numerator, const vector<denom_t>& gen_degrees) const {
180
    map<long, denom_t> other_denom;
181
    size_t i, s = gen_degrees.size();
182
    for (i=0; i<s; ++i) {
183
        assert(gen_degrees[i]>0);
184
        other_denom[gen_degrees[i]]++;
185
    }
186
    // convert numerator to mpz
187
    vector<mpz_class> other_num(numerator.size());
188
    convert(other_num, numerator);
189
    performAdd(other_num, other_denom);
190
}
191

192
//modifies other_num!!
193
void HilbertSeries::performAdd(vector<mpz_class>& other_num, const map<long, denom_t>& oth_denom) const {
194
    map<long, denom_t> other_denom(oth_denom);  //TODO redesign, dont change other_denom
195
    // adjust denominators
196
    denom_t diff;
197
    map<long, denom_t>::iterator it;
198
    for (it = denom.begin(); it != denom.end(); ++it) {  // augment other
199
        denom_t& ref = other_denom[it->first];
200
        diff = it->second - ref;
201
        if (diff > 0) {
202
            ref += diff;
203
            poly_mult_to(other_num, it->first, diff);
204
        }
205
    }
206
    for (it = other_denom.begin(); it != other_denom.end(); ++it) {  // augment this
207
        denom_t& ref = denom[it->first];
208
        diff = it->second - ref;
209
        if (diff > 0) {
210
            ref += diff;
211
            poly_mult_to(num, it->first, diff);
212
        }
213
    }
214
    assert (denom == other_denom);
215

216
    // now just add the numerators
217
    poly_add_to(num,other_num);
218
    remove_zeros(num);
219
    is_simplified = false;
220
}
221

222
void HilbertSeries::collectData() const {
223
    if (denom_classes.empty()) return;
224
	if (verbose) verboseOutput() << "Adding " << denom_classes.size() << " denominator classes..." << flush;
225
    map< vector<denom_t>, vector<num_t> >::iterator it;
226
    for (it = denom_classes.begin(); it != denom_classes.end(); ++it) {
227
        
228
        INTERRUPT_COMPUTATION_BY_EXCEPTION
229
        
230
        performAdd(it->second, it->first);
231
    }
232
    denom_classes.clear();
233
	if (verbose) verboseOutput() << " done." << endl;
234
}
235

236
// simplify, see class description
237
void HilbertSeries::simplify() const {
238
    if (is_simplified)
239
        return;
240
    collectData();
241
/*    if (verbose) {
242
        verboseOutput() << "Hilbert series before simplification: "<< endl << *this;
243
    }*/
244
    vector<mpz_class> q, r, poly; //polynomials
245
    // In denom_cyclo we collect cyclotomic polynomials in the denominator.
246
    // During this method the Hilbert series is given by num/(denom*cdenom)
247
    // where denom | cdenom are exponent vectors of (1-t^i) | i-th cyclotminc poly.
248
    map<long, denom_t> cdenom;
249

250
    map<long, denom_t> save_denom=denom;
251
    map<long, denom_t>::reverse_iterator rit;
252
    long i;
253
    for (rit = denom.rbegin(); rit != denom.rend(); ++rit) {
254
        // check if we can divide the numerator by (1-t^i)
255
        i = rit->first;
256
        denom_t& denom_i = rit->second;
257
        poly = coeff_vector<mpz_class>(i);
258
        while (denom_i > 0) {
259
            poly_div(q, r, num, poly);
260
            if (r.size() == 0) { // numerator is divisable by poly
261
                num = q;
262
                denom_i--;
263
            }
264
            else {
265
                break;
266
            }
267
        }
268
        if (denom_i == 0)
269
            continue;
270

271
        // decompose (1-t^i) into cyclotomic polynomial
272
        for(long d=1; d<=i/2; ++d) {
273
            if (i % d == 0)
274
                cdenom[d] += denom_i;
275
        }
276
        cdenom[i] += denom_i;
277
        // the product of the cyclo. is t^i-1 = -(1-t^i)
278
        if (denom_i%2 == 1)
279
            v_scalar_multiplication(num,mpz_class(-1));
280
    } // end for
281
    denom.clear();
282
 
283
    map<long, denom_t>::iterator it = cdenom.begin(); 
284
    while (it != cdenom.end()) {
285
        // check if we can divide the numerator by i-th cyclotomic polynomial
286
        
287
        INTERRUPT_COMPUTATION_BY_EXCEPTION
288
            
289
        i = it->first;
290
        denom_t& cyclo_i = it->second;
291
        poly = cyclotomicPoly<mpz_class>(i);
292
        while (cyclo_i > 0) {
293
            poly_div(q, r, num, poly);
294
            if (r.empty()) { // numerator is divisable by poly
295
                num = q;
296
                cyclo_i--;
297
            }
298
            else {
299
               break;
300
            }
301
        }
302

303
        if (cyclo_i == 0) {
304
            cdenom.erase(it++);
305
        } else {
306
            ++it;
307
        }
308
    }
309
    // done with canceling
310
    // save this representation
311
    cyclo_num = num;
312
    cyclo_denom = cdenom;
313

314
    // now collect the cyclotomic polynomials in (1-t^i) factors
315
    it = cdenom.find(1);
316
    if (it != cdenom.end())
317
        dim = it->second;
318
    else
319
        dim = 0;
320
    period = lcm_of_keys(cdenom);
321
    if (period_bounded && period > 10*PERIOD_BOUND) {
322
        if (verbose) {
323
            errorOutput() << "WARNING: Period is too big, the representation of the Hilbert series may have more than dimension many factors in the denominator!" << endl;
324
        }
325
        denom=save_denom;
326
    }
327
    else{
328
        while(true){
329
            //create a (1-t^k) factor in the denominator out of all cyclotomic poly.
330
            
331
            INTERRUPT_COMPUTATION_BY_EXCEPTION
332
            
333
            long k=1;                
334
            bool empty=true;
335
            vector<mpz_class> existing_factor(1,1); //collects the existing cyclotomic gactors in the denom
336
            for(it=cdenom.begin();it!=cdenom.end();++it){          // with multiplicvity 1
337
                if(it-> second>0){
338
                    empty=false;
339
                    k=libnormaliz::lcm(k,it->first);
340
                    existing_factor=poly_mult(existing_factor,cyclotomicPoly<mpz_class>(it->first));
341
                    it->second--;
342
                }     
343
            }
344
            if(empty)
345
                break;
346
            denom[k]++;
347
            vector<mpz_class> new_factor=coeff_vector<mpz_class>(k);
348
            vector<mpz_class> quotient, dummy;
349
            poly_div(quotient,dummy,new_factor,existing_factor);
350
            assert(dummy.empty()); //assert remainder r is 0
351
            num=poly_mult(num,quotient);
352
        }
353
    }
354

355
/*    if (verbose) {
356
        verboseOutput() << "Simplified Hilbert series: " << endl << *this;
357
    }*/
358
    if (!hsop_denom.empty()){
359
        compute_hsop_num();
360
    }
361
    is_simplified = true;
362
    computeDegreeAsRationalFunction();
363
    quasi_poly.clear();
364
}
365

366
void HilbertSeries::computeDegreeAsRationalFunction() const {
367
    simplify();
368
    long num_deg = num.size() - 1 + shift;
369
    long denom_deg = 0;
370
    for (auto it = denom.begin(); it != denom.end(); ++it) {
371
            denom_deg += it->first * it->second;
372
    }
373
    degree = num_deg - denom_deg;
374
}
375

376
long HilbertSeries::getDegreeAsRationalFunction() const {
377
    simplify();
378
    return degree;
379
}
380

381
long HilbertSeries::getPeriod() const {
382
    simplify();
383
    return period;
384
}
385

386
bool HilbertSeries::isHilbertQuasiPolynomialComputed() const {
387
    return is_simplified && !quasi_poly.empty();
388
}
389

390
vector< vector<mpz_class> > HilbertSeries::getHilbertQuasiPolynomial() const {
391
    computeHilbertQuasiPolynomial();
392
    if (quasi_poly.empty()) throw NotComputableException("HilbertQuasiPolynomial");
393
    return quasi_poly;
394
}
395

396
mpz_class HilbertSeries::getHilbertQuasiPolynomialDenom() const {
397
    computeHilbertQuasiPolynomial();
398
    if (quasi_poly.empty()) throw NotComputableException("HilbertQuasiPolynomial");
399
    return quasi_denom;
400
}
401

402
void HilbertSeries::computeHilbertQuasiPolynomial() const {
403

404
    if (isHilbertQuasiPolynomialComputed() || nr_coeff_quasipol==0) return;
405
    simplify();
406
    
407
    omp_set_nested(0);
408
    
409
    vector<long> denom_vec=to_vector(denom);    
410
    if(nr_coeff_quasipol > (long) denom_vec.size()){
411
        if(verbose)
412
            verboseOutput() << "Number of coeff of quasipol too large. Reset to deault value." << endl;
413
        nr_coeff_quasipol=-1;
414
    }
415
        
416
    if (period_bounded &&  period > PERIOD_BOUND) {
417
        if (verbose) {
418
            errorOutput()<<"WARNING: We skip the computation of the Hilbert-quasi-polynomial because the period "<< period <<" is too big!" <<endl;
419
        }
420
        return;
421
    }
422
    if (verbose && period > 1) {
423
        verboseOutput() << "Computing Hilbert quasipolynomial of period "
424
                        << period <<" ..." << flush;
425
    }
426
    long i,j;
427
    //period und dim encode the denominator
428
    //now adjust the numerator
429
    long num_size = num.size();
430
    vector<mpz_class> norm_num(num_size);  //normalized numerator
431
    for (i = 0;  i < num_size;  ++i) {
432
        norm_num[i] = num[i];
433
    }
434
    map<long, denom_t>::reverse_iterator rit;
435
    long d;
436
    vector<mpz_class> r;
437
    for (rit = denom.rbegin(); rit != denom.rend(); ++rit) {
438
        
439
        INTERRUPT_COMPUTATION_BY_EXCEPTION
440
            
441
        d = rit->first;
442
        //nothing to do if it already has the correct t-power
443
        if (d != period) {
444
            //norm_num *= (1-t^p / 1-t^d)^denom[d]
445
            //first by multiply: norm_num *= (1-t^p)^denom[d]
446
            poly_mult_to(norm_num, period, rit->second);
447

448
            //then divide: norm_num /= (1-t^d)^denom[d]
449
            for (i=0; i < rit->second; ++i) {
450
                poly_div(norm_num, r, norm_num, coeff_vector<mpz_class>(d));
451
                assert(r.empty()); //assert remainder r is 0
452
            }
453
        }
454
    }
455
    
456
    // determine the common period of the coefficients that will be computed and printed
457
    long reduced_period;
458
    if(nr_coeff_quasipol>=0){
459
        reduced_period=1;
460
        for(long j=0;j< nr_coeff_quasipol;++j)
461
            reduced_period=lcm(reduced_period,denom_vec[j]);
462
    }
463
    else
464
        reduced_period=period;
465
    
466
    //cut numerator into period many pieces and apply standard method
467
    // we make only reduced_period many components
468
    quasi_poly = vector< vector<mpz_class> >(reduced_period);
469
    long nn_size = norm_num.size();
470
    for (j=0; j<reduced_period; ++j) {
471
        quasi_poly[j].reserve(dim);
472
    }
473
    for (i=0; i<nn_size; ++i) {
474
        if(i%period<reduced_period)
475
            quasi_poly[i%period].push_back(norm_num[i]);
476
    }
477

478
    #pragma omp parallel for
479
    for (j=0; j<reduced_period; ++j) {
480
        
481
        INTERRUPT_COMPUTATION_BY_EXCEPTION
482
        
483
        quasi_poly[j] = compute_polynomial(quasi_poly[j], dim);
484
    }
485
    
486
    //substitute t by t/period:
487
    //dividing by period^dim and multipling the coeff with powers of period
488
    mpz_class pp=1;
489
    for (i = dim-2; i >= 0; --i) {
490
        pp *= period; //p^i   ok, it is p^(dim-1-i)
491
        for (j=0; j<reduced_period; ++j) {
492
            quasi_poly[j][i] *= pp;
493
        }
494
    } //at the end pp=p^dim-1
495
    //the common denominator for all coefficients, dim! * pp
496
    quasi_denom = permutations<mpz_class>(1,dim) * pp;
497
    //substitute t by t-j
498
    for (j=0; j<reduced_period; ++j) {
499
        // X |--> X - (j + shift)
500
        linear_substitution<mpz_class>(quasi_poly[j], j + shift); // replaces quasi_poly[j]
501
    }
502
    //divide by gcd //TODO operate directly on vector
503
    Matrix<mpz_class> QP(quasi_poly);
504
    mpz_class g = QP.matrix_gcd();
505
    g = libnormaliz::gcd(g,quasi_denom);
506
    quasi_denom /= g;
507
    QP.scalar_division(g);
508
    //we use a normed shift, so that the cylcic shift % period always yields a non-negative integer
509
    long normed_shift = -shift;
510
    while (normed_shift < 0) normed_shift += reduced_period;
511
    for (j=0; j<reduced_period; ++j) {
512
        quasi_poly[j] = QP[(j+normed_shift)%reduced_period]; // QP[ (j - shift) % p ]
513
    }
514
    
515
    long delete_coeff=0;
516
    if(nr_coeff_quasipol>=0)
517
        delete_coeff=(long) quasi_poly[0].size()-nr_coeff_quasipol;
518
    
519
    for(size_t i=0;i<quasi_poly.size();++i) // delete coefficients that have not been computed completely
520
        for(long j=0;j <delete_coeff;++j)
521
            quasi_poly[i][j]=0;
522
        
523
    if (verbose && period > 1) {
524
        verboseOutput() << " done." << endl;
525
    }
526
}
527

528
// returns the numerator, repr. as vector of coefficients, the h-vector
529
const vector<mpz_class>& HilbertSeries::getNum() const {
530
    simplify();
531
    return num;
532
}
533
// returns the denominator, repr. as a map of the exponents of (1-t^i)^e
534
const map<long, denom_t>& HilbertSeries::getDenom() const {
535
    simplify();
536
    return denom;
537
}
538

539
// returns the numerator, repr. as vector of coefficients
540
const vector<mpz_class>& HilbertSeries::getCyclotomicNum() const {
541
    simplify();
542
    return cyclo_num;
543
}
544
// returns the denominator, repr. as a map of the exponents of (1-t^i)^e
545
const map<long, denom_t>& HilbertSeries::getCyclotomicDenom() const {
546
    simplify();
547
    return cyclo_denom;
548
}
549

550
const map<long, denom_t>& HilbertSeries::getHSOPDenom() const {
551
    simplify();
552
    return hsop_denom;
553
}
554

555
const vector<mpz_class>& HilbertSeries::getHSOPNum() const {
556
    simplify();
557
    assert(v_is_nonnegative(hsop_num));
558
    return hsop_num;
559
}
560

561
// shift
562
void HilbertSeries::setShift(long s) {
563
    if (shift != s) {
564
        is_simplified = false;
565
        // remove quasi-poly //TODO could also be adjusted
566
        quasi_poly.clear();
567
        quasi_denom = 1;
568
        shift = s;
569
    }
570
}
571

572
void HilbertSeries::setHSOPDenom(vector<denom_t> new_denom){
573
    hsop_denom=count_in_map<long,denom_t>(new_denom);
574
}
575

576
void HilbertSeries::setHSOPDenom(map<long,denom_t> new_denom){
577
    hsop_denom=new_denom;
578
}
579

580
long HilbertSeries::getShift() const {
581
    return shift;
582
}
583

584
void HilbertSeries::adjustShift() {
585
    collectData();
586
    size_t adj = 0; // adjust shift by
587
    while (adj < num.size() && num[adj] == 0) adj++;
588
    if (adj > 0) {
589
        shift += adj;
590
        num.erase(num.begin(),num.begin()+adj);
591
        if (cyclo_num.size() != 0) {
592
            assert (cyclo_num.size() >= adj);
593
            cyclo_num.erase(cyclo_num.begin(),cyclo_num.begin()+adj);
594
        }
595
    }
596
}
597

598
// methods for textual transfer of a Hilbert Series
599
string HilbertSeries::to_string_rep() const {
600

601
    collectData();
602
    ostringstream s;
603

604
    s << num.size() << " ";
605
    s << num;
606
    vector<denom_t> denom_vector(to_vector(denom));
607
    s << denom_vector.size() << " ";
608
    s << denom_vector;
609
    return s.str();
610
}
611

612
void HilbertSeries::from_string_rep(const string& input) {
613

614
    istringstream s(input);
615
    long i,size;
616

617
    s >> size;
618
    num.resize(size);
619
    for (i = 0; i < size; ++i) {
620
        s >> num[i];
621
    }
622

623
    vector<denom_t> denom_vector;
624
    s >> size;
625
    denom_vector.resize(size);
626
    for (i = 0; i < size; ++i) {
627
        s >> denom_vector[i];
628
    }
629

630
    denom = count_in_map<long,denom_t>(denom_vector);
631
    is_simplified = false;
632
}
633

634
// writes in a human readable format
635
ostream& operator<< (ostream& out, const HilbertSeries& HS) {
636
    HS.collectData();
637
    out << "(";
638
    // i == 0
639
    if (HS.num.size()>0) out << " " << HS.num[0];
640
    if (HS.shift != 0)   out << "*t^" << HS.shift;
641
    for (size_t i=1; i<HS.num.size(); ++i) {
642
             if ( HS.num[i]== 1 ) out << " +t^"<< i + HS.shift;
643
        else if ( HS.num[i]==-1 ) out << " -t^"<< i + HS.shift;
644
        else if ( HS.num[i] > 0 ) out << " +"  << HS.num[i] << "*t^" << i + HS.shift;
645
        else if ( HS.num[i] < 0 ) out << " -"  <<-HS.num[i] << "*t^" << i + HS.shift;
646
    }
647
    out << " ) / (";
648
    if (HS.denom.empty()) {
649
        out << " 1";
650
    }
651
    map<long, denom_t>::const_iterator it;
652
    for (it = HS.denom.begin(); it != HS.denom.end(); ++it) { 
653
        if ( it->second != 0 ) out << " (1-t^"<< it->first <<")^" << it->second;
654
    }
655
    out << " )" << std::endl;
656
    return out;
657
}
658

659

660

661
//---------------------------------------------------------------------------
662
// polynomial operations, for polynomials repr. as vector of coefficients
663
//---------------------------------------------------------------------------
664

665
// returns the coefficient vector of 1-t^i
666
template<typename Integer>
667
vector<Integer> coeff_vector(size_t i) {
668
    vector<Integer> p(i+1,0);
669
    p[0] =  1;
670
    p[i] = -1;
671
    return p;
672
}
673

674
template<typename Integer>
675
void remove_zeros(vector<Integer>& a) {
676
    size_t i=a.size();
677
    while ( i>0 && a[i-1]==0 ) --i;
678

679
    if (i < a.size()) {
680
        a.resize(i);
681
    }
682
}
683

684
// a += b  (also possible to define the += op for vector)
685
template<typename Integer>
686
void poly_add_to (vector<Integer>& a, const vector<Integer>& b) {
687
    size_t b_size = b.size();
688
    if (a.size() < b_size) {
689
        a.resize(b_size);
690
    }
691
    for (size_t i=0; i<b_size; ++i) {
692
        a[i]+=b[i];
693
    }
694
    remove_zeros(a);
695
}
696

697
// a += b*t^m
698
template<typename Integer>
699
void poly_add_to_tm (vector<Integer>& a, const vector<Integer>& b,long m) {
700
    size_t b_size=b.size();
701
    size_t b_m = b_size+m;
702
    if (a.size() < b_m) {
703
        a.resize(b_m);
704
    }
705
    for (size_t i=0; i<b_size; ++i) {
706
        a[i+m]+=b[i];
707
    }
708
    remove_zeros(a);
709
}
710

711
// a -= b  (also possible to define the -= op for vector)
712
template<typename Integer>
713
void poly_sub_to (vector<Integer>& a, const vector<Integer>& b) {
714
    size_t b_size = b.size();
715
    if (a.size() < b_size) {
716
        a.resize(b_size);
717
    }
718
    for (size_t i=0; i<b_size; ++i) {
719
        a[i]-=b[i];
720
    }
721
    remove_zeros(a);
722
}
723

724
// a *= t^m
725
template<typename Integer>
726
void poly_mult_by_tm(vector<Integer>& a, long m){
727
        long a_ori_size=a.size();
728
        a.resize(a_ori_size+m);
729
        for(long i=a_ori_size-1; i>=0;--i)
730
            a[i+m]=a[i];
731
        for(long i=0;i<m;++i)
732
            a[i]=0; 
733
}
734

735
// a * b
736

737
/* template<typename Integer>
738
vector<Integer> old_poly_mult(const vector<Integer>& a, const vector<Integer>& b) {
739
    size_t a_size = a.size();
740
    size_t b_size = b.size();
741

742
    
743
    vector<Integer> p( a_size + b_size - 1 );
744
    size_t i,j;
745
    for (i=0; i<a_size; ++i) {
746
        if (a[i] == 0) continue;
747
        for (j=0; j<b_size; ++j) {
748
            if (b[j] == 0) continue;
749
            p[i+j] += a[i]*b[j];
750
        }
751
    }
752
    return p;
753
}*/
754

755
template<typename Integer>
756
vector<Integer> karatsubamult(const vector<Integer>& a, const vector<Integer>& b) {
757

758
    size_t a_size = a.size();
759
    size_t b_size = b.size();
760
    if(a_size*b_size<=1000 || a_size <=10 || b_size<=10){
761
        return poly_mult(a,b);
762
    }
763

764
    size_t m=(a_size+1)/2;
765
    if(2*m<(b_size+1)){
766
        m=(b_size+1)/2;
767
    }
768
    
769
    vector<Integer> f0(m),f1(m),g0(m),g1(m);
770
    for(size_t i=0;i<m && i<a_size;++i)
771
        f0[i]=a[i];
772
    for(size_t i=m;i<a_size;++i)
773
        f1[i-m]=a[i];
774
    for(size_t i=0;i<m && i<b_size;++i)
775
        g0[i]=b[i];
776
    for(size_t i=m;i< b_size;++i)
777
        g1[i-m]=b[i];
778
    remove_zeros(f0);
779
    remove_zeros(f1);
780
    remove_zeros(g0);
781
    remove_zeros(g1);
782
    
783
    vector<Integer> sf=f0;
784
    vector<Integer> sg=g0;
785
    
786
    vector<Integer> mix;
787
    vector<Integer> h00;
788
    vector<Integer> h11;
789

790
    #pragma omp parallel // num_threads(3)
791
    {
792
    
793
    #pragma omp single nowait
794
    {
795
    h00=karatsubamult(f0,g0); // h00 = f0 * g0
796
    }
797
    
798
   #pragma omp single nowait
799
    {  
800
    h11=karatsubamult(f1,g1); // h11 = f1 * g1
801
    }
802

803
    #pragma omp single nowait
804
    {
805
    poly_add_to(sf,f1); // f0+f1
806
    poly_add_to(sg,g1); // g0 + g1
807
    mix=karatsubamult(sf,sg); // (f0 + f1)*(g0 + g1)
808
    }
809
    
810
    } // parallel
811
    
812
    f0.clear();
813
    g0.clear();
814
    f1.clear();
815
    g1.clear();
816
    
817
    poly_sub_to(mix,h00);  // mix = mix - f0*g0
818
    poly_sub_to(mix,h11);  // mix = mix - f1*g1
819
    
820
    poly_add_to_tm(h00,mix,m);
821
    poly_add_to_tm(h00,h11,2*m);
822

823
    return h00;
824
}
825

826
template<typename Integer>
827
vector<Integer> poly_mult(const vector<Integer>& a, const vector<Integer>& b) {
828
    size_t a_size = a.size();
829
    size_t b_size = b.size();
830
    
831
    if(a_size*b_size>1000 && a_size >10 && b_size>10){
832
        // omp_set_nested(1);
833
        return karatsubamult(a,b);
834
        // omp_set_nested(0);
835
    }
836
    
837
    vector<Integer> p( a_size + b_size - 1 );
838
    size_t i,j;
839
    for (i=0; i<a_size; ++i) {
840
        if (a[i] == 0) continue;
841
        for (j=0; j<b_size; ++j) {
842
            if (b[j] == 0) continue;
843
            p[i+j] += a[i]*b[j];
844
        }
845
    }
846
    return p;
847
}
848

849
// a *= (1-t^d)^e
850
template<typename Integer>
851
void poly_mult_to(vector<Integer>& a, long d, long e) {
852
    assert(d>0);
853
    assert(e>=0);
854
    long i;
855
    a.reserve(a.size() + d*e);
856
    while (e>0) {
857
        a.resize(a.size() + d);
858
        for (i=a.size()-1; i>=d; --i) {
859
            a[i] -= a[i-d];
860
        }
861
        e--;
862
    }
863
}
864

865
// division with remainder, a = q * b + r, deg(r) < deg(b), needs |leadcoef(b)| = 1
866
template<typename Integer>
867
void poly_div(vector<Integer>& q, vector<Integer>& r, const vector<Integer>& a, const vector<Integer>&b) {
868
    assert(b.back()!=0); // no unneeded zeros
869
    assert(b.back()==1 || b.back()==-1); // then division is always possible
870
    r = a;
871
    remove_zeros(r);
872
    size_t b_size = b.size();
873
    int degdiff = r.size()-b_size; // degree differenz
874
    if (r.size() < b_size) {
875
        q = vector<Integer>();
876
    } else {
877
        q = vector<Integer>(degdiff+1);
878
    }
879
    Integer divisor;
880
    size_t i=0;
881

882
    while (r.size() >= b_size) {
883
        
884
        divisor = r.back()/b.back();
885
        q[degdiff] = divisor;
886
        // r -= divisor * t^degdiff * b
887
        for (i=0; i<b_size; ++i) {
888
            r[i+degdiff] -= divisor * b[i];
889
        }
890
        remove_zeros(r);
891
        degdiff = r.size()-b_size;
892
    }
893

894
    return;
895
}
896

897
template<typename Integer>
898
vector<Integer> cyclotomicPoly(long n) {
899
    
900
    // the static variable is initialized only once and then stored
901
    static map<long, vector<Integer> > CyclotomicPoly = map<long, vector<Integer> >();
902
    if (CyclotomicPoly.count(n) == 0) { //it was not computed so far
903
        vector<Integer> poly, q, r;
904
        for (long i = 1; i <= n; ++i) {
905
            // compute needed and uncomputed factors
906
            if( n % i == 0 && CyclotomicPoly.count(i) == 0) {
907
                // compute the i-th poly by dividing X^i-1 by the 
908
                // d-th cycl.poly. with d divides i
909
                poly = vector<Integer>(i+1);
910
                poly[0] = -1; poly[i] = 1;  // X^i - 1
911
                for (long d = 1; d < i; ++d) { // <= i/2 should be ok
912
                    if( i % d == 0) {
913
                        poly_div(q, r, poly, CyclotomicPoly[d]);
914
                        assert(r.empty());
915
                        poly = q;
916
                    }
917
                }
918
                CyclotomicPoly[i] = poly;
919
                //cout << i << "-th cycl. pol.: " << CyclotomicPoly[i];
920
            }
921
        }
922
    }
923
    assert(CyclotomicPoly.count(n)>0);
924
    return CyclotomicPoly[n];
925
}
926

927

928

929
//---------------------------------------------------------------------------
930
// computing the Hilbert polynomial from h-vector
931
//---------------------------------------------------------------------------
932

933
// The algorithm follows "Cohen-Macaulay rings", 4.1.5 and 4.1.9.
934
// The E_vector is the vector of higher multiplicities.
935
// It is assumed that (d-1)! is used as a common denominator in the calling routine.
936

937
template<typename Integer>
938
vector<Integer> compute_e_vector(vector<Integer> Q, int dim){
939
    size_t j;
940
    int i;
941
    vector <Integer> E_Vector(dim,0); 
942
    // cout << "QQQ " << Q;  
943
    // Q.resize(dim+1);
944
    int bound=Q.size();
945
    if(bound>dim)
946
        bound=dim;  
947
    for (i = 0; i <bound; i++) {
948
        for (j = 0; j < Q.size()-i; j++) {  
949
            E_Vector[i] += Q[j];
950
        }
951
        E_Vector[i]/=permutations<Integer>(1,i);
952
        for (j = 1; j <Q.size()-i; j++) {
953
            Q[j-1]=j*Q[j];
954
        }
955
    }
956
    return E_Vector;
957
}
958

959
//---------------------------------------------------------------------------
960

961
template<typename Integer>
962
vector<Integer> compute_polynomial(vector<Integer> h_vector, int dim) {
963
    // handle dimension 0
964
    if (dim == 0)
965
        return vector<Integer>(dim);
966

967
    vector<Integer> Hilbert_Polynomial = vector<Integer>(dim);
968
    int i,j;
969
    
970
    Integer mult_factor;
971
    vector <Integer> E_Vector=compute_e_vector(h_vector, dim);
972
    vector <Integer> C(dim,0);
973
    C[0]=1;
974
    for (i = 0; i <dim; i++) {
975
        mult_factor=permutations<Integer>(i,dim);
976
        if (((dim-1-i)%2)==0) {
977
            for (j = 0; j <dim; j++) {
978
                Hilbert_Polynomial[j]+=mult_factor*E_Vector[dim-1-i]*C[j];
979
            }
980
        }
981
        else {
982
            for (j = 0; j <dim; j++) {
983
                Hilbert_Polynomial[j]-=mult_factor*E_Vector[dim-1-i]*C[j];
984
            }
985
        }
986
        for (j = dim-1; 0 <j; j--) {
987
            C[j]=(unsigned long)(i+1)*C[j]+C[j-1];
988
        }
989
        C[0]=permutations<Integer>(1,i+1);
990
    }
991

992
    return Hilbert_Polynomial;
993
}
994

995
//---------------------------------------------------------------------------
996

997
// substitutes t by (t-a), overwrites the polynomial!
998
template<typename Integer>
999
void linear_substitution(vector<Integer>& poly, const Integer& a) {
1000
    long deg = poly.size()-1;
1001
    // Iterated division by (t+a)
1002
    for (long step=0; step<deg; ++step) {
1003
        for (long i = deg-1; i >= step; --i) {
1004
            poly[i] -= a * poly[i+1];
1005
        }
1006
        //the remainders are the coefficients of the transformed polynomial
1007
    }
1008
}
1009

1010
//---------------------------------------------------------------------------
1011
IntegrationData::IntegrationData(){
1012
}
1013

1014
void IntegrationData::set_nr_coeff_quasipol(long nr_coeff){
1015
    weighted_Ehrhart_series.first.set_nr_coeff_quasipol(nr_coeff);
1016
}
1017

1018
string IntegrationData::getPolynomial() const{
1019
    return polynomial;
1020
}
1021

1022
long IntegrationData::getDegreeOfPolynomial() const{
1023
    return degree_of_polynomial; 
1024
}
1025

1026
void IntegrationData::setDegreeOfPolynomial(const long d){
1027
    degree_of_polynomial=d; 
1028
}
1029

1030
IntegrationData::IntegrationData(const string& poly){
1031
    polynomial=poly;
1032
    polynomial_is_homogeneous=false; // to be on the safe side
1033
}
1034

1035
bool IntegrationData::isWeightedEhrhartQuasiPolynomialComputed() const{
1036
    return weighted_Ehrhart_series.first.isHilbertQuasiPolynomialComputed();
1037
}
1038

1039
vector< vector<mpz_class> > IntegrationData::getWeightedEhrhartQuasiPolynomial() const{
1040
    return weighted_Ehrhart_series.first.getHilbertQuasiPolynomial();
1041
}
1042

1043
void IntegrationData::computeWeightedEhrhartQuasiPolynomial(){
1044
    weighted_Ehrhart_series.first.computeHilbertQuasiPolynomial();
1045
}
1046

1047

1048
mpz_class IntegrationData::getWeightedEhrhartQuasiPolynomialDenom() const{
1049
    return weighted_Ehrhart_series.first.getHilbertQuasiPolynomialDenom()*weighted_Ehrhart_series.second;         
1050
}
1051

1052
const vector<mpz_class>& IntegrationData::getNum_ZZ() const{
1053
    return weighted_Ehrhart_series.first.getNum();
1054
}
1055

1056
const map<long, denom_t>& IntegrationData::getDenom() const{
1057
    return weighted_Ehrhart_series.first.getDenom();
1058
}
1059

1060
const vector<mpz_class>& IntegrationData::getCyclotomicNum_ZZ() const{
1061
    return weighted_Ehrhart_series.first.getCyclotomicNum();   
1062
}
1063

1064
const map<long, denom_t>& IntegrationData::getCyclotomicDenom() const{
1065
    return weighted_Ehrhart_series.first.getCyclotomicDenom();  
1066
}
1067

1068
const pair<HilbertSeries, mpz_class>& IntegrationData::getWeightedEhrhartSeries() const{
1069
    return weighted_Ehrhart_series;
1070
}
1071

1072
mpq_class IntegrationData::getIntegral() const{
1073
    return integral;
1074
}
1075

1076
mpz_class IntegrationData::getNumeratorCommonDenom() const{
1077
    return weighted_Ehrhart_series.second;
1078
}
1079

1080
mpq_class IntegrationData::getVirtualMultiplicity() const{
1081
        return virtual_multiplicity;
1082
}
1083

1084
void IntegrationData::setIntegral(const mpq_class I){
1085
    integral=I;
1086
}
1087

1088
void IntegrationData::setVirtualMultiplicity(const mpq_class I){
1089
        virtual_multiplicity=I;
1090
}
1091

1092
void IntegrationData::setWeightedEhrhartSeries(const pair<HilbertSeries, mpz_class>& E){
1093
    weighted_Ehrhart_series=E;
1094
    weighted_Ehrhart_series.first.adjustShift();
1095
}
1096

1097
void IntegrationData::setHomogeneity(const bool hom){
1098
    polynomial_is_homogeneous=hom;
1099
}
1100

1101

1102
bool IntegrationData::isPolynomialHomogeneous() const{
1103
    return polynomial_is_homogeneous;
1104
}
1105

1106
} //end namespace libnormaliz
1107

1108
#ifdef NMZ_MIC_OFFLOAD
1109
#pragma offload_attribute (pop)
1110
#endif
1111

1112