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

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

2
#include "gmp.h"

3
/* #include "gmp-impl.h" */

4
/**************************************************************************\

5
@---------------------------------------------------------------------------

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

7
@ FILE: MP_gauss.c

8
@---------------------------------------------------------------------------

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

10
@

11
\**************************************************************************/

12


13


14
/************************************************************************\

15
@-------------------------------------------------------------------------

16
@ int MP_trf_gauss(M, Trf, rows, cols)

17
@ MP_INT **M, **Trf;

18
@ int rows, cols;

19
@ The Matrix M is transformed to a matrix M such that

20
@ M_new = Trf * M_old is Gauss-reduced, with Trf integral with

21
@ determinant +-1.

22
@ CAUTION: The entries of the matrix M are changed.

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

24
\************************************************************************/

25


26
int MP_trf_gauss(M, Trf, rows, cols)

27
MP_INT **M, **Trf;

28
int rows, cols;

29
{

30
  int i,j,k;

31
  int n;

32
  int step;

33
  MP_INT a1,a2,gcd, *v, f;

34
  int tester;

35
  int spos = 0;

36
  MP_INT x1,x2,y1,y2;

37
  MP_INT Mi, Ms;

38


39
  n = rows;

40
  mpz_init(&f); mpz_init(&a1); mpz_init(&a2); mpz_init(&gcd);

41
  mpz_init(&x1); mpz_init(&x2);

42
  mpz_init(&y1); mpz_init(&y2);

43
  mpz_init(&Ms); mpz_init(&Mi);

44
  for(i=0;i<n;i++)

45
   for(j=0;j<n; j++)

46
      mpz_set_si(&Trf[i][j], 0);

47
  for(i=0;i<n;i++)

48
      mpz_set_si(&Trf[i][i], 1);

49
  for(step = 0; step < n; step++)

50
  {

51
    tester = FALSE;

52
    while(tester == FALSE)

53
    {

54
       i = step;

55
       while(i<n && mpz_cmp_si(&M[i][spos], 0) == 0)

56
          i++;

57
       if(i<n)

58
       {

59
         tester = TRUE;

60
         if(i != step)

61
         {

62
           v = M[i];

63
           M[i] = M[step];

64
           M[step] = v;

65
           v = Trf[i];

66
           Trf[i] = Trf[step];

67
           Trf[step] = v;

68
         }

69
       }

70
       else

71
        spos++;

72
       if(spos == cols)

73
       {

74
         mpz_clear(&f); mpz_clear(&a1); mpz_clear(&a2); mpz_clear(&gcd);

75
         mpz_clear(&x1); mpz_clear(&x2);

76
         mpz_clear(&y1); mpz_clear(&y2);

77
         mpz_clear(&Ms); mpz_clear(&Mi);

78
         return(step);

79
       }

80
    }

81
    for(i=step+1;i<n;i++)

82
    {

83
      if(mpz_cmp_si(&M[i][spos], 0) != 0)

84
      {

85
        mpz_gcdext(&gcd, &a1, &a2, &M[step][spos], &M[i][spos]);

86
       if(mpz_cmp_si(&a1, 0) == 0)

87
       {

88
         v = M[step];

89
         M[step] = M[i];

90
         M[i] = v;

91
         v = Trf[step];

92
         Trf[step] = Trf[i];

93
         Trf[i] = v;

94
         mpz_set(&f, &a1); mpz_set(&a1, &a2); mpz_set(&a2, &f);

95
       }

96
       mpz_div(&Ms, &M[step][spos], &gcd);

97
       mpz_div(&Mi, &M[i][spos], &gcd);

98
       mpz_set(&M[step][spos], &Ms);

99
       mpz_set(&M[i][spos], &Mi);

100
       for(j=0;j<n;j++)

101
       {

102
         mpz_mul(&x1, &Trf[step][j], &a1);

103
         mpz_mul(&x2, &Trf[i][j], &a2);

104
         mpz_mul(&y1, &Trf[i][j], &Ms);

105
         mpz_mul(&y2, &Trf[step][j], &Mi);

106
         mpz_add(&Trf[step][j], &x1, &x2);

107
         mpz_sub(&Trf[i][j], &y1, &y2);

108
/**************************

109
         x = Trf->array.SZ[step][j];

110
         y = Trf->array.SZ[i][j];

111
         Trf->array.SZ[step][j] = a1 * x + a2 * y;

112
         Trf->array.SZ[i][j] = Ms * y - Mi * x;

113
**************************/

114
       }

115
       for(j=spos + 1; j<cols; j++)

116
       {

117
         mpz_mul(&x1, &M[step][j], &a1);

118
         mpz_mul(&x2, &M[i][j], &a2);

119
         mpz_mul(&y1, &M[i][j], &Ms);

120
         mpz_mul(&y2, &M[step][j], &Mi);

121
         mpz_add(&M[step][j], &x1, &x2);

122
         mpz_sub(&M[i][j], &y1, &y2);

123
/****************************

124
         x = M->array.SZ[step][j];

125
         y = M->array.SZ[i][j];

126
         M->array.SZ[step][j] = a1 * x + a2 * y;

127
         M->array.SZ[i][j] = Ms * y - Mi * x;

128
**************************/

129
       }

130
       mpz_set(&M[step][spos], &gcd);

131
       mpz_set_si(&M[i][spos], 0);

132
      }

133
    }

134
  }

135
  mpz_clear(&f); mpz_clear(&a1); mpz_clear(&a2); mpz_clear(&gcd);

136
  mpz_clear(&x1); mpz_clear(&x2);

137
  mpz_clear(&y1); mpz_clear(&y2);

138
  mpz_clear(&Ms); mpz_clear(&Mi);

139
  return(n);

140
}

