CoCalc -- MP

GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
gap4r8 / pkg / carat / carat-2.1b1 / functions / Longtools / MP_hnf.c
⁵⁶³⁶⁰³ views
1
#include "typedef.h"
2
#include "gmp.h"
3
/* #include "gmp-impl.h" */
4

5
/**************************************************************************\
6
@---------------------------------------------------------------------------
7
@---------------------------------------------------------------------------
8
@ FILE: MP_hnf.c
9
@---------------------------------------------------------------------------
10
@---------------------------------------------------------------------------
11
@
12
@
13
@ The algorithms do the same as the algorithms in the
14
@ file MP_gauss.c, the only difference the the matrix is in
15
@ Hermite-normal-form after the algorithm
16
@
17
\**************************************************************************/
18

19
/**************************************************************************\
20
@---------------------------------------------------------------------------
21
@ int MP_trf_hnf(M, Trf, rows, cols)
22
@ MP_INT **M, **Trf;
23
@ int rows, cols;
24
@
25
@---------------------------------------------------------------------------
26
\**************************************************************************/
27
/************************************************************************\
28
| The Matrix M is transformed to a matrix M' such that
29
| M' = Trf * M is Gauss-reduced, with Trf integral with determinante +-1.
30
| CAUTION: The entries of the matrix M are changed.
31
\************************************************************************/
32

