CoCalc -- add

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

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

14
/**************************************************************************\
15
@---------------------------------------------------------------------------
16
@ matrix_TYP *imat_add (L_mat, R_mat, Lc, Rc)
17
@ matrix_TYP *L_mat, *R_mat ;
18
@ int Lc, Rc;
19
@ 
20
@  calculates Lc * L_mat + Rc * R_mat
21
@---------------------------------------------------------------------------
22
@
23
\**************************************************************************/
24
/*{{{}}}*/
25
/*{{{  imat_add, exported*/
26
matrix_TYP *imat_add (L_mat, R_mat, Lc, Rc)
27
matrix_TYP *L_mat, *R_mat ;
28
int Lc, Rc;
29
{
30
matrix_TYP *S_mat;
31
int **L, **R, **S, rS, cS;
32
int i, j;
33

34
  rS = L_mat->rows;
35
  cS = L_mat->cols;
36
  S_mat = init_mat(rS,cS,"");
37
  
38
  S_mat->flags.Symmetric=(R_mat->flags.Symmetric&&L_mat->flags.Symmetric);
39
  S_mat->flags.Diagonal =(R_mat->flags.Diagonal && L_mat->flags.Diagonal);
40
  S_mat->flags.Scalar   =(R_mat->flags.Scalar   && L_mat->flags.Scalar);
41
  S_mat->flags.Integral =(R_mat->flags.Integral && L_mat->flags.Integral);
42

43
  /*
44
   *  ovfl_mul() is an expensive routine defined in ../tools/ovfl_mul.c
45
   *  It checks an integer-overflow using a division :(.
46
   *  But as the kgv should tend to grow fast, this seems to be saver
47
   */
48
  if(L_mat->kgv != 1) {
49
    Rc = ovfl_mul( Rc, L_mat->kgv );
50
  }
51
  if(R_mat->kgv != 1) {
52
    Lc = ovfl_mul( Lc, R_mat->kgv );
53
  }
54
  S_mat->kgv = ovfl_mul( L_mat->kgv , R_mat->kgv );
55
  
56
  S = S_mat->array.SZ;
57
  R = R_mat->array.SZ;
58
  L = L_mat->array.SZ;
59
#ifdef TEST_OVFL
60
  for ( i  = 0 ; i < rS; i++) {
61
    if(S_mat->flags.Diagonal) {
62
      S[i][i] = ovfl_mul( Lc , L[i][i] ) + ovfl_mul( Rc , R[i][i] );
63
    } else {
64
      for ( j  = 0 ; j  < cS  ; j ++  ) {
65
        S[i][j] = ovfl_mul( Lc , L[i][j] ) + ovfl_mul( Rc , R[i][j] );
66
      }
67
    }
68
  } 
69
#else
70
  for ( i  = 0 ; i < rS; i++) {
71
    if(S_mat->flags.Diagonal) {
72
      S[i][i] = Lc * L[i][i] + Rc * R[i][i];
73
    } else {
74
      for ( j  = 0 ; j  < cS  ; j ++  ) {
75
        S[i][j] = Lc * L[i][j] + Rc * R[i][j];
76
      }
77
    }
78
  }  
79
#endif
80
  Check_mat(S_mat);
81
  return(S_mat);
82
}
83

84
/*}}}  */
85
/*{{{  imat_addeq, exported*/
86

