CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In

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

| Download

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

Views: 418386
Tweet
Share
Share
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 <fstream>

27
#include <set>

28
#include <algorithm>

29
#include <math.h>

30


31
#include "libQnormaliz/Qmatrix.h"

32
#include "libQnormaliz/Qvector_operations.h"

33
#include "libQnormaliz/Qnormaliz_exception.h"

34
#include "libQnormaliz/Qsublattice_representation.h"

35


36
//---------------------------------------------------------------------------

37


38
namespace libQnormaliz {

39
using namespace std;

40


41
//---------------------------------------------------------------------------

42
//Public

43
//---------------------------------------------------------------------------

44


45
template<typename Number>

46
Matrix<Number>::Matrix(){

47
    nr=0;

48
    nc=0;

49
}

50


51
//---------------------------------------------------------------------------

52


53
template<typename Number>

54
Matrix<Number>::Matrix(size_t dim){

55
    nr=dim;

56
    nc=dim;

57
    elem = vector< vector<Number> >(dim, vector<Number>(dim));

58
    for (size_t i = 0; i < dim; i++) {

59
        elem[i][i]=1;

60
    }

61
}

62


63
//---------------------------------------------------------------------------

64


65
template<typename Number>

66
Matrix<Number>::Matrix(size_t row, size_t col){

67
    nr=row;

68
    nc=col;

69
    elem = vector< vector<Number> >(row, vector<Number>(col));

70
}

71


72
//---------------------------------------------------------------------------

73


74
template<typename Number>

75
Matrix<Number>::Matrix(size_t row, size_t col, Number value){

76
    nr=row;

77
    nc=col;

78
    elem = vector< vector<Number> > (row, vector<Number>(col,value));

79
}

80


81
//---------------------------------------------------------------------------

82


83
template<typename Number>

84
Matrix<Number>::Matrix(const vector< vector<Number> >& new_elem){

85
    nr=new_elem.size();

86
    if (nr>0) {

87
        nc=new_elem[0].size();

88
        elem=new_elem;

89
        //check if all rows have the same length

90
        for (size_t i=1; i<nr; i++) {

91
            if (elem[i].size() != nc) {

92
                throw BadInputException("Inconsistent lengths of rows in matrix!");

93
            }

94
        }

95
    } else {

96
        nc=0;

97
    }

98
}

99


100
//---------------------------------------------------------------------------

101


102
template<typename Number>

103
Matrix<Number>::Matrix(const list< vector<Number> >& new_elem){

104
    nr = new_elem.size();

105
    elem = vector< vector<Number> > (nr);

106
    nc = 0;

107
    size_t i=0;

108
    typename list< vector<Number> >::const_iterator it=new_elem.begin();

109
    for(; it!=new_elem.end(); ++it, ++i) {

110
        if(i == 0) {

111
            nc = (*it).size();

112
        } else {

113
            if ((*it).size() != nc) {

114
                throw BadInputException("Inconsistent lengths of rows in matrix!");

115
            }

116
        }

117
        elem[i]=(*it);

118
    }

119
}

120


121
//---------------------------------------------------------------------------

122
/*

123
template<typename Number>

124
void Matrix<Number>::write(istream& in){

125
    size_t i,j;

126
    for(i=0; i<nr; i++){

127
        for(j=0; j<nc; j++) {

128
            in >> elem[i][j];

129
        }

130
    }

131
}

132
*/

133
//---------------------------------------------------------------------------

134


135
template<typename Number>

136
void Matrix<Number>::write_column(size_t col, const vector<Number>& data){

137
    assert(col >= 0);

138
    assert(col < nc); 

139
    assert(nr == data.size());

140


141
    for (size_t i = 0; i < nr; i++) {

142
        elem[i][col]=data[i];

143
    }

144
}

145


146
//---------------------------------------------------------------------------

147


148
template<typename Number>

149
void Matrix<Number>::print(const string& name,const string& suffix) const{

150
    string file_name = name+"."+suffix;

151
    const char* file = file_name.c_str();

152
    ofstream out(file);

153
    print(out);

154
    out.close();

155
}

156


157
//---------------------------------------------------------------------------

158


159
template<typename Number>

160
void Matrix<Number>::print_append(const string& name,const string& suffix) const{

161
    string file_name = name+"."+suffix;

162
    const char* file = file_name.c_str();

163
    ofstream out(file,ios_base::app);

164
    print(out);

165
    out.close();

166
}

167


168
//---------------------------------------------------------------------------

169


170
template<typename Number>

171
void Matrix<Number>::print(ostream& out) const{

172
    size_t i,j;

173
    out<<nr<<endl<<nc<<endl;

174
    for (i = 0; i < nr; i++) {

175
        for (j = 0; j < nc; j++) {

176
            out<<elem[i][j]<<" ";

177
        }

178
        out<<endl;

179
    }

180
}

181


182
//---------------------------------------------------------------------------

183


184
template<typename Number>

185
void Matrix<Number>::pretty_print(ostream& out, bool with_row_nr) const{

186
    if(nr>1000000 && !with_row_nr){

187
        print(out);

188
        return;

189
    }

190
    size_t i,j,k;

191
    vector<size_t> max_length = maximal_decimal_length_columnwise();

192
    size_t max_index_length = decimal_length(nr);

193
    for (i = 0; i < nr; i++) {

194
        if (with_row_nr) {

195
            for (k= 0; k <= max_index_length - decimal_length(i); k++) {

196
                out<<" ";

197
            }

198
            out << i << ": ";

199
        }

200
        for (j = 0; j < nc; j++) {

201
            ostringstream to_print;

202
            to_print << elem[i][j];

203
            for (k= 0; k <= max_length[j] - to_print.str().size(); k++) {

204
                out<<" ";

205
            }

206
            out<< to_print.str();

207
        }

208
        out<<endl;

209
    }

210
}

211


212
/*

213
 * string to_print;

214
            ostringstream(to_print) << elem[i][j];

215
            cout << elem[i][j] << " S " << to_print << " L " << decimal_length(elem[i][j]) << endl;

216
            for (k= 0; k <= max_length[j] - to_print.size(); k++) {

217
                out<<" ";

218
            }

219
            out << to_print;

220
*/

221
//---------------------------------------------------------------------------

222


223
template<typename Number>

224
size_t Matrix<Number>::nr_of_rows () const{

225
    return nr;

226
}

227


228
//---------------------------------------------------------------------------

229


230
template<typename Number>

231
size_t Matrix<Number>::nr_of_columns () const{

232
    return nc;

233
}

234


235
//---------------------------------------------------------------------------

236


237
template<typename Number>

238
void Matrix<Number>::set_nr_of_columns(size_t c){

239
    nc=c;

240
}

241


242
//---------------------------------------------------------------------------

243


244
template<typename Number>

245
void Matrix<Number>::random (int mod) {

246
    size_t i,j;

247
    int k;

248
    for (i = 0; i < nr; i++) {

249
        for (j = 0; j < nc; j++) {

250
            k = rand();

251
            elem[i][j]=k % mod;

252
        }

253
    }

254
}

255
//---------------------------------------------------------------------------

256
  

257
template<typename Number>

258
void Matrix<Number>::set_zero() {

259
    size_t i,j;

260
    for (i = 0; i < nr; i++) {

261
        for (j = 0; j < nc; j++) {

262
            elem[i][j] = 0;

263
        }

264
    }

265
}

266


267
//---------------------------------------------------------------------------

268


269
template<typename Number>

270
void Matrix<Number>::select_submatrix(const Matrix<Number>& mother, const vector<key_t>& rows){

271


272
    assert(nr>=rows.size());

273
    assert(nc>=mother.nc);

274
    

275
    size_t size=rows.size(), j;

276
    for (size_t i=0; i < size; i++) {

277
        j=rows[i];

278
        for(size_t k=0;k<mother.nc;++k)

279
            elem[i][k]=mother[j][k];

280
    }

281
}

282


283
//---------------------------------------------------------------------------

284


285
template<typename Number>

286
void Matrix<Number>::select_submatrix_trans(const Matrix<Number>& mother, const vector<key_t>& rows){

287


288
    assert(nc>=rows.size());

289
    assert(nr>=mother.nc);

290
    

291
    size_t size=rows.size(), j;

292
    for (size_t i=0; i < size; i++) {

293
        j=rows[i];

294
        for(size_t k=0;k<mother.nc;++k)

295
            elem[k][i]=mother[j][k];

296
    }

297
}

298


299
//---------------------------------------------------------------------------

300


301
template<typename Number>

302
Matrix<Number> Matrix<Number>::submatrix(const vector<key_t>& rows) const{

303
    size_t size=rows.size(), j;

304
    Matrix<Number> M(size, nc);

305
    for (size_t i=0; i < size; i++) {

306
        j=rows[i];

307
        assert(j >= 0);

308
        assert(j < nr);

309
        M.elem[i]=elem[j];

310
    }

311
    return M;

312
}

313


314
//---------------------------------------------------------------------------

315


316
template<typename Number>

317
Matrix<Number> Matrix<Number>::submatrix(const vector<int>& rows) const{

318
    size_t size=rows.size(), j;

319
    Matrix<Number> M(size, nc);

320
    for (size_t i=0; i < size; i++) {

321
        j=rows[i];

322
        assert(j >= 0);

323
        assert(j < nr);

324
        M.elem[i]=elem[j];

325
    }

326
    return M;

327
}

328


329
//---------------------------------------------------------------------------

330


331
template<typename Number>

332
Matrix<Number> Matrix<Number>::submatrix(const vector<bool>& rows) const{

333
    assert(rows.size() == nr);

334
    size_t size=0;

335
    for (size_t i = 0; i <rows.size(); i++) {

336
        if (rows[i]) {

337
            size++;

338
        }

339
    }

340
    Matrix<Number> M(size, nc);

341
    size_t j = 0;

342
    for (size_t i = 0; i < nr; i++) {

343
        if (rows[i]) {

344
            M.elem[j++] = elem[i];

345
        }

346
    }

347
    return M;

348
}

349


350
//---------------------------------------------------------------------------

351


352
template<typename Number>

353
Matrix<Number>& Matrix<Number>::remove_zero_rows() {

354
    size_t from = 0, to = 0; // maintain to <= from

355
    while (from < nr && v_is_zero(elem[from])) from++; //skip zero rows

356
    while (from < nr) {  // go over matrix

357
        // now from is a non-zero row

358
        if (to != from) elem[to].swap(elem[from]);

359
        ++to; ++from;

360
        while (from < nr && v_is_zero(elem[from])) from++; //skip zero rows

361
    }

362
    nr = to;

363
    elem.resize(nr);

364
    return *this;

365
}

366


367
//---------------------------------------------------------------------------

368


369
template<typename Number>

370
void Matrix<Number>::swap(Matrix<Number>& x) {

371
    size_t tmp = nr; nr = x.nr; x.nr = tmp;

372
    tmp = nc; nc = x.nc; x.nc = tmp;

373
    elem.swap(x.elem);

374
}

375


376
//---------------------------------------------------------------------------

377


378
template<typename Number>

379
void Matrix<Number>::resize(size_t nr_rows, size_t nr_cols) {

380
    nc = nr_cols; //for adding new rows with the right length

381
    resize(nr_rows);

382
    resize_columns(nr_cols);

383
}

384


385
template<typename Number>

386
void Matrix<Number>::resize(size_t nr_rows) {

387
    if (nr_rows > elem.size()) {

388
        elem.resize(nr_rows, vector<Number>(nc));

389
    }

390
    nr = nr_rows;

391
}

392


393
template<typename Number>

394
void Matrix<Number>::resize_columns(size_t nr_cols) {

395
    for (size_t i=0; i<nr; i++) {

396
        elem[i].resize(nr_cols);

397
    }

398
    nc = nr_cols;

399
}

400


401
//---------------------------------------------------------------------------

402


403
template<typename Number>

404
vector<Number> Matrix<Number>::diagonal() const{

405
    assert(nr == nc); 

406
    vector<Number> diag(nr);

407
    for(size_t i=0; i<nr;i++){

408
        diag[i]=elem[i][i];

409
    }

410
    return diag;

411
}

412


413
//---------------------------------------------------------------------------

414


415
template<typename Number>

416
size_t Matrix<Number>::maximal_decimal_length() const{

417
    size_t i,maxim=0;

418
    vector<size_t> maxim_col;

419
    maxim_col=maximal_decimal_length_columnwise();

420
    for (i = 0; i <nr; i++)

421
        maxim=max(maxim,maxim_col[i]);

422
    return maxim;

423
}

424


425
//---------------------------------------------------------------------------

426


427
template<typename Number>

428
vector<size_t> Matrix<Number>::maximal_decimal_length_columnwise() const{

429
    size_t i,j=0;

430
    vector<size_t> maxim(nc,0);

431
    for (i = 0; i <nr; i++) {

432
        for (j = 0; j <nc; j++) {

433
            maxim[j]=max(maxim[j],decimal_length(elem[i][j]));

434
/*            if(elem[i][j]<0){

435
                if(elem[i][j]<neg_max[j])

436
                    neg_max[j]=elem[i][j];

437
                continue;

438
            }

439
            if(elem[i][j]>pos_max[j])

440
                pos_max[j]=elem[i][j];

441
*/

442
        }

443
    }

444
    /* for(size_t j=0;j<nc;++j)

445
        maxim[j]=max(decimal_length(neg_max[j]),decimal_length(pos_max[j])); */

446
    return maxim;

447
}

448


449
//---------------------------------------------------------------------------

450


451
template<typename Number>

452
void Matrix<Number>::append(const Matrix<Number>& M) {

453
    assert (nc == M.nc);

454
    elem.reserve(nr+M.nr);

455
    for (size_t i=0; i<M.nr; i++) {

456
        elem.push_back(M.elem[i]);

457
    }

458
    nr += M.nr;

459
}

460


461
//---------------------------------------------------------------------------

462


463
template<typename Number>

464
void Matrix<Number>::append(const vector<vector<Number> >& M) {

465
    if(M.size()==0)

466
        return;

467
    assert (nc == M[0].size());

468
    elem.reserve(nr+M.size());

469
    for (size_t i=0; i<M.size(); i++) {

470
        elem.push_back(M[i]);

471
    }

472
    nr += M.size();

473
}

474


475
//---------------------------------------------------------------------------

476


477
template<typename Number>

478
void Matrix<Number>::append(const vector<Number>& V) {

479
    assert (nc == V.size());

480
    elem.push_back(V);

481
    nr++;

482
}

483


484
//---------------------------------------------------------------------------

485


486
template<typename Number>

487
void Matrix<Number>::append_column(const vector<Number>& v) {

488
    assert (nr == v.size());

489
    for (size_t i=0; i<nr; i++) {

490
        elem[i].resize(nc+1);

491
        elem[i][nc] = v[i];

492
    }

493
    nc++;

494
}

495


496
//---------------------------------------------------------------------------

497


498
template<typename Number>

499
void Matrix<Number>::remove_row(const vector<Number>& row) {

500
    size_t tmp_nr = nr;

501
    for (size_t i = 1; i <= tmp_nr; ++i) {

502
        if (elem[tmp_nr-i] == row) {

503
            elem.erase(elem.begin()+(tmp_nr-i));

504
            nr--;

505
        }

506
    }

507
}

508


509
//---------------------------------------------------------------------------

510


511
template<typename Number>

512
void Matrix<Number>::remove_duplicate_and_zero_rows() {

513
    bool remove_some = false;

514
    vector<bool> key(nr, true);

515


516
    set<vector<Number> > SortedRows;

517
    SortedRows.insert( vector<Number>(nc,0) );

518
    typename set<vector<Number> >::iterator found;

519
    for (size_t i = 0; i<nr; i++) {

520
        found = SortedRows.find(elem[i]);

521
        if (found != SortedRows.end()) {

522
            key[i] = false;

523
            remove_some = true;

524
        }

525
        else

526
            SortedRows.insert(found,elem[i]);

527
    }

528


529
    if (remove_some) {

530
        *this = submatrix(key);

531
    }

532
}

533


534
//---------------------------------------------------------------------------

535


536
template<typename Number>

537
Matrix<Number> Matrix<Number>::add(const Matrix<Number>& A) const{

538
    assert (nr == A.nr);

539
    assert (nc == A.nc);

540
    

541
    Matrix<Number> B(nr,nc);

542
    size_t i,j;

543
    for(i=0; i<nr;i++){

544
        for(j=0; j<nc; j++){

545
            B.elem[i][j]=elem[i][j]+A.elem[i][j];

546
        }

547
    }

548
    return B;

549
}

550


551
//---------------------------------------------------------------------------

552


553
template<typename Number>

554
Matrix<Number> Matrix<Number>::multiplication(const Matrix<Number>& A) const{

555
    assert (nc == A.nr);

556


557
    Matrix<Number> B(nr,A.nc,0);  //initialized with 0

558
    size_t i,j,k;

559
    for(i=0; i<B.nr;i++){

560
        for(j=0; j<B.nc; j++){

561
            for(k=0; k<nc; k++){

562
                B.elem[i][j]=B.elem[i][j]+elem[i][k]*A.elem[k][j];

563
            }

564
        }

565
    }

566
    return B;

567
}

568


569
//---------------------------------------------------------------------------

570


571
template<typename Number>

572
Matrix<Number> Matrix<Number>::multiplication_cut(const Matrix<Number>& A, const size_t& c) const{

573
    assert (nc == A.nr);

574
    assert(c<= A.nc);

575


576
    Matrix<Number> B(nr,c,0);  //initialized with 0

577
    size_t i,j,k;

578
    for(i=0; i<B.nr;i++){

579
        for(j=0; j<c; j++){

580
            for(k=0; k<nc; k++){

581
                B.elem[i][j]=B.elem[i][j]+elem[i][k]*A.elem[k][j];

582
            }

583
        }

584
    }

585
    return B;

586
}

587


588


589
//---------------------------------------------------------------------------

590
/*

591
template<typename Number>

592
Matrix<Number> Matrix<Number>::multiplication(const Matrix<Number>& A, long m) const{

593
    assert (nc == A.nr);

594


595
    Matrix<Number> B(nr,A.nc,0);  //initialized with 0

596
    size_t i,j,k;

597
    for(i=0; i<B.nr;i++){

598
        for(j=0; j<B.nc; j++){

599
                for(k=0; k<nc; k++){

600
                B.elem[i][j]=(B.elem[i][j]+elem[i][k]*A.elem[k][j])%m;

601
                if (B.elem[i][j]<0) {

602
                    B.elem[i][j]=B.elem[i][j]+m;

603
                }

604
            }

605
        }

606
    }

607
    return B;

608
}

609
*/

610


611
//---------------------------------------------------------------------------

612


613
template<typename Number>

614
bool Matrix<Number>::equal(const Matrix<Number>& A) const{

615
    if ((nr!=A.nr)||(nc!=A.nc)){  return false; }

616
    size_t i,j;

617
    for (i=0; i < nr; i++) {

618
        for (j = 0; j < nc; j++) {

619
            if (elem[i][j]!=A.elem[i][j]) {

620
                return false;

621
            }

622
        }

623
    }

624
    return true;

625
}

626


627
//---------------------------------------------------------------------------

628
/*

629
template<typename Number>

630
bool Matrix<Number>::equal(const Matrix<Number>& A, long m) const{

631
    if ((nr!=A.nr)||(nc!=A.nc)){  return false; }

632
    size_t i,j;

633
    for (i=0; i < nr; i++) {

634
        for (j = 0; j < nc; j++) {

635
            if (((elem[i][j]-A.elem[i][j]) % m)!=0) {

636
                return false;

637
            }

638
        }

639
    }

640
    return true;

641
} */

642


643
//---------------------------------------------------------------------------

644


645
template<typename Number>

646
Matrix<Number> Matrix<Number>::transpose()const{

647
    Matrix<Number> B(nc,nr);

648
    size_t i,j;

649
    for(i=0; i<nr;i++){

650
        for(j=0; j<nc; j++){

651
            B.elem[j][i]=elem[i][j];

652
        }

653
    }

654
    return B;

655
}

656


657
//---------------------------------------------------------------------------

658


659
template<typename Number>

660
void Matrix<Number>::scalar_multiplication(const Number& scalar){

661
    size_t i,j;

662
    for(i=0; i<nr;i++){

663
        for(j=0; j<nc; j++){

664
            elem[i][j] *= scalar;

665
        }

666
    }

667
}

668


669
//---------------------------------------------------------------------------

670


671
template<typename Number>

672
void Matrix<Number>::scalar_division(const Number& scalar){

673
    size_t i,j;

674
    assert(scalar != 0);

675
    for(i=0; i<nr;i++){

676
        for(j=0; j<nc; j++){

677
            // assert (elem[i][j]%scalar == 0);

678
            elem[i][j] /= scalar;

679
        }

680
    }

681
}

682


683
//---------------------------------------------------------------------------

684


685
/*

686
template<typename Number>

687
void Matrix<Number>::reduction_modulo(const Number& modulo){

688
    size_t i,j;

689
    for(i=0; i<nr;i++){

690
        for(j=0; j<nc; j++){

691
            elem[i][j] %= modulo;

692
            if (elem[i][j] < 0) {

693
                elem[i][j] += modulo;

694
            }

695
        }

696
    }

697
}

698
*/

699


700


701
//---------------------------------------------------------------------------

702


703
template<typename Number>

704
void Matrix<Number>::simplify_rows() {

705
    // vector<Number> g(nr);

706
    for (size_t i = 0; i <nr; i++) {

707
        v_simplify(elem[i]);

708
    }

709
    // return g;

710
}

711


712
//---------------------------------------------------------------------------

713


714
template<typename Number>

715
Matrix<Number> Matrix<Number>::multiply_rows(const vector<Number>& m) const{  //row i is multiplied by m[i]

716
  Matrix M = Matrix(nr,nc);

717
  size_t i,j;

718
  for (i = 0; i<nr; i++) {

719
     for (j = 0; j<nc; j++) {

720
        M.elem[i][j] = elem[i][j]*m[i];

721
     }

722
  }

723
  return M;

724
}

725


726
//---------------------------------------------------------------------------

727


728
template<typename Number>

729
void Matrix<Number>::MxV(vector<Number>& result, const vector<Number>& v) const{

730
    assert (nc == v.size());

731
    result.resize(nr);

732
    for(size_t i=0; i<nr;i++){

733
        result[i]=v_scalar_product(elem[i],v);

734
    }

735
}

736


737
//---------------------------------------------------------------------------

738


739
template<typename Number>

740
vector<Number> Matrix<Number>::MxV(const vector<Number>& v) const{

741
    vector<Number> w(nr);

742
    MxV(w, v);

743
    return w;

744
}

745


746
//---------------------------------------------------------------------------

747


748
template<typename Number>

749
vector<Number> Matrix<Number>::VxM(const vector<Number>& v) const{

750
    assert (nr == v.size());

751
    vector<Number> w(nc,0);

752
    size_t i,j;

753
    for (i=0; i<nc; i++){

754
        for (j=0; j<nr; j++){

755
            w[i] += v[j]*elem[j][i];

756
        }

757
        if(!check_range(w[i]))

758
            break;

759
    }

760


761
        return w;

762
}

763


764
//---------------------------------------------------------------------------

765


766
template<typename Number>

767
vector<Number> Matrix<Number>::VxM_div(const vector<Number>& v, const Number& divisor, bool& success) const{

768
    assert (nr == v.size());

769
    vector<Number> w(nc,0);

770
    success=true;

771
    size_t i,j;

772
    for (i=0; i<nc; i++){

773
        for (j=0; j<nr; j++){

774
            w[i] += v[j]*elem[j][i];

775
        }

776
        if(!check_range(w[i])){

777
            success=false;

778
            break;

779
        }

780
    }

781


782
    if(success)      

783
        v_scalar_division(w,divisor);  

784
        

785
    return w;

786
}

787


788
//---------------------------------------------------------------------------

789


790
template<typename Number>

791
bool Matrix<Number>::is_diagonal() const{

792


793
    for(size_t i=0;i<nr;++i)

794
        for(size_t j=0;j<nc;++j)

795
            if(i!=j && elem[i][j]!=0)

796
                return false;

797
    return true;

798
}

799


800
//---------------------------------------------------------------------------

801


802
template<typename Number>

803
vector<long> Matrix<Number>::pivot(size_t corner){

804
    assert(corner < nc);

805
    assert(corner < nr);

806
    size_t i,j;

807
    Number help=0;

808
    vector<long> v(2,-1);

809


810
    for (i = corner; i < nr; i++) {

811
        for (j = corner; j < nc; j++) {

812
            if (elem[i][j]!=0) {

813
                if ((help==0)||(Iabs(elem[i][j])<help)) {

814
                    help=Iabs(elem[i][j]);

815
                    v[0]=i;

816
                    v[1]=j;

817
                    if (help == 1) return v;

818
                }

819
            }

820
        }

821
    }

822
    

823
    return v;

824
}

825


826
//---------------------------------------------------------------------------

827


828
template<typename Number>

829
long Matrix<Number>::pivot_column(size_t row,size_t col){

830
    assert(col < nc);

831
    assert(row < nr);

832
    size_t i;

833
    long j=-1;

834
    Number help=0;

835


836
    for (i = row; i < nr; i++) {

837
        if (elem[i][col]!=0) {

838
            if ((help==0)||(Iabs(elem[i][col])<help)) {

839
                help=Iabs(elem[i][col]);

840
                j=i;

841
                if (help == 1) return j;

842
            }

843
        }

844
    }

845


846
    return j;

847
}

848


849
//---------------------------------------------------------------------------

850


851
template<typename Number>

852
long Matrix<Number>::pivot_column(size_t col){

853
    return pivot_column(col,col);

854
}

855


856
//---------------------------------------------------------------------------

857


858
template<typename Number>

859
void Matrix<Number>::exchange_rows(const size_t& row1, const size_t& row2){

860
    if (row1 == row2) return;

861
    assert(row1 < nr);

862
    assert(row2 < nr);

863
    elem[row1].swap(elem[row2]);

864
}

865


866
//---------------------------------------------------------------------------

867


868
template<typename Number>

869
void Matrix<Number>::exchange_columns(const size_t& col1, const size_t& col2){

870
    if (col1 == col2) return;

871
    assert(col1 < nc);

872
    assert(col2 < nc);

873
    for(size_t i=0; i<nr;i++){

874
        std::swap(elem[i][col1], elem[i][col2]);

875
    }

876
}

877


878
//---------------------------------------------------------------------------

879
 

880
template<typename Number>

881
bool Matrix<Number>::reduce_row (size_t row, size_t col) {

882
    assert(col < nc);

883
    assert(row < nr);

884
    size_t i,j;

885
    Number help;

886
    for (i =row+1; i < nr; i++) {

887
        if (elem[i][col]!=0) {

888
            help=elem[i][col] / elem[row][col];

889
            for (j = col; j < nc; j++) {

890
                elem[i][j] -= help*elem[row][j];

891
                if (!check_range(elem[i][j]) ) {

892
                    return false;

893
                }

894
            }

895
            // v_el_trans<Number>(elem[row],elem[i],-help,col);

896
        }

897
    }

898
    return true;

899
}

900


901
//---------------------------------------------------------------------------

902


903
template<typename Number>

904
bool  Matrix<Number>::reduce_row (size_t corner) {

905
    return reduce_row(corner,corner);

906
}

907


908
//---------------------------------------------------------------------------

909
 

910
template<typename Number>

911
bool Matrix<Number>::reduce_rows_upwards () {

912
// assumes that "this" is in row echelon form

913
// and reduces eevery column in which the rank jumps 

914
// by its lowest element

915
//

916
// Aplies v_simplify to make rows "nice"

917
    if(nr==0)

918
        return true;

919


920
    for(size_t row=0;row<nr;++row){

921
        size_t col;

922
        for(col=0;col<nc;++col)

923
            if(elem[row][col]!=0)

924
                break;

925
        if(col==nc) // zero row

926
            continue;

927
        if(elem[row][col]<0)

928
            v_scalar_multiplication<Number>(elem[row],-1); // make corner posizive

929
        

930
        for(long i=row-1;i>=0;--i){

931
            Number quot;            

932
            //minimal_remainder(elem[i][col],elem[row][col],quot,rem);

933
            quot=elem[i][col]/elem[row][col];

934
            elem[i][col]=0; // rem

935
            for(size_t j=col+1;j<nc;++j){

936
                elem[i][j]-=quot* elem[row][j];

937
            }                                           

938
        }

939
    }

940
    

941
    simplify_rows();

942
           

943
    return true;

944
}

945


946
//---------------------------------------------------------------------------

947
 

948
template<typename Number>

949
bool Matrix<Number>::linear_comb_columns(const size_t& col,const size_t& j,

950
            const Number& u,const Number& w,const Number& v,const Number& z){

951
                       

952
    for(size_t i=0;i<nr;++i){

953
        Number rescue=elem[i][col];

954
        elem[i][col]=u*elem[i][col]+v*elem[i][j];

955
        elem[i][j]=w*rescue+z*elem[i][j];

956
        if ( (!check_range(elem[i][col])  || !check_range(elem[i][j]) )) {

957
            return false;

958
        }        

959
    }

960
    return true;

961
}

962


963
//---------------------------------------------------------------------------

964


965
template<typename Number>

966
bool Matrix<Number>::gcd_reduce_column (size_t corner, Matrix<Number>& Right){

967
    assert(corner < nc);

968
    assert(corner < nr);

969
    Number d,u,w,z,v;

970
    for(size_t j=corner+1;j<nc;++j){

971
       d =elem[corner][corner],elem[corner]; // ext_gcd(elem[corner][corner],elem[corner][j],u,v);

972
       u=1;

973
       v=0;

974
       w=-elem[corner][j]/d;

975
       z=elem[corner][corner]/d;

976
       // Now we multiply the submatrix formed by columns "corner" and "j" 

977
       // and rows corner,...,nr from the right by the 2x2 matrix

978
       // | u w |

979
       // | v z |              

980
       if(!linear_comb_columns(corner,j,u,w,v,z))

981
           return false; 

982
       if(!Right.linear_comb_columns(corner,j,u,w,v,z))

983
           return false;  

984
    }   

985
    return true;

986
}

987


988


989
//---------------------------------------------------------------------------

990


991
template<typename Number>

992
bool Matrix<Number>::column_trigonalize(size_t rk, Matrix<Number>& Right) { 

993
    assert(Right.nr == nc);

994
    assert(Right.nc == nc);

995
    vector<long> piv(2,0);       

996
    for(size_t j=0;j<rk;++j){

997
            piv=pivot(j);

998
            assert(piv[0]>=0); // protect against wrong rank

999
            exchange_rows (j,piv[0]);

1000
            exchange_columns (j,piv[1]);

1001
            Right.exchange_columns(j,piv[1]);

1002
            if(!gcd_reduce_column(j, Right))

1003
                return false;

1004
    }

1005
    return true;

1006
}

1007


1008
//---------------------------------------------------------------------------

1009


1010
template<typename Number>

1011
Number Matrix<Number>::compute_vol(bool& success){

1012
        

1013
    assert(nr<=nc);

1014
    

1015
    Number det=1;

1016
    for(size_t i=0;i<nr;++i){

1017
        det*=elem[i][i]; 

1018
        if(!check_range(det)){

1019
            success=false;

1020
            return 0;

1021
        }

1022
    }

1023
            

1024
    det=Iabs(det);

1025
    success=true;

1026
    return det;

1027
}

1028


1029
//---------------------------------------------------------------------------

1030


1031
template<typename Number>

1032
size_t Matrix<Number>::row_echelon_inner_elem(bool& success){

1033


1034
    size_t pc=0;

1035
    long piv=0, rk=0;

1036
    success=true;

1037


1038
    if(nr==0)

1039
        return 0;

1040
    

1041
    for (rk = 0; rk < (long) nr; rk++){

1042
        for(;pc<nc;pc++){

1043
            piv=pivot_column(rk,pc);

1044
            if(piv>=0)

1045
                break;

1046
        }

1047
        if(pc==nc)

1048
            break;

1049
        do{

1050
            exchange_rows (rk,piv);

1051
            if(!reduce_row(rk,pc)){

1052
                success=false;

1053
                return rk;

1054
            }

1055
            piv=pivot_column(rk,pc);

1056
        }while (piv>rk);

1057
    }

1058
                

1059
    return rk;

1060
}

1061


1062
//---------------------------------------------------------------------------

1063


1064
/*

1065
template<typename Number>

1066
size_t Matrix<Number>::row_echelon_inner_bareiss(bool& success, Number& det){

1067
// no overflow checks since this is supposed to be only used with GMP

1068


1069
    success=true;

1070
    if(nr==0)

1071
        return 0;

1072
    assert(using_GMP<Number>());

1073


1074
    size_t pc=0;

1075
    long piv=0, rk=0;

1076
    vector<bool> last_time_mult(nr,false),this_time_mult(nr,false);

1077
    Number last_div=1,this_div=1;

1078
    size_t this_time_exp=0,last_time_exp=0;

1079
    Number det_factor=1;

1080
    

1081
    for (rk = 0; rk < (long) nr; rk++){

1082


1083
        for(;pc<nc;pc++){

1084
            piv=pivot_column(rk,pc);

1085
            if(piv>=0)

1086
                break;

1087
        }

1088
        if(pc==nc)

1089
            break;

1090
                        

1091
        if(!last_time_mult[piv]){

1092
            for(size_t k=rk;k<nr;++k)

1093
                if(elem[k][pc]!=0 && last_time_mult[k]){

1094
                    piv=k;

1095
                    break;                

1096
                }        

1097
        }        

1098
        

1099
        exchange_rows (rk,piv);

1100
        v_bool_entry_swap(last_time_mult,rk,piv);

1101
        

1102
        if(!last_time_mult[rk])

1103
            for(size_t i=0;i<nr;++i)

1104
                last_time_mult[i]=false;

1105
                     

1106
        Number a=elem[rk][pc];

1107
        this_div=Iabs(a);

1108
        this_time_exp=0;

1109
        

1110
        for(size_t i=rk+1;i<nr;++i){

1111
            if(elem[i][pc]==0){

1112
                this_time_mult[i]=false;

1113
                continue;

1114
            }

1115
            

1116
            this_time_exp++;

1117
            this_time_mult[i]=true;

1118
            bool divide=last_time_mult[i] && (last_div!=1);

1119
            if(divide)

1120
                last_time_exp--;

1121
            Number b=elem[i][pc];

1122
            elem[i][pc]=0;

1123
            if(a==1){

1124
                for(size_t j=pc+1;j<nc;++j){

1125
                    elem[i][j]=elem[i][j]-b*elem[rk][j]; 

1126
                    if(divide){

1127
                        elem[i][j]/=last_div;

1128
                    }

1129
                }            

1130
            }

1131
            else{

1132
                if(a==-1){

1133
                    for(size_t j=pc+1;j<nc;++j){

1134
                        elem[i][j]=-elem[i][j]-b*elem[rk][j]; 

1135
                        if(divide){

1136
                            elem[i][j]/=last_div;

1137
                        }            

1138
                    }

1139
                }

1140
                else{            

1141
                    for(size_t j=pc+1;j<nc;++j){

1142
                        elem[i][j]=a*elem[i][j]-b*elem[rk][j]; 

1143
                       if(divide){

1144
                            elem[i][j]/=last_div;

1145
                        }                   

1146
                    }

1147
                }

1148
            }

1149
        }

1150
        

1151
        for(size_t i=0;i<last_time_exp;++i)

1152
            det_factor*=last_div;

1153
        last_time_mult=this_time_mult;

1154
        last_div=this_div;

1155
        last_time_exp=this_time_exp;

1156
    }

1157
    

1158
    det=0;

1159
    if (nr <= nc && rk == (long) nr) { // must allow nonsquare matrices

1160
        det=1;

1161
        for(size_t i=0;i<nr;++i)

1162
            det*=elem[i][i];            

1163
        det=Iabs<Number>(det/det_factor);        

1164
    }

1165
    

1166
    return rk;

1167
}

1168
*/

1169


1170
//---------------------------------------------------------------------------

1171


1172
template<typename Number>

1173
size_t Matrix<Number>::row_echelon_reduce(bool& success){

1174


1175
    size_t rk=row_echelon_inner_elem(success);

1176
    if(success)

1177
        success=reduce_rows_upwards();

1178
    return rk;

1179
}

1180


1181
//---------------------------------------------------------------------------

1182


1183
template<typename Number>

1184
Number Matrix<Number>::full_rank_index(bool& success){

1185


1186
    size_t rk=row_echelon_inner_elem(success);

1187
    Number index=1;

1188
    if(success){

1189
        for(size_t i=0;i<rk;++i){

1190
            index*=elem[i][i];

1191
            if(!check_range(index)){

1192
                success=false;

1193
                index=0;

1194
                return index;

1195
            }

1196
        }

1197
    }

1198
    assert(rk==nc); // must have full rank

1199
    index=Iabs(index);

1200
    return index;

1201
}

1202
//---------------------------------------------------------------------------

1203


1204
template<typename Number>

1205
Matrix<Number> Matrix<Number>::row_column_trigonalize(size_t& rk, bool& success) {

1206


1207
    Matrix<Number> Right(nc);

1208
    rk=row_echelon_reduce(success);

1209
    if(success)

1210
        success=column_trigonalize(rk,Right); 

1211
    return Right; 

1212
} 

1213


1214
//---------------------------------------------------------------------------

1215


1216
template<typename Number>

1217
size_t Matrix<Number>::row_echelon(bool& success, bool do_compute_vol, Number& det){

1218
    

1219
/*    if(using_GMP<Number>()){

1220
        return row_echelon_inner_bareiss(success,det);;

1221
    }

1222
    else{ */

1223
        size_t rk=row_echelon_inner_elem(success);

1224
        if(do_compute_vol)

1225
            det=compute_vol(success);

1226
        return rk;

1227
//    }

1228
}

1229


1230
//---------------------------------------------------------------------------

1231


1232
template<typename Number>

1233
size_t Matrix<Number>::row_echelon(bool& success){

1234
    

1235
    Number dummy;

1236
    return row_echelon(success,false,dummy);

1237
}

1238


1239
//---------------------------------------------------------------------------

1240


1241
template<typename Number>

1242
size_t Matrix<Number>::row_echelon(bool& success, Number& det){

1243
    

1244
    return row_echelon(success,true,det);

1245
}

1246


1247


1248


1249
//---------------------------------------------------------------------------

1250


1251
template<typename Number>

1252
size_t Matrix<Number>::rank_submatrix(const Matrix<Number>& mother, const vector<key_t>& key){

1253


1254
    assert(nc>=mother.nc);

1255
    if(nr<key.size()){

1256
        elem.resize(key.size(),vector<Number>(nc,0));

1257
        nr=key.size();    

1258
    }

1259
    size_t save_nr=nr;

1260
    size_t save_nc=nc;

1261
    nr=key.size();

1262
    nc=mother.nc;

1263


1264
    select_submatrix(mother,key);

1265


1266
    bool success;

1267
    size_t rk=row_echelon(success);

1268
    

1269
    nr=save_nr;

1270
    nc=save_nc;

1271
    return rk;                               

1272
}

1273


1274
//---------------------------------------------------------------------------

1275


1276
template<typename Number>

1277
size_t Matrix<Number>::rank_submatrix(const vector<key_t>& key) const{

1278


1279
    Matrix<Number> work(key.size(),nc);

1280
    return work.rank_submatrix(*this,key);              

1281
}

1282


1283
//---------------------------------------------------------------------------

1284


1285
template<typename Number>

1286
size_t Matrix<Number>::rank() const{

1287
    vector<key_t> key(nr);

1288
    for(size_t i=0;i<nr;++i)

1289
        key[i]=i;

1290
    return rank_submatrix(key);

1291
}

1292


1293
//---------------------------------------------------------------------------

1294


1295
template<typename Number>

1296
Number Matrix<Number>::vol_submatrix(const Matrix<Number>& mother, const vector<key_t>& key){

1297


1298
    assert(nc>=mother.nc);

1299
    if(nr<key.size()){

1300
        elem.resize(key.size(),vector<Number>(nc,0));

1301
        nr=key.size();    

1302
    }

1303
    size_t save_nr=nr;

1304
    size_t save_nc=nc;

1305
    nr=key.size();

1306
    nc=mother.nc;

1307
    

1308
    select_submatrix(mother,key);

1309


1310
    bool success;

1311
    Number det;

1312
    row_echelon(success,det);

1313
    

1314
    nr=save_nr;

1315
    nc=save_nc;

1316
    return det;                               

1317
}

1318
//---------------------------------------------------------------------------

1319


1320
template<typename Number>

1321
Number Matrix<Number>::vol_submatrix(const vector<key_t>& key) const{

1322


1323
    Matrix<Number> work(key.size(),nc);

1324
    return work.vol_submatrix(*this,key);              

1325
}

1326


1327
//---------------------------------------------------------------------------

1328


1329
template<typename Number>

1330
Number Matrix<Number>::vol() const{

1331
    vector<key_t> key(nr);

1332
    for(size_t i=0;i<nr;++i)

1333
        key[i]=i;

1334
    return vol_submatrix(key);

1335
}

1336


1337
//---------------------------------------------------------------------------

1338


1339
template<typename Number>

1340
vector<key_t>  Matrix<Number>::max_rank_submatrix_lex_inner(bool& success) const{

1341


1342
    success=true;

1343
    size_t max_rank=min(nr,nc);

1344
    Matrix<Number> Test(max_rank,nc);

1345
    Test.nr=0;

1346
    vector<key_t> col;

1347
    col.reserve(max_rank);

1348
    vector<key_t> key;

1349
    key.reserve(max_rank);

1350
    size_t rk=0;

1351
    

1352
    vector<vector<bool> > col_done(max_rank,vector<bool>(nc,false));

1353
    

1354
    vector<Number> Test_vec(nc);

1355
     

1356
    for(size_t i=0;i<nr;++i){    

1357
        Test_vec=elem[i];            

1358
        for(size_t k=0;k<rk;++k){

1359
            if(Test_vec[col[k]]==0)

1360
                continue;

1361
            Number a=Test[k][col[k]];

1362
            Number b=Test_vec[col[k]];

1363
            for(size_t j=0;j<nc;++j)

1364
                if(!col_done[k][j]){

1365
                Test_vec[j]=a*Test_vec[j]-b*Test[k][j];

1366
                if (!check_range(Test_vec[j]) ) {

1367
                    success=false;

1368
                    return key;

1369
                }

1370
            }

1371
        }

1372
        

1373
        size_t j=0;

1374
        for(;j<nc;++j)

1375
            if(Test_vec[j]!=0)

1376
                break;

1377
        if(j==nc)     // Test_vec=0

1378
            continue;

1379
            

1380
        col.push_back(j);

1381
        key.push_back(i);

1382
        

1383
        if(rk>0){

1384
            col_done[rk]=col_done[rk-1];

1385
            col_done[rk][col[rk-1]]=true;

1386
        }

1387


1388
        Test.nr++;

1389
        rk++;

1390
        v_simplify(Test_vec);

1391
        Test[rk-1]=Test_vec;

1392
            

1393
        if(rk==max_rank)

1394
            break;   

1395
    }    

1396
    return key;                

1397
}

1398


1399
//---------------------------------------------------------------------------

1400


1401
template<typename Number>

1402
vector<key_t>  Matrix<Number>::max_rank_submatrix_lex() const{

1403
    bool success;

1404
    vector<key_t> key=max_rank_submatrix_lex_inner(success);

1405
    return key;

1406
}

1407


1408
//---------------------------------------------------------------------------

1409


1410
template<typename Number>

1411
bool Matrix<Number>::solve_destructive_inner(bool ZZinvertible,Number& denom) {

1412


1413
    assert(nc>=nr);

1414
    size_t dim=nr;

1415
    bool success;

1416
    

1417
    size_t rk;

1418
    

1419
    if(ZZinvertible){

1420
        rk=row_echelon_inner_elem(success); 

1421
        if(!success)

1422
            return false;        

1423
        assert(rk==nr);

1424
        denom=compute_vol(success);

1425
    }

1426
    else{

1427
        rk=row_echelon(success,denom);

1428
        if(!success)

1429
            return false;    

1430
    }

1431


1432
    if (denom==0) { 

1433
        if(using_GMP<Number>()){

1434
            errorOutput() << "Cannot solve system (denom=0)!" << endl;

1435
            throw ArithmeticException();

1436
        }

1437
        else

1438
            return false;            

1439
    }

1440


1441
    Number S;

1442
    size_t i;

1443
    long j;

1444
    size_t k;

1445
    for (i = nr; i < nc; i++) {

1446
        for (j = dim-1; j >= 0; j--) {

1447
            S=denom*elem[j][i];

1448
            for (k = j+1; k < dim; k++) {

1449
                S-=elem[j][k]*elem[k][i];

1450
            }

1451
            if(!check_range(S))

1452
                return false;

1453
            elem[j][i]=S/elem[j][j];

1454
        }

1455
    }

1456
    return true;

1457
}

1458


1459
//---------------------------------------------------------------------------

1460


1461
template<typename Number>

1462
void Matrix<Number>::customize_solution(size_t dim, Number& denom, size_t red_col, 

1463
                     size_t sign_col, bool make_sol_prime) {

1464
    

1465
    return;

1466
                         

1467
 /*   assert(!(make_sol_prime && (sign_col>0 || red_col>0)));

1468


1469
    for(size_t j=0;j<red_col;++j){  // reduce first red_col columns of solution mod denom

1470
        for(size_t k=0;k<dim;++k){

1471
          elem[k][dim+j]%=denom;

1472
          if(elem[k][dim+j]<0)

1473
              elem[k][dim+j]+=Iabs(denom);

1474
       }

1475
    }

1476
        

1477
    for(size_t j=0;j<sign_col;++j)   // replace entries in the next sign_col columns by their signs

1478
      for(size_t k=0;k<dim;++k){

1479
        if(elem[k][dim+red_col+j]>0){

1480
            elem[k][dim+red_col+j]=1;

1481
            continue;

1482
        } 

1483
        if(elem[k][dim+red_col+j]<0){

1484
            elem[k][dim+red_col+j]=-1;

1485
            continue;

1486
        }       

1487
      }

1488
      

1489
    if(make_sol_prime) // make columns of solution coprime if wanted

1490
        make_cols_prime(dim,nc-1); */

1491
}

1492


1493
//---------------------------------------------------------------------------

1494


1495
template<typename Number>

1496
void Matrix<Number>::solve_system_submatrix_outer(const Matrix<Number>& mother, const vector<key_t>& key, const vector<vector<Number>* >& RS,

1497
        Number& denom, bool ZZ_invertible, bool transpose, size_t red_col, size_t sign_col, 

1498
        bool compute_denom, bool make_sol_prime) {

1499
     

1500
    size_t dim=mother.nc;

1501
    assert(key.size()==dim);

1502
    assert(nr==dim);

1503
    assert(dim+RS.size()<=nc);

1504
    size_t save_nc=nc;

1505
    nc=dim+RS.size();

1506
    

1507
    if(transpose)

1508
       select_submatrix_trans(mother,key);           

1509
    else

1510
       select_submatrix(mother,key);

1511
               

1512
    for(size_t i=0;i<dim;++i)

1513
       for(size_t k=0;k<RS.size();++k)

1514
           elem[i][k+dim]= (*RS[k])[i];

1515
    

1516
    if(solve_destructive_inner(ZZ_invertible,denom)){

1517
        customize_solution(dim, denom,red_col,sign_col,make_sol_prime);        

1518
    } /*

1519
    else{          

1520
       #pragma omp atomic

1521
       GMP_mat++;

1522
    

1523
       Matrix<mpz_class> mpz_this(nr,nc);

1524
       mpz_class mpz_denom;

1525
       if(transpose)

1526
           mpz_submatrix_trans(mpz_this,mother,key);

1527
       else            

1528
           mpz_submatrix(mpz_this,mother,key);

1529
           

1530
       for(size_t i=0;i<dim;++i)

1531
           for(size_t k=0;k<RS.size();++k)

1532
               convert(mpz_this[i][k+dim], (*RS[k])[i]);

1533
       mpz_this.solve_destructive_inner(ZZ_invertible,mpz_denom);

1534
       mpz_this.customize_solution(dim, mpz_denom,red_col,sign_col,make_sol_prime);           

1535
          

1536
       for(size_t i=0;i<dim;++i)  // replace left side by 0, except diagonal if ZZ_invetible

1537
          for(size_t j=0;j<dim;++j){

1538
            if(i!=j || !ZZ_invertible)

1539
                mpz_this[i][j]=0;              

1540
          }

1541
              

1542
       mat_to_Int(mpz_this,*this);

1543
       if(compute_denom)

1544
           convert(denom, mpz_denom);                

1545
    }*/

1546
    nc=save_nc;      

1547
}

1548


1549


1550
//---------------------------------------------------------------------------

1551


1552


1553
template<typename Number>

1554
void Matrix<Number>::solve_system_submatrix(const Matrix<Number>& mother, const vector<key_t>& key, const vector<vector<Number>* >& RS,

1555
         vector< Number >& diagonal, Number& denom, size_t red_col, size_t sign_col) {

1556


1557
    solve_system_submatrix_outer(mother,key,RS,denom,true,false,red_col,sign_col);

1558
    assert(diagonal.size()==nr);

1559
    for(size_t i=0;i<nr;++i)

1560
        diagonal[i]=elem[i][i];

1561
                 

1562
}

1563


1564


1565


1566
//---------------------------------------------------------------------------

1567
// the same without diagonal

1568
template<typename Number>

1569
void Matrix<Number>::solve_system_submatrix(const Matrix<Number>& mother, const vector<key_t>& key, const vector<vector<Number>* >& RS,

1570
         Number& denom, size_t red_col, size_t sign_col, bool compute_denom, bool make_sol_prime) {

1571


1572
    solve_system_submatrix_outer(mother,key,RS,denom,false,false,red_col,sign_col, 

1573
                compute_denom, make_sol_prime);

1574
}

1575


1576
//---------------------------------------------------------------------------

1577


1578
template<typename Number>

1579
void Matrix<Number>::solve_system_submatrix_trans(const Matrix<Number>& mother, const vector<key_t>& key, const vector<vector<Number>* >& RS,

1580
         Number& denom, size_t red_col, size_t sign_col) {

1581
         

1582
    solve_system_submatrix_outer(mother,key,RS,denom,false,true,red_col,sign_col);

1583
}

1584


1585
//---------------------------------------------------------------------------

1586


1587
template<typename Number>

1588
Matrix<Number> Matrix<Number>::extract_solution() const {

1589
    assert(nc>=nr);

1590
    Matrix<Number> Solution(nr,nc-nr); 

1591
    for(size_t i=0;i<nr;++i){

1592
        for(size_t j=0;j<Solution.nc;++j)

1593
            Solution[i][j]=elem[i][j+nr];    

1594
    }

1595
    return Solution;  

1596
}

1597


1598
//---------------------------------------------------------------------------

1599


1600
template<typename Number>

1601
vector<vector<Number>* > Matrix<Number>::row_pointers(){

1602


1603
    vector<vector<Number>* > pointers(nr);

1604
    for(size_t i=0;i<nr;++i)

1605
        pointers[i]=&(elem[i]);

1606
    return pointers;

1607
}

1608


1609
//---------------------------------------------------------------------------

1610


1611
template<typename Number>

1612
vector<vector<Number>* > Matrix<Number>::submatrix_pointers(const vector<key_t>& key){

1613


1614
    vector<vector<Number>* > pointers(key.size());

1615
    for(size_t i=0;i<key.size();++i)

1616
        pointers[i]=&(elem[key[i]]);

1617
    return pointers;

1618
}

1619
//---------------------------------------------------------------------------

1620


1621
template<typename Number>

1622
Matrix<Number> Matrix<Number>::solve(const Matrix<Number>& Right_side,vector<Number>& diagonal,Number& denom) const {

1623


1624
    Matrix<Number> M(nr,nc+Right_side.nc);

1625
    vector<key_t> key=identity_key(nr);

1626
    Matrix<Number> RS_trans=Right_side.transpose();

1627
    vector<vector<Number>* > RS=RS_trans.row_pointers();

1628
    M.solve_system_submatrix(*this,key,RS,diagonal,denom,0,0);

1629
    return M.extract_solution(); 

1630
}

1631


1632
//---------------------------------------------------------------------------

1633


1634
template<typename Number>

1635
Matrix<Number> Matrix<Number>::solve(const Matrix<Number>& Right_side, Number& denom) const {

1636


1637
    Matrix<Number> M(nr,nc+Right_side.nc);

1638
    vector<key_t> key=identity_key(nr);

1639
    Matrix<Number> RS_trans=Right_side.transpose();

1640
    vector<vector<Number>* > RS=RS_trans.row_pointers();

1641
    M.solve_system_submatrix(*this,key,RS,denom,0,0);

1642
    return M.extract_solution(); 

1643
}

1644


1645
//---------------------------------------------------------------------------

1646


1647
template<typename Number>

1648
Matrix<Number> Matrix<Number>::invert(Number& denom) const{

1649
    assert(nr == nc);

1650
    Matrix<Number> Right_side(nr);

1651


1652
    return solve(Right_side,denom);

1653
}

1654


1655
//---------------------------------------------------------------------------

1656


1657
template<typename Number>

1658
Matrix<Number> Matrix<Number>::bundle_matrices(const Matrix<Number>& Right_side) const {

1659


1660
    assert(nr == nc);

1661
    assert(nc == Right_side.nr);

1662
    Matrix<Number> M(nr,nc+Right_side.nc);

1663
    for(size_t i=0;i<nr;++i){

1664
        for(size_t j=0;j<nc;++j)

1665
            M[i][j]=elem[i][j];

1666
        for(size_t j=nc;j<M.nc;++j)

1667
            M[i][j]=Right_side[i][j-nc];

1668
    }

1669
    return M;

1670
}

1671
//---------------------------------------------------------------------------

1672


1673
template<typename Number>

1674
Matrix<Number> Matrix<Number>::invert_unprotected(Number& denom, bool& success) const{

1675
    assert(nr == nc);

1676
    Matrix<Number> Right_side(nr);

1677
    Matrix<Number> M=bundle_matrices(Right_side);

1678
    success=M.solve_destructive_inner(false,denom);

1679
    return M.extract_solution();;

1680
}

1681


1682
//---------------------------------------------------------------------------

1683


1684
template<typename Number>

1685
void Matrix<Number>::invert_submatrix(const vector<key_t>& key, Number& denom, Matrix<Number>& Inv, bool compute_denom, bool make_sol_prime) const{

1686
    assert(key.size() == nc);

1687
    Matrix<Number> unit_mat(key.size());

1688
    Matrix<Number> M(key.size(),2*key.size());        

1689
    vector<vector<Number>* > RS_pointers=unit_mat.row_pointers();

1690
    M.solve_system_submatrix(*this,key,RS_pointers,denom,0,0, compute_denom, make_sol_prime);

1691
    Inv=M.extract_solution();;

1692
}

1693


1694
//---------------------------------------------------------------------------

1695


1696
template<typename Number>

1697
void Matrix<Number>::simplex_data(const vector<key_t>& key, Matrix<Number>& Supp, Number& vol, bool compute_vol) const{

1698
    assert(key.size() == nc);

1699
    invert_submatrix(key,vol,Supp,compute_vol,true);

1700
    Supp=Supp.transpose();

1701
    // Supp.make_prime(); now done internally -- but not in Q !! Therefore

1702
    Supp.simplify_rows();

1703
}

1704
//---------------------------------------------------------------------------

1705


1706
template<typename Number>

1707
vector<Number> Matrix<Number>::solve_rectangular(const vector<Number>& v, Number& denom) const {

1708
    if (nc == 0 || nr == 0) { //return zero-vector as solution

1709
        return vector<Number>(nc,0);

1710
    }

1711
    size_t i;

1712
    vector<key_t>  rows=max_rank_submatrix_lex();

1713
    Matrix<Number> Left_Side=submatrix(rows);

1714
    assert(nc == Left_Side.nr); //otherwise input hadn't full rank //TODO 

1715
    Matrix<Number> Right_Side(v.size(),1);

1716
    Right_Side.write_column(0,v);

1717
    Right_Side = Right_Side.submatrix(rows);

1718
    Matrix<Number> Solution=Left_Side.solve(Right_Side, denom);

1719
    vector<Number> Linear_Form(nc);

1720
    for (i = 0; i <nc; i++) {

1721
        Linear_Form[i] = Solution[i][0];  // the solution vector is called Linear_Form

1722
    }

1723
    vector<Number> test = MxV(Linear_Form); // we have solved the system by taking a square submatrix

1724
                        // now we must test whether the solution satisfies the full system

1725
    for (i = 0; i <nr; i++) {

1726
        if (test[i] != denom * v[i]){

1727
            return vector<Number>();

1728
        }

1729
    }

1730
    Number total_gcd = 1; // libnormaliz::gcd(denom,v_gcd(Linear_Form)); // extract the gcd of denom and solution

1731
    denom/=total_gcd;

1732
    v_scalar_division(Linear_Form,total_gcd);

1733
    return Linear_Form;

1734
}

1735
//---------------------------------------------------------------------------

1736


1737
template<typename Number>

1738
vector<Number> Matrix<Number>::solve_ZZ(const vector<Number>& v) const {

1739


1740
    Number denom;

1741
    vector<Number> result=solve_rectangular(v,denom);

1742
    if(denom!=1)

1743
        result.clear();

1744
    return result;

1745
}

1746
//---------------------------------------------------------------------------

1747


1748
template<typename Number>

1749
vector<Number> Matrix<Number>::find_linear_form() const {

1750


1751
    Number denom;

1752
    vector<Number> result=solve_rectangular(vector<Number>(nr,1),denom);

1753
    return result;

1754
}

1755


1756
//---------------------------------------------------------------------------

1757


1758
template<typename Number>

1759
vector<Number> Matrix<Number>::find_linear_form_low_dim () const{

1760
    size_t rank=(*this).rank();

1761
    if (rank == 0) { //return zero-vector as linear form

1762
        return vector<Number>(nc,0);

1763
    }

1764
    if (rank == nc) { // basis change not necessary

1765
        return (*this).find_linear_form();

1766
    }

1767


1768
    Sublattice_Representation<Number> Basis_Change(*this,true);

1769
    vector<Number> Linear_Form=Basis_Change.to_sublattice(*this).find_linear_form();

1770
    if(Linear_Form.size()!=0)

1771
        Linear_Form=Basis_Change.from_sublattice_dual(Linear_Form);

1772


1773
    return Linear_Form;

1774
}

1775


1776
//---------------------------------------------------------------------------

1777


1778
template<typename Number>

1779
size_t Matrix<Number>::row_echelon_reduce(){

1780


1781
    size_t rk;

1782
    Matrix<Number> Copy(*this);

1783
    bool success;

1784
    rk=row_echelon_reduce(success);

1785


1786
        Shrink_nr_rows(rk);

1787
        return rk;

1788


1789
}

1790
//---------------------------------------------------------------------------

1791


1792
template<typename Number>

1793
Number Matrix<Number>::full_rank_index() const{

1794
    

1795
    Matrix<Number> Copy(*this);

1796
    Number index;

1797
    bool success;

1798
    index=Copy.full_rank_index(success);

1799


1800
        return index;

1801
}

1802


1803
//---------------------------------------------------------------------------

1804


1805
template<typename Number>

1806
size_t Matrix<Number>::row_echelon(){

1807
    

1808
    Matrix<Number> Copy(*this);

1809
    bool success;

1810
    size_t rk;

1811
    rk=row_echelon(success);

1812


1813
        Shrink_nr_rows(rk);

1814
        return rk;

1815


1816
}

1817


1818
//---------------------------------------------------------------------------

1819


1820
template<typename Number>

1821
Matrix<Number> Matrix<Number>::kernel () const{

1822
// computes a ZZ-basis of the solutions of (*this)x=0

1823
// the basis is formed by the rOWS of the returned matrix

1824


1825
    size_t dim=nc;

1826
    if(nr==0)

1827
        return(Matrix<Number>(dim));

1828


1829
    Matrix<Number> Copy(*this);

1830
    size_t rank;

1831
    bool success;

1832
    Matrix<Number> Transf=Copy.row_column_trigonalize(rank,success);

1833
    

1834
    Matrix<Number> ker_basis(dim-rank,dim);

1835
    Matrix<Number> Help =Transf.transpose();

1836
    for (size_t i = rank; i < dim; i++) 

1837
            ker_basis[i-rank]=Help[i];

1838
    ker_basis.row_echelon_reduce();

1839
    return(ker_basis);

1840
}

1841


1842
//---------------------------------------------------------------------------

1843
// Converts "this" into (column almost) Hermite normal form, returns column transformation matrix

1844
/*template<typename Number>

1845
Matrix<Number> Matrix<Number>::AlmostHermite(size_t& rk){

1846


1847
    Matrix<Number> Copy=*this;

1848
    Matrix<Number> Transf;

1849
    bool success;

1850
    Transf=row_column_trigonalize(rk,success);

1851


1852
        return Transf;

1853
} */

1854


1855
//---------------------------------------------------------------------------

1856
// Classless conversion routines

1857
//---------------------------------------------------------------------------

1858


1859
template<typename ToType, typename FromType>

1860
void convert(Matrix<ToType>& to_mat, const Matrix<FromType>& from_mat){

1861
    size_t nrows = from_mat.nr_of_rows();

1862
    size_t ncols = from_mat.nr_of_columns();

1863
    to_mat.resize(nrows, ncols);

1864
    for(size_t i=0; i<nrows; ++i)

1865
        for(size_t j=0; j<ncols; ++j)

1866
            convert(to_mat[i][j], from_mat[i][j]);

1867
}

1868


1869
//---------------------------------------------------------------------------

1870


1871


1872
template<typename Number>

1873
bool weight_lex(const order_helper<Number>& a, const order_helper<Number>& b){

1874
    

1875
        if(a.weight < b.weight)

1876
            return true;

1877
        if(a.weight==b.weight)

1878
            if(*(a.v)< *(b.v))

1879
                return true;

1880
        return false;

1881
}

1882


1883
//---------------------------------------------------------------------------

1884


1885
template<typename Number>

1886
void Matrix<Number>::order_rows_by_perm(const vector<key_t>& perm){

1887
    order_by_perm(elem,perm);    

1888
}

1889


1890
template<typename Number>

1891
Matrix<Number>& Matrix<Number>::sort_by_weights(const Matrix<Number>& Weights, vector<bool> absolute){

1892
    if(nr<=1)

1893
        return *this;

1894
    vector<key_t> perm=perm_by_weights(Weights,absolute);

1895
    order_by_perm(elem,perm);

1896
    return *this;   

1897
}

1898


1899
template<typename Number>

1900
Matrix<Number>& Matrix<Number>::sort_lex(){

1901
    if(nr<=1)

1902
        return *this;

1903
    vector<key_t> perm=perm_by_weights(Matrix<Number>(0,nc),vector<bool>(0));

1904
    order_by_perm(elem,perm);

1905
    return *this;    

1906
}

1907


1908
template<typename Number>

1909
vector<key_t> Matrix<Number>::perm_by_weights(const Matrix<Number>& Weights, vector<bool> absolute){

1910
// the smallest entry is the row with index perm[0], then perm[1] etc.

1911
    

1912
    assert(Weights.nc==nc);

1913
    assert(absolute.size()==Weights.nr);

1914


1915
    list<order_helper<Number> > order;

1916
    order_helper<Number> entry;

1917
    entry.weight.resize(Weights.nr);

1918
    

1919
    for(key_t i=0;i<nr; ++i){

1920
        for(size_t j=0;j<Weights.nr;++j){

1921
            if(absolute[j])

1922
                entry.weight[j]=v_scalar_product(Weights[j],v_abs_value(elem[i]));

1923
            else

1924
                entry.weight[j]=v_scalar_product(Weights[j],elem[i]);                

1925
        }

1926
        entry.index=i;

1927
        entry.v=&(elem[i]);

1928
        order.push_back(entry);        

1929
    }

1930
    order.sort(weight_lex<Number>);

1931
    vector<key_t> perm(nr);

1932
    typename list<order_helper<Number> >::const_iterator ord=order.begin();

1933
    for(key_t i=0;i<nr;++i, ++ord)

1934
        perm[i]=ord->index; 

1935
    

1936
    return perm;

1937
}

1938


1939
//---------------------------------------------------

1940


1941
/* template<typename Number>

1942
Matrix<Number> Matrix<Number>::solve_congruences(bool& zero_modulus) const{

1943
 

1944
    

1945
    zero_modulus=false;

1946
    size_t i,j;

1947
    size_t nr_cong=nr, dim=nc-1;

1948
    if(nr_cong==0)

1949
        return Matrix<Number>(dim); // give back unit matrix

1950
    

1951
    //add slack variables to convert congruences into equaitions

1952
    Matrix<Number> Cong_Slack(nr_cong, dim+nr_cong);

1953
    for (i = 0; i < nr_cong; i++) {

1954
        for (j = 0; j < dim; j++) {

1955
            Cong_Slack[i][j]=elem[i][j];

1956
        }

1957
        Cong_Slack[i][dim+i]=elem[i][dim];

1958
        if(elem[i][dim]==0){

1959
            zero_modulus=true;

1960
            return Matrix<Number>(0,dim);

1961
        }

1962
    }

1963
    

1964
    //compute kernel

1965
    

1966
    Matrix<Number> Help=Cong_Slack.kernel(); // gives the solutions to the the system with slack variables

1967
    Matrix<Number> Ker_Basis(dim,dim);   // must now project to first dim coordinates to get rid of them

1968
    for(size_t i=0;i<dim;++i)

1969
        for(size_t j=0;j<dim;++j)

1970
            Ker_Basis[i][j]=Help[i][j];

1971
    return Ker_Basis;

1972
        

1973
} */

1974


1975
//---------------------------------------------------

1976


1977
template<typename Number>

1978
void Matrix<Number>::saturate(){

1979
    

1980
    // *this=kernel().kernel();

1981
    return;    // no saturation necessary over a field

1982
}

1983


1984
//---------------------------------------------------

1985


1986
template<typename Number>

1987
vector<key_t> Matrix<Number>::max_and_min(const vector<Number>& L, const vector<Number>& norm) const{

1988


1989
    vector<key_t> result(2,0);

1990
    if(nr==0)

1991
        return result;

1992
    key_t maxind=0,minind=0;

1993
    Number maxval=v_scalar_product(L,elem[0]);

1994
    Number maxnorm=1,minnorm=1;

1995
    if(norm.size()>0){

1996
        maxnorm=v_scalar_product(norm,elem[0]);

1997
        minnorm=maxnorm;              

1998
    }

1999
    Number minval=maxval;

2000
    for(key_t i=0;i<nr;++i){

2001
        Number val=v_scalar_product(L,elem[i]);

2002
        if(norm.size()==0){

2003
            if(val>maxval){

2004
                maxind=i;

2005
                maxval=val;            

2006
            }

2007
            if(val<minval){

2008
                minind=i;

2009
                minval=val;            

2010
            }

2011
        }

2012
        else{

2013
            Number nm=v_scalar_product(norm,elem[i]);

2014
            if(maxnorm*val>nm*maxval){

2015
                maxind=i;

2016
                maxval=val;            

2017
            }

2018
            if(minnorm*val<nm*minval){

2019
                minind=i;

2020
                minval=val;            

2021
            }

2022
        }            

2023
    }

2024
    result[0]=maxind;

2025
    result[1]=minind;

2026
    return result;

2027
}

2028


2029
/*

2030
template<typename Number>

2031
size_t Matrix<Number>::extreme_points_first(const vector<Number> norm){

2032
    

2033
    if(nr==0)

2034
        return 1;

2035
    

2036
    vector<long long> norm_copy;

2037
    

2038
    size_t nr_extr=0;

2039
    Matrix<long long> HelpMat(nr,nc);

2040
    try{

2041
        convert(HelpMat,*this);

2042
        convert(norm_copy,norm);

2043
    }

2044
    catch(ArithmeticException){

2045
        return nr_extr;        

2046
    }

2047


2048
    HelpMat.sort_lex();

2049
    

2050
    vector<bool> marked(nr,false);

2051
    size_t no_success=0;

2052
    // size_t nr_attempt=0;

2053
    while(true){

2054
        // nr_attempt++; cout << nr_attempt << endl;

2055
        vector<long long> L=v_random<long long>(nc,10);

2056
        vector<key_t> max_min_ind;

2057
        max_min_ind=HelpMat.max_and_min(L,norm_copy);

2058
            

2059
        if(marked[max_min_ind[0]] && marked[max_min_ind[1]])

2060
            no_success++;

2061
        else

2062
            no_success=0;

2063
        if(no_success > 1000)

2064
            break;

2065
        marked[max_min_ind[0]]=true;

2066
        marked[max_min_ind[1]]=true;

2067
    }

2068
    Matrix<long long> Extr(nr_extr,nc);  // the recognized extreme rays

2069
    Matrix<long long> NonExtr(nr_extr,nc); // the other generators

2070
    size_t j=0;

2071
    vector<key_t> perm(nr);

2072
    for(size_t i=0;i<nr;++i) {

2073
        if(marked[i]){

2074
            perm[j]=i;;

2075
            j++;

2076
        }

2077
    }

2078
    nr_extr=j;

2079
    for(size_t i=0;i<nr;++i) {

2080
        if(!marked[i]){

2081
            perm[j]=i;;

2082
            j++;

2083
        }

2084
    }

2085
    order_rows_by_perm(perm);    

2086
    // cout << nr_extr << "extreme points found"  << endl;

2087
    return nr_extr;

2088
    // exit(0);

2089
}

2090


2091
template<typename Number>

2092
vector<Number> Matrix<Number>::find_inner_point(){

2093
    vector<key_t> simplex=max_rank_submatrix_lex();

2094
    vector<Number> point(nc);

2095
    for(size_t i=0;i<simplex.size();++i)

2096
        point=v_add(point,elem[simplex[i]]);

2097
   return point;    

2098
}

2099
*/

2100


2101
template<typename Number>

2102
void Matrix<Number>::Shrink_nr_rows(size_t new_nr_rows){

2103


2104
    if(new_nr_rows>=nr)

2105
        return;

2106
    nr=new_nr_rows;

2107
    elem.resize(nr);

2108
}

2109


2110
template<typename Number>

2111
Matrix<Number>  readMatrix(const string project){

2112
// reads one matrix from file with name project

2113
// format: nr of rows, nr of colimns, entries

2114
// all separated by white space

2115
    

2116
    string name_in=project;

2117
    const char* file_in=name_in.c_str();

2118
    ifstream in;

2119
    in.open(file_in,ifstream::in);

2120
    if (in.is_open()==false){

2121
        cerr << "Cannot find input file" << endl;

2122
        exit(1);

2123
    }

2124
    

2125
    int nrows,ncols;

2126
    in >> nrows;

2127
    in >> ncols;

2128
    

2129
    if(nrows==0 || ncols==0){

2130
        cerr << "Matrix empty" << endl;

2131
        exit(1);

2132
    }

2133
    

2134
    

2135
    int i,j,entry;

2136
    Matrix<Number> result(nrows,ncols);

2137
    

2138
    for(i=0;i<nrows;++i)

2139
        for(j=0;j<ncols;++j){

2140
            in >> entry;

2141
            result[i][j]=entry;

2142
        }

2143
    return result;

2144
}

2145


2146
/*

2147
#ifndef NMZ_MIC_OFFLOAD  //offload with long is not supported

2148
template Matrix<long>  readMatrix(const string project);

2149
#endif // NMZ_MIC_OFFLOAD

2150
template Matrix<long long>  readMatrix(const string project);

2151
template Matrix<mpz_class>  readMatrix(const string project);

2152
*/

2153


2154
}  // namespace

2155


2156