CoCalc -- mink

GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
gap4r8 / pkg / carat / carat-2.1b1 / functions / Reduction / mink_red.c
⁵⁶³⁶⁶⁵ views
1
#include"typedef.h"
2
#include"matrix.h"
3
#include"symm.h"
4
#include"longtools.h"
5
/**************************************************************************\
6
@---------------------------------------------------------------------------
7
@---------------------------------------------------------------------------
8
@ FILE: mink_red.c
9
@---------------------------------------------------------------------------
10
@---------------------------------------------------------------------------
11
@
12
\**************************************************************************/
13

14

15
static void Gramtrans(G, T, U, Uinv, vec, step)
16
matrix_TYP *G, *T;
17
int **U, **Uinv;
18
int *vec;
19
int step;
20
{
21
  int i,j,k,n;
22
  int **u, **uinv;
23
  int **merk1, **merk2;
24

25
  n = G->cols;
26
  u = (int **)malloc(n *sizeof(int *));
27
  uinv = (int **)malloc(n *sizeof(int *));
28
  merk1 = (int **)malloc(n *sizeof(int *));
29
  merk2 = (int **)malloc(n *sizeof(int *));
30
  for(i=0; i<n;i++)
31
  {
32
   u[i] = (int *)malloc(n *sizeof(int));
33
   uinv[i] = (int *)malloc(n *sizeof(int));
34
   merk1[i] = (int *)malloc(n *sizeof(int));
35
   merk2[i] = (int *)malloc(n *sizeof(int));
36
  }
37
  for(i=0; i<step; i++)
38
  {
39
    for(j=0; j<n; j++)
40
    {
41
      u[i][j] = 0;
42
      uinv[i][j] = 0;
43
    }
44
    u[i][i] = 1;
45
    uinv[i][i] = 1;
46
  }
47
  for(i=step; i<n; i++)
48
  {
49
    for(j=0; j<step; j++)
50
    {
51
      u[i][j] = 0;
52
      u[j][i] = 0;
53
      uinv[i][j] = 0;
54
      uinv[j][i] = 0;
55
    }
56
  }
57
  for(i=step; i<n; i++)
58
  {
59
    for(j=step; j<n;j++)
60
    {
61
       u[i][j] = U[i-step][j-step];
62
       uinv[i][j] = Uinv[i-step][j-step];
63
    }
64
  }
65
  for(i=0; i<(n-step); i++)
66
    for(j=0; j<step; j++)
67
       u[i+step][j] = -(vec[j])* (U[i][0]);
68
  for(i=0; i<step; i++)
69
     uinv[step][i] = vec[i];
70
  for(i=0; i<n; i++)
71
  {
72
     for(j=0; j<n; j++)
73
     {
74
       merk1[i][j] = 0;
75
       for(k=0; k<n; k++)
76
          merk1[i][j] += T->array.SZ[i][k] * u[k][j];
77
     }
78
  }
79
  for(i=0; i<n; i++)
80
    for(j=0;j<n;j++)
81
      T->array.SZ[i][j] = merk1[i][j];
82
  for(i=0; i<n; i++)
83
  {
84
     for(j=0; j<n; j++)
85
     {
86
       merk1[i][j] = 0;
87
       for(k=0; k<n; k++)
88
          merk1[i][j] += uinv[i][k] * G->array.SZ[k][j];
89
     }
90
  }
91
  for(i=0; i<n; i++)
92
  {
93
     for(j=0; j<n; j++)
94
     {
95
        merk2[i][j] = 0;
96
        for(k=0; k<n; k++)
97
          merk2[i][j] += merk1[i][k] * uinv[j][k];
98
     }
99
  }
100
  for(i=0; i<n; i++)
101
    for(j=0;j<n;j++)
102
      G->array.SZ[i][j] = merk2[i][j];
103
  for(i=0; i<n; i++)
104
    { free(u[i]); free(uinv[i]); free(merk1[i]); free(merk2[i]);}
105
  free(u); free(uinv); free(merk1); free(merk2);
106
}
107

108