87
/**************************************************************************\
88
@---------------------------------------------------------------------------
89
@ matrix_TYP *imat_addeq (L_mat, R_mat, Lc, Rc)
90
@ matrix_TYP *L_mat, *R_mat ;
91
@ int Lc, Rc;
92
@ 
93
@  calculates L_mat = Lc*L_mat + Rc*R_mat 
94
@
95
@---------------------------------------------------------------------------
96
@
97
\**************************************************************************/
98
matrix_TYP *imat_addeq (L_mat, R_mat, Lc, Rc)
99
matrix_TYP *L_mat, *R_mat ;
100
int Lc, Rc;
101
{
102

103
int **L, **R, rS, cS;
104
int i, j;
105

106
  rS = L_mat->rows;
107
  cS = L_mat->cols;
108

109
  L_mat->flags.Symmetric=(R_mat->flags.Symmetric&&L_mat->flags.Symmetric);
110
  L_mat->flags.Diagonal =(R_mat->flags.Diagonal && L_mat->flags.Diagonal);
111
  L_mat->flags.Scalar   =(R_mat->flags.Scalar   && L_mat->flags.Scalar);
112
  L_mat->flags.Integral =(R_mat->flags.Integral && L_mat->flags.Integral);
113

114
  /*
115
   *  ovfl_mul() is an expensive routine defined in ../tools/ovfl_mul.c
116
   *  It checks an integer-overflow using a division :(.
117
   *  But as the kgv should tend to grow fast, this seems to be saver
118
   */
119
  if(L_mat->kgv != 1) {
120
    Rc = ovfl_mul( Rc, L_mat->kgv );
121
  }
122
  if(R_mat->kgv != 1) {
123
    Lc = ovfl_mul( Lc, R_mat->kgv );
124
  }
125
  L_mat->kgv = ovfl_mul( L_mat->kgv , R_mat->kgv );
126
  
127
  R = R_mat->array.SZ;
128
  L = L_mat->array.SZ;
129
#ifdef TEST_OVFL
130
  for ( i  = 0 ; i < rS; i++) {
131
    if(L_mat->flags.Diagonal) {
132
      L[i][i] = ovfl_mul( Lc , L[i][i] ) + ovfl_mul( Rc , R[i][i] );
133
    } else {
134
      for ( j  = 0 ; j  < cS  ; j ++  ) {
135
        L[i][j] = ovfl_mul( Lc , L[i][j] ) + ovfl_mul( Rc , R[i][j] );
136
      }
137
    }
138
  } 
139
#else
140
  for ( i  = 0 ; i < rS; i++) {
141
    if(L_mat->flags.Diagonal) {
142
      L[i][i] = Lc * L[i][i] + Rc * R[i][i];
143
    } else {
144
      for ( j  = 0 ; j  < cS  ; j ++  ) {
145
        L[i][j] = Lc * L[i][j] + Rc * R[i][j];
146
      }
147
    }
148
  }  
149
#endif
150
  Check_mat(L_mat);
151
  return(L_mat);
152
}
153

154

155
/**************************************************************************\
156
@---------------------------------------------------------------------------
157
@ matrix_TYP *rmat_add (L_mat, R_mat, L_coeff, R_coeff)
158
@ matrix_TYP *L_mat, *R_mat ;
159
@ rational L_coeff, R_coeff;
160
@
161
@  calculates L_coeff * L_mat + R_coeff * R_mat
162
@---------------------------------------------------------------------------
163
@
164
\**************************************************************************/
165
/*}}}  */
166
/*{{{  rmat_add, exported*/
167
matrix_TYP *rmat_add (L_mat, R_mat, L_coeff, R_coeff)
168
matrix_TYP *L_mat, *R_mat ;
169
rational L_coeff, R_coeff;
170
{
171
matrix_TYP *S_mat;
172
int **S, rS, cS , temp1, temp2;
173
int i, j;
174
rational Lc, Rc;
175

176
  temp1 = temp2 = 0;
177
  Lc.z = L_coeff.z;
178
  Lc.n = L_coeff.n;
179
  Rc.z = R_coeff.z;
180
  Rc.n = R_coeff.n;
181
  Normal( &Lc );
182
  Normal( &Rc );
183

184
  rS = L_mat->rows;
185
  cS = L_mat->cols;
186
  S_mat = init_mat(rS,cS,"il");
187
  
188
  S_mat->flags.Symmetric=(R_mat->flags.Symmetric && L_mat->flags.Symmetric);
189
  S_mat->flags.Diagonal =(R_mat->flags.Diagonal && L_mat->flags.Diagonal);
190
  S_mat->flags.Scalar   =(R_mat->flags.Scalar   && L_mat->flags.Scalar);
191
  S_mat->flags.Integral =(R_mat->flags.Integral && L_mat->flags.Integral);
192
  if(!(S_mat->flags.Integral)) {
193
    temp1= GGT(Lc.z, L_mat->kgv);
194
    temp2= GGT(Rc.z, R_mat->kgv);
195
    Lc.z = Lc.z / temp1 * R_mat->kgv / temp2 * Rc.n;
196
    Normal(&Lc);
197
    Rc.z = Rc.z / temp2 * L_mat->kgv / temp1 * Lc.n;
198
    Normal(&Rc);
199
  
200
    S_mat->kgv = L_mat->kgv * R_mat->kgv * Lc.n * Rc.n;
201
    S_mat->kgv /= temp1;
202
    S_mat->kgv /= temp2;
203
    if (S_mat->kgv == 1) {
204
      S_mat->flags.Integral = TRUE;
205
    }
206
  }
207
  else {
208
    if((Lc.n == 1) && (Rc.n == 1)) {
209
      S_mat->kgv = 1;
210
    } else {
211
      S_mat->kgv = Lc.n * Rc.n ;
212
      S_mat->flags.Integral = FALSE;
213
      Lc.z *= Rc.n;
214
      Rc.z *= Lc.n;
215
    }
216
  }
217
  S= S_mat->array.SZ;
218
  for ( i  = 0 ; i < rS; i++) {
219
    if(S_mat->flags.Diagonal) {
220
      S[i][i] = Lc.z * L_mat->array.SZ[i][i] + Rc.z * R_mat->array.SZ[i][i];
221
    } else {
222
      for ( j  = 0 ; j  < cS  ; j ++  ) {
223
        S[i][j] = Lc.z * L_mat->array.SZ[i][j] + Rc.z * R_mat->array.SZ[i][j];
224
      }
225
    }
226
  }
227

228
  Check_mat(S_mat);
229
  return(S_mat);
230
}
231