33
int MP_trf_hnf(M, Trf, rows, cols)
34
MP_INT **M, **Trf;
35
int rows, cols;
36
{
37
  int i,j,k;
38
  int n;
39
  int step;
40
  MP_INT a1,a2,gcd, *v, f;
41
  int tester;
42
  int spos = 0;
43
  MP_INT x1,x2,y1,y2;
44
  MP_INT Mi, Ms;
45
  int rang = 0;
46

47
  n = rows;
48
  mpz_init(&f); mpz_init(&a1); mpz_init(&a2); mpz_init(&gcd);
49
  mpz_init(&x1); mpz_init(&x2);
50
  mpz_init(&y1); mpz_init(&y2);
51
  mpz_init(&Ms); mpz_init(&Mi);
52
  for(i=0;i<n;i++)
53
   for(j=0;j<n; j++)
54
      mpz_set_si(&Trf[i][j], 0);
55
  for(i=0;i<n;i++)
56
      mpz_set_si(&Trf[i][i], 1);
57
  for(step = 0; step < n && spos != cols; step++)
58
  {
59
    tester = FALSE;
60
    while(tester == FALSE && spos != cols)
61
    {
62
       i = step;
63
       while(i<n && mpz_cmp_si(&M[i][spos], 0) == 0)
64
          i++;
65
       if(i<n)
66
       {
67
         tester = TRUE;
68
         if(i != step)
69
         {
70
           v = M[i];
71
           M[i] = M[step];
72
           M[step] = v;
73
           v = Trf[i];
74
           Trf[i] = Trf[step];
75
           Trf[step] = v;
76
         }
77
       }
78
       else
79
        spos++;
80
    }
81
    for(i=step+1;i<n && spos != cols;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
    if(spos != cols)
135
      rang++;
136
  }
137
  mpz_set_si(&y2, 2);
138
  mpz_set_si(&y1, 1);
139
  for(step = 0; step < rang; step++)
140
  {
141
    for(i=0;i<cols && mpz_cmp_si(&M[step][i], 0) == 0;i++);
142
    if(i != cols && mpz_cmp_si(&M[step][i], 0) < 0)
143
    {
144
       for(j=i;j<cols;j++)
145
         mpz_neg(&M[step][j], &M[step][j]);
146
       for(j=0;j<rows;j++)
147
         mpz_neg(&Trf[step][j], &Trf[step][j]);
148
    }
149
  }
150
  for(step = 0; step <rang-1;step++)
151
  {
152
     for(i=step+1;i<rang;i++)
153
     {
154
       for(spos = 0; spos<cols && mpz_cmp_si(&M[i][spos], 0) == 0;spos++);
155
       if(spos != cols)
156
       {
157
          mpz_set_si(&f, 0);
158
          mpz_div(&a1, &M[i][spos], &y2);
159
          if(mpz_cmp(&a1, &M[step][spos]) < 0)
160
          {
161
             mpz_div(&f, &M[step][spos], &M[i][spos]);
162
             mpz_mul(&x1, &f, &M[i][spos]);
163
             mpz_sub(&x2, &M[step][spos], &x1);
164
             if(mpz_cmp(&a1, &x2) < 0)
165
               mpz_add(&f, &f, &y1);
166
          }
167
          mpz_neg(&a1, &a1);
168
          if(mpz_cmp(&a1, &M[step][spos]) > 0)
169
          {
170
             mpz_div(&f, &M[step][spos], &M[i][spos]);
171
             mpz_mul(&x1, &f, &M[i][spos]);
172
             mpz_sub(&x2, &M[step][spos], &x1);
173
             if(mpz_cmp(&a1, &x2) > 0)
174
               mpz_sub(&f, &f, &y1);
175
             mpz_mod(&x1, &M[i][spos], &y2);
176
             if(mpz_cmp_si(&y2, 0) == 0 && mpz_cmp(&a1, &x2) == 0)
177
               mpz_sub(&f, &f, &y1);
178
          }
179
          if(mpz_cmp_si(&f, 0) != 0)
180
          {
181
             for(j=spos;j<cols;j++)
182
             {
183
                mpz_mul(&x1, &f, &M[i][j]);
184
                mpz_sub(&M[step][j], &M[step][j], &x1);
185
             }
186
             for(j=0;j<rows;j++)
187
             {
188
                mpz_mul(&x1, &f, &Trf[i][j]);
189
                mpz_sub(&Trf[step][j], &Trf[step][j], &x1);
190
             }
191
          }
192
       }
193
     }
194
  }
195
  
196
  mpz_clear(&f); mpz_clear(&a1); mpz_clear(&a2); mpz_clear(&gcd);
197
  mpz_clear(&x1); mpz_clear(&x2);
198
  mpz_clear(&y1); mpz_clear(&y2);
199
  mpz_clear(&Ms); mpz_clear(&Mi);
200
  return(rang);
201
}
202

203

204

205

206

207
/**************************************************************************\
208
@---------------------------------------------------------------------------
209
@ int MP_hnf(M, rows, cols)
210
@ MP_INT **M;
211
@ int rows, cols;
212
@---------------------------------------------------------------------------
213
@
214
\**************************************************************************/
215
int MP_hnf(M, rows, cols)
216
MP_INT **M;
217
int rows, cols;
218
{
219
  int i,j,k;
220
  int n;
221
  int step;
222
  MP_INT a1,a2,gcd, *v, f;
223
  int tester;
224
  int spos = 0;
225
  MP_INT x1,x2,y1,y2;
226
  MP_INT Mi, Ms;
227
  int rang = 0;
228

229
  n = rows;
230
  mpz_init(&f); mpz_init(&a1); mpz_init(&a2); mpz_init(&gcd);
231
  mpz_init(&x1); mpz_init(&x2);
232
  mpz_init(&y1); mpz_init(&y2);
233
  mpz_init(&Ms);
234
  mpz_init(&Mi);
235

236
  for(step = 0; step < n && spos != cols; step++)
237
  {
238
    tester = FALSE;
239
    while(tester == FALSE && spos != cols)
240
    {
241
       i = step;
242
       while(i<n && mpz_cmp_si(&M[i][spos], 0) == 0)
243
          i++;
244
       if(i<n)
245
       {
246
         tester = TRUE;
247
         if(i != step)
248
         {
249
           v = M[i];
250
           M[i] = M[step];
251
           M[step] = v;
252
         }
253
       }
254
       else
255
        spos++;
256
    }
257
    for(i=step+1;i<n && spos != cols;i++)
258
    {
259
      if(mpz_cmp_si(&M[i][spos], 0) != 0)
260
      {
261
        mpz_gcdext(&gcd, &a1, &a2, &M[step][spos], & M[i][spos]);
262
       if(mpz_cmp_si(&a1, 0) == 0)
263
       {
264
         v = M[step];
265
         M[step] = M[i];
266
         M[i] = v;
267
         mpz_set(&f, &a1);
268
         mpz_set(&a1, &a2);
269
         mpz_set(&a2, &f);
270
       }
271
       mpz_div(&Ms, &M[step][spos], &gcd);
272
       mpz_div(&Mi, &M[i][spos], &gcd);
273
       mpz_set(&M[step][spos], &Ms);
274
       mpz_set(&M[i][spos], &Mi);
275
       for(j=spos + 1; j<cols; j++)
276
       {
277
         mpz_mul(&x1, &M[step][j], &a1);
278
         mpz_mul(&x2, &M[i][j], &a2);
279
         mpz_mul(&y1, &M[i][j], &Ms);
280
         mpz_mul(&y2, &M[step][j], &Mi);
281
         mpz_add(&M[step][j], &x1, &x2);
282
         mpz_sub(&M[i][j], &y1, &y2);
283
/****************************
284
         x = M->array.SZ[step][j];
285
         y = M->array.SZ[i][j];
286
         M->array.SZ[step][j] = a1 * x + a2 * y;
287
         M->array.SZ[i][j] = Ms * y - Mi * x;
288
**************************/
289
       }
290
       mpz_set(&M[step][spos], &gcd);
291
       mpz_set_si(&M[i][spos], 0);
292
      }
293
    }
294
    if(spos != cols)
295
      rang++;
296
  }
297
  mpz_set_si(&y2, 2);
298
  mpz_set_si(&y1, 1);
299
  for(step = 0; step < rang; step++)
300
  {
301
    for(i=0;i<cols && mpz_cmp_si(&M[step][i], 0) == 0;i++);
302
    if(i != cols && mpz_cmp_si(&M[step][i], 0) < 0)
303
    {
304
       for(j=i;j<cols;j++)
305
         mpz_neg(&M[step][j], &M[step][j]);
306
    }
307
  }
308
  for(step = 0; step <rang-1;step++)
309
  {
310
     for(i=step+1;i<rang;i++)
311
     {
312
       for(spos = 0; spos<cols && mpz_cmp_si(&M[i][spos], 0) == 0;spos++);
313
       if(spos != cols)
314
       {
315

316
          mpz_div(&f,&M[step][spos],&M[i][spos]);
317

318
          if (mpz_cmp_si(&M[step][spos],0) < 0){
319
             mpz_mod(&x1,&M[step][spos],&M[i][spos]);
320
             if (mpz_cmp_si(&x1,0) != 0)
321
                mpz_sub_ui(&f,&f,1);
322
          }
323

324
          /*
325
          mpz_set_si(&f, 0);
326
          mpz_div(&a1, &M[i][spos], &y2);
327
          if(mpz_cmp(&a1, &M[step][spos]) < 0)
328
          {
329
             mpz_div(&f, &M[step][spos], &M[i][spos]);
330
             mpz_mul(&x1, &f, &M[i][spos]);
331
             mpz_sub(&x2, &M[step][spos], &x1);
332
             if(mpz_cmp(&a1, &x2) < 0)
333
               mpz_add(&f, &f, &y1);
334
          }
335
          mpz_neg(&a1, &a1);
336
          if(mpz_cmp(&a1, &M[step][spos]) >= 0)
337
          {
338
             mpz_div(&f, &M[step][spos], &M[i][spos]);
339
             mpz_mul(&x1, &f, &M[i][spos]);
340
             mpz_sub(&x2, &M[step][spos], &x1);
341
             if(mpz_cmp(&a1, &x2) > 0)
342
               mpz_sub(&f, &f, &y1);
343
             mpz_mod(&x1, &M[i][spos], &y2);
344
             if(mpz_cmp_si(&y2, 0) == 0 && mpz_cmp(&a1, &x2) == 0)
345
               mpz_sub(&f, &f, &y1);
346
          } */
347
          if(mpz_cmp_si(&f, 0) != 0)
348
          {
349
             for(j=spos;j<cols;j++)
350
             {
351
                mpz_mul(&x1, &f, &M[i][j]);
352
                mpz_sub(&M[step][j], &M[step][j], &x1);
353
             }
354
             i--;
355
          }
356
       }
357
     }
358
  }
359
  mpz_clear(&f); mpz_clear(&a1); mpz_clear(&a2); mpz_clear(&gcd);
360
  mpz_clear(&x1); mpz_clear(&x2);
361
  mpz_clear(&y1); mpz_clear(&y2);
362
  mpz_clear(&Ms);
363
  mpz_clear(&Mi);
364

365
  return(rang);
366
}
367

368

369

370
/**************************************************************************\
371
@---------------------------------------------------------------------------
372
@ int MP_hnf_simultaneous(M, rows, cols, B, Bcols)
373
@ MP_INT **M, **B;
374
@ int rows, cols, Bcols;
375
@
376
@---------------------------------------------------------------------------
377
\**************************************************************************/
378
int MP_hnf_simultaneous(M, rows, cols, B, Bcols)
379
MP_INT **M, **B;
380
int rows, cols, Bcols;
381
{
382
  int i,j,k;
383
  int n;
384
  int step;
385
  MP_INT a1,a2,gcd, *v, f;
386
  int tester;
387
  int spos = 0;
388
  MP_INT x1,x2,y1,y2;
389
  MP_INT Mi, Ms;
390
  int rang = 0;
391

392
  n = rows;
393
  mpz_init(&f); mpz_init(&a1); mpz_init(&a2); mpz_init(&gcd);
394
  mpz_init(&x1); mpz_init(&x2);
395
  mpz_init(&y1); mpz_init(&y2);
396
  mpz_init(&Ms); mpz_init(&Mi);
397
  for(step = 0; step < n && spos != cols; step++)
398
  {
399
    tester = FALSE;
400
    while(tester == FALSE && spos != cols)
401
    {
402
       i = step;
403
       while(i<n && mpz_cmp_si(&M[i][spos], 0) == 0)
404
          i++;
405
       if(i<n)
406
       {
407
         tester = TRUE;
408
         if(i != step)
409
         {
410
           v = M[i];
411
           M[i] = M[step];
412
           M[step] = v;
413
           v = B[i];
414
           B[i] = B[step];
415
           B[step] = v;
416
         }
417
       }
418
       else
419
        spos++;
420
    }
421
    for(i=step+1;i<n && spos != cols;i++)
422
    {
423
      if(mpz_cmp_si(&M[i][spos], 0) != 0)
424
      {
425
        mpz_gcdext(&gcd, &a1, &a2, &M[step][spos], & M[i][spos]);
426
       if(mpz_cmp_si(&a1, 0) == 0)
427
       {
428
         v = M[step];
429
         M[step] = M[i];
430
         M[i] = v;
431
         v = B[step];
432
         B[step] = B[i];
433
         B[i] = v;
434
         mpz_set(&f, &a1); mpz_set(&a1, &a2); mpz_set(&a2, &f);
435
       }
436
       mpz_div(&Ms, &M[step][spos], &gcd);
437
       mpz_div(&Mi, &M[i][spos], &gcd);
438
       mpz_set(&M[step][spos], &Ms);
439
       mpz_set(&M[i][spos], &Mi);
440
       for(j=0;j<Bcols;j++)
441
       {
442
         mpz_mul(&x1, &B[step][j], &a1);
443
         mpz_mul(&x2, &B[i][j], &a2);
444
         mpz_mul(&y1, &B[i][j], &Ms);
445
         mpz_mul(&y2, &B[step][j], &Mi);
446
         mpz_add(&B[step][j], &x1, &x2);
447
         mpz_sub(&B[i][j], &y1, &y2);
448
/**************************
449
         x = Trf->array.SZ[step][j];
450
         y = Trf->array.SZ[i][j];
451
         Trf->array.SZ[step][j] = a1 * x + a2 * y;
452
         Trf->array.SZ[i][j] = Ms * y - Mi * x;
453
**************************/
454
       }
455
       for(j=spos + 1; j<cols; j++)
456
       {
457
         mpz_mul(&x1, &M[step][j], &a1);
458
         mpz_mul(&x2, &M[i][j], &a2);
459
         mpz_mul(&y1, &M[i][j], &Ms);
460
         mpz_mul(&y2, &M[step][j], &Mi);
461
         mpz_add(&M[step][j], &x1, &x2);
462
         mpz_sub(&M[i][j], &y1, &y2);
463
/****************************
464
         x = M->array.SZ[step][j];
465
         y = M->array.SZ[i][j];
466
         M->array.SZ[step][j] = a1 * x + a2 * y;
467
         M->array.SZ[i][j] = Ms * y - Mi * x;
468
**************************/
469
       }
470
       mpz_set(&M[step][spos], &gcd);
471
       mpz_set_si(&M[i][spos], 0);
472
      }
473
    }
474
    if(spos != cols)
475
      rang++;
476
  }
477
  mpz_set_si(&y2, 2);
478
  mpz_set_si(&y1, 1);
479
  for(step = 0; step < rang; step++)
480
  {
481
    for(i=0;i<cols && mpz_cmp_si(&M[step][i], 0) == 0;i++);
482
    if(i != cols && mpz_cmp_si(&M[step][i], 0) < 0)
483
    {
484
       for(j=i;j<cols;j++)
485
         mpz_neg(&M[step][j], &M[step][j]);
486
       for(j=0;j<Bcols;j++)
487
         mpz_neg(&B[step][j], &B[step][j]);
488
    }
489
  }
490
  for(step = 0; step <rang-1;step++)
491
  {
492
     for(i=step+1;i<rang;i++)
493
     {
494
       for(spos = 0; spos<cols && mpz_cmp_si(&M[i][spos], 0) == 0;spos++);
495
       if(spos != cols)
496
       {
497
          mpz_set_si(&f, 0);
498
          mpz_div(&a1, &M[i][spos], &y2);
499
          if(mpz_cmp(&a1, &M[step][spos]) < 0)
500
          {
501
             mpz_div(&f, &M[step][spos], &M[i][spos]);
502
             mpz_mul(&x1, &f, &M[i][spos]);
503
             mpz_sub(&x2, &M[step][spos], &x1);
504
             if(mpz_cmp(&a1, &x2) < 0)
505
               mpz_add(&f, &f, &y1);
506
          }
507
          mpz_neg(&a1, &a1);
508
          if(mpz_cmp(&a1, &M[step][spos]) > 0)
509
          {
510
             mpz_div(&f, &M[step][spos], &M[i][spos]);
511
             mpz_mul(&x1, &f, &M[i][spos]);
512
             mpz_sub(&x2, &M[step][spos], &x1);
513
             if(mpz_cmp(&a1, &x2) > 0)
514
               mpz_sub(&f, &f, &y1);
515
             mpz_mod(&x1, &M[i][spos], &y2);
516
             if(mpz_cmp_si(&y2, 0) == 0 && mpz_cmp(&a1, &x2) == 0)
517
               mpz_sub(&f, &f, &y1);
518
          }
519
          if(mpz_cmp_si(&f, 0) != 0)
520
          {
521
             for(j=spos;j<cols;j++)
522
             {
523
                mpz_mul(&x1, &f, &M[i][j]);
524
                mpz_sub(&M[step][j], &M[step][j], &x1);
525
             }
526
             for(j=0;j<Bcols;j++)
527
             {
528
                mpz_mul(&x1, &f, &B[i][j]);
529
                mpz_sub(&B[step][j], &B[step][j], &x1);
530
             }
531
          }
532
       }
533
     }
534
  }
535
  
536
  mpz_clear(&f); mpz_clear(&a1); mpz_clear(&a2); mpz_clear(&gcd);
537
  mpz_clear(&x1); mpz_clear(&x2);
538
  mpz_clear(&y1); mpz_clear(&y2);
539
  mpz_clear(&Ms); mpz_clear(&Mi);
540
  return(rang);
541
}
542

543
Product

Resources

Company
