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
//---------------------------------------------------------------------------
25

26
#include <algorithm>
27
#include <sstream>
28
#include <math.h>
29
#include "libnormaliz/integer.h"
30
#include "libnormaliz/vector_operations.h"
31

32
//---------------------------------------------------------------------------
33

34
namespace libnormaliz {
35
using namespace std;
36

37
bool try_convert(long& ret, const long long& val) {
38
    if (fits_long_range(val)) {
39
        ret = val;
40
        return true;
41
    }
42
    return false;
43
}
44

45
bool try_convert(long& ret, const mpz_class& val) {
46
    if (!val.fits_slong_p()) {
47
        return false;
48
    }
49
    ret = val.get_si();
50
    return true;
51
}
52

53
bool try_convert(long long& ret, const mpz_class& val) {
54
    if (val.fits_slong_p()) {
55
        ret = val.get_si();
56
        return true;
57
    }
58
    if (sizeof(long long)==sizeof(long)) {
59
        return false;
60
    }
61
    mpz_class quot;
62
    ret = mpz_fdiv_q_ui(quot.get_mpz_t(), val.get_mpz_t(), LONG_MAX); //returns remainder
63
    if (!quot.fits_slong_p()){
64
        return false;
65
    }
66
    ret += ((long long) quot.get_si())*((long long) LONG_MAX);
67
    return true;
68
}
69

70
bool try_convert(mpz_class& ret, const long long& val) {
71
    if (fits_long_range(val)) {
72
        ret = mpz_class(long(val));
73
    } else {
74
        ret = mpz_class(long (val % LONG_MAX)) + mpz_class(LONG_MAX) * mpz_class(long(val/LONG_MAX));
75
    }
76
    return true;
77
}
78

79
bool fits_long_range(long long a) {
80
    return sizeof(long long) == sizeof(long) || (a <= LONG_MAX && a >= LONG_MIN);
81
}
82

83
template<typename FromType>
84
bool try_convert(nmz_float& ret, const FromType& val){
85
    ret= (nmz_float) val;
86
    return true;
87
}
88

89
bool try_convert(nmz_float& ret, const mpz_class& val){    
90
    ret=val.get_d();
91
    return true;
92
}
93

94
template<typename ToType>
95
bool try_convert(ToType& ret, const nmz_float& val){
96
    mpz_class bridge;
97
    if(!try_convert(bridge,val))
98
        return false;
99
    return try_convert(ret,bridge);
100
}
101

102
bool try_convert(mpz_class& ret, const nmz_float& val){    
103
    ret=mpz_class(val);
104
    return true;
105
}
106

107
//---------------------------------------------------------------------------
108

109
template <typename Integer>
110
Integer gcd(const Integer& a, const Integer& b){
111
    if (a==0) {
112
        return Iabs<Integer>(b);
113
    }
114
    if (b==0) {
115
        return Iabs<Integer>(a);
116
    }
117
    Integer q0,q1,r;
118
    q0=Iabs<Integer>(a);
119
    r=Iabs<Integer>(b);
120
    do {
121
        q1=r;
122
        r=q0%q1;
123
        q0=q1;
124
    } while (r!=0);
125
    return q1;
126
}
127

128
template<> mpz_class gcd<mpz_class>(const mpz_class& a, const mpz_class& b) {
129
    mpz_class g;
130
    mpz_gcd (g.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
131
    return g;
132
}
133

134
template void sign_adjust_and_minimize<long long>(const long long& a, const long long& b, long long& d, long long& u, long long&v);
135
template long long ext_gcd<long long>(const long long& a, const long long& b, long long& u, long long&v);
136

137

138
template <typename Integer>
139
void sign_adjust_and_minimize(const Integer& a, const Integer& b, Integer& d, Integer& u, Integer&v){
140
    if(d<0){
141
        d=-d;
142
        u=-u;
143
        v=-v;    
144
    }
145
    // cout << u << " " << v << endl;
146
    if(b==0)
147
        return;
148
        
149
    Integer sign=1;
150
    if(a<0)
151
        sign=-1;
152
    Integer u1= (sign*u) % (Iabs(b)/d);
153
    if(u1==0)
154
        u1+=Iabs(b)/d;
155
    u=sign*u1;
156
    v=(d-a*u)/b;
157
}
158

159

160
template <typename Integer>
161
Integer ext_gcd(const Integer& a, const Integer& b, Integer& u, Integer&v){
162

163
    u=1;
164
    v=0;
165
    Integer d=a;
166
    if (b==0) {
167
        sign_adjust_and_minimize(a,b,d,u,v);
168
        return(d);
169
    }
170
    Integer v1=0;
171
    Integer v3=b;
172
    Integer q,t1,t3;
173
    while(v3!=0){
174
        q=d/v3;
175
        t3=d-q*v3;
176
        t1=u-q*v1;
177
        u=v1;
178
        d=v3;
179
        v1=t1;
180
        v3=t3;
181
    }
182
    sign_adjust_and_minimize(a,b,d,u,v);
183
    return(d);
184
}
185

186
//---------------------------------------------------------------------------
187

188
template <typename Integer>
189
Integer lcm(const Integer& a, const Integer& b){
190
    if ((a==0)||(b==0)) {
191
        return 0;
192
    }
193
    else
194
        return Iabs<Integer>(a*b/gcd<Integer>(a,b));
195
}
196

197
template<> mpz_class lcm<mpz_class>(const mpz_class& a, const mpz_class& b) {
198
    mpz_class g;
199
    mpz_lcm (g.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
200
    return g;
201
}
202

203
//---------------------------------------------------------------------------
204

205
template <typename Integer>
206
size_t decimal_length(Integer a){
207
    ostringstream test;
208
    test << a;
209
    return test.str().size();
210
}
211

212
//---------------------------------------------------------------------------
213

214
template <typename Integer>
215
Integer permutations(const size_t& a, const size_t& b){
216
    unsigned long i;
217
    Integer P=1;
218
    for (i = a+1; i <= b; i++) {
219
        P*=i;
220
    }
221
    return P;
222
}
223

224
//---------------------------------------------------------------------------
225

226
template<typename Integer> 
227
Integer permutations_modulo(const size_t& a, const size_t& b, long m) {
228
    unsigned long i;
229
    Integer P=1;
230
    for (i = a+1; i <= b; i++) {
231
        P*=i; P%=m;
232
    }
233
    return P;
234
}
235
//---------------------------------------------------------------------------
236

237
template<typename Integer>
238
Integer int_max_value_dual(){
239
    Integer k=sizeof(Integer)*8-2;  // number of bytes convetred to number of bits
240
    Integer test=1;
241
    test = test << k;  // (maximal positive number)/2^k
242
    return test;
243
}
244

245
//---------------------------------------------------------------------------
246

247
template<typename Integer>
248
Integer int_max_value_primary(){
249
    Integer k=sizeof(Integer)*8-12;  // number of bytes convetred to number of bits
250
    Integer test=1;
251
    test = test << k;  // (maximal positive number)/2^k
252
    // test=0; // 10000;
253
    return test;
254
}
255

256
//---------------------------------------------------------------------------
257

258
template<>
259
mpz_class int_max_value_dual<mpz_class>(){
260
    assert(false);
261
    return 0;
262
}
263

264
//---------------------------------------------------------------------------
265

266
template<>
267
mpz_class int_max_value_primary<mpz_class>(){
268
    assert(false);
269
    return 0;
270
}
271

272

273
//---------------------------------------------------------------------------
274

275
template<typename Integer>
276
void check_range_list(const CandidateList<Integer>& ll){
277
    check_range_list(ll.Candidates);
278
}
279

280
//---------------------------------------------------------------------------
281

282
template<typename Integer>
283
void check_range_list(const std::list<Candidate<Integer> >& ll){
284

285
    if (using_GMP<Integer>())
286
        return;
287
        
288
    typename list<Candidate<Integer> >::const_iterator v=ll.begin();
289
    
290
    Integer test=int_max_value_dual<Integer>();
291
    // cout << "test " << test << endl;
292
    
293
    for(;v!=ll.end();++v){
294
        for(size_t i=0;i<v->values.size();++i)
295
            if(Iabs(v->values[i])>= test){
296
            // cout << *v;
297
            // cout << "i " << i << " " << Iabs((*v)[i]) << endl;
298
                throw ArithmeticException("Vector entry out of range. Imminent danger of arithmetic overflow.");
299
            }
300
                    
301
    }
302
    
303

304
}
305

306
//---------------------------------------------------------------------------
307
 template<typename Integer>
308
void minimal_remainder(const Integer& a, const Integer&b, Integer& quot, Integer& rem) {
309

310
    quot=a/b;
311
    rem=a-quot*b;
312
    if(rem==0)
313
        return;
314
    if(2*Iabs(rem)>Iabs(b)){
315
        if((rem<0 && b>0) || (rem >0 && b<0)){                
316
            rem+=b;
317
            quot--;
318
        }
319
        else{
320
            rem-=b;
321
            quot++;                
322
        }
323
    }
324
}
325

326
//---------------------------------------------------------------------------
327

328
mpq_class dec_fraction_to_mpq(string s){
329
    
330
    // cout << "s " << s <<endl;
331
    
332
    mpz_class sign=1;
333
    if(s[0]=='+')
334
        s=s.substr(1);
335
    else
336
        if(s[0]=='-'){
337
            s=s.substr(1);
338
            sign=-1;
339
    }
340
    
341
    if(s[0]=='+' || s[0]=='-')
342
            throw BadInputException("Error in decimal fraction");
343
    
344
    string int_string,frac_string,exp_string;
345
    size_t frac_part_length=0;
346
    size_t pos_point=s.find(".");
347
    size_t pos_E=s.find("e");
348
    if(pos_point!=string::npos){
349
        int_string=s.substr(0,pos_point);
350
        if(pos_E!=string::npos){
351
            frac_part_length=pos_E-(pos_point+1);
352
        }
353
        else
354
            frac_part_length=s.size()-(pos_point+1);
355
        frac_string=s.substr(pos_point+1,frac_part_length);
356
        if(frac_string[0]=='+' || frac_string[0]=='-')
357
            throw BadInputException("Error in decimal fraction");
358
    }
359
    else
360
        int_string=s.substr(0,pos_E);
361
    if(pos_E!=string::npos)
362
        exp_string=s.substr(pos_E+1,s.size()-(pos_E+1));
363
    
364
    /* cout << "int  " << int_string << endl;
365
    cout << "frac " << frac_string << endl;
366
    cout << "exp  " << exp_string << endl; */
367
    
368
    // remove leading 0
369
    while(int_string[0]=='0')
370
        int_string=int_string.substr(1);
371
    while(frac_string[0]=='0')
372
        frac_string=frac_string.substr(1);
373
    
374
    mpq_class int_part, frac_part, exp_part;
375
    if(!int_string.empty())
376
        int_part=mpz_class(int_string);
377

378
    if(pos_E==0)
379
        int_part=1;
380
    
381
    // cout << "int_part " << int_part << endl;
382
    
383
    mpz_class den=1;
384
    if(!frac_string.empty()){
385
        frac_part=mpz_class(frac_string);
386
        for(size_t i=0;i<frac_part_length;++i)
387
            den*=10;        
388
    }
389
    // cout << "frac_part " << frac_part << endl;
390
    mpq_class result=int_part;
391
    if(frac_part!=0)
392
        result+=frac_part/den;
393
    if(!exp_string.empty()){
394
        long expo=stol(exp_string);
395
        long abs_expo=Iabs(expo);
396
        // cout << "expo " << expo << endl;
397
        mpz_class factor=1;
398
        for(long i=0;i< abs_expo;++i)
399
            factor*=10;
400
        if(expo>=0)
401
            result*=factor;
402
        else
403
            result/=factor;
404
    }
405
    /* cout <<" result " << sign*result << endl;
406
    cout << "==========" << endl; */
407
    return sign*result;
408
}
409

410
//----------------------------------------------------------------------
411
// the next function produce an integer quotient and determine whether
412
// there is a remainder
413

414
bool int_quotient(long& Quot, const long long& Num, const long& Den){
415
    
416
    Quot=Iabs(Num)/Iabs(Den);
417
    return Quot*Iabs(Den)!=Iabs(Num);    
418
}
419

420
bool int_quotient(long long& Quot, const long& Num, const long& Den){
421
    
422
    Quot=Iabs(Num)/Iabs(Den);
423
    return Quot*Iabs(Den)!=Iabs(Num);    
424
}
425

426

427
bool int_quotient(mpz_class& Quot, const mpz_class& Num, const mpz_class& Den){
428
    
429
    Quot=Iabs(Num)/Iabs(Den);
430
    return Quot*Iabs(Den)!=Iabs(Num);    
431
}
432

433
template<typename IntegerRet>
434
bool int_quotient(IntegerRet& Quot, const nmz_float& Num, const nmz_float& Den){
435
   
436
    nmz_float FloatQuot=Iabs(Num)/Iabs(Den); // cout << "FF " << FloatQuot << endl;
437
    nmz_float IntQuot=trunc(FloatQuot+nmz_epsilon);      // cout << "II " << IntQuot << endl;
438
    Quot=convertTo<IntegerRet>(IntQuot);     // cout << "QQ " <<  Quot << endl;
439
    return FloatQuot-IntQuot > nmz_epsilon;    
440
}
441

442

443
} //end namespace libnormaliz
444

445