CoCalc -- calculate_representatives.c

GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
gap4r8 / pkg / carat / carat-2.1b1 / functions / Graph / calculate_representatives.c
⁵⁶³⁵⁸⁰ views
1
/* last change: 24.05.2002 by Oliver Heidbuechel */
2

3

4
#include <ZZ.h>
5
#include <typedef.h>
6
#include <presentation.h>
7
#include <matrix.h>
8
#include <bravais.h>
9
#include <base.h>
10
#include <datei.h>
11
#include <graph.h>
12
#include <gmp.h>
13
#include <zass.h>
14
#include <tools.h>
15
#include <longtools.h>
16
#include <sort.h>
17

18

19
/* ================================================================= */
20
/* Der Quellcode wurde normalop.c (identify) entnommen!!!            */
21
/* Stelle G->normal auf H^1(G, Q^n/Z^n) dar.                         */
22
/* ================================================================= */
23
/* G: Gruppe							     */
24
/* H_G_Z: H^1(G, Q^n/Z^n)                                            */
25
/* ================================================================= */
26
static matrix_TYP **stelle_normalisator_dar(bravais_TYP *G,
27
                                            matrix_TYP **H_G_Z)
28
{
29
  matrix_TYP **N;
30

31
  int i;
32

33
  /* put_bravais(G, 0, 0); */
34

35
  N = (matrix_TYP **)calloc(G->normal_no, sizeof(matrix_TYP *));
36
  for (i = 0; i < G->normal_no; i++){
37
     N[i] = normalop(H_G_Z[0], H_G_Z[1], H_G_Z[2], G, G->normal[i], i == (G->normal_no - 1));
38
     if (INFO_LEVEL & 16){
39
        put_mat(N[i], NULL, "N[i]", 2);
40
     }
41
  }
42

43
  return(N);
44
}
45

46
/*
47
Uebersetzungstabelle:
48

49
cozykle:	H_G_Z[0]
50
D:		H_G_Z[1]
51
R:		H_G_Z[2]
52

53
*/
54

55

56

57
/* ================================================================= */
58
/* Der Quellcode wurde dem Programm Extensions entnommen!!!          */
59
/* Berechne die Ordnung von H^1(G,Q^n/Z^n)!                          */
60
/* ================================================================= */
61
/* H_G_Z: H^1(G, Q^n/Z^n)                                            */
62
/* ================================================================= */
63
static int CohomSize(matrix_TYP **H_G_Z)
64
{
65
  int n;
66

67
  MP_INT cohom_size;
68

69

70
  if (H_G_Z[0]->cols <1){
71
     return(1);
72
  }
73
  else {
74
     cohom_size = cohomology_size(H_G_Z[1]);
75
     n = mpz_get_ui(&cohom_size);
76
     mpz_clear(&cohom_size);
77
     return(n);
78
  }
79
}
80

81

82

83

84
/* ================================================================= */
85
/* Berechne den Stabilisator von elem + Ker phi                      */
86
/* ================================================================= */
87
/* H1_mod_ker: Daten ueber H^1 / ker (phi)                           */
88
/* S1: Stabilisator                                                  */
89
/* new_rep: representation from S1 on H^1(G,Q^n/lattice) / Ker(phi)  */
90
/* S1_word_no: Anzahl der Elemente von new_rep                       */
91
/* elem: Repraesentant von elem + Ker phi                            */
92
/* H_GL_Z: H^1(GL, Q^n/Z^n) isom. zu H^1(G, Q^n/L)                   */
93
/* counter: hier wird die Anzahl der Stabilisatorerzeuger gespeichert*/
94
/* ================================================================= */
95
static matrix_TYP **stab_preimage_Ker(H1_mod_ker_TYP H1_mod_ker,
96
                                      matrix_TYP **S1,
97
                                      matrix_TYP **new_rep,
98
				      int S1_word_no,
99
				      matrix_TYP *elem,
100
				      matrix_TYP **H_GL_Z,
101
				      int *counter)