232
/**************************************************************************\
233
@---------------------------------------------------------------------------
234
@ matrix_TYP *rmat_addeq(L_mat, R_mat, L_coeff, R_coeff)
235
@ matrix_TYP *L_mat, *R_mat ;
236
@ rational L_coeff, R_coeff;
237
@
238
@  calculates L_coeff = L_coeff * L_mat + R_coeff * R_mat
239
@---------------------------------------------------------------------------
240
@
241
\**************************************************************************/
242
/*}}}  */
243
/*{{{  rmat_addeq, exported*/
244
matrix_TYP *rmat_addeq (L_mat, R_mat, L_coeff, R_coeff)
245
matrix_TYP *L_mat, *R_mat ;
246
rational L_coeff, R_coeff;
247
{
248
int **L, **R, rS, cS , temp1, Rk;
249
int i, j;
250
rational Rc, Lc;
251

252
  Rk    =
253
  temp1 = 0;
254
  rS = L_mat->rows;
255
  cS = L_mat->cols;
256
  Lc.z= L_coeff.z;
257
  Lc.n= L_coeff.n;
258
  Rc.z= R_coeff.z;
259
  Rc.n= R_coeff.n;
260
  Normal( &Lc );
261
  Normal( &Rc );
262

263
  L_mat->flags.Symmetric=(R_mat->flags.Symmetric && L_mat->flags.Symmetric);
264
  L_mat->flags.Diagonal =(R_mat->flags.Diagonal && L_mat->flags.Diagonal);
265
  L_mat->flags.Scalar   =(R_mat->flags.Scalar   && L_mat->flags.Scalar);
266
  L_mat->flags.Integral =(R_mat->flags.Integral && L_mat->flags.Integral);
267
  if(!(L_mat->flags.Integral)) {
268
    Lc.n= L_coeff.n * L_mat->kgv; 
269
    Normal(&Lc);
270
    Rc.n= R_coeff.n * R_mat->kgv; 
271
    Normal(&Rc);
272
    temp1= GGT( Lc.n, Rc.n);
273
    Lc.z= Lc.z * Rc.n / temp1;
274
    Rc.z= Rc.z * Lc.n / temp1;
275
    L_mat->kgv= Lc.n *  Rc.n / temp1;
276
    if (L_mat->kgv == 1) {
277
      L_mat->flags.Integral = TRUE;
278
    }
279
  } else {
280
    if((L_coeff.n == 1) && (R_coeff.n == 1)) {
281
      L_mat->kgv = 1;
282
    } else {
283
      temp1= GGT(Lc.n,Rc.n);
284
      L_mat->kgv = Lc.n * Rc.n / temp1;
285
      L_mat->flags.Integral = FALSE;
286
      Lc.z = Lc.z * Rc.n / temp1;
287
      Rc.z = Rc.z * Lc.n /temp1;
288
    }
289
  }
290
  L = L_mat->array.SZ;
291
  R = R_mat->array.SZ;
292
  for ( i  = 0 ; i < rS; i++) {
293
    if(L_mat->flags.Diagonal) {
294
      L[i][i] = L[i][i] * Lc.z + Rc.z * R[i][i];
295
    } else {
296
      for ( j  = 0 ; j  < cS  ; j ++  ) {
297
        L[i][j] = L[i][j] * Lc.z + Rc.z*R[i][j];
298
      }
299
    }
300
  }
301
  Check_mat(L_mat);
302
  return(L_mat);
303
}
304

305
/*}}}  */
306
/*{{{  pmat_add, exported*/
307