141


142


143
/**************************************************************************\

144
@---------------------------------------------------------------------------

145
@ int MP_row_gauss(M, rows, cols)

146
@ MP_INT **M;

147
@ int rows, cols;

148
@

149
@ The same as MP_trf_gauss without calulating the transformation

150
@ matrix.

151
@---------------------------------------------------------------------------

152
@

153
\**************************************************************************/

154


155
int MP_row_gauss(M, rows, cols)

156
MP_INT **M;

157
int rows, cols;

158
{

159
  int i,j,k;

160
  int n;

161
  int step;

162
  MP_INT a1,a2,gcd, *v, f;

163
  int tester;

164
  int spos = 0;

165
  MP_INT x1,x2,y1,y2;

166
  MP_INT Mi, Ms;

167


168
  n = rows;

169
  mpz_init(&f); mpz_init(&a1); mpz_init(&a2); mpz_init(&gcd);

170
  mpz_init(&x1); mpz_init(&x2);

171
  mpz_init(&y1); mpz_init(&y2);

172
  mpz_init(&Ms); mpz_init(&Mi);

173
  for(step = 0; step < n; step++)

174
  {

175
    tester = FALSE;

176
    while(tester == FALSE)

177
    {

178
       i = step;

179
       while(i<n && mpz_cmp_si(&M[i][spos], 0) == 0)

180
          i++;

181
       if(i<n)

182
       {

183
         tester = TRUE;

184
         if(i != step)

185
         {

186
           v = M[i];

187
           M[i] = M[step];

188
           M[step] = v;

189
         }

190
       }

191
       else

192
        spos++;

193
       if(spos == cols)

194
       {

195
         mpz_clear(&f); mpz_clear(&a1); mpz_clear(&a2); mpz_clear(&gcd);

196
         mpz_clear(&x1); mpz_clear(&x2);

197
         mpz_clear(&y1); mpz_clear(&y2);

198
         mpz_clear(&Ms); mpz_clear(&Mi);

199
         return(step);

200
       }

201
    }

202
    for(i=step+1;i<n;i++)

203
    {

204
      if(mpz_cmp_si(&M[i][spos], 0) != 0)

205
      {

206
        mpz_gcdext(&gcd, &a1, &a2, &M[step][spos], & M[i][spos]);

207
       if(mpz_cmp_si(&a1, 0) == 0)

208
       {

209
         v = M[step];

210
         M[step] = M[i];

211
         M[i] = v;

212
         mpz_set(&f, &a1); mpz_set(&a1, &a2); mpz_set(&a2, &f);

213
       }

214
       mpz_div(&Ms, &M[step][spos], &gcd);

215
       mpz_div(&Mi, &M[i][spos], &gcd);

216
       mpz_set(&M[step][spos], &Ms);

217
       mpz_set(&M[i][spos], &Mi);

218
       for(j=spos + 1; j<cols; j++)

219
       {

220
         mpz_mul(&x1, &M[step][j], &a1);

221
         mpz_mul(&x2, &M[i][j], &a2);

222
         mpz_mul(&y1, &M[i][j], &Ms);

223
         mpz_mul(&y2, &M[step][j], &Mi);

224
         mpz_add(&M[step][j], &x1, &x2);

225
         mpz_sub(&M[i][j], &y1, &y2);

226
/****************************

227
         x = M->array.SZ[step][j];

228
         y = M->array.SZ[i][j];

229
         M->array.SZ[step][j] = a1 * x + a2 * y;

230
         M->array.SZ[i][j] = Ms * y - Mi * x;

231
**************************/

232
       }

233
       mpz_set(&M[step][spos], &gcd);

234
       mpz_set_si(&M[i][spos], 0);

235
      }

236
    }

237
  }

238
  mpz_clear(&f); mpz_clear(&a1); mpz_clear(&a2); mpz_clear(&gcd);

239
  mpz_clear(&x1); mpz_clear(&x2);

240
  mpz_clear(&y1); mpz_clear(&y2);

241
  mpz_clear(&Ms); mpz_clear(&Mi);

242
  return(n);

243
}

