1 #include "vectorsF2.h"
2 #include <stdio.h>
3 #include <stdlib.h>
4
5 #define WL vectorsF2_WL
6
7 #define MC 0x80000000UL /* permet de diagonaliser la matrice dans Diag() */
8
9 unsigned long MMC[WL] =
10 { MC, MC >> 1, MC >> 2, MC >> 3, MC >> 4, MC >> 5, MC >> 6, MC >> 7,
11 MC >> 8, MC >> 9, MC >> 10,
12 MC >> 11, MC >> 12, MC >> 13, MC >> 14, MC >> 15, MC >> 16, MC >> 17,
13 MC >> 18, MC >> 19, MC >> 20,
14 MC >> 21, MC >> 22, MC >> 23, MC >> 24, MC >> 25, MC >> 26, MC >> 27,
15 MC >> 28, MC >> 29, MC >> 30, MC >> 31
16 };
17
InverseMatrix(Matrix * InvM,Matrix * M)18 lebool InverseMatrix (Matrix * InvM, Matrix * M)
19 {
20
21 Matrix Temp;
22 int j, rang;
23 if (M->nblignes != M->l) {
24 printf ("Matrix M is not square!\n");
25 exit (1);
26 }
27 AllocMat (&Temp, M->nblignes, M->l, 2);
28 for (j = 0; j < M->l; j++)
29 CopyBV (&(Temp.lignes[j][0]), &(M->lignes[j][0]));
30 for (j = 0; j < M->l; j++) {
31 BVCanonic (&(Temp.lignes[j][1]), j);
32 }
33 /* DispMat(&Temp,2,M->l,M->nblignes,0); */
34 rang = CompleteElimination (&Temp, M->nblignes, M->l, 2);
35 /* DispMat(&Temp,2,M->l,M->nblignes,0); */
36 /* printf("rang=%d",rang); */
37 for (j = 0; j < M->l; j++)
38 CopyBV (&(InvM->lignes[j][0]), &(Temp.lignes[j][1]));
39 return (rang == M->l);
40 FreeMat (&Temp);
41 }
42
43
44 /* ********************************************************************** */
45 /* lebool Diag( Matrix m, int kg, */
46 /* int t, int l, int *gr ) */
47 /* Evalue si la matrice de travail m sur kg lignes est de plein rang en */
48 /* la diagonalisant. On procede sur t BitVect en considerant les l */
49 /* premiers bits de chacun. La fonction retourne TRUE si la matrice m */
50 /* est de plein rang t*l et *gr est inchange. La fonction retourne */
51 /* FALSE sinon et *gr prend pour valeur le numero du BitVect ou il y a */
52 /* eu echec moins un ( = dimension pour laquelle on a resolution l ). */
53 /* ********************************************************************** */
Diag(Matrix * m,int kg,int t,int l,int * gr)54 lebool Diag (Matrix * m, int kg, int t, int l, int *gr)
55 {
56 int i, j, cl, rang;
57
58 rang = 0;
59
60 /* On diagonalise la matrice des entrees (i,j) sur l bits */
61 /* avec 0 <= i < kg et 0 <= j < t . */
62
63 for (j = 0; j < t; j++) {
64 cl = 1;
65 while (cl <= l) {
66 /* On cherche dans la j-eme colonne, commencant a la */
67 /* rang-eme ligne, la premiere entree dont le bit le plus */
68 /* significatif est non nul. Bref, on cherche un pivot. */
69 i = rang;
70
71 while ((i < kg)
72 && (m->lignes[i][j].vect[(cl - 1) / WL] < MMC[(cl - 1) % WL]))
73 i++;
74 if (i < kg) { /* pivot trouve ... */
75 ExchangeVect (m, rang, i);
76 for (i = rang + 1; i < kg; i++) {
77 if (m->lignes[i][j].vect[(cl - 1) / WL] & MMC[(cl - 1) % WL])
78 XorVect (m, i, rang, j, m->t);
79 }
80 rang++;
81 } else { /* pas de pivot trouve ... ==> pas de plein
82 rang ... */
83
84 *gr = j; /* no de groupe ou il y a echec moins un */
85 return FALSE; /* c'est j car on indexe a partir de 0 . */
86 }
87 cl++;
88 }
89 }
90 return TRUE; /* on a trouve tous les pivots ==> plein
91 rang ! */
92 }
93
94
CompleteElimination(Matrix * m,int nblignes,int l,int t)95 int CompleteElimination (Matrix * m, int nblignes, int l, int t)
96 {
97 int i, j, cl, rang;
98
99 rang = 0;
100
101 j = 0;
102 while (j < t) {
103 cl = 0;
104 while (cl < l) {
105 /* On cherche dans la j-eme colonne, commencant a la */
106 /* rang-eme ligne, la premiere entree dont le bit le plus */
107 /* significatif est non nul. Bref, on cherche un pivot. */
108 i = rang;
109 while ((i < nblignes)
110 && ((((m->lignes)[i])[j]).vect[(cl) / WL] < MMC[(cl) % WL]))
111 i++;
112 if (i < nblignes) { /* pivot trouve ... */
113 ExchangeVect (m, rang, i);
114 for (i = 0; i < nblignes; i++)
115 if (i != rang)
116 if ((((m->lignes)[i])[j]).vect[(cl) / WL] & MMC[(cl) % WL])
117 XorVect (m, i, rang, 0, m->t);
118
119 rang++;
120 if (rang == nblignes)
121 return rang;
122 } else
123 return rang;
124 cl++;
125 }
126 j++;
127 }
128 return rang;
129 }
130
131
GaussianElimination(Matrix * m,int nblignes,int l,int t)132 int GaussianElimination (Matrix * m, int nblignes, int l, int t)
133 {
134 int i, j, cl, rang;
135
136 rang = 0;
137
138 j = 0;
139 while (j < t) {
140 cl = 0;
141 while (cl < l) {
142 /* On cherche dans la j-eme colonne, commencant a la */
143 /* rang-eme ligne, la premiere entree dont le bit le plus */
144 /* significatif est non nul. Bref, on cherche un pivot. */
145 i = rang;
146 while ((i < nblignes)
147 && ((((m->lignes)[i])[j]).vect[(cl) / WL] < MMC[(cl) % WL]))
148 i++;
149 if (i < nblignes) { /* pivot trouve ... */
150 ExchangeVect (m, rang, i);
151 for (i = rang + 1; i < nblignes; i++)
152 if ((((m->lignes)[i])[j]).vect[(cl) / WL] & MMC[(cl) % WL])
153 XorVect (m, i, rang, 0, m->t);
154
155 rang++;
156 if (rang == nblignes)
157 return rang;
158 }
159 cl++;
160 }
161 j++;
162 }
163 return rang;
164 }
165
166
SpecialGaussianElimination(Matrix * m,int nblignes,int l,int t,int * indices)167 int SpecialGaussianElimination (Matrix * m, int nblignes, int l, int t,
168 int *indices)
169 {
170 int i, j, cl, rang;
171
172 rang = 0;
173
174 j = 0;
175 while (j < t) {
176 cl = 0;
177 while (cl < l) {
178 /* On cherche dans la j-eme colonne, commencant a la */
179 /* rang-eme ligne, la premiere entree dont le bit le plus */
180 /* significatif est non nul. Bref, on cherche un pivot. */
181 i = rang;
182 while ((i < nblignes)
183 && ((((m->lignes)[i])[indices[j]]).vect[(cl) / WL] <
184 MMC[(cl) % WL]))
185 i++;
186 if (i < nblignes) { /* pivot trouve ... */
187 ExchangeVect (m, rang, i);
188 for (i = rang + 1; i < nblignes; i++)
189 if ((((m->lignes)[i])[indices[j]]).vect[(cl) /
190 WL] & MMC[(cl) % WL])
191 XorVect (m, i, rang, 0, m->t);
192
193 rang++;
194 if (rang == nblignes)
195 return rang;
196 }
197 cl++;
198 }
199 j++;
200 }
201 return rang;
202 }
203
204
MultMatrixByBV(BitVect * A,Matrix * M,BitVect * B)205 void MultMatrixByBV (BitVect * A, Matrix * M, BitVect * B)
206 {
207 int i, j, res;
208 if (M->l < B->n * WL) {
209 printf ("Error in MultMatrixByBV(): sizes do not match\n");
210 exit (1);
211 }
212 if (A->n * WL < M->nblignes) {
213 printf ("Error in MultMatrixByBV(): sizes do not match\n");
214 exit (1);
215 }
216 if (M->t != 1) {
217 printf ("Error in MultMatrixByBV(): Not implemented for M->t > 1\n");
218 exit (1);
219 }
220 PutBVToZero (A);
221 for (i = 0; i < M->nblignes; i++) {
222 res = 0;
223 for (j = 0; j < M->l; j++)
224 res += ValBitBV (&(M->lignes[i][0]), j) * ValBitBV (B, j);
225 res %= 2;
226 PutBitBV (A, i, res);
227 }
228 }
229
230
231 /* ********************************************************************** */
232 /* int GetBitBV (BitVect A, int noBit) */
233 /* Fonction qui permet de prendre la valeur du noBit-ieme bit */
234 /* (l'indexation commence au bit 0) */
235 /* ********************************************************************** */
ValBitBV(BitVect * A,int noBit)236 int ValBitBV (BitVect * A, int noBit)
237 {
238 int k;
239 unsigned long mask;
240 k = noBit / WL;
241 mask = 0x80000000UL >> (noBit - k * WL);
242 if (A->vect[k] & mask)
243 return 1;
244 else
245 return 0;
246 }
247
248
249 /* ********************************************************************** */
250 /* void SetBitBV (BitVect A, int noBit, int valBit) */
251 /* Fonction qui permet de rendre la valeur du noBit-ieme bit `a valBit */
252 /* (l'indexation commence au bit 0) (valBit=1 ou valBit=0) */
253 /* ********************************************************************** */
PutBitBV(BitVect * A,int noBit,int valBit)254 void PutBitBV (BitVect * A, int noBit, int valBit)
255 {
256 int k;
257 unsigned long mask;
258
259 k = noBit / WL;
260 if (valBit == 1) {
261 mask = 0x80000000UL >> (noBit - k * WL);
262 A->vect[k] |= mask;
263 } else {
264 mask = 0xffffffffUL ^ (0x80000000UL >> (noBit - k * WL));
265 A->vect[k] &= mask;
266 }
267 }
268
269
270 /* ********************************************************************** */
271 /* SetBVToZero( BitVect A) */
272 /* Initialise le vecteur de bit a zero */
273 /* ********************************************************************** */
PutBVToZero(BitVect * A)274 void PutBVToZero (BitVect * A)
275 {
276 int i;
277 for (i = 0; i < A->n; i++)
278 A->vect[i] = 0UL;
279 }
280
281
282 /* ********************************************************************** */
283 /* void CopyBV(BitVect A, BitVect B) */
284 /* Copie le contenu de B dans A (A=B) */
285 /* ********************************************************************** */
CopyBV(BitVect * A,BitVect * B)286 void CopyBV (BitVect * A, BitVect * B)
287 {
288 int i;
289
290 if (A->n != B->n) {
291 printf
292 ("Error in CopyBV(): vectors of different dimensions! (%d and %d bits)\n",
293 A->n * WL, B->n * WL);
294 exit (1);
295 }
296
297 if (B->n == 0) {
298 printf ("Nothing to copy!\n");
299 exit (1);
300 }
301 for (i = 0; i < B->n; i++)
302 A->vect[i] = B->vect[i];
303 }
CopyBVPart(BitVect * A,BitVect * B,int l)304 void CopyBVPart (BitVect * A, BitVect * B, int l)
305 {
306
307 int i, n;
308
309 n = (l - 1) / WL + 1;
310
311 if (A->n < n) {
312 printf ("Error in CopyBVPart() : The vector A is not large enough!\n");
313 exit (1);
314 }
315 if (B->n == 0) {
316 printf ("Nothing to copy!\n");
317 exit (1);
318 }
319
320 for (i = 0; i < n; i++)
321 A->vect[i] = B->vect[i];
322
323 if (l % WL) {
324 BitVect m;
325 AllocBV (&m, A->n * WL);
326 Mask (&m, l);
327 ANDBVSelf (A, &m);
328 FreeBV (&m);
329 }
330 }
331
332
333 /* ********************************************************************** */
334 /* void EgalBV(BitVect A, BitVect B) */
335 /* Compare le contenu de B avec celui de A. Retourne TRUE si les deux */
336 /* contiennent la meme information. */
337 /* ********************************************************************** */
CompareBV(BitVect * A,BitVect * B)338 lebool CompareBV (BitVect * A, BitVect * B)
339 {
340
341 int i;
342
343 if (A->n != B->n) {
344 printf ("Error in EgalBV(): Vectors of different sizes\n");
345 exit (1);
346 }
347
348 for (i = 0; i < A->n; i++)
349 if (A->vect[i] != B->vect[i])
350 return FALSE;
351 return TRUE;
352 }
353
354
BVisZero(BitVect * A)355 lebool BVisZero (BitVect * A)
356 {
357 int j = 0;
358 while (j < A->n)
359 if (A->vect[j++] != 0UL)
360 return FALSE;
361 return TRUE;
362 }
363
364
365 /* ********************************************************************** */
366 /* void XORBV(BitVect A, BitVect B, BitVect C) */
367 /* Cette fonction effectue A = B ^ C */
368 /* ********************************************************************** */
XORBV(BitVect * A,BitVect * B,BitVect * C)369 void XORBV (BitVect * A, BitVect * B, BitVect * C)
370 {
371
372 int i;
373
374 if ((A->n != B->n) || (B->n != C->n)) {
375 printf ("Error in XORBV(): Vectors of different sizes\n");
376 exit (1);
377 }
378
379 for (i = 0; i < B->n; i++)
380 A->vect[i] = B->vect[i] ^ C->vect[i];
381 }
382
383
384 /* ********************************************************************** */
385 /* void XOR2BV(BitVect A, BitVect B, BitVect C, BitVect D) */
386 /* Cette fonction effectue A = B ^ C ^ D */
387 /* ********************************************************************** */
XOR2BV(BitVect * A,BitVect * B,BitVect * C,BitVect * D)388 void XOR2BV (BitVect * A, BitVect * B, BitVect * C, BitVect * D)
389 {
390
391 int i;
392
393 if ((A->n != B->n) || (B->n != C->n) || (C->n != D->n)) {
394 printf ("Error in XOR2BV(): Vectors of different sizes\n");
395 exit (1);
396 }
397
398 for (i = 0; i < B->n; i++)
399 A->vect[i] = B->vect[i] ^ C->vect[i] ^ D->vect[i];
400 }
401
402
403 /* ********************************************************************** */
404 /* void ANDBV(BitVect A, BitVect B, BitVect C) */
405 /* Cette fonction effectue A = B & C */
406 /* ********************************************************************** */
ANDBV(BitVect * A,BitVect * B,BitVect * C)407 void ANDBV (BitVect * A, BitVect * B, BitVect * C)
408 {
409
410 int i;
411
412 if ((A->n != B->n) || (B->n != C->n)) {
413 printf ("Error in ANDBV(): Vectors of different sizes\n");
414 exit (1);
415 }
416
417 for (i = 0; i < B->n; i++)
418 A->vect[i] = B->vect[i] & C->vect[i];
419 }
420
421
ANDBVMask(BitVect * A,BitVect * B,int t)422 void ANDBVMask (BitVect * A, BitVect * B, int t)
423 {
424 int n, m, j;
425 if (A->n != B->n) {
426 printf ("Error in ANDBVMask(): Vectors of different sizes\n");
427 exit (1);
428 }
429
430 if (t > B->n * WL)
431 CopyBV (A, B);
432 else if (t == 0)
433 PutBVToZero (A);
434 else {
435 n = t / WL;
436 m = t - n * WL;
437 for (j = 0; j < n; j++) {
438 A->vect[j] = B->vect[j];
439
440 }
441 if (m != 0) {
442 A->vect[j] = B->vect[j] & (0xffffffffUL << (WL - m));
443 j++;
444 }
445 /* printf("n=%d j=%d %d m=%d ",n,j,A->n,m); */
446 for (; j < A->n; j++) {
447 A->vect[j] = 0UL;
448 }
449 }
450 }
451
452
ANDBVInvMask(BitVect * A,BitVect * B,int t)453 void ANDBVInvMask (BitVect * A, BitVect * B, int t)
454 {
455 int n, m, j;
456 if (A->n != B->n) {
457 printf ("Error in ANDBV(): Vectors of different sizes\n");
458 exit (1);
459 }
460 if (t > B->n * WL)
461 PutBVToZero (A);
462 else if (t == 0)
463 CopyBV (A, B);
464 else {
465 n = t / WL;
466 m = t - n * WL;
467 for (j = 0; j < n; j++)
468 A->vect[j] = 0UL;
469 if (m == 0)
470 A->vect[j] = B->vect[j];
471 else {
472 A->vect[j] = B->vect[j] & (0xffffffffUL >> m);
473
474 }
475 j++;
476 for (; j < A->n; j++)
477 A->vect[j] = B->vect[j];
478 }
479 }
480
481
482
483
484 /* ********************************************************************** */
485 /* void ANDBVSelf(BitVect A, BitVect B) */
486 /* Cette fonction effectue A &=B */
487 /* ********************************************************************** */
ANDBVSelf(BitVect * A,BitVect * B)488 void ANDBVSelf (BitVect * A, BitVect * B)
489 {
490 int i;
491
492 if ((A->n != B->n)) {
493 printf ("Error in ANDBVSelf(): Vectors of different sizes\n");
494 exit (1);
495 }
496
497 for (i = 0; i < B->n; i++)
498 A->vect[i] &= B->vect[i];
499 }
500
501
502 /* ********************************************************************** */
503 /* void XORBVSelf(BitVect A, BitVect B) */
504 /* Cette fonction effectue A ^=B */
505 /* ********************************************************************** */
XORBVSelf(BitVect * A,BitVect * B)506 void XORBVSelf (BitVect * A, BitVect * B)
507 {
508 int i;
509
510 if ((A->n != B->n)) {
511 printf ("Error in XORBVSelf(): Vectors of different sizes\n");
512 exit (1);
513 }
514
515 for (i = 0; i < B->n; i++)
516 A->vect[i] ^= B->vect[i];
517 }
518
519
520 /* ********************************************************************** */
521 /* void BVLShift( BitVect R, BitVect A, int n ) */
522 /* Effectue : R = A << n ; */
523 /* ********************************************************************** */
BVLShift(BitVect * R,BitVect * A,int n)524 void BVLShift (BitVect * R, BitVect * A, int n)
525 {
526 int i;
527 int WLmn;
528 unsigned long temp;
529
530 if ((R->n != A->n)) {
531 printf ("Error in BVLShift(): Vectors of different sizes\n");
532 exit (1);
533 }
534
535 for (i = 0; i < A->n; i++)
536 R->vect[i] = A->vect[i];
537 while (n >= 32) {
538 for (i = 1; i < A->n; i++)
539 R->vect[i - 1] = R->vect[i];
540 R->vect[A->n - 1] = 0UL;
541 n -= 32;
542 }
543 if (n > 0) {
544 WLmn = WL - n;
545 R->vect[0] <<= n;
546 for (i = 1; i < A->n; i++) {
547 temp = R->vect[i] >> WLmn;
548 R->vect[i - 1] |= temp;
549 R->vect[i] <<= n;
550 }
551 }
552
553 }
554
555
556 /* ********************************************************************** */
557 /* void BVRShift( BitVect R, BitVect A, int n ) */
558 /* Effectue : R = A >> n ; */
559 /* ********************************************************************** */
BVRShift(BitVect * R,BitVect * A,int n)560 void BVRShift (BitVect * R, BitVect * A, int n)
561 {
562 int i;
563 int WLmn;
564 unsigned long temp;
565 if ((R->n != A->n)) {
566 printf ("Error in BVRShift(): Vectors of different sizes\n");
567 exit (1);
568 }
569
570 for (i = 0; i < A->n; i++)
571 R->vect[i] = A->vect[i];
572 while (n >= 32) {
573 for (i = A->n; i > 1; i--)
574 R->vect[i - 1] = R->vect[i - 2];
575 R->vect[0] = 0UL;
576 n -= 32;
577 }
578 if (n > 0) {
579 WLmn = WL - n;
580 R->vect[A->n - 1] >>= n;
581 for (i = A->n - 2; i >= 0; i--) {
582 temp = R->vect[i] << WLmn;
583 R->vect[i + 1] |= temp;
584 R->vect[i] >>= n;
585 }
586 }
587 }
588
589
590 /* ********************************************************************** */
591 /* void BVLShiftSelf( BitVect R, int n ) */
592 /* Effectue : R <<= n ; */
593 /* ********************************************************************** */
BVLShiftSelf(BitVect * R,int n)594 void BVLShiftSelf (BitVect * R, int n)
595 {
596 int i;
597 int WLmn;
598 unsigned long temp;
599
600 while (n >= 32) {
601 for (i = 1; i < R->n; i++)
602 R->vect[i - 1] = R->vect[i];
603 R->vect[R->n - 1] = 0UL;
604 n -= 32;
605 }
606 if (n > 0) {
607 WLmn = WL - n;
608 R->vect[0] <<= n;
609 for (i = 1; i < R->n; i++) {
610 temp = R->vect[i] >> WLmn;
611 R->vect[i - 1] |= temp;
612 R->vect[i] <<= n;
613 }
614 }
615 }
616
617
618 /* ********************************************************************** */
619 /* void BVLS1Self( BitVect R ) */
620 /* Effectue : R <<= 1 ; */
621 /* Version specialisee de BVLShiftSelf pour utilisation frequente. */
622 /* ********************************************************************** */
BVLS1Self(BitVect * R)623 void BVLS1Self (BitVect * R)
624 {
625 int i;
626 R->vect[0] <<= 1;
627 for (i = 1; i < R->n; i++) {
628 if (R->vect[i] & MC)
629 R->vect[i - 1] |= 0x1UL;
630 R->vect[i] <<= 1;
631 }
632 }
633
634
635 /* ********************************************************************** */
636 /* void BVRShiftSelf( BitVect R, int n ) */
637 /* Effectue : R >>= n ; */
638 /* ********************************************************************** */
BVRShiftSelf(BitVect * R,int n)639 void BVRShiftSelf (BitVect * R, int n)
640 {
641 int i;
642 int WLmn;
643 unsigned long temp;
644
645 while (n >= 32) {
646 for (i = R->n - 1; i > 0; i--)
647 R->vect[i] = R->vect[i - 1];
648 R->vect[0] = 0UL;
649 n -= 32;
650 }
651 if (n > 0) {
652 WLmn = WL - n;
653 R->vect[R->n - 1] >>= n;
654 for (i = R->n - 2; i >= 0; i--) {
655 temp = R->vect[i] << WLmn;
656 R->vect[i + 1] |= temp;
657 R->vect[i] >>= n;
658 }
659 }
660 }
661
662
663 /* ********************************************************************** */
664 /* void invertBV(Bitvect A) */
665 /* fait A ~=A */
666 /* ********************************************************************** */
InverseBV(BitVect * A)667 void InverseBV (BitVect * A)
668 {
669 int i;
670 for (i = 0; i < A->n; i++)
671 A->vect[i] = ~A->vect[i];
672 }
673
674
675 /* ********************************************************************** */
676 /* lebool CheckCD( BitVect ds1, BitVect ds2 ) */
677 /* Verifie si les ensembles ds1 et ds2 ont des bits communs. */
678 /* ********************************************************************** */
VerifBitsCommuns(BitVect * ds1,BitVect * ds2)679 lebool VerifBitsCommuns (BitVect * ds1, BitVect * ds2)
680 {
681 int i;
682 unsigned long temp = 0UL;
683 if ((ds1->n != ds2->n)) {
684 printf ("Error in VerifBitsCommuns(): Vectors of different sizes\n");
685 exit (1);
686 }
687 for (i = 0; i < ds1->n; i++)
688 temp |= (ds1->vect[i] & ds2->vect[i]);
689 if (temp)
690 return TRUE;
691 else
692 return FALSE;
693 }
694
695
696 /* -------------------------------------------- */
697 /* Fonctions pour la manipulation des matrices. */
698 /* -------------------------------------------- */
BVCanonic(BitVect * A,int l)699 void BVCanonic (BitVect * A, int l)
700 {
701 int n;
702 PutBVToZero (A);
703 n = l / WL;
704 if (n > A->n) {
705 printf
706 ("Error in BVCanonic(): vector A is not long enough to store BVCanonic[%d].\n",
707 l);
708 exit (1);
709 }
710 A->vect[n] = 0x80000000UL >> (l - n * WL);
711 }
712
Mask(BitVect * A,int l)713 void Mask (BitVect * A, int l)
714 {
715
716 InvMask (A, l);
717 InverseBV (A);
718 }
719
InvMask(BitVect * A,int l)720 void InvMask (BitVect * A, int l)
721 {
722 AllOnes (A);
723 BVRShiftSelf (A, l);
724 }
725
AllOnes(BitVect * A)726 void AllOnes (BitVect * A)
727 {
728 int i;
729 for (i = 0; i < A->n; i++)
730 A->vect[i] = 0xffffffffUL;
731 }
732
AllocBV(BitVect * A,int l)733 void AllocBV (BitVect * A, int l)
734 {
735 int n;
736 n = (l - 1) / WL + 1;
737 A->vect = (unsigned long *) calloc ((size_t) n, sizeof (unsigned long));
738 A->n = n;
739 }
740
FreeBV(BitVect * A)741 void FreeBV (BitVect * A)
742 {
743 if (A->vect != NULL) {
744 free (A->vect);
745 }
746 A->n = 0;
747 }
748
AllocMat(Matrix * m,int nblines,int l,int t)749 void AllocMat (Matrix * m, int nblines, int l, int t)
750 {
751 int i, j;
752 m->lignes = (BitVect **) calloc ((size_t) nblines, sizeof (BitVect *));
753 for (i = 0; i < nblines; i++) {
754 if (!(m->lignes[i] = (BitVect *) calloc ((size_t) t, sizeof (BitVect)))) {
755 printf ("\n*** Memoire insuffisante pour AllocMat() ! nl=%d***\n",
756 nblines);
757 exit (1);
758 }
759 for (j = 0; j < t; j++)
760 AllocBV (&(m->lignes[i][j]), l);
761
762 }
763 m->nblignes = nblines;
764 m->t = t;
765 m->l = l;
766
767 }
768
769
770 /* ********************************************************************** */
771 /* void FreeSpace( Matrix m, int nl ) */
772 /* Libere l'espace des nl vecteurs de la matrice m. */
773 /* ********************************************************************** */
FreeMat(Matrix * m)774 void FreeMat (Matrix * m)
775 {
776 int i, j;
777 for (i = 0; i < m->nblignes; i++) {
778 for (j = 0; j < m->t; j++)
779 FreeBV (&(m->lignes[i][j]));
780 free (m->lignes[i]);
781 }
782 free (m->lignes);
783 m->nblignes = 0;
784 m->l = 0;
785 m->t = 0;
786
787 }
788
789
790 /* ********************************************************************** */
791 /* void CopyMat( Matrix m, Matrix ms, int nl, int t ) */
792 /* Cette procedure sert a copier les t premiers BitVect des nl premieres */
793 /* lignes de la matrice ms dans la matrice de travail m. */
794 /* ********************************************************************** */
CopyMat(Matrix * m,Matrix * ms,int nl,int t)795 void CopyMat (Matrix * m, Matrix * ms, int nl, int t)
796 {
797 int i, j;
798 if (m == NULL) {
799 AllocMat (m, ms->nblignes, ms->l, ms->t);
800 } else if ((ms->nblignes < nl) || (ms->t < t)) {
801 printf ("Error in CopyMat(): source matrix too small %d\n",
802 ms->nblignes / ms->t);
803 exit (1);
804 } else if ((m->nblignes < nl) || (m->t < t)) {
805 printf ("Error in CopyMat(): destination matrix too small\n");
806 exit (1);
807 }
808 for (i = 0; i < nl; i++)
809 for (j = 0; j < t; j++) {
810 CopyBV (&(m->lignes[i][j]), &(ms->lignes[i][j]));
811 }
812 }
813
814
815 /* ********************************************************************** */
816 /* void CopyNTupleMat( Matrix m, Matrix ms, int nl, */
817 /* int *colonnes, int t ) */
818 /* Cette procedure sert a copier les t-1 BitVect indiqu'es par le vecteur */
819 /* *colonnes plus la colonne 0 des nl premieres lignes de la matrice ms */
820 /* dans la matrice de travail m. */
821 /* ********************************************************************** */
CopyNTupleMat(Matrix * m,Matrix * ms,int nl,int * colonnes,int t)822 void CopyNTupleMat (Matrix * m, Matrix * ms, int nl, int *colonnes, int t)
823 {
824
825 int i, j, k, n;
826
827 if (m == NULL)
828 AllocMat (m, ms->nblignes, ms->l, t);
829 else {
830 if ((ms->nblignes != m->nblignes) || (ms->l != m->l))
831 printf ("Error in CopieNTupleMat(): matrices of different sizes\n");
832 }
833 n = (ms->l - 1) / WL;
834 for (i = 0; i < nl; i++) {
835 for (k = 0; k <= n; k++)
836 (m->lignes[i])[0].vect[k] = (ms->lignes[i])[0].vect[k];
837 for (j = 1; j < t; j++)
838 for (k = 0; k <= n; k++)
839 (m->lignes[i])[j].vect[k] =
840 (ms->lignes[i])[colonnes[j - 1]].vect[k];
841
842 }
843 }
844
845
846 /* ********************************************************************** */
847 /* void SwapVect( Matrix m, int i, int j ) */
848 /* Pour interchanger les lignes i et j de la matrice de travail m. */
849 /* ********************************************************************** */
ExchangeVect(Matrix * m,int i,int j)850 void ExchangeVect (Matrix * m, int i, int j)
851 {
852 BitVect *p;
853 if (i != j) {
854 p = m->lignes[i];
855 m->lignes[i] = m->lignes[j];
856 m->lignes[j] = p;
857 }
858 }
859
860
TransposeMatrices(Matrix * T,Matrix * M,int mmax,int smax,int L)861 void TransposeMatrices (Matrix * T, Matrix * M, int mmax, int smax, int L)
862 {
863
864 int s, l, m;
865
866 for (s = 0; s < smax; s++)
867 for (l = 0; l < L; l++) {
868 PutBVToZero (&T->lignes[l][s]);
869 for (m = 0; m < mmax; m++) {
870 /* printf("m=%d l=%d s=%d\n",m,l,s);fflush(stdout); */
871 if (M->lignes[m][s].vect[0] & (0x80000000UL >> l)) {
872 T->lignes[l][s].vect[0] |= (0x80000000UL >> m);
873 }
874 }
875 }
876 }
877
878
879 /* ********************************************************************** */
880 /* void XorVect( Matrix m, */
881 /* int r, int s, int min, int max ) */
882 /* Effectue un Xor entre la s-eme et la r-eme ligne de la matrice de */
883 /* travail m pour les colonnes (BitVect) min a max-1 seulement. */
884 /* Le resultat est mis dans la r-eme ligne. ( m[r] ^= m[s] ) */
885 /* ********************************************************************** */
XorVect(Matrix * m,int r,int s,int min,int max)886 void XorVect (Matrix * m, int r, int s, int min, int max)
887 {
888 int j;
889 for (j = min; j < max; j++)
890 XORBVSelf (&(m->lignes[r][j]), &(m->lignes[s][j]));
891 }
892
893
894 /* ********************************************************************** */
895 /* void displaymat(Matrix m, int t, int l, int kg) */
896 /* Affiche la matrice m sur kg lignes par t x l colonnes */
897 /* ********************************************************************** */
DispMat(Matrix * m,int t,int l,int kg,lebool mathematica)898 void DispMat (Matrix * m, int t, int l, int kg, lebool mathematica)
899 {
900 int i, j;
901
902 i = kg;
903
904 printf ("\n");
905 if (mathematica)
906 printf ("{");
907 for (i = 0; i < kg; i++) {
908 if (!mathematica)
909 printf ("[");
910 for (j = 0; j < t; j++) {
911 DispBitVect (&(m->lignes[i][j]), l, mathematica);
912 }
913 if (mathematica) {
914 if (i != kg - 1)
915 printf (",\n");
916 else
917 printf ("}\n");
918 } else
919 printf ("]\n");
920 }
921 printf ("\n\n");
922
923 }
924
925
926 /* ********************************************************************** */
927 /* void displaybitvect(BitVect A, int l) */
928 /* Affiche le BitVect A sur l bits seulement */
929 /* ********************************************************************** */
DispBitVect(BitVect * A,int l,int mathematica)930 void DispBitVect (BitVect * A, int l, int mathematica)
931 {
932 int j;
933 unsigned Un;
934 Un = 1UL;
935 j = 0;
936 if (mathematica) {
937 printf ("{");
938 while (j < l - 1) {
939 printf ("%ld,",
940 (A->vect[j / WL] >> (((WL * A->n) - j - 1) % WL)) & Un);
941 j++;
942 }
943 printf ("%ld}", (A->vect[j / WL] >> (((WL * A->n) - j - 1) % WL)) & Un);
944 } else
945 while (j < l) {
946 printf ("%ld",
947 (A->vect[j / WL] >> (((WL * A->n) - j - 1) % WL)) & Un);
948 j++;
949 }
950 }
951
952
953 /* ********************************************************************** */
954 /* MultMatrixByMatrix ( Matrix *A, Matrix *B, Matrix *C) */
955 /* Fait la multiplication matricielle A = B x C */
956 /* ********************************************************************** */
957
MultMatrixByMatrix(Matrix * A,Matrix * B,Matrix * C)958 void MultMatrixByMatrix (Matrix * A, Matrix * B, Matrix * C)
959 {
960 int i, j;
961 if (B->l != C->nblignes) {
962 printf ("Tailles de matrices non-compatibles, kaput.\n");
963 exit (1);
964 }
965
966 if (A->nblignes != B->nblignes || A->l != C->l) {
967 printf ("Matrice preallouee de mauvaise taille.\n");
968 exit (1);
969 }
970
971 for (i = 0; i < A->nblignes; i++)
972 PutBVToZero (A->lignes[i]);
973
974 for (i = 0; i < B->nblignes; i++)
975 for (j = 0; j < B->l; j++) {
976 if (ValBitBV (B->lignes[i], j))
977 XORBVSelf (A->lignes[i], C->lignes[j]);
978 }
979 }
980
981
982 /* ********************************************************************** */
983 /* MatrixTwoPow ( Matrix *A, Matrix *B, unsigned int e) */
984 /* Fait l'exponentiation : A = B^(2^e) */
985 /* ********************************************************************** */
986
MatrixTwoPow(Matrix * A,Matrix * B,unsigned int e)987 void MatrixTwoPow (Matrix * A, Matrix * B, unsigned int e)
988 {
989 unsigned int i;
990 Matrix tempMatrix;
991 Matrix *AA = &tempMatrix;
992
993 if (B->nblignes != B->l) {
994 printf ("Matrice non carree.\n");
995 exit (1);
996 }
997 if (A->nblignes != B->nblignes || A->l != B->l) {
998 printf ("Matrice preallouee de mauvaise taille.\n");
999 exit (1);
1000 }
1001
1002 AllocMat (AA, B->nblignes, B->l, 1);
1003
1004 if (e == 0) {
1005 CopyMat (A, B, B->nblignes, 1);
1006 return;
1007 }
1008 /* A = B^(2^1) */
1009 MultMatrixByMatrix (A, B, B);
1010
1011 for (i = 1; i < e - 1; i += 2) {
1012 /* AA = A * A */
1013 MultMatrixByMatrix (AA, A, A);
1014
1015 /* A = AA * AA */
1016 MultMatrixByMatrix (A, AA, AA);
1017 }
1018
1019 if (i == e - 1) {
1020 /* AA = A * A */
1021 MultMatrixByMatrix (AA, A, A);
1022 CopyMat (A, AA, AA->nblignes, 1);
1023 }
1024
1025 FreeMat (AA);
1026 }
1027
1028 /* ********************************************************************** */
1029 /* MatrixPow ( Matrix *A, Matrix *B, unsigned int e) */
1030 /* Fait l'exponentiation : A = B^e */
1031 /* ********************************************************************** */
1032
1033 #ifdef USE_LONGLONG
MatrixPow(Matrix * A,Matrix * B,longlong e)1034 void MatrixPow (Matrix * A, Matrix * B, longlong e)
1035 #else
1036 void MatrixPow (Matrix * A, Matrix * B, long e)
1037 #endif
1038 {
1039 int i;
1040 Matrix C;
1041 Matrix D;
1042
1043 if (B->nblignes != B->l) {
1044 printf ("Matrice non carree.\n");
1045 exit (1);
1046 }
1047
1048 if (A->nblignes != B->nblignes || A->l != B->l) {
1049 printf ("Matrice preallouee de mauvaise taille.\n");
1050 exit (1);
1051 }
1052
1053 AllocMat (&C, B->nblignes, B->l, 1);
1054
1055 if (e < 0) {
1056 InverseMatrix (&C, B);
1057 MatrixPow (A, &C, -e);
1058 FreeMat (&C);
1059 return;
1060 }
1061 AllocMat (&D, B->nblignes, B->l, 1);
1062
1063 /* A = I */
1064 for (i = 0; i < A->nblignes; i++)
1065 BVCanonic (A->lignes[i], i);
1066
1067 /* C = B^1 */
1068 CopyMat (&C, B, B->nblignes, 1);
1069
1070 while (e) {
1071 if (e & 1) {
1072 CopyMat (&D, A, B->nblignes, 1);
1073 MultMatrixByMatrix (A, &D, &C);
1074 }
1075
1076 e >>= 1;
1077
1078 if (e) {
1079 CopyMat (&D, &C, B->nblignes, 1);
1080 MultMatrixByMatrix (&C, &D, &D);
1081 }
1082 }
1083
1084 FreeMat (&C);
1085 FreeMat (&D);
1086 }
1087
1088
1089 /* FIN vectorsF2.c */
1090