102
{
103
   int reps_no, *length, m, k,
104
       ***WORDS, *WORDS_no;
105

106
   matrix_TYP ***reps, **S_ksi, *temp, **S1_inv, *std;
107

108

109

110
   /* Berechne Standardform von elem als Element von H^1/Ker */
111
   /*if (elem == NULL){
112
      temp = init_mat(H1_mod_ker.i->rows, 1, "");
113
   }
114
   else{*/
115
      temp = mat_mul(H1_mod_ker.i, elem);
116
   /*}*/
117
   for (k = 0; k < temp->rows; k++){
118
      temp->array.SZ[k][0] %= H1_mod_ker.D->array.SZ[k][k];
119
      if (temp->array.SZ[k][0] < 0)
120
         temp->array.SZ[k][0] += H1_mod_ker.D->array.SZ[k][k];
121
   }
122
   std = init_mat(H1_mod_ker.erz_no, 1, "");
123
   for (k = 0; k < H1_mod_ker.erz_no; k++){
124
      std->array.SZ[k][0] = temp->array.SZ[k + H1_mod_ker.D_first][0];
125
   }
126
   free_mat(temp);
127

128
   /* Bahn von elem sowie dessen Stabilisator berechnen */
129
   reps = H1_mod_ker_orbit_alg(H1_mod_ker, new_rep, S1_word_no, &reps_no,
130
                              &length, &WORDS, &WORDS_no, std);
131
   free_mat(std);
132
   if (reps_no != 2 || length[1] < 1){
133
      fprintf(stderr, "ERROR 1 in stab_preimage_Ker!\n");
134
      exit(78);
135
   }
136

137
   /* Berechne Stabilisator */
138
   S1_inv = (matrix_TYP **)calloc(S1_word_no, sizeof(matrix_TYP *));
139
   counter[0] = 0;
140
   S_ksi = (matrix_TYP **)calloc(WORDS_no[1], sizeof(matrix_TYP *));
141
   for (m = 0; m < WORDS_no[1]; m++){
142
      /* simplify the words */
143
      normalize_word(WORDS[1][m]);
144
      if (WORDS[1][m][0] == 1 && WORDS[1][m][1] < 0)
145
         WORDS[1][m][1] *= (-1);
146
      if (word_already_there(WORDS[1], m) == 0){
147
         S_ksi[counter[0]] = graph_mapped_word(WORDS[1][m], S1, S1_inv, H_GL_Z[1]);
148
         if (mat_search(S_ksi[counter[0]], S_ksi, counter[0], mat_comp) != -1){
149
            free_mat(S_ksi[counter[0]]);
150
         }
151
         else{
152
            mat_quicksort(S_ksi, 0, counter[0], mat_comp);
153
            counter[0]++;
154
         }
155
      }
156
   }
157

158
   /* clean */
159
   /* ===== */
160
   for (m = 0; m < S1_word_no; m++){
161
      if (S1_inv[m] != NULL)
162
         free_mat(S1_inv[m]);
163
   }
164
   free(S1_inv);
165
   for (m = 0; m < WORDS_no[1]; m++){
166
      free(WORDS[1][m]);
167
   }
168
   for (m = 0; m < length[1]; m++){
169
      free_mat(reps[1][m]);
170
   }
171
   free(reps[1]);
172
   free(WORDS[1]);
173
   free(reps);
174
   free(WORDS);
175
   free(WORDS_no);
176
   free(length);
177

178
   return(S_ksi);
179
}
180

181

182

183