109
static void kurvectrans(SV, U, no, step)
110
matrix_TYP *SV;
111
int **U;
112
int no, step;
113
{
114
  int i,j, k;
115
  int n;
116
  int *vec;
117

118
  n = SV->cols;
119
  vec = (int *)malloc(n *sizeof(int));
120
  for(i=0;i<no;i++)
121
  {
122
    for(j=0; j<step; j++)
123
      vec[j] = SV->array.SZ[i][j];
124
    for(j=step;j<n; j++)
125
    {
126
      vec[j] = 0;
127
      for(k=step; k<n;k++)
128
        vec[j] += SV->array.SZ[i][k] * U[k-step][j-step];	
129
    }
130
    for(j=0; j<step; j++)
131
      vec[j] -= vec[step] * SV->array.SZ[no][j];
132
    for(j=0; j<n;j++)
133
      SV->array.SZ[i][j] = vec[j];
134
  }
135
  for(i=no+1;i<SV->rows;i++)
136
  {
137
    for(j=0; j<step; j++)
138
      vec[j] = SV->array.SZ[i][j];
139
    for(j=step;j<n; j++)
140
    {
141
      vec[j] = 0;
142
      for(k=step; k<n;k++)
143
        vec[j] += SV->array.SZ[i][k] * U[k-step][j-step];	
144
    }
145
    for(j=0; j<step; j++)
146
      vec[j] -= vec[step] * SV->array.SZ[no][j];
147
    for(j=0; j<n;j++)
148
      SV->array.SZ[i][j] = vec[j];
149
  }
150
  for(i=0; i<n; i++)
151
    SV->array.SZ[no][i] = 0;
152
  SV->array.SZ[no][step] = 1;
153
  free(vec);
154
}
155

156
static int **trafoinverse(vec, step, dim)
157
int *vec;
158
int step, dim;
159
{
160
  int i,j, s, k, v, w, udet, merk;
161
  int **U;
162
  int n;
163

164
  n = dim - step;
165
  U = (int **)malloc(n *sizeof(int *));
166
  for(i=0; i<n; i++)
167
   U[i] = (int *)malloc(n *sizeof(int));
168

169
  for(i=0; i<n;i++)
170
   U[0][i] = vec[i+step];
171
  for(k=0;k<n && U[0][k] == 0; k++);
172
  for(i=1; i<=k;i++)
173
   for(j=0;j<n;j++)
174
     U[i][j] = 0;
175
  for(i=1; i<=k;i++)
176
     U[i][i-1] = 1;
177
  udet = U[0][k];
178
  for(s=k+1; s<n;s++)
179
  {
180
    if(U[0][s] == 0)
181
    {
182
       for(i=0; i<s; i++)
183
        U[s][i] = 0;
184
       for(i=s+1; i<n; i++)
185
        U[s][i] = 0;
186
       U[s][s] = 1;
187
    }
188
    else
189
    {
190
       gcd_darstell(udet, U[0][s], &v, &w, &merk);
191
       U[s][s] = v;
192
       for(i=0; i<s; i++)
193
         U[s][i] = -(U[0][i]/udet) * w;
194
       udet = merk;
195
       for(i=s+1; i<n; i++)
196
        U[s][i] = 0;
197
    }
198
  }
199
  return(U);
200
}
201

202
static int wechsel(vec, step, udim)
203
int *vec;
204
int step, udim;
205
{
206
  int i, j;
207
  int a,b,c;
208
  i=step;
209
  while(i<udim && vec[i] == 0)
210
     i++;
211
  if(i== udim)
212
    return(FALSE);
213
  a = vec[i];
214
  for(j=i+1; j<udim && a != 1 && a!= -1; j++)
215
  {
216
    if(vec[j] != 0)
217
    {
218
      b = vec[j];
219
      while ((c=a%b)!=0){a=b; b=c;}
220
      a = b;
221
    }
222
  }
223
  if(a == 1 || a == -1)
224
    return(TRUE);
225
  return(FALSE);
226
}
227

228

229

