CoCalc -- p_solve

GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
gap4r8 / pkg / carat / carat-2.1b1 / functions / Matrix / p_solve_mat.c
⁵⁶³⁶⁸⁰ views
1
#include "typedef.h"
2
#include "tools.h"
3
#include "matrix.h"
4

5
/**************************************************************************\
6
@---------------------------------------------------------------------------
7
@---------------------------------------------------------------------------
8
@ FILE: p_solve.c
9
@---------------------------------------------------------------------------
10
@---------------------------------------------------------------------------
11
@
12
\**************************************************************************/
13

14
/**************************************************************************\
15
@---------------------------------------------------------------------------
16
@ matrix_TYP **p_solve (anz, L_mat, R_mat, option)
17
@ int *anz;
18
@ matrix_TYP **L_mat, **R_mat ;
19
@ int option;
20
@
21
@ Solves L_mat[i] * X = X * R_mat[i] modulo act_prime simultaneous
22
@  for 0 <= i<= anz.
23
@ The return is a basis of the space of solution.
24
@ the dimension of this space is returned by the pointer anz.
25
@ 
26
@---------------------------------------------------------------------------
27
@
28
\**************************************************************************/
29

30
/*{{{}}}*/
31
/*{{{  p_solve*/
32
matrix_TYP **p_solve (anz, L_mat, R_mat, option)
33
int *anz;
34
matrix_TYP **L_mat, **R_mat ;
35
int option;
36
{
37
/*{{{  variables*/
38
int **M, **N, *v, f, help, flag ;
39
int *M_j, *M_act, numax;
40
matrix_TYP **E;
41
boolean L_null, R_null ;
42
int *ss, *nuin;
43
int i, j, jj, k, kk, l, cR, cL, cM, rR, rL, rM, rN, rg = 0, act, p;
44
int **TS;
45

46
/*}}}  */
47
  /*{{{  var init & error-check*/
48

49
  p = act_prime;
50
  L_null = (L_mat == NULL);
51
  R_null = (R_mat == NULL);
52

53
  cL = L_null ? 1 : L_mat[0]->cols;
54
  rL = L_null ? 1 : L_mat[0]->rows;
55
  rR = R_null ? 1 : R_mat[0]->rows;
56
  cR = R_null ? 1 : R_mat[0]->cols;
57

58
  if(!(L_null || R_null)) {
59
    if((cL != rL) || (rR != cR)) {
60
      fprintf(stderr,"Error in input: Matrizes of wrong dimension\n");
61
      exit(3);
62
    }
63
  }
64

65
  cM = cL * rR;
66
  rM = *anz*rL*cR;
67
  rN = rL * cR;
68
  rg = *anz*cR;    
69

70
  /*}}}  */
71
  /*{{{  alloc memory*/
72
  N    = (int **)malloc(rN*sizeof(int *));
73
  M    = (int **)malloc(rM*sizeof(int *));
74
  nuin = (int *) malloc((cM+1)*sizeof(int));
75
  ss   = (int *) malloc(cM*sizeof(int));
76

77
  /*}}}  */
78
  
79
  for ( i = 0  ; i < *anz; i++) {
80
  /*{{{  */
81
    if ( !L_null ) {
82
      if ( (L_mat[i]->cols != cL) || (L_mat[i]->rows != rL)) {
83
        fprintf (stderr, "Error in input\n");
84
        exit (3);
85
      }
86
    }
87
    if ( !R_null ) {
88
      if (( R_mat[i]->rows != rR ) || ( R_mat[i]->cols != cR )) {
89
        fprintf (stderr, "Error in input\n");
90
        exit (3);
91
      }
92
    }
93

94
    for ( j = 0; j < rN; j++ ) {
95
      N[j] = (int *)calloc(cM,sizeof(int));
96
    }
97

98
    if ( !L_null ) {
99
      /*{{{  */
100
      flag = L_mat[i]->prime != 0;
101
      TS = L_mat[i]->array.SZ;
102
      for ( j = 0, jj = 0; j < rL; j++, jj += rR) {
103
        for ( k = 0, kk = 0; k < cL; k++, jj -= rR) {
104
          if(flag) {
105
            help = TS[j][k];
106
          } else if ((help = TS[j][k] % p) < 0) {
107
            help += p;
108
          }
109
          if(help) {
110
            for( l = 0; l < rR; l++,jj++, kk++) {
111
              N[jj][kk] = help;                 
112
            }
113
          } else {
114
            jj += rR;
115
            kk += rR;
116
          }
117
        } 
118
      }
119
      /*}}}  */
120
    }
121
    if ( !R_null ) {
122
      /*{{{  */
123
      flag = R_mat[i]->prime != 0;
124
      TS = R_mat[i]->array.SZ;
125
      for ( j = 0; j < cR; j++) {
126
        for ( k = 0; k < rR; k++) {
127
          if(flag) {
128
            help = TS[k][j];
129
          } else if ((help = TS[k][j] % p) < 0) {
130
            help += p;
131
          }
132
          if(help) {
133
            for(l=0,jj=j, kk=k; l<cL; l++, jj+=cR, kk+=rR) {
134
              N[jj][kk] = S(N[jj][kk],-help);
135
            }
136
          }
137
        }
138
      }
139
      /*}}}  */
140
    }
141

142
    /* permutierte Eingabe in Liste */
143
    for(j = 0; j < rL; j++) {
144
      for(k = 0; k < cR; k++) {
145
        M[j*rg+i*cR+k] = N[(rL-1-j)*cR+k];
146
      }
147
    }
148
  /*}}}  */
149
  }
150
  free(N);
151
  rg = 0;
152
    /*------------------------------------------------------------*\
153
    | Gauss Algorithm            |
154
    \*------------------------------------------------------------*/
155
  for ( i = 0; i < cM; i++ ) {
156
#if 0
157
  /*{{{  auskommentiert*/
158
    for(pri = 0; pri < rM; pri++) {
159
      for(prj = 0; prj < cM; prj++)
160
        if(M[pri][prj] != 0)
161
          printf("%d ",M[pri][prj]);
162
        else
163
          printf("  ");
164
        printf("\n");
165
      }
166
        printf(" i = %d \n", i);
167
        /* ende der Einfuegung */
168
  /*}}}  */
169
#endif
170
    ss[i] = -1;
171
    /*---------------------------------------------------------*\
172
    | Find the row with non-zero entry in i-th column           |
173
    \*---------------------------------------------------------*/
174
    for ( j = rg; (j < rM) && (M[j][i] == 0); j++);
175
    if ( j == rM ) {
176
      continue;     
177
    }
178
    act = j;
179
    /*---------------------------------------------------------*\
180
    | Normalize act-th row and mark the non-zero entries        |
181
    \*---------------------------------------------------------*/
182
    nuin[0] = 1;
183
    f = P(1,-M[j][i]);
184
    for ( k = i ; k < cM; k++ ) {
185
      if((help = M[act][k])) {
186
        M[act][k] = P(help,f);
187
        nuin[nuin[0]++] = k;
188
      }                         
189
    }
190
    /*---------------------------------------------------------*\
191
    | Swap act-th row and rg-th row                  |
192
    \*---------------------------------------------------------*/
193
    v      = M[rg];
194
    M[rg]  = M[act];
195
    M[act] = v;
196
    ss[rg] = i;
197
    act    = rg ++;
198
    /*---------------------------------------------------------*\
199
    | Clear i-th column  downwards                |
200
    \*---------------------------------------------------------*/
201
    M_act = M[act];
202
    for ( j = act+1; j < rM; j++) {
203
      if ( (f = S(0,-M[j][i])) != 0 ) {
204
        M_j = M[j];
205
        M_j[i] = 0;
206
        numax = nuin[0];
207
        if( f == 1) {
208
          for(k=2; k < numax; ++k ) {
209
	    kk = nuin[k];
210
            M_j[kk] = S(M_j[kk],M_act[kk]);
211
	  }
212
        } else {
213
          for(k=2; k < numax; ++k ) {
214
	    kk = nuin[k];
215
            M_j[kk] = S(M_j[kk],P(M_act[kk],f));           
216
          }
217
        }
218
      }
219
    }
220
  }
221
  
222
  if ( (*anz = cM - rg) != 0 ) {
223
  /*{{{  */
224
    /*========================================================*\
225
    || Matrix has not full rank: clear it upwards             ||
226
    \*=========================================================*/
227
    for (i = rg-1; i > 0; i--) {
228
      /*{{{  */
229
      nuin[0] = 2;
230
      nuin[1] = ss[i];
231
      M_act = M[i];
232
      for(j = ss[i]+1; j < cM; j++) {
233
        if(M_act[j]) {
234
          nuin[nuin[0]++] = j;
235
        }
236
      }
237
      if(nuin[0] == 2) {
238
        j = ss[i];
239
        for (k = i-1; k > 0; k--) {
240
          M[k][j] = 0;            
241
        }
242
      } else {
243
        for ( j = i-1; j >= 0; j--) {
244
          if ( (f = S(0,-M[j][ss[i]])) != 0 ) {
245
            M_j = M[j];
246
            M_j[ss[i]] = 0;
247
            numax = nuin[0];
248
            if( f == 1) {
249
              for(k=2; k < numax; ++k ) {
250
		kk = nuin[k];
251
                M_j[kk] = S(M_j[kk],M_act[kk]);                
252
              }
253
            } else {
254
              for(k=2; k < numax; ++k ) {
255
		kk = nuin[k];
256
                M_j[kk] = S(M_j[kk],P(M_act[kk],f));           
257
              }
258
            }
259
          }
260
        }
261
      }
262
      /*}}}  */
263
    }
264
    if((option & 2) || (!L_null && !R_null) ) {
265
      /*{{{  */
266
      E = (matrix_TYP **)malloc(*anz*sizeof(matrix_TYP *));
267
      for ( i = 0 , k = 0, l = 0; i < cM; i++ ) {
268
        if ( i == ss[k] ) {
269
          k ++;
270
          continue;
271
        }
272
        E[l] = init_mat(cL, rR, "ip");
273
        E[l]->prime = act_prime;
274
        N = E[l]->array.SZ;
275
        for ( j = 0; j < k; j ++) {
276
          if ( M[j][i] != 0 ) {
277
            N[ss[j]/rR][ss[j]%rR] = p-M[j][i];
278
          }
279
        }
280
        N[i/rR][i%rR] = 1;
281
        l ++;
282
      }
283
      /*}}}  */
284
    } else {
285
      /*{{{  */
286
      E = (matrix_TYP **)malloc(1*sizeof(matrix_TYP *));
287
      E[0] = init_mat(*anz, cM,"ip");
288
      E[0]->prime = act_prime;
289
      N = E[0]->array.SZ;
290
      for ( i = 0 , k = 0, l = 0; i < cM; i++ ) {
291
        if ( i == ss[k] ) {
292
          k ++;
293
          continue;
294
        }
295
        for ( j = 0; j < k; j ++) {
296
          if ( M[j][i] != 0 ) {
297
            N[l][ss[j]] = p-M[j][i];
298
          }
299
        }
300
        N[l][i] = 1;
301
        l ++;
302
      }
303
      *anz = 1;
304
      /*}}}  */
305
    }
306
  /*}}}  */
307
  } else {
308
    E = NULL;
309
  }
310
  
311
  /*{{{  free memory*/
312
  for (i = 0; i < rM; i++ ) {
313
    free(M[i]);
314
  }
315
  free(M);
316
  free(ss);
317
  free(nuin);
318
  /*}}}  */
319
  
320
  return (E);
321
}
322
/*}}}  */
323

324
Product

Resources

Company