184
/* ================================================================= */
185
/* Der Quellcode wurde subgroupgraph.c entnommen!!!                  */
186
/* Berechne Vertreter der Bahnen der Untergruppen unter dem          */
187
/* Stab des Gitters geschnitten mit dem Stabilisator des Cocykels    */
188
/* als Element von H^1(G, Q^n/L)!!!                                  */
189
/* ================================================================= */
190
/* G: Punktruppe                                                     */
191
/* GL: G konjugiert mit Gitter lattice                               */
192
/* H_G_Z: H^1(G, Q^n/Z^n)                                            */
193
/* H_GL_Z: H^1(GL, Q^n/Z^n) isom. zu H^1(G, Q^n/L)                   */
194
/* lattice: Gitter                                                   */
195
/* phi: H^1(G, Q^n/L) isom. zu H^1(GL,Q^n,Z^n) -> H^1(G,Q^n/Z^n)     */
196
/* H_1_mod_ker: Daten ueber H^1 / ker (phi)                          */
197
/* preimage: Element aus H^1(GL, Q^n/Z^n), welches auf Cozykel zur   */
198
/*           Raumgruppe abgebildet wird.                             */
199
/* kernel_mat: Matrix, die den Kern von phi beschreibt               */
200
/* orbit_no: speichere die Anzahl der Vertreter hier                 */
201
/* orbit_length: speichere die Laengen der Bahnen hier               */
202
/* ================================================================= */
203
matrix_TYP **calculate_representatives(bravais_TYP *G,
204
                                       bravais_TYP *GL,
205
				       matrix_TYP **H_G_Z,
206
				       matrix_TYP **H_GL_Z,
207
				       matrix_TYP *lattice,
208
                                       matrix_TYP *phi,
209
                                       H1_mod_ker_TYP H1_mod_ker,
210
				       matrix_TYP *preimage,
211
				       matrix_TYP *kernel_mat,
212
				       int *orbit_no,
213
				       int **orbit_length)