244


245


246
/**************************************************************************\

247
@---------------------------------------------------------------------------

248
@ int MP_row_gauss_simultaneous(M, rows, cols, B, Bcols)

249
@ MP_INT **M, **B;

250
@ int rows, cols, Bcols;

251
@

252
@ applies an integral gaussian algorithm (on the rows) to to 2-dimensional

253
@ array 'M' and does the simultaneous row operations on B.

254
@---------------------------------------------------------------------------

255
@

256
\**************************************************************************/

257
int MP_row_gauss_simultaneous(M, rows, cols, B, Bcols)

258
MP_INT **M, **B;

259
int rows, cols, Bcols;

260
{

261
  int i,j,k;

262
  int n;

263
  int step;

264
  MP_INT a1,a2,gcd, *v, f;

265
  int tester;

266
  int spos = 0;

267
  MP_INT x1,x2,y1,y2;

268
  MP_INT Mi, Ms;

269


270
  n = rows;

271
  mpz_init(&f); mpz_init(&a1); mpz_init(&a2); mpz_init(&gcd);

272
  mpz_init(&x1); mpz_init(&x2);

273
  mpz_init(&y1); mpz_init(&y2);

274
  mpz_init(&Ms); mpz_init(&Mi);

275
  for(step = 0; step < n; step++)

276
  {

277
    tester = FALSE;

278
    while(tester == FALSE)

279
    {

280
       i = step;

281
       while(i<n && mpz_cmp_si(&M[i][spos], 0) == 0)

282
          i++;

283
       if(i<n)

284
       {

285
         tester = TRUE;

286
         if(i != step)

287
         {

288
           v = M[i];

289
           M[i] = M[step];

290
           M[step] = v;

291
           v = B[i];

292
           B[i] = B[step];

293
           B[step] = v;

294
         }

295
       }

296
       else

297
        spos++;

298
       if(spos == cols)

299
       {

300
         mpz_clear(&f); mpz_clear(&a1); mpz_clear(&a2); mpz_clear(&gcd);

301
         mpz_clear(&x1); mpz_clear(&x2);

302
         mpz_clear(&y1); mpz_clear(&y2);

303
         mpz_clear(&Ms); mpz_clear(&Mi);

304
         return(step);

305
       }

306
    }

307
    for(i=step+1;i<n;i++)

308
    {

309
      if(mpz_cmp_si(&M[i][spos], 0) != 0)

310
      {

311
        mpz_gcdext(&gcd, &a1, &a2, &M[step][spos], & M[i][spos]);

312
       if(mpz_cmp_si(&a1, 0) == 0)

313
       {

314
         v = M[step];

315
         M[step] = M[i];

316
         M[i] = v;

317
         v = B[step];

318
         B[step] = B[i];

319
         B[i] = v;

320
         mpz_set(&f, &a1); mpz_set(&a1, &a2); mpz_set(&a2, &f);

321
       }

322
       mpz_div(&Ms, &M[step][spos], &gcd);

323
       mpz_div(&Mi, &M[i][spos], &gcd);

324
       mpz_set(&M[step][spos], &Ms);

325
       mpz_set(&M[i][spos], &Mi);

326
       for(j=0;j<Bcols;j++)

327
       {

328
         mpz_mul(&x1, &B[step][j], &a1);

329
         mpz_mul(&x2, &B[i][j], &a2);

330
         mpz_mul(&y1, &B[i][j], &Ms);

331
         mpz_mul(&y2, &B[step][j], &Mi);

332
         mpz_add(&B[step][j], &x1, &x2);

333
         mpz_sub(&B[i][j], &y1, &y2);

334
/**************************

335
         x = Trf->array.SZ[step][j];

336
         y = Trf->array.SZ[i][j];

337
         Trf->array.SZ[step][j] = a1 * x + a2 * y;

338
         Trf->array.SZ[i][j] = Ms * y - Mi * x;

339
**************************/

340
       }

341
       for(j=spos + 1; j<cols; j++)

342
       {

343
         mpz_mul(&x1, &M[step][j], &a1);

344
         mpz_mul(&x2, &M[i][j], &a2);

345
         mpz_mul(&y1, &M[i][j], &Ms);

346
         mpz_mul(&y2, &M[step][j], &Mi);

347
         mpz_add(&M[step][j], &x1, &x2);

348
         mpz_sub(&M[i][j], &y1, &y2);

349
/****************************

350
         x = M->array.SZ[step][j];

351
         y = M->array.SZ[i][j];

352
         M->array.SZ[step][j] = a1 * x + a2 * y;

353
         M->array.SZ[i][j] = Ms * y - Mi * x;

354
**************************/

355
       }

356
       mpz_set(&M[step][spos], &gcd);

357
       mpz_set_si(&M[i][spos], 0);

358
      }

359
    }

360
  }

361
  mpz_clear(&f); mpz_clear(&a1); mpz_clear(&a2); mpz_clear(&gcd);

362
  mpz_clear(&x1); mpz_clear(&x2);

363
  mpz_clear(&y1); mpz_clear(&y2);

364
  mpz_clear(&Ms); mpz_clear(&Mi);

365
  return(n);

366
}

