Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download

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

563603 views
Tweet
Share
Share
1
#include "typedef.h"

2
#include "gmp.h"

3
#include "longtools.h"

4
#include "matrix.h"

5
/**************************************************************************\

6
@---------------------------------------------------------------------------

7
@---------------------------------------------------------------------------

8
@ FILE: long_rein_mat.c

9
@---------------------------------------------------------------------------

10
@---------------------------------------------------------------------------

11
@

12
\**************************************************************************/

13


14


15
/**************************************************************************\

16
@---------------------------------------------------------------------------

17
@ matrix_TYP *long_rein_mat(Mat)

18
@ matrix_TYP *Mat;

19
@

20
@ long_rein_mat calculates a matrix M such that the rows of M

21
@ form a Z-basis of the lattice of all integral points in the

22
@ the vectorspace generated by the rows of Mat.

23
@---------------------------------------------------------------------------

24
@

25
\**************************************************************************/

26
matrix_TYP *long_rein_mat(Mat)

27
matrix_TYP *Mat;

28
{

29


30
   int i,j,k, step, stepclear;

31
   int cols, rows, rang;

32
   MP_INT **trf, **M, **Mt, **M1, **M2, **M3, **M4, **T1, *merkpointer;

33
   matrix_TYP *erg;

34
   MP_INT a1, a2, x1, x2, y1, y2, merk, g, f;

35
 

36
   cols = Mat->cols;

37
   rows = Mat->rows;

38
   mpz_init(&a1), mpz_init(&a2);

39
   mpz_init(&x1), mpz_init(&x2);

40
   mpz_init(&y1), mpz_init(&y2);

41
   mpz_init(&f), mpz_init(&g);

42
   mpz_init(&merk);

43
   /***************************************************************\

44
   | Set Mt= Mat^{tr} transform Mt to Hermite normal form.

45
   \***************************************************************/

46
   if((Mt = (MP_INT **)malloc(cols *sizeof(MP_INT *))) == NULL)

47
   {

48
     printf("malloc of 'Mt' in 'long_elt_mat' failed\n");

49
     exit(2);

50
   }

51
   for(i=0;i<cols;i++)

52
   {

53
     if((Mt[i] = (MP_INT *)malloc(rows *sizeof(MP_INT))) == NULL)

54
     {

55
       printf("malloc of 'Mt[%d]' in 'long_elt_mat' failed\n", i);

56
       exit(2);

57
     }

58
     for(j=0;j<rows;j++)

59
       mpz_init_set_si(&Mt[i][j], Mat->array.SZ[j][i]);

60
   }

61


62
   rang = MP_hnf(Mt, cols, rows);

63
   /***************************************************************\

64
   | Set trf =  left_tans

65
   \***************************************************************/

66
     if((trf = (MP_INT **)malloc(rows *sizeof(MP_INT *))) == NULL)

67
     {

68
       printf("malloc of 'trf' in 'long_elt_mat' failed\n");

69
       exit(2);

70
     }

71
     for(i=0;i<rows;i++)

72
     {

73
       if((trf[i] = (MP_INT *)malloc(rows *sizeof(MP_INT))) == NULL)

74
       {

75
         printf("malloc of 'trf[%d]' in 'long_elt_mat' failed\n", i);

76
         exit(2);

77
       }

78
       for(j=0;j<rows;j++)

79
         mpz_init_set_si(&trf[i][j], 0);

80
       mpz_set_si(&trf[i][i], 1);

81
     }

82
   /***************************************************************\

83
   | Set M= Mt^{tr} transform Mt to Hermite normal form.

84
   \***************************************************************/

85
   if((M = (MP_INT **)malloc(rows *sizeof(MP_INT *))) == NULL)

86
   {

87
     printf("malloc of 'M' in 'long_elt_mat' failed\n");

88
     exit(2);

89
   }

90
   for(i=0;i<rows;i++)

91
   {

92
     if((M[i] = (MP_INT *)malloc(cols *sizeof(MP_INT))) == NULL)

93
     {

94
       printf("malloc of 'M[%d]' in 'long_elt_mat' failed\n", i);

95
       exit(2);

96
     }

97
     for(j=0;j<cols;j++)

98
       mpz_init_set(&M[i][j], &Mt[j][i]);

99
   }

100


101
   /* output for debugging purposes

102
   dump_MP_mat(trf,rows,rows,"trf");

103
   dump_MP_mat(M,rows,rows,"M"); */

104


105
   rang = MP_hnf_simultaneous(M, rows, cols, trf, rows);

106


107
   /* output for debugging purposes

108
   dump_MP_mat(trf,rows,rows,"trf");

109
   dump_MP_mat(M,rows,rows,"M"); */

110


111
   /***************************************************************\

112
   | Clear the space allocated for 'Mt'

113
   \***************************************************************/

114
   for(i=0;i<cols;i++)

115
   {

116
     for(j=0;j<rows;j++)

117
       mpz_clear(&Mt[i][j]);

118
     free(Mt[i]);

119
   }

120
   free(Mt);

121


122
   /***************************************************************\

123
   | Now the elementary divisor algorithm starts

124
   \***************************************************************/

125
   for(step = 0; step < rang;step++)

126
   {

127
      do

128
      {

129
         /*------------------------------------------------------*\

130
         | Clear the 'step'^th row of M

131
         \*------------------------------------------------------*/

132
          for(i=step;i<rang && mpz_cmp_si(&M[step][i], 0) == 0; i++);

133
          if(i!=step)

134
          {

135


136
            /* swap the step-th col with the i-th col */

137
            for(j=step;j<rang;j++)

138
            {

139
              mpz_set(&merk, &M[j][step]);

140
              mpz_set(&M[j][step], &M[j][i]);

141
              mpz_set(&M[j][i], &merk);

142
            }

143
          }

144
          if(mpz_cmp_si(&M[step][step], 0) < 0)

145
          {

146
            for(i=step;i<rows;i++)

147
              mpz_neg(&M[i][step], &M[i][step]);

148
          }

149
          for(i=step+1;i<rang;i++)

150
          {

151
            if(mpz_cmp_si(&M[step][i], 0) != 0)

152
            {

153
              if(mpz_cmp_si(&M[step][step], 1) == 0)

154
              {

155
                mpz_set(&f, &M[step][i]);

156
                for(j=step+1;j<rang;j++)

157
                {

158
                  mpz_mul(&merk, &M[j][step], &f);

159
                  mpz_sub(&M[j][i], &M[j][i], &merk);

160
                }

161
                mpz_set_si(&M[step][i], 0);

162


163
                   /*

164
                   for(j=0;j<rows;j++)

165
                   {

166
                     mpz_mul(&merk, &trf[j][step], &f);

167
                     mpz_sub(&trf[j][i], &trf[j][i], &merk);

168
                   }*/

169
              }

170
              else

171
              {

172
                 /* changed 11/06/99 (tilman) from

173
                 mpz_gcdext(&g, &a1, &a2, &M[step][step], &M[step][i]); to */

174
                 mpz_abs( &g,&M[step][i]);

175
                 if (mpz_cmp(&M[step][step],&g) == 0){

176
                    mpz_set_si(&a1,1);

177
                    mpz_set_si(&a2,0);

178
                 }

179
                 else{

180
                    mpz_gcdext(&g, &a1, &a2, &M[step][step], &M[step][i]);

181
                 }

182


183


184
                 mpz_div(&x1, &M[step][i], &g);

185
                 mpz_div(&x2, &M[step][step], &g);

186
                 mpz_neg(&x2, &x2);

187


188
                 for(j=step+1; j<rang;j++)

189
                 {

190
                    mpz_mul(&y1, &a1, &M[j][step]);

191
                    mpz_mul(&merk, &a2, &M[j][i]);

192
                    mpz_add(&y1, &y1, &merk);

193


194
                    mpz_mul(&y2, &x1, &M[j][step]);

195
                    mpz_mul(&merk, &x2, &M[j][i]);

196
                    mpz_add(&y2, &y2, &merk);

197


198
                    mpz_set(&M[j][step], &y1);

199
                    mpz_set(&M[j][i], &y2);

200
                 }

201
                 mpz_set(&M[step][step], &g);

202
                 mpz_set_si(&M[step][i], 0);

203


204
                 /*

205
                   for(j=0; j<rows;j++)

206
                   {

207
                      mpz_mul(&y1, &a1, &trf[j][step]);

208
                      mpz_mul(&merk, &a2, &trf[j][i]);

209
                      mpz_add(&y1, &y1, &merk);

210


211
                      mpz_mul(&y2, &x1, &trf[j][step]);

212
                      mpz_mul(&merk, &x2, &trf[j][i]);

213
                      mpz_add(&y2, &y2, &merk);

214


215
                      mpz_set(&trf[j][step], &y1);

216
                      mpz_set(&trf[j][i], &y2);

217
                   }*/

218
              }

219
            }

220
          }

221
         /*------------------------------------------------------*\

222
         | Clear the 'step'^th column of M

223
         \*------------------------------------------------------*/

224
          for(i=step;i<rang && mpz_cmp_si(&M[i][step], 0) == 0; i++);

225
          if(i!=step)

226
          {

227
            merkpointer = M[step];

228
            M[step] = M[i];

229
            M[i] = merkpointer;

230


231
              merkpointer = trf[step];

232
              trf[step] = trf[i];

233
              trf[i] = merkpointer;

234
          }

235
          if(mpz_cmp_si(&M[step][step], 0) < 0)

236
          {

237
            for(i=step;i<rows;i++)

238
              mpz_neg(&M[i][step], &M[i][step]);

239
          }

240
          for(i=step+1;i<rang;i++)

241
          {

242
            if(mpz_cmp_si(&M[i][step], 0) != 0)

243
            {

244
              if(mpz_cmp_si(&M[step][step], 1) == 0)

245
              {

246
                mpz_set(&f, &M[i][step]);

247
                for(j=step+1;j<rang;j++)

248
                {

249
                  mpz_mul(&merk, &M[step][j], &f);

250
                  mpz_sub(&M[i][j], &M[i][j], &merk);

251
                }

252
                mpz_set_si(&M[i][step], 0);

253


254
                /* inserted the next 4 lines: tilman 4/12/96 */

255
                for (j=0;j<rang;j++){

256
                  mpz_mul(&merk, &trf[step][j], &f);

257
                  mpz_sub(&trf[i][j], &trf[i][j], &merk);

258
                }

259


260
              }

261
              else

262
              {

263
                 /* changed on 8/1/97 tilman from

264
                 mpz_gcdext(&g, &a1, &a2, &M[step][step], &M[i][step]);

265
                 to : */

266
                 mpz_abs( &g,&M[i][step]);

267
                 if (mpz_cmp(&M[step][step],&g) == 0){

268
                    mpz_set_si(&a1,1);

269
                    mpz_set_si(&a2,0);

270
                 }

271
                 else{

272
                    mpz_gcdext(&g, &a1, &a2, &M[step][step], &M[i][step]);

273
                 }

274


275
                 mpz_div(&x1, &M[i][step], &g);

276
                 mpz_div(&x2, &M[step][step], &g);

277
                 mpz_neg(&x2, &x2);

278


279
                 for(j=step+1; j<rang;j++)

280
                 {

281
                    mpz_mul(&y1, &a1, &M[step][j]);

282
                    mpz_mul(&merk, &a2, &M[i][j]);

283
                    mpz_add(&y1, &y1, &merk);

284


285
                    mpz_mul(&y2, &x1, &M[step][j]);

286
                    mpz_mul(&merk, &x2, &M[i][j]);

287
                    mpz_add(&y2, &y2, &merk);

288


289
                    mpz_set(&M[step][j], &y1);

290
                    mpz_set(&M[i][j], &y2);

291
                 }

292


293
                 /* inserted the next 12 lines: 4/12/96 tilman */

294
                 for (j=0;j<rows;j++){

295
                    mpz_mul(&y1, &a1, &trf[step][j]);

296
                    mpz_mul(&merk, &a2, &trf[i][j]);

297
                    mpz_add(&y1, &y1, &merk);

298


299
                    mpz_mul(&y2, &x1, &trf[step][j]);

300
                    mpz_mul(&merk, &x2, &trf[i][j]);

301
                    mpz_add(&y2, &y2, &merk);

302


303
                    mpz_set(&trf[step][j], &y1);

304
                    mpz_set(&trf[i][j], &y2);

305
                 } 

306


307
                 mpz_set(&M[step][step], &g);

308
                 mpz_set_si(&M[i][step], 0);

309
              }

310
            }

311
          }

312
         /*------------------------------------------------------*\

313
         | Check whether the column and row are both clear

314
         \*------------------------------------------------------*/

315
         stepclear = TRUE;

316
         for(i=step+1;i<cols && stepclear == TRUE;i++)

317
         {

318
           if(mpz_cmp_si(&M[step][i], 0) != 0)

319
             stepclear = FALSE;

320
         }

321
         for(i=step+1;i<rows && stepclear == TRUE;i++)

322
         {

323
           if(mpz_cmp_si(&M[i][step], 0) != 0)

324
             stepclear = FALSE;

325
         }

326
      }while(stepclear == FALSE);

327
      if(mpz_cmp_si(&M[step][step], 0) < 0)

328
         mpz_neg(&M[step][step], &M[step][step]);

329


330
      /* output for debugging purposes

331
      printf("step = %d\n",step);

332
      dump_MP_mat(trf,rows,rows,"trf");

333
      dump_MP_mat(M,rows,rows,"M"); */

334


335
   }

336


337
   /* output for debugging purposes

338
   dump_MP_mat(trf,rows,rows,"trf");

339
   dump_MP_mat(M,rows,rows,"M"); */

340


341
   /*-------------------------------------------------------------------*\

342
   |  Multiply trf with Mat

343
   \*-------------------------------------------------------------------*/

344
   M1 = matrix_to_MP_mat(Mat);

345
   if( (M2 = (MP_INT **)malloc(rang *sizeof(MP_INT *))) == NULL)

346
   {

347
      printf("malloc of 'M2' in 'long_rein_mat'failed\n");

348
      exit(2);

349
   }

350
   for(i=0;i<rang;i++)

351
   {

352
     if( (M2[i] = (MP_INT *)malloc(cols *sizeof(MP_INT))) == NULL)

353
     {

354
        printf("malloc of 'M2[%d]' in 'long_rein_mat'failed\n", i);

355
        exit(2);

356
     }

357
     for(j=0;j<cols;j++)

358
        mpz_init(&M2[i][j]);

359
   }

360
   for(i=0;i<rang;i++)

361
     for(j=0;j<cols;j++)

362
     {

363
        for(k=0;k<rows;k++)

364
        {

365
          mpz_mul(&merk, &trf[i][k], &M1[k][j]);

366
          mpz_add(&M2[i][j], &M2[i][j], &merk);

367
        }

368


369
        /* changed 4/12/96 tilman from:

370
        mpz_div(&M2[i][j], &M2[i][i], &M[i][i]);

371
        to : */

372
        mpz_div(&M2[i][j], &M2[i][j], &M[i][i]);

373


374
     }

375


376
   /* output for debugging purposes

377
   dump_MP_mat(M2,rows,rows,"M2"); */

378


379
   /***************************************************************\

380
   | make a pair_reduction to the rows of M2

381
   \***************************************************************/

382
   M3 = init_MP_mat(rang, rang);

383
   for(i=0;i<rang;i++)

384
     for(j=0;j<=i;j++)

385
     {

386
       for(k=0;k<cols;k++)

387
       {

388
        mpz_mul(&merk, &M2[i][k], &M2[j][k]);

389
        mpz_add(&M3[i][j], &merk, &M3[i][j]);

390
       }

391
     }

392
   for(i=0;i<rang;i++)

393
     for(j=0;j<i;j++)

394
       mpz_set(&M3[j][i], &M3[i][j]);

395


396
   /* output for debugging purposes

397
   dump_MP_mat(M3,rows,rows,"M3"); */

398


399
   T1 = init_MP_mat(rang, rang);

400
   for(i=0;i<rang;i++)

401
     mpz_set_si(&T1[i][i], 1);

402
   MP_pair_red(M3, T1, rang);

403


404
   /* output for debugging purposes

405
   dump_MP_mat(T1,rows,rows,"T1"); */

406


407
   M4 = init_MP_mat(rang, cols);

408
   for(i=0;i<rang;i++)

409
    for(j=0;j<cols;j++)

410
      for(k=0;k<rang;k++)

411
      {

412
         mpz_mul(&merk, &T1[i][k], &M2[k][j]);

413
         mpz_add(&M4[i][j], &M4[i][j], &merk);

414
      }

415
   

416
   /* output for debugging purposes

417
   dump_MP_mat(M4,rows,rows,"M4"); */

418


419
   erg = MP_mat_to_matrix(M4, rang, cols);

420


421
   free_MP_mat(M, rows, cols); free(M);

422
   free_MP_mat(M1, rows, cols); free(M1);

423
   free_MP_mat(M2, rang, cols); free(M2);

424
   free_MP_mat(M3, rang, rang); free(M3);

425
   free_MP_mat(M4, rang, cols); free(M4);

426
   free_MP_mat(trf, rows,rows); free(trf);

427
   /* changed 11/3/97 tilman from

428
   free_MP_mat(T1, rang, rang); free(M2);

429
   to: */

430
   free_MP_mat(T1, rang, rang); free(T1);

431
   merkpointer = NULL;

432
   mpz_clear(&a1), mpz_clear(&a2);

433
   mpz_clear(&x1), mpz_clear(&x2);

434
   mpz_clear(&y1), mpz_clear(&y2);

435
   mpz_clear(&f), mpz_clear(&g);

436
   mpz_clear(&merk);

437
   return(erg);

438
}