214
{
215
   matrix_TYP **S1,
216
              **S1_inv,
217
              **lNli,
218
	      **Ninv,
219
	      **new_rep,
220
	      *invlattice,
221
	      **N_H_GL_Z,
222
	      *id,
223
	      **S_ksi,
224
              **kernel_gen,
225
              **kernel_elements,
226
	      **orbit_rep,
227
	      **rep;
228

229
   int l,
230
       norm_no, S1_word_no, cohomSize, 
231
       kernel_order, S_ksi_no,
232
       *kernel_list;
233

234
   rational eins;
235

236

237

238
   eins.z = eins.n = 1;
239

240

241
   /* calculate S1 = Stab_N(L) */
242
   /* ======================== */
243
   if (phi->cols > 0){
244
      norm_no = GL->normal_no;
245
      invlattice = mat_inv(lattice);
246
      id = init_mat(G->dim, G->dim, "1");
247
      Ninv = (matrix_TYP **)calloc(norm_no, sizeof(matrix_TYP *));
248
      lNli = (matrix_TYP **)calloc(norm_no, sizeof(matrix_TYP *));
249
      for (l = 0; l < norm_no; l++){
250
         lNli[l] = mat_kon(lattice, GL->normal[l], invlattice);
251
      }
252
      free_mat(invlattice);
253
      N_H_GL_Z = stelle_normalisator_dar(GL, H_GL_Z);
254
      S1 = calculate_S1(id, lNli, norm_no, &S1_word_no, N_H_GL_Z, Ninv, H_GL_Z[1]);
255

256
      S1_inv = (matrix_TYP **)calloc(S1_word_no, sizeof(matrix_TYP *));
257
      for (l = 0; l < norm_no; l++){
258
         free_mat(lNli[l]);
259
	 free_mat(N_H_GL_Z[l]);
260
	 if (Ninv[l] != NULL)
261
	    free_mat(Ninv[l]);
262
      }
263
      free_mat(id);
264
      free(lNli);
265
      free(Ninv);
266
      free(N_H_GL_Z);
267
   }
268
   else{
269
      /* trivial part */
270
      S1_word_no = 0;
271
   }
272

273
   cohomSize = CohomSize(H_GL_Z);
274

275
   /* Kern */
276
   kernel_gen = col_to_list(kernel_mat);
277
   kernel_elements = (matrix_TYP **)calloc(cohomSize, sizeof(matrix_TYP *));
278
   kernel_list = aufspannen(cohomSize, kernel_elements, kernel_gen,
279
                            kernel_mat->cols, H_GL_Z[1], &kernel_order);
280

281

282
   if (H1_mod_ker.flag == 0 && H1_mod_ker.erz_no > 0){
283
      if (preimage == NULL){
284
         rep = orbit_ker(kernel_elements, kernel_order, H_GL_Z[1], S1,
285
                         S1_word_no, cohomSize, kernel_list, orbit_no, orbit_length);
286
      }
287
      else{
288
         /* representation from S1 on H^1(G,Q^n/lattice) / Ker(phi) */
289
         new_rep = new_representation(S1, S1_word_no, H1_mod_ker, H_GL_Z[1]);
290

291
         /* affine representation of the kernel */
292
         kernel_elements_2_affine(kernel_elements, cohomSize);
293

294
         /* calculate Stab_S1 (preimage + Kern) */
295
         S_ksi = stab_preimage_Ker(H1_mod_ker, S1, new_rep, S1_word_no, preimage, H_GL_Z, &S_ksi_no);
296

297
         /* calculate representatives */
298
         orbit_rep = orbit_ksi_plus_ker(preimage, kernel_elements, kernel_order,
299
                                        H_GL_Z[1], S_ksi, S_ksi_no, cohomSize,
300
	   			        kernel_list, orbit_no, orbit_length);
301
         rep = (matrix_TYP **)calloc(orbit_no[0], sizeof(matrix_TYP *));
302
         for (l = 0; l < orbit_no[0]; l++){
303
            rep[l] = mat_add(orbit_rep[l], preimage, eins, eins);
304
	    free_mat(orbit_rep[l]);
305
	 }
306
	 free(orbit_rep);
307
         for (l = 0; l < S_ksi_no; l++)
308
            free_mat(S_ksi[l]);
309
         free(S_ksi);
310
         for (l = 0; l < S1_word_no; l++)
311
            free_mat(new_rep[l]);
312
         free(new_rep);
313
      }
314
   }
315
   else{
316
      switch (H1_mod_ker.flag){
317
         case 0:
318
         case 2:
319
            /* orbit on 0 + Ker(phi) */
320
            rep = orbit_ker(kernel_elements, kernel_order, H_GL_Z[1], S1,
321
                            S1_word_no, cohomSize, kernel_list, orbit_no, orbit_length);
322
            break;
323
         case 1:
324
         case 3:
325
            rep = (matrix_TYP **)calloc(1, sizeof(matrix_TYP *));
326
	    rep[0] = init_mat(0, 0, "");
327
	    orbit_length[0] = (int *)calloc(1, sizeof(int));
328
	    orbit_length[0][0] = 1;
329
	    orbit_no[0] = 1;
330
            break;
331
	 default:
332
            fprintf(stderr, "ERROR 1 in cacluate_representative!\n");
333
            exit(55);
334
      }
335
   }
336

337

338
   /* clean */
339
   /* ===== */
340
   if (phi->cols > 0){
341
      for (l = 0; l < S1_word_no; l++){
342
         free_mat(S1[l]);
343
         if (S1_inv[l] != NULL)
344
            free_mat(S1_inv[l]);
345
       }
346
       free(S1);
347
       free(S1_inv);
348
   }
349
   for (l = 0; l < cohomSize; l++){
350
      if (kernel_elements[l] != NULL)
351
         free_mat(kernel_elements[l]);
352
   }
353
   free(kernel_elements);
354
   free(kernel_list);
355
   for (l = 0; l < kernel_mat->cols; l++){
356
      free_mat(kernel_gen[l]);
357
   }
358
   free(kernel_gen);
359

360

361
   return(rep);
362
}
363

364

365

366

367
/*
368
Uebersetzungstabelle:
369
=====================
370

371
data->Z[j]		GL
372
data->Z[i]		G
373
data->X[j]		H_GL_Z
374
data->X[i]		H_G_Z
375
data->N[j]		N_H_GL_Z
376

377

378
*/
379

380

381
Product

Resources

Company
