CoCalc -- Qvector_operations.cpp

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.

GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
Project: cocalc-sagemath-dev-slelievre
Path: gap4r8 / pkg / NormalizInterface-1.0.2 / Normaliz.git / Qsource / libQnormaliz / Qvector_operations.cpp
Views: ⁴¹⁸³⁸⁶
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 <iostream>
27
#include <string>
28
#include <algorithm>
29
#include <list>
30

31
#include "libQnormaliz/Qinteger.h"
32
#include "libQnormaliz/Qvector_operations.h"
33
#include "libQnormaliz/Qmatrix.h"
34

35
//---------------------------------------------------------------------------
36

37
namespace libQnormaliz {
38
using namespace std;
39

40
//---------------------------------------------------------------------------
41

42
template<typename Number>
43
Number v_scalar_product(const vector<Number>& av,const vector<Number>& bv){
44
    //loop stretching ; brings some small speed improvement
45

46
    Number ans = 0;
47
    size_t i,n=av.size();
48

49
    typename vector<Number>::const_iterator a=av.begin(), b=bv.begin();
50

51
    if( n >= 16 )
52
    {
53
        for( i = 0; i < ( n >> 4 ); ++i, a += 16, b +=16 ){
54
            ans += a[0] * b[0];
55
            ans += a[1] * b[1];
56
            ans += a[2] * b[2];
57
            ans += a[3] * b[3];
58
            ans += a[4] * b[4];
59
            ans += a[5] * b[5];
60
            ans += a[6] * b[6];
61
            ans += a[7] * b[7];
62
            ans += a[8] * b[8];
63
            ans += a[9] * b[9];
64
            ans += a[10] * b[10];
65
            ans += a[11] * b[11];
66
            ans += a[12] * b[12];
67
            ans += a[13] * b[13];
68
            ans += a[14] * b[14];
69
            ans += a[15] * b[15];
70
        }
71

72
        n -= i<<4;
73
    }
74

75
    if( n >= 8)
76
    {
77
        ans += a[0] * b[0];
78
        ans += a[1] * b[1];
79
        ans += a[2] * b[2];
80
        ans += a[3] * b[3];
81
        ans += a[4] * b[4];
82
        ans += a[5] * b[5];
83
        ans += a[6] * b[6];
84
        ans += a[7] * b[7];
85

86
        n -= 8;
87
        a += 8;
88
        b += 8;
89
    }
90

91
    if( n >= 4)
92
    {
93
        ans += a[0] * b[0];
94
        ans += a[1] * b[1];
95
        ans += a[2] * b[2];
96
        ans += a[3] * b[3];
97

98
        n -= 4;
99
        a += 4;
100
        b += 4;
101
    }
102

103
    if( n >= 2)
104
    {
105
        ans += a[0] * b[0];
106
        ans += a[1] * b[1];
107

108
        n -= 2;
109
        a += 2;
110
        b += 2;
111
    }
112

113
    if(n>0)
114
        ans += a[0]*b[0];
115
        
116
    return ans;
117
}
118

119
//---------------------------------------------------------------------------
120

121
template<typename Number>
122
Number v_scalar_product_unequal_vectors_end(const vector<Number>& a,const vector<Number>& b){
123
    Number ans = 0;
124
    size_t i,n=a.size(),m=b.size();
125
    for (i = 1; i <= n; i++) {
126
        ans+=a[n-i]*b[m-i];
127
    }
128
    return ans;
129
}
130

131
//---------------------------------------------------------------------------
132

133
template<typename Number>
134
vector<Number> v_add(const vector<Number>& a,const vector<Number>& b){
135
   assert(a.size() == b.size());
136
    size_t i,s=a.size();
137
    vector<Number> d(s);
138
    for (i = 0; i <s; i++) {
139
        d[i]=a[i]+b[i];
140
    }
141
    return d;
142
}
143

144
//---------------------------------------------------------------------------
145

146
template<typename Number>
147
void v_add_result(vector<Number>& result, const size_t s, const vector<Number>& a,const vector<Number>& b){
148
   assert(a.size() == b.size() && a.size() == result.size());
149
    size_t i;
150
    // vector<Number> d(s);
151
    for (i = 0; i <s; i++) {
152
        result[i]=a[i]+b[i];
153
    }
154
    // return d;
155
}
156

157
//---------------------------------------------------------------------------
158

159
template<typename Number>
160
vector<Number>& v_abs(vector<Number>& v){
161
    size_t i, size=v.size();
162
    for (i = 0; i < size; i++) {
163
        if (v[i]<0) v[i] = Iabs(v[i]);
164
    }
165
    return v;
166
}
167

168
//---------------------------------------------------------------------------
169

170
template<typename Number>
171
vector<Number> v_abs_value(vector<Number>& v){
172
    size_t i, size=v.size();
173
    vector<Number> w=v;
174
    for (i = 0; i < size; i++) {
175
        if (v[i]<0) w[i] = Iabs(v[i]);
176
    }
177
    return w;
178
}
179

180
//---------------------------------------------------------------------------
181

182
// the following function removes the denominators and then extracts the Gcd of the numerators
183
mpq_class v_simplify(vector<mpq_class>& v){
184
    size_t size=v.size();
185
    mpz_class d=1;
186
    for (size_t i = 0; i < size; i++)
187
        //d=lcm(d,v[i].get_den());  // GMP C++ function only available in GMP >= 6.1
188
        mpz_lcm(d.get_mpz_t(), d.get_mpz_t(), v[i].get_den().get_mpz_t());
189
    for (size_t i = 0; i < size; i++)
190
        v[i]*=d;
191
    mpz_class g=0;
192
    for (size_t i = 0; i < size; i++)
193
        //g=gcd(g,v[i].get_num());  //  GMP C++ function only available in GMP >= 6.1
194
        mpz_gcd(g.get_mpz_t(), g.get_mpz_t(), v[i].get_num().get_mpz_t());
195
    if (g==0)
196
        return 0;
197
    for (size_t i = 0; i < size; i++)
198
        v[i]/=g;
199
    return 1;
200
}
201

202
//---------------------------------------------------------------------------
203

204
template<typename T>
205
vector<T> v_merge(const vector<T>& a, const T& b) {
206
    size_t s=a.size();
207
    vector<T> c(s+1);
208
    for (size_t i = 0; i < s; i++) {
209
        c[i]=a[i];
210
    }
211
    c[s] = b;
212
    return c;
213
}
214

215
//---------------------------------------------------------------------------
216

217
template<typename T>
218
vector<T> v_merge(const vector<T>& a,const vector<T>& b){
219
    size_t s1=a.size(), s2=b.size(), i;
220
    vector<T> c(s1+s2);
221
    for (i = 0; i < s1; i++) {
222
        c[i]=a[i];
223
    }
224
    for (i = 0; i < s2; i++) {
225
        c[s1+i]=b[i];
226
    }
227
    return c;
228
}
229
//---------------------------------------------------------------------------
230

231
template<typename T>
232
vector<T> v_cut_front(const vector<T>& v, size_t size){
233
    size_t s,k;
234
    vector<T> tmp(size);
235
    s=v.size()-size;
236
    for (k = 0; k < size; k++) {
237
        tmp[k]=v[s+k];
238
    }
239
    return tmp;
240
}
241

242
//---------------------------------------------------------------------------
243

244
template<typename Number>
245
vector<key_t> v_non_zero_pos(const vector<Number>& v){
246
    vector<key_t> key;
247
    size_t size=v.size();
248
    key.reserve(size);
249
    for (key_t i = 0; i <size; i++) {
250
        if (v[i]!=0) {
251
            key.push_back(i);
252
        }
253
    }
254
    return key;
255
}
256

257
//---------------------------------------------------------------------------
258

259
template<typename Number>
260
bool v_is_zero(const vector<Number>& v) {
261
    for (size_t i = 0; i < v.size(); ++i) {
262
        if (v[i] != 0) return false;
263
    }
264
    return true;
265
}
266

267
//---------------------------------------------------------------------------
268

269
template<typename Number>
270
bool v_is_symmetric(const vector<Number>& v) {
271
    for (size_t i = 0; i < v.size()/2; ++i) {
272
        if (v[i] != v[v.size()-1-i]) return false;
273
    }
274
    return true;
275
}
276

277
//---------------------------------------------------------------------------
278

279
template<typename Number>
280
bool v_is_nonnegative(const vector<Number>& v) {
281
    for (size_t i = 0; i < v.size(); ++i) {
282
        if (v[i] <0) return false;
283
    }
284
    return true;
285
}
286

287

288
//---------------------------------------------------------------------------
289

290
template<typename Number>
291
void v_el_trans(const vector<Number>& av,vector<Number>& bv, const Number& F, const size_t& start){
292

293
    size_t i,n=av.size();
294

295
    typename vector<Number>::const_iterator a=av.begin();
296
    typename vector<Number>::iterator b=bv.begin();
297

298
    a += start;
299
    b += start;
300
    n -= start;
301

302

303
    if( n >= 8 )
304
    {
305
        for( i = 0; i < ( n >> 3 ); ++i, a += 8, b += 8 ){
306
            b[0] += F*a[0];
307
            b[1] += F*a[1];
308
            b[2] += F*a[2];
309
            b[3] += F*a[3];
310
            b[4] += F*a[4];
311
            b[5] += F*a[5];
312
            b[6] += F*a[6];
313
            b[7] += F*a[7];
314
        }
315
        n -= i << 3;
316
    }
317

318
    if( n >= 4)
319
    {
320
        b[0] += F*a[0];
321
        b[1] += F*a[1];
322
        b[2] += F*a[2];
323
        b[3] += F*a[3];
324

325
        n -=4;
326
        a +=4;
327
        b +=4;
328
    }
329

330
    if( n >= 2)
331
    {
332
        b[0] += F*a[0];
333
        b[1] += F*a[1];
334

335
        n -=2;
336
        a +=2;
337
        b +=2;
338
    }
339

340
    if(n>0)
341
        b[0] += F*a[0];
342
}
343

344
//---------------------------------------------------------------
345

346
vector<bool> v_bool_andnot(const vector<bool>& a, const vector<bool>& b) {
347
    assert(a.size() == b.size());
348
    vector<bool> result(a);
349
    for (size_t i=0; i<b.size(); ++i) {
350
        if (b[i])
351
            result[i]=false;
352
    }
353
    return result;
354
}
355

356
// swaps entry i and j of the vector<bool> v
357
void v_bool_entry_swap(vector<bool>& v, size_t i, size_t j) {
358
    if (v[i] != v[j]) {
359
        v[i].flip();
360
        v[j].flip();
361
    }
362
}
363

364

365
vector<key_t> identity_key(size_t n){
366
    vector<key_t> key(n);
367
    for(size_t k=0;k<n;++k)
368
        key[k]=k;
369
    return key;
370
}
371

372
//---------------------------------------------------------------
373
// Sorting
374

375
template <typename T>
376
void order_by_perm(vector<T>& v, const vector<key_t>& permfix){
377
    
378
    vector<key_t> perm=permfix; // we may want to use permfix a second time
379
    vector<key_t> inv(perm.size());
380
    for(key_t i=0;i<perm.size();++i)
381
        inv[perm[i]]=i;
382
    for(key_t i=0;i<perm.size();++i){
383
        key_t j=perm[i];
384
        swap(v[i],v[perm[i]]);        
385
        swap(perm[i],perm[inv[i]]);        
386
        swap(inv[i],inv[j]);                
387
    }
388
}
389

390
// vector<bool> is special
391
template <>
392
void order_by_perm(vector<bool>& v, const vector<key_t>& permfix){
393
    
394
    vector<key_t> perm=permfix; // we may want to use permfix a second time
395
    vector<key_t> inv(perm.size());
396
    for(key_t i=0;i<perm.size();++i)
397
        inv[perm[i]]=i;
398
    for(key_t i=0;i<perm.size();++i){
399
        key_t j=perm[i];
400
        // v.swap(v[i],v[perm[i]]);
401
        v_bool_entry_swap(v,i,perm[i]);
402
        swap(perm[i],perm[inv[i]]);        
403
        swap(inv[i],inv[j]);                
404
    }
405
}
406

407
// make random vector of length n with entries between -m and m
408
template <typename Number>
409
vector<Number> v_random(size_t n, long m){
410
    vector<Number> result(n);
411
    for(size_t i=0;i<n;++i)
412
        result[i]=rand()%(2*m+1)-m;
413
    return result;    
414
}
415

416
template bool v_is_nonnegative<long>(const vector<long>&);
417
template bool v_is_nonnegative<long long>(const vector<long long>&);
418
template bool v_is_nonnegative<mpz_class>(const vector<mpz_class>&);
419

420
template bool v_is_symmetric<long>(const vector<long>&);
421
template bool v_is_symmetric<long long>(const vector<long long>&);
422
template bool v_is_symmetric<mpz_class>(const vector<mpz_class>&);
423
// 
424

425
template void v_add_result<long     >(vector<long     >&, size_t, const vector<long     >&, const vector<long     >&);
426
template void v_add_result<long long>(vector<long long>&, size_t, const vector<long long>&, const vector<long long>&);
427
template void v_add_result<mpz_class>(vector<mpz_class>&, size_t, const vector<mpz_class>&, const vector<mpz_class>&);
428

429
template long v_scalar_product(const vector<long>& a,const vector<long>& b);
430
template long long v_scalar_product(const vector<long long>& a,const vector<long long>& b);
431
template mpz_class v_scalar_product(const vector<mpz_class>& a,const vector<mpz_class>& b);
432

433
} // end namespace libQnormaliz
434

435
Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.

Product

Resources

Company

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more, all in one place.

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.