439


440
/*************************************************************************

441
@

442
@------------------------------------------------------------------------

443
@

444
@ void long_rein_formspace(matrix_TYP **forms,int number,int option)

445
@

446
@ option:       drives the behaviour of the function, for

447
@               option == 0 there is no assumption on the symetry

448
@                           of these matrices made.

449
@               option == 1 the matrices are assumed to be symetric.

450
@               option == 2 the matrices are assumed to be skew symetric.

451
@

452
@------------------------------------------------------------------------

453
@

454
**************************************************************************/

455
int long_rein_formspace(matrix_TYP **forms,int number,int option)

456
{

457
   int i,

458
       j,

459
       k,

460
       l,

461
       rank,

462
       dim = forms[0]->cols;

463


464
   matrix_TYP *rein,

465
              *res;

466


467
   if (option == 0){

468
      rein = init_mat(number,dim*dim,"");

469
   }

470
   else if (option == 1){

471
      rein = init_mat(number,(dim*(dim+1))/2,"");

472
   }

473
   else if (option == 2){

474
      rein = init_mat(number,(dim*(dim+2))/2,"");

475
   }

476
   else{

477
      fprintf(stderr,"you must specify one of the options 0,1,2 in");

478
      fprintf(stderr," long_rein_formspace\n");

479
      exit(1);

480
   }

481


482
   for (i=0;i<number;i++){

483
      l=0;

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

485
         if (option == 0){

486
            for (k=0;k<dim;k++){

487
               rein->array.SZ[i][l] = forms[i]->array.SZ[j][k];

488
               l++;

489
            }

490
         }

491
         else if (option == 1){

492
            for (k=0;k<=j;k++){

493
               rein->array.SZ[i][l] = forms[i]->array.SZ[j][k];

494
               l++;

495
            }

496
         }

497
         else{

498
            for (k=0;k<j;k++){

499
               rein->array.SZ[i][l] = forms[i]->array.SZ[j][k];

500
               l++;

501
            }

502
         }

503
      }

504
   } 

505


506
   res = long_rein_mat(rein);

507


508
   for (i=0;i<res->rows;i++){

509
      l=0;

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

511
         if (option == 0){

512
            for (k=0;k<dim;k++){

513
               forms[i]->array.SZ[j][k] = res->array.SZ[i][l];

514
               l++;

515
            }

516
         }

517
         else if (option == 1){

518
            for (k=0;k<=j;k++){

519
               forms[i]->array.SZ[j][k] = res->array.SZ[i][l];

520
               forms[i]->array.SZ[k][j] = res->array.SZ[i][l];

521
               l++;

522
            }

523
         }

524
         else{

525
            for (k=0;k<j;k++){

526
               forms[i]->array.SZ[j][k] = res->array.SZ[i][l];

527
               forms[i]->array.SZ[k][j] = -res->array.SZ[i][l];

528
               l++;

529
            }

530
            forms[i]->array.SZ[j][j] = 0;

531
         }

532
      }

533
   } 

534


535
   rank = res->rows;

536


537
   for (i=0;i<number;i++) Check_mat(forms[i]);

538


539
   free_mat(rein);

540
   free_mat(res);

541


542
   return rank;

543
}

544


545