367


368


369
/*************************************************************************

370
@

371
@------------------------------------------------------------------------

372
@

373
@ void MP_row_gauss_reverse(MP_INT **A,int rows,int cols,int option)

374
@

375
@ Performs a manhattan style gauss algorithm on the MP_INT matrix

376
@ M, which has to be gauss reduced by MP_row_gauss before.

377
@ The algorithm is not done complete for sake of speed!

378
@

379
@ MP_INT **A   : The matrix in question. It will be changed!

380
@ int rows     : the rows of A. It is suffisient to feed the rank in here.

381
@ int cols     : the number of columns of A.

382
@ int option   : this flag determines the behaviour of the function:

383
@                for option == 0 it will only perform Z elementary

384
@                tranformations, for option == 1 it will divide every

385
@                row by the gcd of it's entries (so do it Q style)!

386
@------------------------------------------------------------------------

387
@

388
*************************************************************************/

389
void MP_row_gauss_reverse(MP_INT **A,int rows,int cols,int option)

390
{

391
   int i,

392
       j,

393
       k,

394
       actcol;

395


396
   MP_INT x,

397
          y;

398


399
   mpz_init(&x);

400
   mpz_init(&y);

401


402
   for (k=rows-1;k>=0;k--){

403


404
      /* find the first non zero entry in the k-th row of A */

405
      for (actcol=0;actcol<cols && (mpz_cmp_si(A[k]+actcol,0) ==0);actcol++);

406


407
      /* make the first entry in this row positive for nice purposes */

408
      if (mpz_cmp_si(A[k]+actcol,0)<0)

409
        for (j=actcol;j<cols;j++) mpz_neg(A[k]+j,A[k]+j);

410


411
      /* bare in mind: option == 1 means we are doing gauss over Q, so

412
         divide the row by the gcd of it's entries */

413
      if (option == 1){

414
         mpz_set(&x,A[k]+actcol);

415
         for (i=actcol+1;i<cols;i++) mpz_gcd(&x,&x,A[k]+i);

416
         mpz_abs(&x,&x);

417
         if (mpz_cmp_si(&x,1) != 0)

418
            for (i=actcol;i<cols;i++) mpz_div(A[k]+i,A[k]+i,&x);

419
      }

420


421
      for (i=0;i<k;i++){

422
         mpz_div(&x,A[i]+actcol,A[k]+actcol);

423
         for (j=actcol;j<cols;j++){

424
            mpz_mul(&y,&x,A[k]+j);

425
            mpz_sub(A[i]+j,A[i]+j,&y);

426
         }

427
      }

428
   }

429


430
   mpz_clear(&x);

431
   mpz_clear(&y);

432


433
   return;

434
}

435


436


437