230
/**************************************************************************\
231
@---------------------------------------------------------------------------
232
@ matrix_TYP *mink_red(G, Trf)
233
@ matrix_TYP *G, *Trf;
234
@
235
@ calculates a matrices A and Trf such that A = Trf * G * Trf^{tr}
236
@ is Minkowski_reduced
237
@---------------------------------------------------------------------------
238
@
239
\**************************************************************************/
240
matrix_TYP *mink_red(G, Trf)
241
matrix_TYP *G, *Trf;
242
{
243
  int i, j, k, l, step, anz;
244
  int n, bound, merk;
245
  int **U, **Uinv;
246
  int ergmax;
247
  matrix_TYP *erg;
248
  matrix_TYP *kurvecs;
249
  matrix_TYP *UU, *UUi;
250
  matrix_TYP *T, *Ti;
251

252
  extern int **trafoinverse();
253
  extern int wechsel();
254
  extern void Gramtrans();
255
  extern void kurvectrans();
256
  extern int matrixinverses();
257

258
  n = G->cols;
259

260
  /* ------------------------------------------------------------------- *\
261
  |     setze T = I_n  und erg = G                                        |
262
  \* ------------------------------------------------------------------- */
263
  T = init_mat(n,n,"1");
264
  erg = init_mat(G->rows, G->cols, "");
265
  for(i=0;i<G->rows;i++)
266
    for(j=0;j<G->cols;j++)
267
      erg->array.SZ[i][j] = G->array.SZ[i][j];
268
  /* ------------------------------------------------------------------- *\
269
  |     ordne erg tolal                                                     |
270
  \* ------------------------------------------------------------------- */
271
  for(i=0; i<n; i++)
272
  {
273
    k=i;
274
    for(j=i+1; j<n; j++)
275
    {
276
      if(erg->array.SZ[j][j] < erg->array.SZ[k][k])
277
        k=j;
278
    }
279
    if(k!=i)
280
    {
281
       for(j=0; j<n; j++)
282
       {
283
          l=erg->array.SZ[i][j];
284
          erg->array.SZ[i][j] = erg->array.SZ[k][j];
285
          erg->array.SZ[k][j] = l;
286
       }
287
       for(j=0; j<n; j++)
288
       {
289
          l=erg->array.SZ[j][i];
290
          erg->array.SZ[j][i] = erg->array.SZ[j][k];
291
          erg->array.SZ[j][k] = l;
292
          l=T->array.SZ[j][i];
293
          T->array.SZ[j][i] = T->array.SZ[j][k];
294
          T->array.SZ[j][k] = l;
295
       }
296
    }
297
  }
298

299
  /* ------------------------------------------------------------------- *\
300
  |     Allocieren und Initialisierung des Speicherplatzes                |
301
  \* ------------------------------------------------------------------- */
302
  erg->flags.Symmetric = TRUE;
303
  if(n == 0 || n == 1)
304
     return(erg);
305

306
  if((UUi = (matrix_TYP *)malloc(sizeof(matrix_TYP))) == NULL){
307
    printf("malloc of 'UUi' in 'mink_red' failed\n");
308
    exit(2);
309
  }
310
  UUi->kgv = 1;
311
  ergmax = erg->array.SZ[n-1][n-1];
312
  kurvecs = short_vectors(erg, ergmax, 0, 0,0,&anz);
313

314
  for(step = 0; step < n;step++)
315
  {
316
    bound = erg->array.SZ[step][step];
317
    for(i=0; i<kurvecs->rows; i++)
318
    {
319
     if(kurvecs->array.SZ[i][n] < bound && wechsel(kurvecs->array.SZ[i],step,n)== 1)
320
     {
321
        Uinv = trafoinverse(kurvecs->array.SZ[i], step, n);
322
        UUi->cols = UUi->rows = n-step;
323
        UUi->array.SZ = Uinv;
324
        UU = long_mat_inv(UUi);
325
        U = UU->array.SZ;
326
        Gramtrans(erg, T, U, Uinv, kurvecs->array.SZ[i], step);
327
        kurvectrans(kurvecs, U, i, step);
328
        bound = erg->array.SZ[step][step];
329
        for(j=0; j<(n-step); j++)
330
          free(Uinv[j]);
331
        free_mat(UU); free(Uinv);
332
     }
333
    }
334
  }
335

336
  /* printf("step = %d\n", step);
337
  put_mat(erg, NULL, "erg in mink_red", 2); */
338

339
  for(i=1; i<n;i++)
340
  {
341
    if(erg->array.SZ[i][i-1] < 0)
342
    {
343
      for(j=0; j<n; j++)
344
      {
345
         erg->array.SZ[i][j] = - erg->array.SZ[i][j];
346
         erg->array.SZ[j][i] = - erg->array.SZ[j][i];
347
         T->array.SZ[j][i] = - T->array.SZ[j][i];
348
      }
349
    }
350
  }
351
  erg->kgv =  G->kgv;
352
  if(erg->kgv != 1)
353
    erg->flags.Integral = FALSE;
354

355
  /* changed 16/1/97 tilman from
356
  free(kurvecs);
357
  to: */
358
  free_mat(kurvecs);
359

360
  free(UUi);
361
  Ti = long_mat_inv(T);
362
  free_mat(T);
363
  for(i=0;i<n;i++)
364
    for(j=0;j<n;j++)
365
       Trf->array.SZ[i][j] = Ti->array.SZ[i][j];
366
  free_mat(Ti);
367

368
  /* put_mat(erg, NULL, "erg in mink_red", 2); */
369

370
  return(erg);
371
}
372

373
Product

Resources

Company