308
/**************************************************************************\
309
@---------------------------------------------------------------------------
310
@ matrix_TYP *pmat_add (L_mat, R_mat, L_coeff, R_coeff)
311
@ matrix_TYP *L_mat, *R_mat ;
312
@ int L_coeff, R_coeff;
313
@
314
@ calculates L_coeff * L_mat + R_coeff * R_mat modulo L_mat->prime
315
@---------------------------------------------------------------------------
316
@
317
\**************************************************************************/
318
matrix_TYP *pmat_add (L_mat, R_mat, L_coeff, R_coeff)
319
matrix_TYP *L_mat, *R_mat ;
320
int L_coeff, R_coeff;
321
{
322
matrix_TYP *S_mat;
323
int **L, **R, **A, rS, cS ;
324
int i, j;
325

326
  rS = L_mat->rows;
327
  cS = L_mat->cols;
328
  S_mat = init_mat(rS,cS,"p");
329
  S_mat->prime = L_mat->prime;
330
  S_mat->flags.Symmetric = (R_mat->flags.Symmetric&&L_mat->flags.Symmetric);
331
  S_mat->flags.Diagonal  = (R_mat->flags.Diagonal && L_mat->flags.Diagonal);
332
  S_mat->flags.Scalar    = (R_mat->flags.Scalar   && L_mat->flags.Scalar);
333
  A= S_mat->array.SZ;
334
  L= L_mat->array.SZ;
335
  R= R_mat->array.SZ;
336
  for ( i  = 0 ; i < rS; i++) {
337
    if(S_mat->flags.Diagonal) {
338
      A[i][i] = S(P(L_coeff,L[i][i]), P(R_coeff,R[i][i]));
339
    } else {
340
      for ( j  = 0 ; j  < cS  ; j ++  ) {
341
        A[i][j] = S(P(L_coeff,L[i][j]), P(R_coeff,R[i][j]));
342
      }
343
    }
344
  }
345
  Check_mat (S_mat);
346
  return(S_mat);
347
}
348

349
/*}}}  */
350
/*{{{  pmat_addeq, exported*/
351
/**************************************************************************\
352
@---------------------------------------------------------------------------
353
@ matrix_TYP *pmat_addeq (L_mat, R_mat, L_coeff, R_coeff)
354
@ matrix_TYP *L_mat, *R_mat ;
355
@ int L_coeff, R_coeff;
356
@
357
@ calculates L_mat = L_coeff * L_mat + R_coeff * R_mat modulo L_mat->prime
358
@---------------------------------------------------------------------------
359
@
360
\**************************************************************************/
361
matrix_TYP *pmat_addeq (L_mat, R_mat, L_coeff, R_coeff)
362
matrix_TYP *L_mat, *R_mat ;
363
int L_coeff, R_coeff;
364
{
365
int **L, **R, rS, cS ;
366
int i, j;
367

368
  rS = L_mat->rows;
369
  cS = L_mat->cols;
370
  L_mat->flags.Symmetric=(R_mat->flags.Symmetric&&L_mat->flags.Symmetric);
371
  L_mat->flags.Diagonal =(R_mat->flags.Diagonal && L_mat->flags.Diagonal);
372
  L_mat->flags.Scalar   =(R_mat->flags.Scalar   && L_mat->flags.Scalar);
373
  L= L_mat->array.SZ;
374
  R= R_mat->array.SZ;
375
  for ( i  = 0 ; i < rS; i++) {
376
    if(L_mat->flags.Diagonal) {
377
      L[i][i] = S(P(L_coeff,L[i][i]), P(R_coeff,R[i][i]));
378
    } else {
379
      for ( j  = 0 ; j  < cS  ; j ++  ) {
380
        L[i][j] = S(P(L_coeff,L[i][j]), P(R_coeff,R[i][j]));
381
      }
382
    }
383
  }
384
  Check_mat (L_mat);
385
  return(L_mat);
386
}
387

388
/*}}}  */
389
/*{{{  mat_add, exported*/
390

391
/**************************************************************************\
392
@---------------------------------------------------------------------------
393
@ matrix_TYP *mat_add (L_mat, R_mat, L_coeff, R_coeff)
394
@ matrix_TYP *L_mat, *R_mat ;
395
@ rational L_coeff, R_coeff;
396
@
397
@ calculates L_coeff * L_mat + R_coeff * R_mat
398
@---------------------------------------------------------------------------
399
@
400
\**************************************************************************/
401
matrix_TYP *mat_add (L_mat, R_mat, L_coeff, R_coeff)
402
matrix_TYP *L_mat, *R_mat ;
403
rational L_coeff, R_coeff;
404
{
405
int Lc, Rc;
406
int lp, help;
407

408
  help = 0;
409
  Lc = Rc = 1;
410
  if (( L_mat->cols != R_mat->cols) || (L_mat->rows != R_mat->rows )) {
411
    fprintf (stderr, "Can't add %dx%d with %dx%d\n",
412
                      L_mat->rows, L_mat->cols, R_mat->rows, R_mat->cols );
413
    exit (3);
414
  }
415
  if( L_mat->prime != 0 ) {
416
    if (L_mat->prime != R_mat->prime) {
417
      fprintf(stderr, "Error: Different primes in mat_add.\n");
418
      exit (3);
419
    }
420
    /*============================================================*\
421
    ||                                                            ||
422
    || The prime-field case:                                      ||
423
    || Initialize actuell prime and make sure that the            ||
424
    || Coeffizients are element of Fp.                            ||
425
    ||                                                            ||
426
    \*============================================================*/
427
    init_prime(L_mat->prime);
428
    lp = act_prime;
429
    if(L_coeff.n != 1) Lc = P(1,-L_coeff.n%lp);
430
    if(L_coeff.z != 1) Lc = P(Lc,L_coeff.z%lp);
431
    if(R_coeff.n != 1) Rc = P(1,-R_coeff.n%lp);
432
    if(R_coeff.z != 1) Rc = P(Rc,R_coeff.z%lp);
433
    return  pmat_add(L_mat,R_mat,Lc,Rc);
434
  } else {
435
    if ( L_mat->array.N != NULL ) {
436
      rat2kgv( L_mat );
437
    }    
438
    if ( R_mat->array.N != NULL ) {
439
      rat2kgv( R_mat );
440
    }
441
    return rmat_add(L_mat,R_mat,L_coeff, R_coeff);
442
  }
443
}
444

445
/*}}}  */
446
/*{{{  mat_addeq, exported*/
447
/**************************************************************************\
448
@---------------------------------------------------------------------------
449
@ matrix_TYP *mat_addeq (L_mat, R_mat, L_coeff, R_coeff)
450
@ matrix_TYP *L_mat, *R_mat ;
451
@ rational L_coeff, R_coeff;
452
@
453
@ calculates L_mat = L_coeff * L_mat + R_coeff * R_mat
454
@---------------------------------------------------------------------------
455
@
456
\**************************************************************************/
457
matrix_TYP *mat_addeq (L_mat, R_mat, L_coeff, R_coeff)
458
matrix_TYP *L_mat, *R_mat ;
459
rational L_coeff, R_coeff;
460
{
461
int Lc, Rc;
462
int lp, help;
463

464
  help = 0;
465
  Lc = Rc = 1;
466
  if (( L_mat->cols != R_mat->cols) || (L_mat->rows != R_mat->rows )) {
467
    fprintf (stderr, "Can't add %dx%d with %dx%d\n",
468
                      L_mat->rows, L_mat->cols, R_mat->rows, R_mat->cols );
469
    exit (3);
470
  }
471
  if( L_mat->prime != 0 ) {
472
    if (L_mat->prime != R_mat->prime) {
473
      fprintf(stderr, "Error: Different primes in mat_add.\n");
474
      exit (3);
475
    }
476
    /*============================================================*\
477
    ||                                                            ||
478
    || The prime-field case:                                      ||
479
    || Initialize actuell prime and make sure that the            ||
480
    || Coeffizients are element of Fp.                            ||
481
    ||                                                            ||
482
    \*============================================================*/
483
    init_prime(L_mat->prime);
484
    lp = act_prime;
485
    if(L_coeff.n != 1) Lc = P(1,-L_coeff.n%lp);
486
    if(L_coeff.z != 1) Lc = P(Lc,L_coeff.z%lp);
487
    if(R_coeff.n != 1) Rc = P(1,-R_coeff.n%lp);
488
    if(R_coeff.z != 1) Rc = P(Rc,R_coeff.z%lp);
489
    return pmat_addeq(L_mat,R_mat,Lc,Rc);
490
  } else {
491
    if ( L_mat->array.N != NULL ) {
492
      rat2kgv( L_mat );
493
    }    
494
    if ( R_mat->array.N != NULL ) {
495
      rat2kgv( R_mat );
496
    }
497
    return rmat_addeq(L_mat,R_mat,L_coeff, R_coeff);
498
  }
499
}
500

501
/*}}}  */
502

503
Product

Resources

Company
